757. 变异注释VEP/ANNOVAR对比¶

变异检测（GATK等）给你的是：
  chr1  12345  A  G  → 只知道位置和碱基变化

变异注释要回答：
  1. 在哪个基因上？→ TP53
  2. 在基因的什么区域？→ 外显子3
  3. 对蛋白有什么影响？→ p.Arg175His（错义突变）
  4. 功能预测？→ SIFT=Damaging, PolyPhen=Probably_damaging
  5. 临床意义？→ ClinVar: Pathogenic
  6. 人群频率？→ gnomAD: 0.0001（罕见）

类比：
  变异检测 = 找到嫌疑人
  变异注释 = 查嫌疑人的身份证、犯罪记录、家庭背景

# ===== 安装VEP =====
# 方法一：conda安装（推荐）
conda install -c bioconda ensembl-vep  # conda安装VEP

# 方法二：Docker
docker pull ensemblorg/ensembl-vep  # 拉取Docker镜像

# 方法三：源码安装
git clone https://github.com/Ensembl/ensembl-vep.git  # 克隆源码
cd ensembl-vep  # 进入目录
perl INSTALL.pl  # 运行安装脚本
# 安装时会提示下载缓存文件（homo_sapiens_vep_xxx_GRCh38.tar.gz）
# 强烈建议下载，否则每次运行都要联网查询

# 验证安装
vep --help  # 查看帮助信息

# ===== VEP基本使用 =====
vep \
  --input_file variants.vcf \           # 输入VCF文件
  --output_file annotated.vcf \         # 输出文件
  --vcf \                               # 输出VCF格式（也可选tab/json）
  --cache \                             # 使用本地缓存（不联网）
  --dir_cache ~/.vep \                  # 缓存目录
  --assembly GRCh38 \                   # 参考基因组版本
  --species homo_sapiens \              # 物种
  --force_overwrite \                   # 覆盖已有输出
  --offline \                           # 完全离线模式
  --fork 4                              # 使用4个线程

# ===== 常用功能选项 =====
vep \
  --input_file variants.vcf \
  --output_file annotated_full.vcf \
  --vcf \
  --cache \
  --offline \
  --assembly GRCh38 \
  --everything \                        # 启用所有注释（包含以下所有）
  # --sift b \                          # SIFT预测(b=both:预测+分数)
  # --polyphen b \                      # PolyPhen预测
  # --ccds \                            # CCDS转录本ID
  # --hgvs \                            # HGVS命名法
  # --symbol \                          # 基因symbol
  # --numbers \                         # 外显子/内含子编号
  # --domains \                         # 蛋白域
  # --regulatory \                      # 调控区域注释
  # --canonical \                       # 标记标准转录本
  # --protein \                         # 蛋白序列
  # --biotype \                         # 转录本类型
  # --af \                              # 全球等位基因频率
  # --af_gnomade \                      # gnomAD外显子频率
  # --max_af \                          # 最大等位基因频率
  # --pubmed \                          # PubMed文献
  # --variant_class \                   # 变异分类
  --fork 8                              # 8线程并行

# ===== VEP插件：CADD分数 =====
vep \
  --input_file variants.vcf \
  --output_file annotated_cadd.vcf \
  --vcf --cache --offline --assembly GRCh38 \
  --plugin CADD,/path/to/whole_genome_SNVs.tsv.gz \  # CADD分数插件
  --plugin LOFTEE \                     # LoF变异评估插件
  --plugin SpliceAI,snv=/path/to/spliceai_scores.raw.snv.hg38.vcf.gz \  # 剪接AI预测
  --fork 8

# ===== VEP过滤注释结果 =====
# 用filter_vep筛选高影响变异
filter_vep \
  --input_file annotated.vcf \          # 输入注释后的VCF
  --output_file high_impact.vcf \       # 输出筛选结果
  --filter "IMPACT is HIGH" \           # 只保留高影响变异
  --vcf_info_field ANN                  # VCF INFO字段名

# 复合筛选
filter_vep \
  --input_file annotated.vcf \
  --output_file filtered.vcf \
  --filter "IMPACT is HIGH or (SIFT is deleterious and PolyPhen is probably_damaging)" \
  --filter "MAX_AF < 0.01 or not MAX_AF"  # 罕见变异（MAF<1%）

