遗传变异致病性预测¶

# 方法1: 使用CADD在线服务
# https://cadd.gs.washington.edu/score
# 上传VCF文件获取评分

# 方法2: 本地安装CADD
# 下载预计算分数（全基因组SNV约80GB）
wget https://krishna.gs.washington.edu/download/CADD/v1.7/GRCh38/whole_genome_SNVs.tsv.gz
wget https://krishna.gs.washington.edu/download/CADD/v1.7/GRCh38/whole_genome_SNVs.tsv.gz.tbi

# 查询特定变异的CADD分数
tabix whole_genome_SNVs.tsv.gz 1:942451-942451
# 输出: chr, pos, ref, alt, RawScore, PHRED

# 方法3: 使用CADD脚本评分自定义变异
git clone https://github.com/kircherlab/CADD-scripts.git
cd CADD-scripts

# 安装依赖和下载数据
./install.sh

# 对VCF文件评分
./CADD.sh -g GRCh38 -v v1.7 -o scored_variants.tsv.gz input.vcf

# 解析CADD分数
zcat scored_variants.tsv.gz | awk '$6 >= 20' > cadd_deleterious.tsv
# PHRED >= 20: top 1% most deleterious
# PHRED >= 30: top 0.1%

# 方法4: 使用Python/R注释
# Python
import pandas as pd
from pysam import TabixFile

cadd = TabixFile("whole_genome_SNVs.tsv.gz")
records = cadd.fetch("1", 942450, 942451)
for r in records:
    fields = r.split("\t")
    print(f"Chr:{fields[0]} Pos:{fields[1]} Ref:{fields[2]} Alt:{fields[3]} PHRED:{fields[5]}")

# 下载REVEL预计算分数
wget https://sites.google.com/site/revelgenomics/downloads/revel-v1.3_all_chromosomes.zip
unzip revel-v1.3_all_chromosomes.zip

# 格式: chr, hg19_pos, grch38_pos, ref, alt, aaref, aaalt, REVEL
# 查询特定变异
zcat revel_all_chr.csv.gz | awk -F',' '$1=="1" && $3==942451'

# 批量注释VCF文件（使用VEP + REVEL插件）
vep --input_file variants.vcf \
    --output_file annotated.vcf \
    --format vcf \
    --vcf \
    --assembly GRCh38 \
    --cache \
    --plugin REVEL,revel_scores.tsv.gz \
    --plugin CADD,whole_genome_SNVs.tsv.gz \
    --plugin SpliceAI,snv=spliceai_scores.raw.snv.hg38.vcf.gz,indel=spliceai_scores.raw.indel.hg38.vcf.gz \
    --everything \
    --fork 8

# 解析VEP输出，提取高分变异
grep -v "^#" annotated.vcf | \
    awk '{if ($0 ~ /REVEL=/) {match($0, /REVEL=([0-9.]+)/, a); if (a[1] > 0.5) print}}' \
    > high_revel_variants.vcf

# R中使用REVEL注释
library(tidyverse)

# 读取REVEL数据
revel <- read_csv("revel_all_chr.csv.gz", 
    col_names=c("chr","hg19_pos","grch38_pos","ref","alt",
                "aaref","aaalt","REVEL"))

# 读取候选变异
variants <- read_tsv("candidate_variants.tsv")

# 合并注释
annotated <- variants %>%
    left_join(revel, by=c("chr","pos"="grch38_pos","ref","alt"))

# 筛选高致病性变异
# ClinGen校准阈值（Pejaver et al., 2022, AJHG）:
#   PP3_Supporting: >= 0.644
#   PP3_Moderate:   >= 0.773
#   PP3_Strong:     >= 0.932
#   BP4_Moderate:   <= 0.183
pathogenic <- annotated %>%
    filter(REVEL >= 0.644) %>%  # ClinGen PP3_Supporting阈值
    arrange(desc(REVEL))

cat("REVEL >= 0.644 (PP3 Supporting) 的变异数:", nrow(pathogenic), "\n")

# 下载AlphaMissense预计算分数
# https://zenodo.org/record/8208688
# 文件: AlphaMissense_hg38.tsv.gz (约4.5GB)

import pandas as pd

# 读取AlphaMissense分数
am_scores = pd.read_csv("AlphaMissense_hg38.tsv.gz", 
    sep="\t", comment="#",
    names=["CHROM","POS","REF","ALT","genome","uniprot_id",
           "transcript_id","protein_variant","am_pathogenicity","am_class"])

