628_拷贝数变异分析CNV¶

一句话概述： 拷贝数变异（CNV）是基因组某段DNA的拷贝数增加或减少，CNVkit是WES/Panel数据的首选工具，GATK CNV Pipeline适合大规模标准化分析，Sequenza专门用于肿瘤纯度和倍性估计——三者互补覆盖了从种系到体细胞CNV的全部场景。

核心知识点速查表¶

概念	白话解释
Copy Number Variation (CNV)	拷贝数变异——某段DNA多了或少了几个拷贝
Gain/Amplification	扩增——拷贝数>2（正常二倍体是2）
Loss/Deletion	缺失——拷贝数<2（杂合缺失=1，纯合缺失=0）
LOH	杂合性缺失——一个等位基因丢失，只剩一个
Tumor purity	肿瘤纯度——样本中肿瘤细胞的比例
Ploidy	倍性——细胞中染色体的套数
Log2 ratio	log2比值——肿瘤/正常的覆盖度比值取log2
B-allele frequency (BAF)	B等位基因频率——杂合位点的等位基因比例
Segmentation	分段——将连续的CNV信号分成离散的区段
Focal CNV	局灶性CNV——影响一个或几个基因的小范围CNV

一、安装¶

# === CNVkit安装 ===
conda install -c bioconda cnvkit    # conda安装
cnvkit.py version                   # 验证

# === GATK安装（含CNV流程）===
conda install -c bioconda gatk4     # GATK4安装
gatk --version                      # 验证

# === Sequenza安装 ===
# R包安装
Rscript -e 'install.packages("sequenza")'
# Python工具安装
pip install sequenza-utils          # Sequenza Python工具

二、CNVkit分析（WES/Panel首选）¶

2.1 完整流程¶

# === CNVkit一条命令完成全部分析 ===
cnvkit.py batch \
    tumor.bam \                    # 肿瘤BAM
    --normal normal.bam \          # 配对正常BAM
    --targets targets.bed \        # 捕获区域（exome/panel的target bed）
    --fasta reference.fa \         # 参考基因组
    --access access.hg38.bed \     # 可访问区域
    --output-reference normal_ref.cnn \  # 参考模型输出
    --output-dir cnvkit_results/ \ # 输出目录
    --diagram \                    # 生成染色体图
    --scatter \                    # 生成散点图
    -p 16                          # 16线程

# === 输出文件说明 ===
# tumor.targetcoverage.cnn  — 目标区域覆盖度
# tumor.antitargetcoverage.cnn — 非目标区域覆盖度
# tumor.cnr               — 每个bin的log2比值
# tumor.cns               — 分段后的CNV（关键结果）
# tumor-diagram.pdf        — 染色体图
# tumor-scatter.pdf        — 散点图

2.2 分步运行（更灵活）¶

# === 第1步：计算目标和非目标区域 ===
cnvkit.py target targets.bed \
    --annotate refFlat.txt \       # 基因注释
    --split \                      # 拆分大区间
    -o targets_split.bed           # 输出

cnvkit.py antitarget targets_split.bed \
    --access access.hg38.bed \     # 可访问区域
    -o antitargets.bed             # 输出

# === 第2步：计算覆盖度 ===
cnvkit.py coverage tumor.bam targets_split.bed -o tumor.targetcoverage.cnn
cnvkit.py coverage tumor.bam antitargets.bed -o tumor.antitargetcoverage.cnn
cnvkit.py coverage normal.bam targets_split.bed -o normal.targetcoverage.cnn
cnvkit.py coverage normal.bam antitargets.bed -o normal.antitargetcoverage.cnn

# === 第3步：构建参考 ===
cnvkit.py reference \
    normal.targetcoverage.cnn \
    normal.antitargetcoverage.cnn \
    --fasta reference.fa \
    -o normal_reference.cnn