# ===== 安装ANNOVAR =====
# ANNOVAR需要注册下载：https://annovar.openbioinformatics.org/
# 解压后即可使用
tar xvfz annovar.latest.tar.gz  # 解压
cd annovar  # 进入目录

# ===== 下载注释数据库 =====
# 基因注释（必需）
perl annotate_variation.pl \
  -buildver hg38 \                      # 基因组版本
  -downdb -webfrom annovar \            # 从ANNOVAR服务器下载
  refGene \                             # RefSeq基因注释
  humandb/                              # 保存目录

# 下载其他常用数据库
perl annotate_variation.pl -buildver hg38 -downdb -webfrom annovar ensGene humandb/      # Ensembl基因
perl annotate_variation.pl -buildver hg38 -downdb -webfrom annovar clinvar_20230416 humandb/  # ClinVar
perl annotate_variation.pl -buildver hg38 -downdb -webfrom annovar gnomad40_exome humandb/   # gnomAD
perl annotate_variation.pl -buildver hg38 -downdb -webfrom annovar dbnsfp42a humandb/       # dbNSFP(含SIFT/PolyPhen/CADD等)
perl annotate_variation.pl -buildver hg38 -downdb -webfrom annovar cosmic70 humandb/        # COSMIC
perl annotate_variation.pl -buildver hg38 -downdb -webfrom annovar avsnp150 humandb/        # dbSNP

# ===== 模式一：基因注释(gene-based) =====
# 判断变异在基因的哪个区域
perl annotate_variation.pl \
  -geneanno \                           # 基因注释模式
  -dbtype refGene \                     # 使用RefSeq基因模型
  -buildver hg38 \                      # 基因组版本
  input.avinput \                       # 输入文件（ANNOVAR格式）
  humandb/                              # 数据库目录

# ===== 模式二：区域注释(region-based) =====
# 变异是否落在特定区域（如保守区域、调控元件）
perl annotate_variation.pl \
  -regionanno \                         # 区域注释模式
  -dbtype phastConsElements46way \      # 46种物种保守性
  -buildver hg38 \
  input.avinput \
  humandb/

# ===== 模式三：过滤注释(filter-based) =====
# 用数据库信息过滤变异（频率、致病性等）
perl annotate_variation.pl \
  -filter \                             # 过滤注释模式
  -dbtype gnomad40_exome \              # gnomAD频率
  -buildver hg38 \
  input.avinput \
  humandb/

# ===== table_annovar：一次完成所有注释 =====
perl table_annovar.pl \
  input.vcf \                           # 输入VCF文件
  humandb/ \                            # 数据库目录
  -buildver hg38 \                      # 基因组版本
  -out annotated \                      # 输出前缀
  -remove \                             # 删除临时文件
  -protocol refGene,ensGene,clinvar_20230416,gnomad40_exome,dbnsfp42a,avsnp150 \  # 注释数据库列表
  -operation g,g,f,f,f,f \              # 每个数据库的注释模式(g=gene,f=filter)
  -nastring . \                         # 缺失值用"."表示
  -vcfinput \                           # 输入是VCF格式
  -polish \                             # 美化输出
  -thread 8                             # 8线程

# 输出文件：
# annotated.hg38_multianno.txt  → 制表符分隔的注释结果
# annotated.hg38_multianno.vcf  → VCF格式的注释结果

# ===== VCF转ANNOVAR输入格式 =====
# ANNOVAR需要特定格式：chr start end ref alt
perl convert2annovar.pl \
  -format vcf4 \                        # VCF4格式
  -allsample \                          # 处理所有样本
  -withfreq \                           # 包含频率信息
  input.vcf \                           # 输入VCF
  > input.avinput                       # 输出ANNOVAR格式

# 或者直接用table_annovar.pl的-vcfinput选项（推荐）

# ===== Python解析ANNOVAR结果 =====
import pandas as pd  # 导入pandas

