罕见病致病变异筛选¶

WES/WGS测序 → 变异检测(GATK) → 频率过滤 → 功能注释 → 遗传模式筛选 → ACMG分类 → 候选变异

# 假设已完成 GATK Best Practices 变异检测
# 得到过滤后的 VCF 文件

# 用 ANNOVAR 做全面注释
perl table_annovar.pl \
    proband.vcf \                     # 先证者VCF
    humandb/ \                         # 数据库目录
    -buildver hg38 \                  # 基因组版本
    -out proband_annotated \          # 输出前缀
    -remove \                          # 删除临时文件
    -protocol refGene,gnomad40_genome,gnomad40_exome,clinvar_20240917,dbnsfp42a,dbscsnv11,intervar_20180118 \
    -operation g,f,f,f,f,f,f \        # 操作类型
    -nastring . \                      # 缺失值标记
    -vcfinput                          # 输入是VCF格式

# 或用 VEP (Variant Effect Predictor)
vep --input_file proband.vcf \
    --output_file proband_vep.vcf \
    --cache \                          # 使用本地缓存
    --dir_cache ~/.vep \              # 缓存目录
    --assembly GRCh38 \               # 基因组版本
    --everything \                     # 输出所有注释
    --vcf \                            # 输出VCF格式
    --force_overwrite                  # 覆盖已有文件

# ========== Python脚本：变异过滤 ==========
import pandas as pd

# 读取 ANNOVAR 注释结果
variants = pd.read_csv("proband_annotated.hg38_multianno.txt",
                        sep="\t")                       # 读取注释结果

print(f"原始变异数: {len(variants)}")

# 第 1 步：频率过滤（核心！）
# 罕见病的致病变异在健康人群中应该极罕见
# 常染色体显性：gnomAD AF < 0.0001 (万分之一)
# 常染色体隐性：gnomAD AF < 0.01 (百分之一)

# 过滤 gnomAD 频率
for col in ['gnomAD_genome_ALL', 'gnomAD_exome_ALL']:  # 检查两个数据库
    variants[col] = pd.to_numeric(variants[col],
                                   errors='coerce')     # 转为数值
    variants[col] = variants[col].fillna(0)              # 缺失值设为0

# 保留罕见变异（AF < 0.01）
rare = variants[
    (variants['gnomAD_genome_ALL'] < 0.01) &             # 基因组频率<1%
    (variants['gnomAD_exome_ALL'] < 0.01)                # 外显子组频率<1%
]
print(f"频率过滤后: {len(rare)}")

# 第 2 步：功能过滤
# 只保留可能影响蛋白功能的变异
functional = rare[
    rare['ExonicFunc.refGene'].isin([
        'nonsynonymous SNV',              # 错义突变
        'stopgain',                        # 获得终止密码子
        'stoploss',                        # 丢失终止密码子
        'frameshift deletion',             # 移码缺失
        'frameshift insertion',            # 移码插入
        'nonframeshift deletion',          # 非移码缺失
        'nonframeshift insertion',         # 非移码插入
        'splicing',                        # 剪接位点变异
    ]) |
    rare['Func.refGene'].str.contains('splicing', na=False)  # 剪接区域
]
print(f"功能过滤后: {len(functional)}")

# 第 3 步：质量过滤
# 深度≥10, 基因型质量≥20
quality = functional[
    (functional['QUAL'] > 30)                             # 变异质量>30
]
print(f"质量过滤后: {len(quality)}")

# 保存候选变异
quality.to_csv("candidate_variants.tsv", sep="\t", index=False)

# Trio 分析是罕见病诊断的"杀手锏"
# 先证者 + 父亲 + 母亲 一起测序分析

# GATK Joint Genotyping
gatk GenotypeGVCFs \
    -R reference.fa \
    -V proband.g.vcf.gz \              # 先证者
    -V father.g.vcf.gz \               # 父亲
    -V mother.g.vcf.gz \               # 母亲
    -O trio.vcf.gz                      # 联合基因型

