肿瘤突变负荷 TMB 计算¶

# === TMB 计算完整流程 ===

# 1. 先完成体细胞突变检测（接 319 的流程）
# 假设已有 filtered_somatic.vcf.gz

# 2. 提取 PASS 的非同义突变
bcftools view -f "PASS" \              # 只保留PASS变异
  filtered_somatic.vcf.gz | \          # 输入VCF
  bcftools query \                     # 查询字段
  -f '%CHROM\t%POS\t%REF\t%ALT\t%INFO/CSQ\n' | \  # 提取注释信息
  grep -E "missense|nonsense|frameshift|splice" \    # 保留非同义变异
  > nonsynonymous_variants.txt         # 输出非同义变异列表

# 3. 计算靶区域大小（BED文件）
# WES 外显子靶区域
awk '{sum += $3 - $2} END {print sum / 1e6, "Mb"}' \
  target_regions.bed                   # 计算BED覆盖总长度(Mb)

# 4. 计算 TMB
VARIANT_COUNT=$(wc -l < nonsynonymous_variants.txt)  # 变异数
TARGET_SIZE=$(awk '{sum+=$3-$2} END {print sum/1e6}' target_regions.bed)  # 区域大小Mb
TMB=$(echo "scale=2; $VARIANT_COUNT / $TARGET_SIZE" | bc)  # TMB计算
echo "TMB = $TMB mutations/Mb"         # 打印TMB值

# 用 Python 进行更精细的 TMB 计算
import cyvcf2                          # 高效VCF解析库
import pandas as pd                    # 数据处理

def calculate_tmb(vcf_path, bed_path, min_vaf=0.05, min_depth=20):
    """
    计算肿瘤突变负荷
    vcf_path: 过滤后的体细胞突变VCF
    bed_path: 靶区域BED文件
    min_vaf: 最小VAF阈值
    min_depth: 最小深度阈值
    """
    # 计算靶区域大小
    bed = pd.read_csv(bed_path, sep="\t", header=None,  # 读取BED
                      names=["chrom", "start", "end"])
    target_size_mb = (bed["end"] - bed["start"]).sum() / 1e6  # 总大小Mb

    # 遍历VCF统计非同义突变
    vcf = cyvcf2.VCF(vcf_path)         # 打开VCF
    mutation_count = 0                  # 突变计数器
    excluded_drivers = set()            # 排除的已知驱动突变

    for variant in vcf:                 # 遍历每个变异
        if variant.FILTER is not None:  # 跳过非PASS变异
            continue

        # 检查深度和VAF
        tumor_dp = variant.format("DP")[0][0]  # 肿瘤样本深度
        tumor_ad = variant.format("AD")[0]     # 等位基因深度
        if tumor_dp < min_depth:        # 深度不够跳过
            continue
        vaf = tumor_ad[1] / tumor_dp if tumor_dp > 0 else 0  # 计算VAF
        if vaf < min_vaf:               # VAF太低跳过
            continue

        # 检查是否为非同义突变（通过VEP注释）
        info = variant.INFO              # 获取INFO字段
        csq = info.get("CSQ", "")       # 获取VEP注释
        if any(term in csq for term in  # 检查功能影响
               ["missense", "nonsense", "frameshift",
                "stop_gained", "splice_donor", "splice_acceptor"]):
            mutation_count += 1          # 计数+1

    # 计算TMB
    tmb = mutation_count / target_size_mb  # 突变数 / 区域大小
    return {
        "mutation_count": mutation_count,   # 突变总数
        "target_size_mb": round(target_size_mb, 2),  # 靶区域大小
        "tmb": round(tmb, 2),               # TMB值
        "tmb_high": tmb >= 10               # 是否TMB-High
    }

# 使用示例
result = calculate_tmb("final_somatic.vcf.gz", "target.bed")
print(f"突变数: {result['mutation_count']}")
print(f"靶区域: {result['target_size_mb']} Mb")
print(f"TMB: {result['tmb']} mut/Mb")
print(f"TMB-High: {'是' if result['tmb_high'] else '否'}")