SnpEff变异注释（SnpEff/SnpSift/ANNOVAR/VEP变异功能注释与过滤）¶

# ===== SnpEff安装 =====
conda install -c bioconda snpeff snpsift

# 或手动安装
wget https://snpeff.blob.core.windows.net/versions/snpEff_latest_core.zip
unzip snpEff_latest_core.zip

# ===== 下载注释数据库 =====
# 查看可用数据库
snpEff databases | grep "GRCh38"

# 下载人类GRCh38数据库
snpEff download -v GRCh38.105

# ===== 运行SnpEff注释 =====
snpEff ann \
  -v \
  GRCh38.105 \
  input_variants.vcf.gz \
  > annotated.vcf

# 带统计报告
snpEff ann \
  -v \
  -stats snpeff_stats.html \
  -csvStats snpeff_stats.csv \
  GRCh38.105 \
  input_variants.vcf.gz | \
  bgzip > annotated.vcf.gz
tabix -p vcf annotated.vcf.gz

# ===== 高级参数 =====
snpEff ann \
  -v \
  -canon \                    # 只注释canonical转录本
  -nodownstream \             # 不注释下游区域
  -noUpstream \               # 不注释上游区域
  -noIntergenic \             # 不注释基因间区
  -no-intron \                # 不注释内含子
  -noLog \
  -lof \                      # 预测Loss of Function
  -hgvs \                     # 使用HGVS命名法
  GRCh38.105 \
  input.vcf.gz > annotated.vcf

# ===== 查看注释结果 =====
# SnpEff在VCF的INFO字段添加ANN注释
# 格式：ANN=Allele|Annotation|Impact|Gene_Name|Gene_ID|Feature_Type|Feature_ID|Transcript_BioType|Rank|HGVS.c|HGVS.p|cDNA.pos/cDNA.length|CDS.pos/CDS.length|AA.pos/AA.length|Distance|ERRORS_WARNINGS_INFO

# 示例：
# ANN=T|missense_variant|MODERATE|BRCA1|ENSG00000012048|transcript|ENST00000357654.9|protein_coding|10/23|c.1129G>A|p.Gly377Ser|1234/5592|1129/5451|377/1817||

# 查看HTML报告
# 包含变异效应统计、各类变异数量、Ts/Tv比率等

# ===== SnpSift数据库注释 =====

# 1. dbSNP注释（添加rsID）
snpSift annotate \
  -v \
  dbSnp156.vcf.gz \
  annotated.vcf.gz > annotated_dbsnp.vcf

# 2. ClinVar注释
snpSift annotate \
  -v \
  -info CLNSIG,CLNDN,CLNREVSTAT \
  clinvar.vcf.gz \
  annotated_dbsnp.vcf > annotated_clinvar.vcf

# 3. dbNSFP注释（包含SIFT, PolyPhen2, CADD等多种预测）
# 下载dbNSFP（4.x系列最新为4.9a，5.x系列最新为5.3）
# 下载地址参见 https://sites.google.com/site/jpopgen/dbNSFP 或 https://www.dbnsfp.org/
wget https://dbnsfp.s3.amazonaws.com/dbNSFP4.9a.zip
# 准备dbNSFP for SnpSift
zcat dbNSFP4.9a_variant.chr*.gz | \
  awk 'NR==1 || $1!="#"' | bgzip > dbNSFP4.9a.txt.gz
tabix -s 1 -b 2 -e 2 dbNSFP4.9a.txt.gz

snpSift dbnsfp \
  -v \
  -db dbNSFP4.9a.txt.gz \
  -f SIFT_pred,Polyphen2_HDIV_pred,Polyphen2_HVAR_pred,\
CADD_phred,REVEL_score,MutationTaster_pred,\
GERP_RS,phyloP100way_vertebrate,gnomAD_genomes_AF \
  annotated_clinvar.vcf > annotated_full.vcf

# 4. gnomAD人群频率注释（推荐使用v4.1，v3.1.2的AF计算有bug已在v4.1修复）
snpSift annotate \
  -v \
  -info AF,AF_eas,AF_nfe,AF_afr \
  gnomad.genomes.v4.1.sites.vcf.gz \
  annotated_full.vcf > annotated_gnomad.vcf

# ===== SnpSift过滤 =====

# 基于Impact过滤（只保留HIGH和MODERATE影响）
snpSift filter \
  "(ANN[*].IMPACT = 'HIGH') | (ANN[*].IMPACT = 'MODERATE')" \
  annotated_full.vcf > high_moderate.vcf