# === 第4步：修正和分段 ===
cnvkit.py fix \
    tumor.targetcoverage.cnn \     # 肿瘤目标覆盖
    tumor.antitargetcoverage.cnn \ # 肿瘤非目标覆盖
    normal_reference.cnn \         # 正常参考
    -o tumor.cnr                   # 输出log2比值

cnvkit.py segment tumor.cnr \
    -o tumor.cns                   # 输出分段结果

# === 第5步：可视化 ===
cnvkit.py scatter tumor.cnr -s tumor.cns \
    -o tumor_scatter.pdf           # 散点图

cnvkit.py diagram tumor.cnr -s tumor.cns \
    -o tumor_diagram.pdf           # 染色体图

2.3 无配对正常样本分析¶

# === 使用Panel of Normals (PoN) ===
# 当没有配对正常样本时，用一组正常样本构建参考

# 第1步：处理多个正常样本
for normal in normals/*.bam; do
    name=$(basename "$normal" .bam)
    cnvkit.py coverage ${normal} targets_split.bed -o ${name}.target.cnn
    cnvkit.py coverage ${normal} antitargets.bed -o ${name}.antitarget.cnn
done

# 第2步：构建PoN参考
cnvkit.py reference \
    normals/*.target.cnn \         # 所有正常样本的覆盖度
    --fasta reference.fa \
    -o pon_reference.cnn           # PoN参考模型

# 第3步：用PoN分析肿瘤样本
cnvkit.py fix tumor.target.cnn tumor.antitarget.cnn pon_reference.cnn -o tumor.cnr
cnvkit.py segment tumor.cnr -o tumor.cns

2.4 CNVkit结果解读¶

# === Python解读CNVkit结果 ===
import pandas as pd  # 数据处理

# 读取分段结果
cns = pd.read_csv("tumor.cns", sep="\t")  # 读取CNS文件

# 关键列说明：
# chromosome — 染色体
# start, end — 区段起止位置
# log2 — log2(肿瘤/正常)比值
# cn — 估计的拷贝数
# probes — 该区段包含的探针/bin数

# 筛选显著CNV
gains = cns[cns['log2'] > 0.3]   # 扩增（log2>0.3 约对应CN>2.5）
losses = cns[cns['log2'] < -0.3]  # 缺失（log2<-0.3 约对应CN<1.5）

print(f"扩增区段: {len(gains)}")
print(f"缺失区段: {len(losses)}")

# log2比值与拷贝数的关系（假设纯度100%，倍性2）：
# log2 = 0   → CN = 2（正常）
# log2 = 0.58 → CN = 3（单拷贝增加）
# log2 = 1.0  → CN = 4（加倍）
# log2 = -1.0 → CN = 1（杂合缺失）
# log2 = -inf → CN = 0（纯合缺失）

三、GATK CNV Pipeline¶

3.1 体细胞CNV（GATK推荐流程）¶

# === GATK4 CNV Pipeline ===

# 第1步：准备intervals
gatk PreprocessIntervals \
    -R reference.fa \              # 参考基因组
    -L targets.interval_list \     # 捕获区域
    --bin-length 0 \               # 用原始interval（WES设0）
    --padding 250 \                # 边界扩展250bp
    -O preprocessed.interval_list  # 输出

# 第2步：收集覆盖度计数
gatk CollectReadCounts \
    -I tumor.bam \                 # 输入BAM
    -L preprocessed.interval_list \  # 区间
    --interval-merging-rule OVERLAPPING_ONLY \
    -O tumor.counts.hdf5           # 输出HDF5格式

# 对正常样本也做同样操作
gatk CollectReadCounts \
    -I normal.bam \
    -L preprocessed.interval_list \
    --interval-merging-rule OVERLAPPING_ONLY \
    -O normal.counts.hdf5

# 第3步：构建Panel of Normals
gatk CreateReadCountPanelOfNormals \
    -I normal1.counts.hdf5 \       # 多个正常样本
    -I normal2.counts.hdf5 \
    -I normal3.counts.hdf5 \
    -O pon.hdf5                    # 输出PoN

# 第4步：去噪和标准化
gatk DenoiseReadCounts \
    -I tumor.counts.hdf5 \         # 肿瘤计数
    --count-panel-of-normals pon.hdf5 \  # PoN
    --standardized-copy-ratios tumor.standardized.tsv \
    --denoised-copy-ratios tumor.denoised.tsv

# 第5步：收集等位基因计数
gatk CollectAllelicCounts \
    -I tumor.bam \                 # 肿瘤BAM
    -R reference.fa \              # 参考基因组
    -L common_snps.interval_list \ # 常见SNP位点
    -O tumor.allelicCounts.tsv     # 等位基因计数

# 第6步：分段
gatk ModelSegments \
    --denoised-copy-ratios tumor.denoised.tsv \
    --allelic-counts tumor.allelicCounts.tsv \
    --output segments_output \     # 输出目录
    --output-prefix tumor          # 输出前缀

# 第7步：确定CNV状态
gatk CallCopyRatioSegments \
    -I segments_output/tumor.cr.seg \  # 分段结果
    -O tumor.called.seg            # 确定CNV状态（+/-/0）

3.2 GATK CNV可视化¶

# === GATK自带可视化 ===
gatk PlotDenoisedCopyRatios \
    --standardized-copy-ratios tumor.standardized.tsv \
    --denoised-copy-ratios tumor.denoised.tsv \
    --sequence-dictionary reference.dict \
    --output-prefix tumor_plot \
    -O plots_dir/

gatk PlotModeledSegments \
    --denoised-copy-ratios tumor.denoised.tsv \
    --allelic-counts segments_output/tumor.hets.tsv \
    --segments segments_output/tumor.modelFinal.seg \
    --sequence-dictionary reference.dict \
    --output-prefix tumor_modeled \
    -O plots_dir/

四、Sequenza分析（肿瘤纯度+倍性）¶

4.1 数据预处理¶

# === Sequenza-utils预处理 ===

# 第1步：处理参考基因组GC含量
sequenza-utils gc_wiggle \
    --fasta reference.fa \         # 参考基因组
    -o reference.gc50.wig.gz       # GC含量文件

# 第2步：处理BAM文件
sequenza-utils bam2seqz \
    --normal normal.bam \          # 正常BAM
    --tumor tumor.bam \            # 肿瘤BAM
    --fasta reference.fa \         # 参考基因组
    -gc reference.gc50.wig.gz \    # GC含量
    --output tumor.seqz.gz         # 输出seqz文件

# 第3步：合并小区间
sequenza-utils seqz_binning \
    --seqz tumor.seqz.gz \        # 输入seqz
    --window 50 \                  # 窗口大小50bp
    -o tumor.small.seqz.gz        # 输出压缩后的seqz

4.2 R分析¶

# === Sequenza R分析 ===
library(sequenza)  # 加载sequenza

# 读取seqz文件
seqz_data <- sequenza.extract(
    "tumor.small.seqz.gz",        # 预处理后的seqz文件
    chromosome.list = paste0("chr", c(1:22, "X"))  # 分析的染色体
)

# 拟合模型——估计纯度和倍性
CP <- sequenza.fit(seqz_data)  # 拟合

# 查看最佳纯度和倍性估计
cat("最佳肿瘤纯度:", CP$max.cellularity, "\n")
cat("最佳倍性:", CP$max.ploidy, "\n")

# 生成完整结果
sequenza.results(
    sequenza.extract = seqz_data,
    cp.table = CP,
    sample.id = "tumor_sample",    # 样本ID
    out.dir = "sequenza_results/"  # 输出目录
)

# === 输出文件 ===
# tumor_sample_CP_contours.pdf    — 纯度/倍性等高线图
# tumor_sample_genome_view.pdf    — 全基因组CNV视图
# tumor_sample_chromosome_view.pdf — 各染色体CNV视图
# tumor_sample_segments.txt       — 分段结果（含绝对拷贝数）
# tumor_sample_confints_CP.txt    — 纯度/倍性置信区间

五、三工具对比¶

特性	CNVkit	GATK CNV	Sequenza
适用数据	WES/Panel/WGS	WES/WGS	WGS/WES
需要配对正常	推荐但非必须	需要PoN	必须
肿瘤纯度估计	有限	有限	核心功能
倍性估计	基本支持	不直接支持	核心功能
等位基因特异	支持	支持	支持（BAF）
可视化	好（内置多种图）	好（自带脚本）	好（自动生成）
Pipeline集成	灵活	标准化好	独立
速度	快	中等	中等
推荐场景	WES/Panel首选	大规模标准化	肿瘤纯度估计

六、常见报错与解决¶

报错信息	原因	解决方案
`No targets in BED`	BED文件格式错误	检查BED文件（tab分隔，无header）
`Too few bins`	捕获区域太小/太少	增加bin-length或换WGS模式
`NaN in log2 ratio`	覆盖度为0的区域	过滤低覆盖区域
`Segmentation failed`	数据点太少	减少分段最小长度参数
`Cellularity estimation failed`	肿瘤纯度太低	纯度<20%时估计不可靠
`GATK HDF5 error`	HDF5库版本问题	确认GATK版本和依赖兼容

七、面试高频题¶

Q1：CNV检测的基本原理是什么？¶

答：两个核心信号：(1) 读深度（Read Depth）——拷贝数增加的区域会有更多reads覆盖，拷贝数减少的区域reads更少。用肿瘤和正常样本的覆盖度比值（log2 ratio）来检测CNV。(2) B等位基因频率（BAF）——在杂合SNP位点，正常情况两个等位基因约各占50%（BAF≈0.5）。杂合缺失时只剩一个等位基因（BAF→0或1），拷贝数增加可能导致BAF偏移。结合log2 ratio和BAF可以区分不同类型的CNV（如CN-LOH: log2=0但BAF偏移）。

Q2：CNVkit为什么适合WES/Panel数据？¶

答： WES/Panel数据只捕获了基因组的一小部分（外显子区域），传统WGS-CNV方法（基于均匀覆盖假设）不适用。CNVkit的创新在于：(1) 同时利用on-target（捕获区域）和off-target（非捕获区域的少量reads）信号，提高分辨率；(2) 用GC校正和重复区域校正消除WES数据特有的偏差；(3) 支持从多个正常样本构建PoN（Panel of Normals），比单配对更稳定。缺点是对WES数据，CNV断点定位不如WGS精确（因为数据不连续）。

Q3：肿瘤纯度和倍性为什么影响CNV分析？¶

答：实际肿瘤样本是肿瘤细胞和正常细胞的混合物。如果肿瘤纯度只有50%，一个真正的纯合缺失（CN=0）在混合样本中看起来只是CN≈1（因为50%的正常细胞还有正常拷贝数）。同理，肿瘤如果是四倍体（倍性=4），CN=4实际上是"正常"的，不算扩增。不校正纯度和倍性会导致：(1) 低估CNV的幅度；(2) 错误判断CN-LOH；(3) 无法计算绝对拷贝数。Sequenza通过联合拟合log2 ratio和BAF来同时估计纯度和倍性。

Q4：CNV分析在肿瘤研究中有什么应用？¶

答： (1) 发现驱动基因——常见的扩增癌基因（如ERBB2/HER2乳腺癌、MYC多种肿瘤）和缺失抑癌基因（如TP53、RB1）；(2) 指导靶向治疗——HER2扩增→曲妥珠单抗，MET扩增→MET抑制剂；(3) 肿瘤分型——GBM根据CNV分子分型；(4) 肿瘤进化分析——追踪治疗前后CNV变化了解耐药机制；(5) 预后评估——基因组不稳定性（大量CNV）通常预后较差。

参考资料：CNVkit: Talevich et al., PLoS Comput Biol 2016 | GATK CNV: GATK Best Practices | Sequenza: Favero et al., Ann Oncol 2015 | CNV综述: Zhao et al., BMC Bioinformatics 2013