348_WES三人组分析Trio¶

# ========== 三人组GATK联合分析完整流程 ==========

# 假设已完成：比对→排序→标记重复→BQSR
# 样本命名：CHILD.recal.bam, FATHER.recal.bam, MOTHER.recal.bam

# 步骤1：各样本生成GVCF
for SAMPLE in CHILD FATHER MOTHER; do
    gatk HaplotypeCaller \
        -R hg38.fa \                         # 参考基因组
        -I ${SAMPLE}.recal.bam \             # 校正后BAM
        -O ${SAMPLE}.g.vcf.gz \              # 输出GVCF
        -ERC GVCF \                          # 生成GVCF格式
        --dbsnp dbsnp_146.hg38.vcf.gz \      # dbSNP数据库
        -L exome_targets.bed \               # 外显子靶向区域
        -G StandardAnnotation \              # 标准注释集
        --native-pair-hmm-threads 4          # 并行线程
    echo "完成 ${SAMPLE} 的GVCF生成"
done

# 步骤2：三人组联合变异分型
# 先合并GVCF到GenomicsDB（推荐用于大规模分析）
gatk GenomicsDBImport \
    -V CHILD.g.vcf.gz \         # 患者GVCF
    -V FATHER.g.vcf.gz \        # 父亲GVCF
    -V MOTHER.g.vcf.gz \        # 母亲GVCF
    --genomicsdb-workspace-path trio_db \  # 数据库目录
    -L exome_targets.bed \       # 目标区域
    --reader-threads 4           # 读取线程数

# 联合基因分型
gatk GenotypeGVCFs \
    -R hg38.fa \
    -V gendb://trio_db \         # 从数据库读取
    -O trio_raw.vcf.gz \         # 原始变异VCF
    -G StandardAnnotation \      # 标准注释
    --dbsnp dbsnp_146.hg38.vcf.gz

echo "联合基因分型完成"

# 步骤3：变异质量分层校正（VQSR）
# SNV的VQSR
gatk VariantRecalibrator \
    -R hg38.fa \
    -V trio_raw.vcf.gz \
    -O snv_recal \               # 校正模型输出
    --tranches-file snv.tranches \
    -mode SNP \                  # SNV模式
    --resource:hapmap,known=false,training=true,truth=true,prior=15.0 hapmap.vcf.gz \
    --resource:omni,known=false,training=true,truth=true,prior=12.0 omni.vcf.gz \
    --resource:1000G,known=false,training=true,truth=false,prior=10.0 1000G.vcf.gz \
    --resource:dbsnp,known=true,training=false,truth=false,prior=7.0 dbsnp.vcf.gz \
    -an QD -an MQ -an MQRankSum -an ReadPosRankSum -an FS -an SOR \
    --max-gaussians 4            # 高斯模型数量

gatk ApplyVQSR \
    -V trio_raw.vcf.gz \
    --recal-file snv_recal \
    --tranches-file snv.tranches \
    --truth-sensitivity-filter-level 99.5 \  # 99.5%敏感度阈值
    -mode SNP \
    -O trio_vqsr_snv.vcf.gz      # SNV过滤后VCF

# 步骤4：de novo变异检测
# 方法1：GATK PhaseByTransmission（传统方法）

# 方法2：DeepTrio（推荐，更高精度）
# DeepTrio是DeepVariant的三人组扩展版本
BIN_VERSION="1.6.1"

docker run \
    -v "${PWD}:/data" \                     # 挂载工作目录
    google/deepvariant:${BIN_VERSION} \
    /opt/deepvariant/bin/deeptrio/run_deeptrio \
    --model_type WES \                       # WES模式
    --ref /data/hg38.fa \                   # 参考基因组
    --reads_child /data/CHILD.recal.bam \   # 患者BAM
    --reads_parent1 /data/FATHER.recal.bam \ # 父亲BAM
    --reads_parent2 /data/MOTHER.recal.bam \ # 母亲BAM
    --sample_name_child CHILD \             # 样本名称
    --sample_name_parent1 FATHER \
    --sample_name_parent2 MOTHER \
    --output_vcf_child /data/child.vcf.gz \ # 各样本输出VCF
    --output_vcf_parent1 /data/father.vcf.gz \
    --output_vcf_parent2 /data/mother.vcf.gz \
    --num_shards 16 \                        # 并行分片数
    --regions /data/exome_targets.bed        # 目标区域

echo "DeepTrio分析完成"

# 步骤5：GEMINI三人组遗传模式查询
import subprocess   # 调用命令行工具
import pandas as pd # 数据处理

# 导入三人组变异到GEMINI数据库
subprocess.run([
    "gemini", "load",
    "-v", "trio_annotated.vcf.gz",  # 注释后VCF
    "-t", "snpEff",                  # 注释工具
    "-p", "family.ped",              # 家系PED文件
    "--cores", "8",                  # 并行核数
    "trio.db"                        # 输出数据库
])

# 查询de novo变异
result = subprocess.run(
    ["gemini", "de_novo",
     "--columns", "chrom,start,end,ref,alt,gene,impact,cadd_phred,aaf_gnomad_all",
     "--filter", "cadd_phred > 20 and impact_severity != 'LOW'",  # 高功能影响
     "trio.db"],
    capture_output=True, text=True
)

# 解析输出
lines = result.stdout.strip().split('\n')
if lines and lines[0]:
    import io
    df_denovo = pd.read_csv(io.StringIO(result.stdout), sep='\t')
    print(f"de novo候选变异: {len(df_denovo)}")
    print(df_denovo[['gene', 'ref', 'alt', 'impact', 'cadd_phred']].to_string())
else:
    print("未发现de novo候选变异")