# 读取ANNOVAR的multianno结果
df = pd.read_csv(
    "annotated.hg38_multianno.txt",  # ANNOVAR输出文件
    sep="\t",                          # 制表符分隔
    low_memory=False                   # 避免dtype警告
)

# 查看列名
print(df.columns.tolist())
# ['Chr', 'Start', 'End', 'Ref', 'Alt', 'Func.refGene', 'Gene.refGene',
#  'ExonicFunc.refGene', 'AAChange.refGene', 'CLNSIG', 'gnomAD_exome_AF',
#  'SIFT_pred', 'Polyphen2_HDIV_pred', 'CADD_phred', ...]

# 筛选致病变异
pathogenic = df[
    (df["CLNSIG"].str.contains("Pathogenic", na=False)) &  # ClinVar致病
    (df["Func.refGene"] == "exonic") &                      # 外显子区
    (df["ExonicFunc.refGene"] != "synonymous SNV")          # 非同义突变
]
print(f"致病外显子变异数: {len(pathogenic)}")

# 筛选罕见有害变异
rare_damaging = df[
    (df["gnomAD_exome_AF"].astype(float, errors="ignore") < 0.01) &  # MAF<1%
    (df["SIFT_pred"] == "D") &                                        # SIFT有害
    (df["Polyphen2_HDIV_pred"] == "D")                                # PolyPhen有害
]
print(f"罕见有害变异数: {len(rare_damaging)}")

# ===== 解析VEP的VCF输出 =====
import re  # 导入正则表达式

# VEP注释在VCF的INFO字段的CSQ中
with open("annotated_vep.vcf") as f:
    for line in f:
        if line.startswith("##INFO=<ID=CSQ"):  # 找到CSQ格式定义
            # 提取字段名
            fields = re.search(r'Format: (.+)"', line).group(1).split("|")
            print(f"VEP注释字段: {fields[:10]}...")  # 打印前10个字段
            break

# ===== 实战：同一批变异用两种工具注释 =====

# VEP注释
vep --input_file sample.vcf \
    --output_file sample_vep.vcf \
    --vcf --cache --offline --assembly GRCh38 \
    --everything --canonical --fork 8

# ANNOVAR注释
perl table_annovar.pl sample.vcf humandb/ \
    -buildver hg38 -out sample_annovar \
    -protocol refGene,clinvar_20230416,gnomad40_exome,dbnsfp42a \
    -operation g,f,f,f -vcfinput -polish

# 对比关键指标：
# 1. 注释了多少变异？两者应该基本一致
# 2. 外显子变异的功能分类是否一致？约85-90%一致
# 3. LoF变异是否一致？差异最大的类别（VEP有LOFTEE）
# 4. 临床变异分类？取决于ClinVar版本

场景选择指南：

临床/精准医学 → VEP
  理由：ClinVar/COSMIC内置，LOFTEE插件，REST API支持

大规模队列研究 → ANNOVAR
  理由：速度快，table_annovar一键完成，输出表格易处理

最佳实践 → 两者都用 + 交叉验证
  理由：2025年研究显示两工具在LoF变异上只有~65%一致
  特别是剪接变异差异最大

教学/入门 → VEP（在线版）
  理由：无需安装，网页界面，REST API秒出结果

# ===== 变异注释速查 =====

# VEP基本用法
vep --input_file in.vcf --output_file out.vcf \
    --vcf --cache --offline --assembly GRCh38 \
    --everything --canonical --fork 8

# VEP过滤
filter_vep -i out.vcf -o filtered.vcf \
    --filter "IMPACT is HIGH"

# ANNOVAR一键注释
perl table_annovar.pl in.vcf humandb/ \
    -buildver hg38 -out result \
    -protocol refGene,clinvar,gnomad,dbnsfp \
    -operation g,f,f,f -vcfinput

# 工具选择：
# 临床精准医学 → VEP（深度注释+插件）
# 大规模研究 → ANNOVAR（速度+表格输出）
# 最佳实践 → 两者交叉验证

# 关键阈值：
# CADD ≥ 15 → 可能有害
# CADD ≥ 25 → 很可能有害
# gnomAD AF < 0.01 → 罕见变异
# SIFT = Deleterious + PolyPhen = Probably_damaging → 有害