伴随诊断 CDx 生信分析¶

# === CDx 报告生成：从 VCF 到临床可操作变异 ===
import pandas as pd                    # 数据处理

# 1. 定义临床可操作变异库（简化版）
actionable_db = {
    # 基因: [(突变, 药物, 证据等级, 癌种)]
    "EGFR": [
        ("L858R", "Osimertinib", "Level 1", "NSCLC"),
        ("exon19del", "Osimertinib", "Level 1", "NSCLC"),
        ("T790M", "Osimertinib", "Level 1", "NSCLC"),
    ],
    "BRAF": [
        ("V600E", "Dabrafenib+Trametinib", "Level 1", "Melanoma"),
        ("V600E", "Encorafenib+Cetuximab", "Level 1", "CRC"),
    ],
    "ALK": [
        ("fusion", "Alectinib", "Level 1", "NSCLC"),
    ],
    "BRCA1": [
        ("pathogenic", "Olaparib", "Level 1", "Ovarian"),
    ],
    "BRCA2": [
        ("pathogenic", "Olaparib", "Level 1", "Ovarian"),
    ],
}

def generate_cdx_report(vcf_annotated_df, cancer_type, tmb, msi_status):
    """生成 CDx 风格的临床报告"""
    report = []                        # 报告条目

    # 检查每个变异是否可操作
    for _, var in vcf_annotated_df.iterrows():
        gene = var["gene"]             # 基因名
        protein = var["protein_change"] # 蛋白变化
        vaf = var["vaf"]               # VAF

        if gene in actionable_db:      # 如果在可操作库中
            for mutation, drug, level, indication in actionable_db[gene]:
                if mutation in protein or mutation == "pathogenic":
                    report.append({
                        "标志物": f"{gene} {protein}",
                        "药物": drug,
                        "证据等级": level,
                        "适应症": indication,
                        "VAF": f"{vaf:.1%}",
                        "状态": "检出"
                    })

    # 检查 TMB
    if tmb >= 10:                      # TMB-High
        report.append({
            "标志物": f"TMB-H ({tmb:.1f} mut/Mb)",
            "药物": "Pembrolizumab",
            "证据等级": "Level 1",
            "适应症": "泛癌种",
            "VAF": "N/A",
            "状态": "阳性"
        })

    # 检查 MSI
    if msi_status == "MSI-H":          # MSI-High
        report.append({
            "标志物": "MSI-H/dMMR",
            "药物": "Pembrolizumab",
            "证据等级": "Level 1",
            "适应症": "泛癌种",
            "VAF": "N/A",
            "状态": "阳性"
        })

    # 生成报告表
    report_df = pd.DataFrame(report)
    print("=" * 60)
    print("     伴随诊断检测报告 (CDx Report)")
    print("=" * 60)
    print(f"癌种: {cancer_type}")
    print(f"TMB: {tmb:.1f} mut/Mb")
    print(f"MSI: {msi_status}")
    print("-" * 60)
    if len(report_df) > 0:
        print(report_df.to_string(index=False))
    else:
        print("未检出临床可操作变异")
    print("=" * 60)
    return report_df

# 使用示例
# generate_cdx_report(mutations, "NSCLC", tmb=12.5, msi_status="MSS")

# === CDx Panel 数据分析质控 ===

# 1. 检查靶区域覆盖度（CDx要求每个靶位点深度达标）
mosdepth --by target_regions.bed \     # 按靶区域统计
  --thresholds 100,250,500 \           # 覆盖度阈值
  --no-per-base \                      # 不输出逐碱基
  sample_qc \                          # 输出前缀
  aligned.bam                          # 输入BAM

# 2. 检查每个 hotspot 位点的深度
samtools depth -a \                    # 包含零深度位点
  -b hotspot_positions.bed \           # hotspot位点
  aligned.bam | \                      # 输入BAM
  awk '$3 < 500 {print $0, "LOW_COVERAGE"}' \  # 标记低覆盖
  > hotspot_coverage.txt               # 输出

# 3. 统计 on-target 率
picard CollectHsMetrics \
  I=aligned.bam \                      # 输入BAM
  O=hs_metrics.txt \                   # 输出
  R=hg38.fa \                          # 参考基因组
  BAIT_INTERVALS=bait.interval_list \  # 探针区域
  TARGET_INTERVALS=target.interval_list  # 靶区域