# 基于人群频率过滤（稀有变异：MAF < 1%）
snpSift filter \
  "(!(exists gnomAD_genomes_AF)) | (gnomAD_genomes_AF < 0.01)" \
  high_moderate.vcf > rare_damaging.vcf

# 基于致病性预测过滤
snpSift filter \
  "(SIFT_pred = 'D') & (Polyphen2_HDIV_pred = 'D') & (CADD_phred >= 20)" \
  rare_damaging.vcf > predicted_pathogenic.vcf

# 基于ClinVar过滤（已知致病）
snpSift filter \
  "(CLNSIG = 'Pathogenic') | (CLNSIG = 'Likely_pathogenic')" \
  annotated_clinvar.vcf > clinvar_pathogenic.vcf

# 组合复杂过滤
snpSift filter \
  "((ANN[*].IMPACT = 'HIGH') | \
    ((ANN[*].IMPACT = 'MODERATE') & (CADD_phred >= 20))) & \
   (!(exists gnomAD_genomes_AF) | (gnomAD_genomes_AF < 0.01))" \
  annotated_full.vcf > candidate_variants.vcf

# ===== 提取特定字段为表格 =====
snpSift extractFields \
  annotated_full.vcf \
  CHROM POS ID REF ALT QUAL FILTER \
  "ANN[0].GENE" "ANN[0].IMPACT" "ANN[0].EFFECT" \
  "ANN[0].HGVS_C" "ANN[0].HGVS_P" \
  "SIFT_pred" "Polyphen2_HDIV_pred" "CADD_phred" \
  "gnomAD_genomes_AF" "CLNSIG" \
  > variants_table.tsv

# ===== ANNOVAR安装和数据库下载 =====
# 从ANNOVAR官网注册下载
# https://annovar.openbioinformatics.org/

# 下载注释数据库
# 基因注释
perl annotate_variation.pl -buildver hg38 -downdb -webfrom annovar refGene humandb/
perl annotate_variation.pl -buildver hg38 -downdb -webfrom annovar ensGene humandb/

# 人群频率
perl annotate_variation.pl -buildver hg38 -downdb -webfrom annovar gnomad40_genome humandb/
perl annotate_variation.pl -buildver hg38 -downdb -webfrom annovar exac03 humandb/
perl annotate_variation.pl -buildver hg38 -downdb -webfrom annovar 1000g2015aug humandb/

# 致病性预测
perl annotate_variation.pl -buildver hg38 -downdb -webfrom annovar dbnsfp47a humandb/
perl annotate_variation.pl -buildver hg38 -downdb -webfrom annovar clinvar_20231231 humandb/
perl annotate_variation.pl -buildver hg38 -downdb -webfrom annovar avsnp150 humandb/
perl annotate_variation.pl -buildver hg38 -downdb -webfrom annovar cosmic97 humandb/

# ===== VCF转ANNOVAR格式 =====
perl convert2annovar.pl -format vcf4 input.vcf > input.avinput

# ===== ANNOVAR注释（一步完成多数据库注释） =====
perl table_annovar.pl \
  input.vcf \
  humandb/ \
  -buildver hg38 \
  -out annotated \
  -remove \
  -protocol refGene,ensGene,gnomad40_genome,exac03,avsnp150,dbnsfp47a,clinvar_20231231 \
  -operation g,g,f,f,f,f,f \
  -nastring . \
  -vcfinput \
  -polish \
  -thread 8

# 输出文件：
# annotated.hg38_multianno.txt   # 主要注释表格
# annotated.hg38_multianno.vcf   # 注释后的VCF

# ===== 查看结果 =====
head annotated.hg38_multianno.txt
# 包含列：Chr, Start, End, Ref, Alt, Func.refGene, Gene.refGene,
#         GeneDetail.refGene, ExonicFunc.refGene, AAChange.refGene,
#         gnomad40_genome_AF, SIFT_pred, Polyphen2_HDIV_pred,
#         CADD_phred, CLNSIG, ...

# ===== ANNOVAR过滤变异 =====
# 在R或Python中过滤（ANNOVAR输出为制表符分隔文本，易于处理）

# R中处理ANNOVAR输出
library(dplyr)

# 读取ANNOVAR注释结果
anno <- read.table("annotated.hg38_multianno.txt", header = TRUE, sep = "\t",
                    quote = "", fill = TRUE, comment.char = "")