# am_class: 'likely_pathogenic', 'likely_benign', 'ambiguous'
# am_pathogenicity: 0-1 连续分数
#
# 重要阈值说明：
# - 0.564 = 原始论文默认阈值（90%精度），但ClinGen校准表明此阈值
#   未达到PP3 Supporting的证据强度，不建议直接用于临床分类
# - ClinGen校准阈值（Bergquist et al., Genetics in Medicine, 2025）:
#   PP3_Supporting: >= 0.719
#   PP3_Moderate:   >= 0.898
#   PP3_Strong:     >= 0.978
#   BP4_Moderate:   <= 0.070（较少变异达到此区间）
# - 临床变异分类应使用ClinGen校准阈值，而非原始论文阈值

# 查询特定变异
query = am_scores[(am_scores["CHROM"]=="chr1") & 
                   (am_scores["POS"]==942451)]
print(query)

# 统计全基因组分类分布
print(am_scores["am_class"].value_counts())
# likely_benign: ~50M
# likely_pathogenic: ~18M  
# ambiguous: ~22M

# 筛选致病候选
pathogenic_am = am_scores[am_scores["am_class"] == "likely_pathogenic"]

# 与VCF文件合并
import cyvcf2

vcf = cyvcf2.VCF("patient_variants.vcf.gz")
results = []
for variant in vcf:
    chrom = f"chr{variant.CHROM}"
    pos = variant.POS
    ref = variant.REF
    alt = variant.ALT[0]

    match = am_scores[(am_scores["CHROM"]==chrom) & 
                       (am_scores["POS"]==pos) & 
                       (am_scores["REF"]==ref) & 
                       (am_scores["ALT"]==alt)]
    if not match.empty:
        results.append({
            "variant": f"{chrom}:{pos}{ref}>{alt}",
            "am_score": match.iloc[0]["am_pathogenicity"],
            "am_class": match.iloc[0]["am_class"],
            "protein_change": match.iloc[0]["protein_variant"]
        })

results_df = pd.DataFrame(results)
print(results_df.sort_values("am_score", ascending=False).head(20))

# 安装SpliceAI
pip install spliceai

# 使用预计算分数（推荐，大约30GB）
# 下载地址: https://basespace.illumina.com/s/otSPW8hnhaZR
# spliceai_scores.raw.snv.hg38.vcf.gz
# spliceai_scores.raw.indel.hg38.vcf.gz

# 查询预计算分数
tabix spliceai_scores.raw.snv.hg38.vcf.gz chr1:942451-942451

# 或实时计算（适合少量变异）
spliceai -I input.vcf \
    -O output.vcf \
    -R /ref/hg38.fa \
    -A grch38 \
    -D 500    # 考虑距离剪接位点500bp范围

# 输出VCF中的SpliceAI注释:
# SpliceAI=ALLELE|SYMBOL|DS_AG|DS_AL|DS_DG|DS_DL|DP_AG|DP_AL|DP_DG|DP_DL
# DS: Delta Score (0-1)  DP: Delta Position (距离)
# AG: Acceptor Gain, AL: Acceptor Loss
# DG: Donor Gain, DL: Donor Loss

# 解读:
# DS >= 0.2: 可能影响剪接
# DS >= 0.5: 很可能影响剪接
# DS >= 0.8: 高置信度影响剪接

# Python中使用SpliceAI
from keras.models import load_model
import numpy as np

# 或更实用地解析预计算结果
import cyvcf2

vcf = cyvcf2.VCF("output_spliceai.vcf.gz")
splice_variants = []

for variant in vcf:
    info = variant.INFO.get("SpliceAI")
    if info:
        fields = info.split("|")
        allele, gene = fields[0], fields[1]
        ds_ag, ds_al, ds_dg, ds_dl = [float(x) for x in fields[2:6]]
        max_ds = max(ds_ag, ds_al, ds_dg, ds_dl)

        if max_ds >= 0.2:
            splice_variants.append({
                "chrom": variant.CHROM,
                "pos": variant.POS,
                "ref": variant.REF,
                "alt": variant.ALT[0],
                "gene": gene,
                "max_delta_score": max_ds,
                "acceptor_gain": ds_ag,
                "acceptor_loss": ds_al,
                "donor_gain": ds_dg,
                "donor_loss": ds_dl
            })

splice_df = pd.DataFrame(splice_variants)
print(f"影响剪接的变异: {len(splice_df)}")
print(splice_df.sort_values("max_delta_score", ascending=False).head(10))

# 综合变异注释与分级
import pandas as pd

