344_遗传病致病变异筛选流程¶

# ========== 第一步：变异注释（ANNOVAR）==========

# 下载ANNOVAR数据库（需要注册）
# perl annotate_variation.pl -buildver hg38 -downdb -webfrom annovar \
#     refGene humandb/

# 运行ANNOVAR多库注释
perl table_annovar.pl \
    input.vcf \                      # 输入VCF文件
    humandb/ \                       # 数据库目录
    -buildver hg38 \                 # 参考基因组版本
    -out annotated_result \          # 输出前缀
    -remove \                        # 删除中间文件
    -protocol refGene,gnomad41_genome,clinvar_20240611,dbnsfp47a \  # 注释数据库
    -operation gx,f,f,f \           # 操作类型：基因/频率/频率/频率
    -nastring . \                    # 缺失值填充
    -vcfinput                        # 输入为VCF格式

# ========== 第二步：频率过滤 + 功能过滤 ==========
import pandas as pd   # 数据框处理

# 读取ANNOVAR注释结果
df = pd.read_csv("annotated_result.hg38_multianno.txt", sep='\t', low_memory=False)

print(f"原始变异数量: {len(df)}")  # 初始变异总数

# 步骤1：保留外显子区域变异
df = df[df['Func.refGene'].isin(['exonic', 'splicing'])]
print(f"外显子/剪接变异: {len(df)}")  # 过滤后数量

# 步骤2：过滤同义突变（不改变氨基酸的变异）
df = df[~df['ExonicFunc.refGene'].isin(['synonymous SNV'])]
print(f"非同义变异: {len(df)}")

# 步骤3：人群频率过滤（MAF < 0.001）
# 将频率列转为数值型，'.'转为NaN
df['AF_gnomad'] = pd.to_numeric(df['AF'], errors='coerce')
# 保留低频或未见变异（gnomAD）
df = df[(df['AF_gnomad'] < 0.001) | df['AF_gnomad'].isna()]
print(f"低频变异: {len(df)}")

# 步骤4：功能预测有害（CADD > 20）
df['CADD'] = pd.to_numeric(df['CADD_phred'], errors='coerce')
df_candidate = df[df['CADD'] >= 20]
print(f"CADD>=20候选变异: {len(df_candidate)}")  # 最终候选数

# ========== 第三步：ClinVar致病性注释 ==========
import re   # 正则表达式

def parse_clinvar(clnsig):
    """解析ClinVar临床意义字段"""
    if pd.isna(clnsig) or clnsig == '.':
        return 'Unknown'  # 无记录

    # 标记致病变异
    pathogenic_terms = ['Pathogenic', 'Likely_pathogenic', 
                        'Pathogenic/Likely_pathogenic']
    for term in pathogenic_terms:
        if term in str(clnsig):
            return 'Pathogenic'  # 致病或可能致病

    return 'Other'  # 其他分类

# 应用ClinVar解析
df_candidate['ClinVar_class'] = df_candidate['CLNSIG'].apply(parse_clinvar)

# 优先展示ClinVar已知致病变异
priority = df_candidate[df_candidate['ClinVar_class'] == 'Pathogenic']
print(f"ClinVar致病变异: {len(priority)}")

# ========== 第四步：Exomiser表型驱动排序 ==========

# 准备Exomiser配置文件
cat > exomiser_config.yml << 'EOF'
analysis:
    vcf: patient.vcf.gz           # 输入VCF
    proband: PATIENT_ID           # 先证者样本ID
    ped: family.ped               # 家系系谱文件（三人组用）
    hpoIds:                       # HPO表型编码（从patient症状得到）
        - HP:0001250  # 癫痫
        - HP:0000252  # 头围减小
        - HP:0001263  # 整体发育迟缓
    analysisMode: FULL            # 全面分析模式
    inheritanceModes:
        AUTOSOMAL_DOMINANT: 0.001   # 显性遗传频率阈值
        AUTOSOMAL_RECESSIVE: 0.001  # 隐性遗传频率阈值
outputOptions:
    outputFormats: [TSV_GENE, TSV_VARIANT, HTML]  # 输出格式
    numGenes: 20                  # 输出前20个候选基因
EOF

# 运行Exomiser
java -jar exomiser-cli-14.0.0.jar \
    --analysis exomiser_config.yml \  # 配置文件
    --output results/                  # 结果目录