# 过滤策略：
# 1. 功能性变异（外显子区域非同义突变）
exonic <- anno %>% filter(Func.refGene %in% c("exonic", "splicing"))
nonsynonymous <- exonic %>% filter(!ExonicFunc.refGene %in% c("synonymous SNV", "."))

# 2. 稀有变异（gnomAD MAF < 1%）
rare <- nonsynonymous %>%
  filter(is.na(gnomad40_genome_AF) | gnomad40_genome_AF < 0.01 | gnomad40_genome_AF == ".")

# 3. 致病性预测（至少一种工具预测有害）
damaging <- rare %>%
  filter(SIFT_pred == "D" | Polyphen2_HDIV_pred == "D" | 
         (!is.na(CADD_phred) & CADD_phred >= 20))

# 4. ClinVar已知致病
clinvar_pathogenic <- anno %>%
  filter(grepl("Pathogenic|Likely_pathogenic", CLNSIG))

# 5. 合并候选变异
candidates <- unique(rbind(damaging, clinvar_pathogenic))
cat("候选致病变异数:", nrow(candidates), "\n")

# 输出结果
write.csv(candidates, "candidate_pathogenic_variants.csv", row.names = FALSE)

# ===== VEP安装 =====
# Docker方式（推荐）
docker pull ensemblorg/ensembl-vep

# Conda方式
conda install -c bioconda ensembl-vep

# 下载缓存数据
vep_install --AUTO cf --ASSEMBLY GRCh38 --SPECIES homo_sapiens --PLUGINS all

# ===== 运行VEP =====
vep \
  --input_file input.vcf.gz \
  --output_file vep_annotated.vcf \
  --vcf \
  --cache \
  --dir_cache $HOME/.vep/ \
  --assembly GRCh38 \
  --species homo_sapiens \
  --force_overwrite \
  --fork 8 \
  --everything \
  --offline \
  --fasta GRCh38.fa

# --everything 等价于同时启用：
# --sift b --polyphen b --ccds --hgvs --symbol --numbers
# --domains --regulatory --canonical --protein --biotype
# --af --af_1kg --af_gnomadg --max_af --pubmed --clinical_signif

# ===== 使用插件增强注释 =====
vep \
  --input_file input.vcf.gz \
  --output_file vep_annotated.vcf \
  --vcf \
  --cache \
  --assembly GRCh38 \
  --fork 8 \
  --force_overwrite \
  --everything \
  --offline \
  --fasta GRCh38.fa \
  --plugin CADD,snv=whole_genome_SNVs.tsv.gz,indels=InDels.tsv.gz \
  --plugin REVEL,file=new_tabbed_revel_grch38.tsv.gz \
  --plugin SpliceAI,snv=spliceai_scores.raw.snv.hg38.vcf.gz,indel=spliceai_scores.raw.indel.hg38.vcf.gz \
  --plugin AlphaMissense,file=AlphaMissense_hg38.tsv.gz \
  --plugin LoFtool \
  --plugin MaxEntScan \
  --custom clinvar.vcf.gz,ClinVar,vcf,exact,0,CLNSIG,CLNDN

# ===== VEP过滤（vep_filter） =====
# 使用filter_vep工具
filter_vep \
  --input_file vep_annotated.vcf \
  --output_file filtered.vcf \
  --filter "IMPACT is HIGH or (IMPACT is MODERATE and CADD_PHRED >= 20)" \
  --filter "MAX_AF < 0.01 or not MAX_AF"

# ===== VEP输出为表格 =====
vep \
  --input_file input.vcf.gz \
  --output_file vep_table.txt \
  --tab \
  --fields "Uploaded_variation,Location,Allele,Gene,Feature,Consequence,\
cDNA_position,CDS_position,Protein_position,Amino_acids,Codons,\
IMPACT,SYMBOL,BIOTYPE,SIFT,PolyPhen,CADD_PHRED,gnomADg_AF,CLIN_SIG" \
  --cache --assembly GRCh38 --fork 8 --everything --offline

# ===== R中综合致病性评估 =====
library(dplyr)

# 读取注释结果
variants <- read.table("annotated_full.tsv", header = TRUE, sep = "\t")