def classify_variant(row):
    """简化的变异分级逻辑"""
    evidence = {"pathogenic": [], "benign": []}

    # 1. 群体频率
    if row.get("gnomAD_AF", 0) > 0.05:
        evidence["benign"].append("BA1")  # 良性, standalone
    elif row.get("gnomAD_AF", 0) < 0.0001:
        evidence["pathogenic"].append("PM2")  # 罕见

    # 2. 计算预测（ClinGen校准阈值，Pejaver et al., 2022）
    # REVEL PP3: >=0.644(Supporting), >=0.773(Moderate), >=0.932(Strong)
    # REVEL BP4: <=0.183(Moderate), <=0.016(Strong)
    if row.get("REVEL", 0) >= 0.932:
        evidence["pathogenic"].append("PP3_Strong")
    elif row.get("REVEL", 0) >= 0.773:
        evidence["pathogenic"].append("PP3_Moderate")
    elif row.get("REVEL", 0) >= 0.644:
        evidence["pathogenic"].append("PP3")  # Supporting
    elif row.get("REVEL", 0) <= 0.183:
        evidence["benign"].append("BP4")  # Moderate

    # 3. 蛋白截断变异
    if row.get("consequence") in ["stop_gained", "frameshift"]:
        if row.get("pLI", 0) > 0.9:  # 基因对LoF不耐受
            evidence["pathogenic"].append("PVS1")

    # 4. ClinVar已知
    if row.get("ClinVar_class") == "Pathogenic":
        evidence["pathogenic"].append("PS1_like")

    # 5. SpliceAI
    if row.get("SpliceAI_max", 0) > 0.5:
        evidence["pathogenic"].append("PP3_splice")

    # 简化分类（实际应遵循ACMG完整规则）
    path_count = len(evidence["pathogenic"])
    benign_count = len(evidence["benign"])

    if "BA1" in evidence["benign"]:
        return "Benign"
    elif "PVS1" in evidence["pathogenic"] and path_count >= 2:
        return "Pathogenic"
    elif path_count >= 3:
        return "Likely_Pathogenic"
    elif benign_count >= 2:
        return "Likely_Benign"
    else:
        return "VUS"

# 应用分级
variants_df["ACMG_class"] = variants_df.apply(classify_variant, axis=1)
print(variants_df["ACMG_class"].value_counts())

# 使用InterVar自动ACMG分类
# pip install intervar
# 或命令行:
python InterVar.py \
    -b hg38 \
    -i input.vcf \
    --input_type=VCF \
    -o intervar_output \
    -d /db/intervar/

# 使用VEP + 多插件一站式注释
vep -i input.vcf -o annotated.vcf --vcf --assembly GRCh38 --cache \
    --plugin CADD,cadd_snvs.tsv.gz \
    --plugin REVEL,revel.tsv.gz \
    --plugin SpliceAI,snv=spliceai_snv.vcf.gz,indel=spliceai_indel.vcf.gz \
    --plugin AlphaMissense,file=AlphaMissense.tsv.gz \
    --plugin gnomADc,gnomad.genomes.vcf.gz \
    --everything --fork 8

#!/bin/bash
# === 变异致病性评估完整流程 ===

VCF="patient_wes.vcf.gz"
REF="/ref/hg38.fa"
THREADS=8

# 1. VEP注释（一站式）
vep -i $VCF -o vep_annotated.vcf --vcf --assembly GRCh38 --cache \
    --plugin CADD,cadd.tsv.gz \
    --plugin REVEL,revel.tsv.gz \
    --plugin SpliceAI,snv=spliceai_snv.vcf.gz,indel=spliceai_indel.vcf.gz \
    --sift b --polyphen b --af_gnomade --everything --fork $THREADS

# 2. 过滤候选致病变异
# 保留: 罕见(AF<0.01) + 蛋白影响 + 高CADD/REVEL
bcftools query -f '%CHROM\t%POS\t%REF\t%ALT\t%INFO/CSQ\n' vep_annotated.vcf | \
    python3 filter_pathogenic.py > candidate_pathogenic.tsv

# 3. InterVar ACMG分级
python InterVar.py -b hg38 -i $VCF --input_type=VCF -o acmg_classified

错义变异: REVEL(>=0.644) 或 AlphaMissense(>=0.719) + CADD(>25)
         注意：ClinGen要求每个实验室/基因预先选定一个工具，不能混用取最高分
截断变异: pLI(>0.9)/LOEUF(<0.35) + NMD预测 + 最后外显子规则
剪接变异: SpliceAI(>0.5) + MaxEntScan + dbscSNV
非编码: CADD(>15) + DANN + RegulomeDB + ENCODE注释
同义变异: SpliceAI + 密码子使用频率变化 + ESE/ESS影响