肿瘤基因组拷贝数变异 CNV¶

# === CNVkit：WES 肿瘤拷贝数分析（最常用） ===

# 1. 构建参考（用正常样本集合）
cnvkit.py batch \
  tumor.bam \                          # 肿瘤BAM
  --normal normal1.bam normal2.bam \   # 正常样本BAM（多个更好）
  --targets target.bed \               # 外显子靶区域
  --fasta hg38.fa \                    # 参考基因组
  --output-reference normal_ref.cnn \  # 输出参考文件
  --output-dir cnvkit_results/ \       # 输出目录
  --diagram --scatter                  # 生成图表

# 分步运行版本：
# 2. 计算覆盖度
cnvkit.py coverage tumor.bam target.bed -o tumor.targetcoverage.cnn   # 靶区域覆盖
cnvkit.py coverage tumor.bam antitarget.bed -o tumor.antitargetcoverage.cnn  # 非靶区域

# 3. 拷贝数分段（segmentation）
cnvkit.py segment tumor.cnr \          # 输入log2比值
  -o tumor.cns \                       # 输出分段结果
  --method cbs                         # CBS分段算法

# 4. 识别基因级别 CNV
cnvkit.py genemetrics tumor.cnr \      # 基因级别分析
  -s tumor.cns \                       # 分段结果
  -t 0.4 \                            # log2阈值（扩增）
  -m 5 \                              # 最少探针数
  > gene_cnv.txt                       # 输出

# 5. 绘制全基因组散点图
cnvkit.py scatter tumor.cnr \          # log2散点图
  -s tumor.cns \                       # 分段叠加
  -o scatter.pdf                       # 输出PDF

# 6. 绘制染色体图
cnvkit.py diagram tumor.cnr \          # 染色体图
  -s tumor.cns \                       # 分段叠加
  -o diagram.pdf                       # 输出PDF

# 7. 导出分段结果（整数拷贝数）
cnvkit.py call tumor.cns \             # 调用整数拷贝数
  --purity 0.7 \                       # 肿瘤纯度（从FACETS等获得）
  --ploidy 2 \                         # 肿瘤倍性
  -o tumor.call.cns                    # 输出调用结果

# 解析 CNVkit 结果，找临床重要基因的 CNV
import pandas as pd                    # 数据处理

# 读取基因级别 CNV 结果
gene_cnv = pd.read_csv("gene_cnv.txt", sep="\t")

# 定义临床重要基因列表
clinical_genes = {
    "amplification": ["ERBB2", "MYC", "MYCN", "EGFR", "CDK4",
                       "MDM2", "FGFR1", "MET", "CCND1"],
    "deletion": ["CDKN2A", "CDKN2B", "RB1", "PTEN", "TP53",
                  "BRCA1", "BRCA2", "ATM"]
}

print("=== 临床重要基因拷贝数变异 ===")

# 检查扩增（log2 > 0.58 约等于拷贝数 >= 3）
for gene in clinical_genes["amplification"]:
    row = gene_cnv[gene_cnv["gene"] == gene]  # 查找基因
    if not row.empty:                  # 如果找到了
        log2 = row.iloc[0]["log2"]     # 获取log2值
        cn = row.iloc[0].get("cn", 2 * 2**log2)  # 拷贝数
        if log2 > 0.58:               # 扩增阈值
            print(f"  ★ {gene} 扩增: log2={log2:.2f}, CN≈{cn:.1f}")

# 检查缺失（log2 < -1.0 约等于纯合缺失）
for gene in clinical_genes["deletion"]:
    row = gene_cnv[gene_cnv["gene"] == gene]
    if not row.empty:
        log2 = row.iloc[0]["log2"]
        if log2 < -1.0:               # 纯合缺失阈值
            print(f"  ★ {gene} 纯合缺失: log2={log2:.2f}")
        elif log2 < -0.4:             # 杂合缺失阈值
            print(f"  △ {gene} 杂合缺失: log2={log2:.2f}")