# 定义致病性评估函数
assess_pathogenicity <- function(sift, polyphen, cadd, revel, gerp) {
  score <- 0
  if (!is.na(sift) && sift == "D") score <- score + 1
  if (!is.na(polyphen) && polyphen %in% c("D", "P")) score <- score + 1
  if (!is.na(cadd) && cadd >= 20) score <- score + 1
  if (!is.na(revel) && revel >= 0.5) score <- score + 1
  if (!is.na(gerp) && gerp >= 2) score <- score + 1

  if (score >= 4) return("Likely Pathogenic")
  if (score >= 3) return("Possibly Pathogenic")
  if (score >= 2) return("VUS")
  return("Likely Benign")
}

# 应用评估
variants$pathogenicity <- mapply(assess_pathogenicity,
                                  variants$SIFT_pred,
                                  variants$Polyphen2_HDIV_pred,
                                  variants$CADD_phred,
                                  variants$REVEL_score,
                                  variants$GERP_RS)

table(variants$pathogenicity)

# ===== 分步过滤pipeline =====

# 步骤1：SnpEff功能注释
snpEff ann -v -canon GRCh38.105 input.vcf.gz | bgzip > step1_annotated.vcf.gz

# 步骤2：保留功能性变异
snpSift filter \
  "(ANN[*].IMPACT = 'HIGH') | (ANN[*].IMPACT = 'MODERATE')" \
  step1_annotated.vcf.gz > step2_functional.vcf

# 步骤3：添加人群频率和致病性注释
snpSift dbnsfp \
  -v \
  -db dbNSFP4.9a.txt.gz \
  -f gnomAD_genomes_AF,SIFT_pred,Polyphen2_HDIV_pred,CADD_phred,REVEL_score \
  step2_functional.vcf > step3_annotated.vcf

# 步骤4：稀有变异过滤（MAF < 1%）
snpSift filter \
  "(!(exists gnomAD_genomes_AF)) | (gnomAD_genomes_AF < 0.01)" \
  step3_annotated.vcf > step4_rare.vcf

# 步骤5：致病性过滤
snpSift filter \
  "(CADD_phred >= 15) | (SIFT_pred = 'D') | (REVEL_score >= 0.5)" \
  step4_rare.vcf > step5_damaging.vcf

# 步骤6：ClinVar注释
snpSift annotate \
  -info CLNSIG,CLNDN \
  clinvar.vcf.gz \
  step5_damaging.vcf > step6_clinvar.vcf

# 步骤7：输出表格
snpSift extractFields \
  step6_clinvar.vcf \
  CHROM POS ID REF ALT QUAL \
  "ANN[0].GENE" "ANN[0].EFFECT" "ANN[0].IMPACT" \
  "ANN[0].HGVS_P" "ANN[0].HGVS_C" \
  gnomAD_genomes_AF SIFT_pred Polyphen2_HDIV_pred CADD_phred REVEL_score \
  CLNSIG CLNDN \
  > final_candidates.tsv

# 统计各步骤变异数
echo "Step 1 (all variants): $(grep -vc '^#' step1_annotated.vcf.gz)"
echo "Step 2 (functional): $(grep -vc '^#' step2_functional.vcf)"
echo "Step 4 (rare): $(grep -vc '^#' step4_rare.vcf)"
echo "Step 5 (damaging): $(grep -vc '^#' step5_damaging.vcf)"
echo "Final candidates: $(wc -l < final_candidates.tsv)"

# ===== 家系变异分析 =====
library(dplyr)

# 读取三人家系注释结果
variants <- read.table("trio_annotated.tsv", header = TRUE, sep = "\t")

# 假设列包含：proband_GT, father_GT, mother_GT
# GT格式：0/0=野生型纯合, 0/1=杂合, 1/1=变异纯合

# === 常染色体显性（de novo）===
de_novo <- variants %>%
  filter(proband_GT == "0/1" &   # 患者杂合
         father_GT == "0/0" &     # 父亲野生型
         mother_GT == "0/0" &     # 母亲野生型
         IMPACT %in% c("HIGH", "MODERATE"))
cat("De novo变异:", nrow(de_novo), "\n")

# === 常染色体隐性（纯合） ===
ar_homozygous <- variants %>%
  filter(proband_GT == "1/1" &   # 患者纯合
         father_GT == "0/1" &     # 父亲携带
         mother_GT == "0/1" &     # 母亲携带
         IMPACT %in% c("HIGH", "MODERATE"))
cat("AR纯合变异:", nrow(ar_homozygous), "\n")