# 检测 de novo 变异（新发突变）
# de novo = 先证者有，父母都没有
bcftools view \
    -i 'GT[0]="0/1" & GT[1]="0/0" & GT[2]="0/0"' \  # 先证者杂合，父母纯合参考
    trio.vcf.gz > denovo_candidates.vcf

# 筛选复合杂合（Compound Heterozygous）
# 同一基因上，一个变异来自父亲，一个来自母亲
# 使用 gemini 工具
gemini load -v trio.vcf.gz -t VEP trio.db  # 加载到gemini数据库
gemini comp_hets \
    --columns "chrom,start,gene,impact" \
    --filter "aaf < 0.01" \              # 频率<1%
    trio.db > compound_hets.txt           # 复合杂合结果

# Exomiser 结合表型信息（HPO）和基因组数据做智能排序
# 安装
wget https://github.com/exomiser/Exomiser/releases/latest/download/exomiser-cli.zip
unzip exomiser-cli.zip

# 准备配置文件
cat << 'EOF' > analysis.yml
analysis:
  vcf: proband.vcf                    # 输入VCF
  ped: family.ped                     # 家系PED文件
  proband: PROBAND_ID                 # 先证者ID
  hpoIds: ['HP:0001250', 'HP:0001249', 'HP:0002069']  # HPO表型术语
  # HP:0001250 = 癫痫
  # HP:0001249 = 智力障碍
  # HP:0002069 = 全面性肌阵挛

  analysisMode: PASS_ONLY             # 只分析PASS变异
  inheritanceModes:                   # 考虑的遗传模式
    AUTOSOMAL_DOMINANT: 0.1           # 常染色体显性
    AUTOSOMAL_RECESSIVE_HOM_ALT: 0.1 # 常染色体隐性纯合
    AUTOSOMAL_RECESSIVE_COMP_HET: 2.0 # 复合杂合
    X_DOMINANT: 0.1                   # X连锁显性
    X_RECESSIVE: 0.1                  # X连锁隐性

  frequencySources:
    - GNOMAD_E_AFR                    # gnomAD外显子组-非洲
    - GNOMAD_E_EAS                    # gnomAD外显子组-东亚
    - GNOMAD_E_SAS                    # gnomAD外显子组-南亚

  pathogenicitySources:
    - POLYPHEN                        # PolyPhen-2预测
    - MUTATION_TASTER                 # MutationTaster预测
    - SIFT                            # SIFT预测
    - CADD                            # CADD评分
    - REVEL                           # REVEL评分

  steps:
    - variantEffectFilter:
        remove: [SYNONYMOUS_VARIANT]  # 移除同义变异
    - frequencyFilter:
        maxFrequency: 1.0             # 最大频率1%
    - pathogenicityFilter:
        keepNonPathogenic: false      # 移除预测为良性的
    - inheritanceFilter: {}           # 遗传模式过滤
    - omimPrioritiser: {}             # OMIM优先排序
    - hiPhivePrioritiser: {}          # 表型相似度排序
EOF

# 运行 Exomiser
java -jar exomiser-cli.jar \
    --analysis analysis.yml            # 运行分析

# 输出：HTML报告，按优先级排序的候选基因和变异

# === 罕见病变异筛选速查 ===

# 1. 变异注释
perl table_annovar.pl input.vcf humandb/ -buildver hg38 \
  -protocol refGene,gnomad40_genome,clinvar,dbnsfp42a -operation g,f,f,f

# 2. 频率过滤阈值：
# 显性遗传: gnomAD AF < 0.0001 (万分之一)
# 隐性遗传: gnomAD AF < 0.01 (百分之一)

# 3. Trio de novo 检测
bcftools view -i 'GT[0]="0/1" & GT[1]="0/0" & GT[2]="0/0"' trio.vcf

# 4. Exomiser（表型驱动排序）
java -jar exomiser-cli.jar --analysis analysis.yml

# 5. HPO 表型术语查询
# https://hpo.jax.org/

# 罕见病诊断率参考：
# Trio WES: ~40%
# Singleton WES: ~25%
# WGS: ~45%
# 重分析(2-3年后): 额外 +10%