变异检测

胚系变异检测与体系变异检测不同，由于胚系突变一般都是较高丰度的突变，因此难度实际上没有体细胞变异高。一些公司会选择使用多套胚系变异检测软件共同进行，尽量降低漏检概率。

在以往，大多数流程会使用的是gatk的HaplotypeCaller进行胚系变异检测，gatk的标准流程推荐使用BQSR，即先对bam文件进行质量校正，以提高变异检测的效果。

但在纯CPU分析里，HaplotypeCaller实在是太慢了，再叠加BQSR后，效率进一步降低（使用GPU流程当我没说）。现在及可见的未来里，一统胚系变异检测江湖的应该是DeepVariant，在我实测中，它的NA12878测试比HaplotypeCaller表现更好，在速度上更是占据绝对优势。

也可以看一看英伟达的Parabricks，有较高的硬件要求。

Call SNP/Indel

在DeepVariant里，暂不知道有没有同步覆盖侧翼位点的参数，因此，这里首先建立一个新的bed。这里是padding了上下游各50bp。

awk 'BEGIN {OFS="\t"} {start=$2; end=$3; new_start=start-50; new_end=end+50; if(new_start<0) new_start=0; print $1,new_start,new_end}' S31285117_Covered.bed > extended.bed

然后，使用侧翼放大后的bed，对bam文件进行变异检测。DeepVariant 1.10.0：

run_deepvariant \
    --model_type WES \
    --ref GRCh38.d1.vd1.fa \
    --reads bam/SRR14724513.markdup.bam \
    --output_vcf vcf/SRR14724513.vcf.gz \
    --output_gvcf vcf/SRR14724513.g.vcf.gz \
    --num_shards 16 \
    --regions extended.bed

这里生成了一个vcf文件和一个gvcf文件。gvcf文件不仅会储存变异位点信息，还会保留野生型的位点信息，这对于往后需要进行家系分析有极大帮助。如果要进行NA12878测试，可使用此时未经处理的vcf.gz。

单倍型计算

由于变异检测是局部进行的，可能会存在一些同一单倍型上的不连续的突变，被拆分为多个位点进行回报（基于HGVS规则，这是合理的），有时对于一些位点，需要进行合并回报，此时则需要确认位点是否属于同一单倍型。

为了提高效率，可以先按染色体拆分VCF，然后并行进行单倍型计算，最后再合并。

# 按染色体拆分VCF
for chrom in $(bcftools view -h vcf/SRR14724513.vcf.gz | grep "^##contig" | sed 's/.*ID=\([^,>]*\).*/\1/' | grep -E "^(chr)?([0-9]+|X|Y|M)$"); do
    bcftools view vcf/SRR14724513.vcf.gz --regions $chrom -O z -o vcf/split/${chrom}.split.vcf.gz
    bcftools index -t vcf/split/${chrom}.split.vcf.gz
done

# 对每个染色体并行进行单倍型计算
for chrom_vcf in vcf/split/*.split.vcf.gz; do
    prefix=$(basename $chrom_vcf .split.vcf.gz)
    whatshap phase \
        --reference=GRCh38.d1.vd1.fa \
        -o vcf/split/${prefix}.phase.vcf \
        $chrom_vcf \
        bam/SRR14724513.markdup.bam
    bgzip -f vcf/split/${prefix}.phase.vcf
    tabix -f -p vcf vcf/split/${prefix}.phase.vcf.gz
done

# 合并所有染色体的结果
bcftools concat -f vcf/split/phase_list.txt -a -O z -o vcf/SRR14724513.phase.vcf.gz
bcftools index -t vcf/SRR14724513.phase.vcf.gz

在下游中，可以进一步使用merge_mnp.py进行合并，但这往往不会达成最佳效果，因此建议是基于现实需求再人工进行合并。

左对齐和位点拆分

将检出的vcf位点进行拆分和左对齐，即将具有多个突变方向的位点，每个突变方向拆分到一行，同时对位点进行对齐。

gatk LeftAlignAndTrimVariants \
    -R GRCh38.d1.vd1.fa \
    -V vcf/SRR14724513.phase.vcf.gz \
    -O vcf/SRR14724513.left.vcf.gz \
    --split-multi-allelics

结果过滤

一般的，需要对结果进行过滤，主要是过滤掉深度、突变丰度不满足要求的位点。按目前的NA12878测试，DeepVariant的F1-score可达99%以上，因此即使位点没有通过我们的过滤阈值，大概率还是真的。因此这里比较保守只过滤不是PASS的位点。

bcftools view -f 'PASS,.' vcf/SRR14724513.left.vcf.gz -Oz -o vcf/SRR14724513.final.vcf.gz
bcftools index -t vcf/SRR14724513.final.vcf.gz

过滤后的VCF同样可以按染色体拆分，用于后续并行化的VEP注释。

注释

使用VEP (Variant Effect Predictor) 进行注释。VEP支持多种自定义数据库和插件，这里使用以下配置：

自定义数据库

ClinVar：临床变异数据库，收录变异的临床意义（致病性、良性等）
cytoBand：细胞遗传学条带，用于定位变异的染色体区域
GnomAD v4.1：大规模人群频率数据库，包含东亚人群频率
Pangolin：剪接位点预测，提供gain/loss分数
EVOScore2：进化保守性评分
AlphaMissense v3：基于AlphaFold的错义突变有害性预测

VEP插件

FlankingSequence：显示变异位点上下游10bp的侧翼序列
MissenseZscoreTranscript：转录本水平的错义突变Z-score

注释命令

# 先按染色体拆分VCF用于并行注释
for chrom in $(bcftools view -h vcf/SRR14724513.final.vcf.gz | grep "^##contig" | sed 's/.*ID=\([^,>]*\).*/\1/' | grep -E "^(chr)?([0-9]+|X|Y|M)$"); do
    bcftools view vcf/SRR14724513.final.vcf.gz --regions $chrom -O z -o vcf/split/${chrom}.split.vcf.gz
    bcftools index -t vcf/split/${chrom}.split.vcf.gz