# === 常染色体隐性（复合杂合） ===
# 在同一基因中找两个不同的杂合变异
ar_het <- variants %>%
  filter(proband_GT == "0/1" & IMPACT %in% c("HIGH", "MODERATE"))

# 找同一基因中有≥2个杂合变异的基因
compound_het_genes <- ar_het %>%
  group_by(GENE) %>%
  filter(n() >= 2) %>%
  ungroup()

# 验证：一个来自父亲，一个来自母亲（trans配置）
# 需要检查父母的基因型

# === X连锁隐性 ===
xlinked <- variants %>%
  filter(CHROM == "chrX" &
         proband_GT %in% c("0/1", "1/1") &  # 男性患者
         mother_GT == "0/1" &                  # 母亲携带
         IMPACT %in% c("HIGH", "MODERATE"))

# ===== ACMG分类标准 =====
# 根据ACMG指南对变异进行5级分类
# Pathogenic, Likely Pathogenic, VUS, Likely Benign, Benign
# 考虑因素：
# - PVS1: 无效变异（Loss of Function）在已知LOF致病基因中
# - PS1: 已知致病的相同氨基酸变化
# - PM2: 稀有变异（人群频率极低）
# - PP3: 多种计算工具预测有害
# - BP1: 错义变异在LOF致病基因中（弱良性证据）

#!/bin/bash
# ============================================
# 变异注释与过滤完整流程
# ============================================

INPUT="raw_variants.vcf.gz"
OUTPUT_DIR="annotation_results"
mkdir -p $OUTPUT_DIR

# 1. SnpEff功能注释
snpEff ann -v -canon -hgvs -lof \
  -stats $OUTPUT_DIR/snpeff_stats.html \
  GRCh38.105 $INPUT | \
  bgzip > $OUTPUT_DIR/snpeff.vcf.gz
tabix -p vcf $OUTPUT_DIR/snpeff.vcf.gz

# 2. dbNSFP致病性注释
snpSift dbnsfp -v \
  -db dbNSFP4.9a.txt.gz \
  -f SIFT_pred,Polyphen2_HDIV_pred,CADD_phred,REVEL_score,\
MutationTaster_pred,GERP_RS,gnomAD_genomes_AF \
  $OUTPUT_DIR/snpeff.vcf.gz > $OUTPUT_DIR/dbnsfp.vcf

# 3. ClinVar注释
snpSift annotate -v \
  -info CLNSIG,CLNDN,CLNREVSTAT \
  clinvar_20240101.vcf.gz \
  $OUTPUT_DIR/dbnsfp.vcf > $OUTPUT_DIR/clinvar.vcf

# 4. 过滤：功能性 + 稀有 + 致病
snpSift filter \
  "((ANN[*].IMPACT = 'HIGH') | \
    ((ANN[*].IMPACT = 'MODERATE') & \
     ((CADD_phred >= 20) | (REVEL_score >= 0.5)))) & \
   (!(exists gnomAD_genomes_AF) | (gnomAD_genomes_AF < 0.01))" \
  $OUTPUT_DIR/clinvar.vcf > $OUTPUT_DIR/filtered.vcf

# 5. 导出表格
snpSift extractFields $OUTPUT_DIR/filtered.vcf \
  CHROM POS ID REF ALT QUAL FILTER \
  "ANN[0].GENE" "ANN[0].EFFECT" "ANN[0].IMPACT" \
  "ANN[0].HGVS_C" "ANN[0].HGVS_P" "ANN[0].BIOTYPE" \
  gnomAD_genomes_AF SIFT_pred Polyphen2_HDIV_pred \
  CADD_phred REVEL_score GERP_RS \
  CLNSIG CLNDN \
  > $OUTPUT_DIR/final_table.tsv

# 统计
echo "Total variants: $(bcftools stats $INPUT | grep 'number of records' | awk '{print $6}')"
echo "After filtering: $(grep -vc '^#' $OUTPUT_DIR/filtered.vcf)"
echo "HIGH impact: $(grep -c 'HIGH' $OUTPUT_DIR/filtered.vcf)"

# 基因组水平
NC_000017.11:g.7675088C>T

# cDNA水平
NM_000546.6:c.743G>A

# 蛋白质水平
NP_000537.3:p.Arg248Gln (或 p.R248Q)

# 规则：
# 编码区SNV: c.位置碱基>碱基
# 缺失: c.位置del
# 插入: c.位置_位置insXXX
# 重复: c.位置_位置dup
# 移码: p.氨基酸位置fs*终止位置