done

# 对每个染色体并行运行VEP
for vcf_file in vcf/split/*.split.vcf.gz; do
    prefix=$(basename $vcf_file .split.vcf.gz)

    vep \
        --offline --cache \
        --dir_cache cache_dir --merged \
        --dir_plugins /opt/vep/.vep/Plugins \
        --force_overwrite --fork 4 \
        -i $vcf_file -o vcf/split/${prefix}.vep.vcf \
        --format vcf --vcf \
        --buffer_size 50000 \
        --fa GRCh38.d1.vd1.fa \
        --shift_3prime 1 --assembly GRCh38 --no_escape --check_existing \
        -exclude_predicted --uploaded_allele --show_ref_allele --numbers --domains \
        --total_length --hgvs --hgvsg --symbol --ccds --uniprot --max_af --pubmed \
        --transcript_filter "stable_id match N[MR]_" \
        --custom file=schema_bundle/hg38_clinvar_20260415.vcf.gz,short_name=ClinVar,format=vcf,type=exact,coords=0,fields=CLNSIG%CLNDN%CLNREVSTAT%CLNSTAR \
        --custom file=schema_bundle/hg38_cytoBand.bed.gz,short_name=cytoBand,format=bed,type=overlap,coords=0 \
        --custom file=schema_bundle/hg38_gnomad.v4.1.filtered.vcf.gz,short_name=GnomAD,format=vcf,type=exact,coords=0,fields=AC_joint%AN_joint%AF_joint%AF_joint_eas%nhomalt_joint_XX%nhomalt_joint_XY \
        --custom file=schema_bundle/hg38_pangolin.vcf.gz,short_name=Pangolin,format=vcf,type=exact,coords=0,fields=gain_score%loss_score \
        --custom file=schema_bundle/hg38_EVOScore2.vcf.gz,short_name=EVOScore2,format=vcf,type=exact,coords=0,fields=EVOScore \
        --custom file=schema_bundle/hg38_AlphaMissense.v3.vcf.gz,short_name=AlphaMissense,format=vcf,type=exact,coords=0,fields=AM%AMC \
        --plugin FlankingSequence,10 \
        --plugin MissenseZscoreTranscript,/opt/vep/.vep/Plugins/missenseByTranscript.hg38.v4.1.bed \
        --fields "Uploaded_variation,Location,REF_ALLELE,Allele,Consequence,IMPACT,DOMAINS,Feature,DISTANCE,EXON,INTRON,SYMBOL,STRAND,HGNC_ID,HGVSc,HGVSp,HGVSg,MAX_AF,Protein_position,Amino_acids,Codons,PUBMED,Existing_variation,cytoBand,ClinVar_CLNSIG,ClinVar_CLNREVSTAT,ClinVar_CLNDN,ClinVar_CLNSTAR,GnomAD_AC_joint,GnomAD_AN_joint,GnomAD_AF_joint,GnomAD_AF_joint_eas,GnomAD_nhomalt_joint_XX,GnomAD_nhomalt_joint_XY,Pangolin_gain_score,Pangolin_loss_score,EVOScore2_EVOScore,AlphaMissense_AM,AlphaMissense_AMC,FlankingSequence,MissenseZscore"
done

# 合并注释结果
bcftools concat -f vcf/split/vep_list.txt -a -O z -o vcf/SRR14724513.vep.vcf.gz
bcftools index -t vcf/SRR14724513.vep.vcf.gz

--transcript_filter "stable_id match N[MR]_" 确保只保留RefSeq转录本（NM_和NR_开头的ID）。

变异报告生成

从VEP注释后的VCF中提取变异信息，生成结构化的报告。支持首选转录本选择、GenCC疾病基因注释、性别感知的杂合性判断：

python3 vep_report.py \
    -i vcf/SRR14724513.vep.vcf.gz \
    -o report/SRR14724513.snv_indel.txt \
    --gencc assets/gencc-submissions.xlsx \
    -t assets/transcripts.json \
    --sex male \
    -n SRR14724513

--sex 参数会影响X和Y染色体上变异的杂合性判断：对于男性，X和Y染色体上的变异被判定为半合子（hemizygous），而非杂合。transcripts.json 提供了每个基因的首选转录本映射，确保报告的一致性。

NA12878

是的，可以进行NA12878的SNP/InDel性能测试了。

性别评估

基于全外显子的VCF进行性别评估，原理是男性的X染色体突变一般都是纯合，因此对若干男性、女性的样本进行统计X的杂合比例，就可获得一个确切的区分阈值，可以使用vcf作为输入文件；在发生X染色体异常时可能不准。