企业官网建设流程全解析

1. 先聊清楚：SV与CNV检测为什么值得单开一条pipeline

做精准医学和基因组学数据分析的人，几乎都绕不开一个现实问题：WES或者WGS数据拿到手后，点突变（SNV/Indel）的处理已经非常成熟，GATK HaplotypeCaller、Mutect2这些工具已经把流程推得相当顺手。可一旦要检测结构变异（SV）和拷贝数变异（CNV），很多团队会下意识地沿用SNP那套pipeline，结果就是掉进坑里：假阳性高得离谱，断点嘈杂，拷贝数分段结果根本没法看。

我自己在被这类问题折磨了很久之后，才理解了一个道理——SV与CNV检测，必须单独设计pipeline，而不是在SNP流程里顺手加几个工具。原因很简单：这两类变异在生物信息学信号上完全不同，点突变只看单个位点的碱基差异，而SV/CNV依赖的是比对信号的结构性改变，比如插入片段长度异常、read方向错乱、覆盖深度波动等。信号源不同，处理策略就必须跟着变。

这篇博客要聊的就是一个可落地的三段式基因组数据处理工程，我叫它pipeline 1.2.3，专门解决SV和CNV检测的问题。整套流程分为三个阶段：信号提取、变异识别与联合分析、过滤注释与结果汇总。它适用于WGS全基因组数据和WES外显子组数据，适合三类人参考：

• 正在搭临床检测流程的实验室生信工程师
• 做肿瘤或罕见病研究、需要批量处理样本的科研人员
• 对SV/CNV检测工具选型感到迷茫的初学者

先说结论：SV与CNV虽然经常被并列讨论，但它们在算法上是两种完全不同的检测任务。SV更依赖双端比对方向和分裂read的断点证据，CNV更依赖read覆盖深度的统计学波动。pipeline 1.2.3的核心思路，就是把这两类任务拆开处理，再在后期统一合并，避免互相干扰。

2. 核心设计思路：为什么pipeline 1.2.3要分三段而不是一把梭

2.1 SV/CNV检测的三个底层信号

在展开pipeline结构之前，必须先把信号源讲透，否则后面理解不了为什么参数要那么调。

做SV/CNV检测时，我们手里的核心输入是BAM文件，也就是测序read比对到参考基因组后的结果。检测SV/CNV，本质上就是从BAM文件里寻找三类异常信号：

第一类：插入片段长度异常。正常双端测序时，一对read的正向和反向端之间的距离大致符合文库插入片段的分布，比如300bp左右的文库，两条read之间的间隔基本稳定在250~350bp。如果某条read对的实际距离远大于预期，或者比对方向出现问题（顺向-顺向、反向-反向），高度提示这个位置存在deletion或duplication。

第二类：分裂read信号。一条read可以被比对上两段不同的基因组区域，这在正常比对中是不允许的。一旦出现这种情况，说明这个read跨越了结构变异的断点。分裂read是定位SV断点最精确的证据，尤其是对于50bp以上的中型变异和更大的结构变异。

第三类：覆盖深度信号。将基因组区域按窗口划分，统计每个窗口的read覆盖深度，拷贝数正常的区域体现为基准深度，拷贝数增加或缺失的区域则表现为明显的深度偏离。CNV检测的核心就是对这个深度信号做归一化、去批次效应、分段估计。

理解了这三个信号，就能明白为什么GATK那套流程没法直接用来做SV/CNV——它根本没有针对上述信号设计统计模型。SV检测工具需要有基于图比对的断点聚类算法，CNV检测工具需要考虑GC偏好、文库制备偏差、区域可比性。这些都不是SNP calling的范畴。

2.2 分段而不是联合建模的三个理由

我在设计pipeline 1.2.3的时候，最初也想过用一套端到端的模型同时输出SV和CNV，比如GRIDSS那样同时利用split read和depth信号。但实际跑下来之后，我逐渐倾向一种组合式设计，把流程拆成三个阶段：

• Stage 1：信号提取与质量控制。此阶段只处理BAM和样本元数据，产出行级别的签名文件，目的是保证后续分析不会把低质量样本误判成拷贝数异常。
• Stage 2：独立检测器调用与联合去重。用多个SV检测器和CNV检测器分别产出初步候选集，再通过共识策略过滤工具特异噪声。
• Stage 3：功能注释与临床关联。将过滤后的变异与人群频率、基因区域、OMIM/ClinVar等注释整合，输出面向临床解读的表格。

这种三段式设计有三个具体好处。

一是便于Debug。如果最终结果出现假阳性，我可以快速定位是哪个信号出了问题。比如BSR（bin-wise scaled read depth）信号波动太大，问题出在Stage 1的样本QC；而如果是某一个工具单独报出的变异过多，问题在Stage 2的参数设置。整条链路的可解释性大幅提升。

二是便于类似样本的批量扩展。精准医学项目通常不是做一个样本，而是几十上百个。分阶段设计配合批量脚本之后，每加一个新样本只需要跑一遍Stage 1，Stage 2和Stage 3可以直接复用之前构建的背景模型。节省的时间非常可观。

三是便于多工具整合。目前没有任何一个单一工具能在SV和CNV两个维度上同时达到高敏感性和低假阳性。DELly对deletion敏感，LUMPY擅长复杂重排，cn.MOPS对WES拷贝数检测稳健，PureCN则考虑肿瘤纯度。pipeline 1.2.3在Stage 2采用多工具并行、共识取交集的策略，让每一个工具的优势互补。

3. 实操前提：BAM质检、参考基因组与样本设计的隐性坑

3.1 进入pipeline之前必须查的几项BAM指标

拿到一个BAM文件，我建议不要直接丢进pipeline。先跑一遍samtools flagstat和picard CollectInsertSizeMetrics，把下面几个指标看一遍：

• 总比对率：人类WGS样本一般要求比对率不低于95%，WES要求比对率不低于85%。如果比对率偏低，后续read depth会表现为全局性下降，容易被误判为拷贝数缺失。
• 重复率：PCR重复率过高（超过15%~20%）会严重抬高read depth的方差，CNV检测的分段边界会变得极其破碎。对于WGS最好查看光学重复与标记重复的比例。
• 插入片段分布：在SV检测里，插入片段是重要信号。如果库检不合格导致插入片段呈双峰分布，SV检测结果会大量出现假阳性断点。我见过不少案例，最后查下来都是文库制备出了问题，而不是算法不行。

这些检查看起来是基础功，但在临床项目里往往是最容易被跳过的一步。我的习惯是：整批样本先统一跑完QC指标生成汇总表，任何偏离阈值的样本单独标注，在进入Stage 2之前直接与测序方沟通是否重新上机或调换样本。这一步虽然增加了流程的前置时间，但能省下后面很多无效分析。

3.2 参考基因组版本与contig命名

SV/CNV检测对参考基因组版本的选择非常敏感，这一点经常被低估。不同版本之间的坐标转换（liftOver）会引入断点重映射误差，尤其是在重复区域和着丝粒附近。所以整个pipeline 1.2.3的输入输出必须锁定同一个参考基因组版本，比如GRCh37或GRCh38，团队内部不能混用。

另外一个容易被忽略的问题：contig命名不一致。有的BAM文件用的是chr1、chr2，有的则用1、2。如果SV检测工具和后期注释数据库的contig命名不统一，结果里会出现大量无法匹配的坐标。我建议在进入pipeline前统一用samtools reheader或者比对阶段直接指定统一的命名方式，最省事的做法是下载参考基因组时就直接选择带chr前缀的版本并保持一致。

3.3 性别、样本关系与批次效应对CNV检测的影响

CNV检测依赖覆盖深度的统计模型，因此性别的影响特别大。男性样本的X染色体覆盖深度大约是女性的一半，Y染色体几乎为零，女性样本的Y染色体则完全不存在。如果性别标注错误，X染色体的拷贝数结果会出现系统性偏差，甚至被模型判定为整条染色体缺失。

对于肿瘤样本，还需要特别考虑纯度和倍性。肿瘤细胞占比低时，拷贝数变化的信号会被大量正常细胞稀释。举例来说，一个纯度为30%的肿瘤样本，如果真实拷贝数是纯合缺失（拷贝数=0），测序深度大约只下降到基准深度的70%，很容易被模型漏检。因此在pipeline 1.2.3中，我强烈建议肿瘤样本搭配PureCN或者FACETS这类能够同时估计肿瘤纯度和倍性的工具，而不是只依赖单纯的深度波动模型。

批次效应也是CNV检测里非常棘手的难题。不同测序批次之间文库制备、上机浓度、测序仪校准状态不同，会导致覆盖深度出现系统性偏差。最实用的做法是：如果样本量够大（比如超过50例），先用主成分分析或者PCA对深度矩阵做批次可视化，发现明显批次聚集时加入Combat-seq或者Peerwise等批次校正步骤，再做后续分段。

4. pipeline 1.2.3核心实操：三个阶段究竟怎么跑

4.1 Stage 1：从BAM到可计算的信号文件

Stage 1在逻辑上解决的是“数据准备”问题，核心目标是把BAM文件转换成一个干净的信号集合，供后续变异检测器消费。

对于SV检测，这一步通常需要：

• 提取双端read信息，输出包含read名称、染色体、比对起止位置、比对链、插入片段大小的文件。
• 提取分裂read信息，用于定位断点。
• 做PCR重复标记，并把重复read排除在SV信号之外，否则容易产生假断点。

对于CNV检测，这一步的产物通常是按固定窗口（比如WGS用1kb、WES按捕获区间）统计的read深度矩阵。与此同时需要计算GC含量、比对复杂度等协变量，为后续归一化做准备。

具体命令行上，我会这样组织：

# 假设输入是hg38比对后的sample.bam # SV信号提取：manta也可以同时产出SV和简单CNV的候选，这里先做SV部分 configManta.py \ --bam sample.bam \ --referenceGenome hg38.fa \ --runDir manta_out python manta_out/runWorkflow.py -m local -j 8 # CNV信号提取：用GATK CollectReadCounts按窗口统计深度 gatk CollectReadCounts \ -I sample.bam \ -L wes_target.preprocessed.interval_list \ --interval-merging-rule OVERLAPPING_ONLY \ -O sample.hdf5

以上只是示例，真正的生产流程会用Snakemake或者Nextflow来编排这些步骤。我的建议是：不要在Stage 1里试图做任何变异判断，只产出信号中间文件即可。这种设计的用意非常明确——pipeline是可重入的，当某个新版本检测器发布时，无需重跑整个流程。

一个值得注意的细节：在WES中做CNV时，窗口划分非常讲究。外显子区域长度差异大、GC含量差异明显，直接按固定大小窗口切分会导致局部窗口深度极低或过高，直接把CNV信号淹没。比较稳妥的方案是用预设的捕获区间作为窗口，同时对探针区域做padding（例如左右扩50bp），再把邻近探针的深度做合并降噪。如果使用GATK的CollectReadCounts，建议配合denoise后的间隔列表使用，而不是原始bed直接硬切。

4.2 Stage 2：多工具联合分析如何协作

Stage 2是整个pipeline最复杂的部分，多个检测器各司其职，再汇总到一个共识层面。这里必须说清楚：为什么我不依赖单一工具。任何一个现有SV/CNV检测工具都有自己偏好的变异类型和偏差来源。比如Manta对deletion和duplication的检测率很高，但它的split read策略在小插入和复杂重排上并不可靠；cn.MOPS在WES CNV上表现优秀，但对肿瘤样本的低纯度容忍度极差。

所以pipeline 1.2.3的策略是：在SV侧同时运行DELly、LUMPY、Manta，在CNV侧同时运行cn.MOPS、ExomeDepth（WES专用）和GATK gCNV。各工具运行完毕后，利用Survivor或者BEDTools对SV结果做合并，利用自建的区间重叠逻辑对CNV结果做一致性评估。

这套策略在业内其实并不新鲜，肿瘤领域的SEQC2项目以及几项大型测序项目都采用类似的一致性方案。我的经验是：共识策略能显著降低假阳性，但会以损失部分灵敏度为代价。所以临床上建议同时输出“至少一个工具检测到”和“两个及以上工具检测到”两层结果，由遗传分析师综合判断。例如肿瘤大panel中如果一个CNV只在Manta中检出而gCNV没有检出，这个位点仍然值得一看，因为很可能是Manta捕获了细小断点而深度信号被周边的CNV背景污染了。

SV和CNV不一致的情形需要单独处理。例如在同一区域，LUMPY给出了一个deletion候选，而cn.MOPS给出拷贝数正常的结论。这类矛盾通常出现在大片段重复区域或者同源序列区域。我在实践中会额外调取BAM文件的可视化截图（用IGV查看覆盖度和split read分布），而不是轻易丢弃任何一个信号。这种人工复核虽然在自动化流程里显得“笨”，却是避免临床漏检的关键。

4.3 Stage 3：过滤、注释与出报告

Stage 3解决的是“变异信号如何变成临床可读信息”的问题。

先把共识集里的候选变异做初步过滤。过滤项包括：

• 与参考基因组中已知的重复区域或片段重叠超过一定比例（比如50%），直接标记为低可信度。
• 在人群数据库中频率过高。SV/CNV同样有频率概念，比如gnomAD SV数据库。如果一个deletion在人群中的等位基因频率超过1%，大概率是中性的多态性事件而不是致病性变异。
• 断点周围是否存在低复杂度序列，这类区域经常产生比对噪音。

接着做功能注释。SV/CNV的注释不仅仅是重叠基因，还需要考虑：

• 基因是常染色体还是性染色体。
• 变异是否落在外显子区域并影响编码序列。
• 是否与已知致病区间重叠，例如ClinVar中的pathogenic CNV、OMIM中记录的剂量敏感基因。
• 是否影响基因融合，特别是肿瘤中的驱动性融合事件。

我通常会使用AnnotSV做注释，因为它整合了多个数据库和基准，输出格式也比较适合后续做Excel筛选。AnnotSV的输出中有一列AnnotSV ranking score，这个分数综合了ACMG指南相关要素、人群频率、基因剂量敏感性、功能影响，非常实用。对于肿瘤样本，还需要补充FusionCatcher或Arriba的融合基因注释结果（如果panel支持转录组信息），用于确定表达的融合是否真的被转录。

最终报告格式我建议采用长表风格，每行一个变异，包含样本ID、染色体、起始、结束、类型、大小、基因、功能影响、注释评分、交叉验证支持的工具有哪些、IGV验证链接。这个表格同时满足生信工程师的后续率值提取和临床遗传学家的审核需求。

4.4 一次真实运行的参数调整记录

这里分享一个我近期跑过的WES样本的案例。样本是某罕见病家系先证者，临床高度怀疑存在小片段缺失导致单倍剂量不足。我按照pipeline 1.2.3执行后，在Stage 2中cn.MOPS报出了一个约130kb的杂合缺失，但ExomeDepth的Z-score为-1.9，刚好低于阈值，没有报出来。两个工具不一致，按照流程进入人工复核环节。

我先查了该区域的GC含量分布，发现该区域平均GC为58%，在WES探针上覆盖度偏低。随后调取IGV截图，确实能够看到覆盖深度略低于相邻区域，但SD并不显著。随后我用GATK gCNV重新跑了这个样本加上另外20个正常对照，发现该区域在gCNV中被判定为纯合缺失趋势。最后综合三个工具的结果和家系共分离数据，这个缺失被报告为“可能致病性”。

这个案例说明：参数阈值不是死的。cn.MOPS对WES的灵敏度较高，ExomeDepth依赖对照样本的同质性和偏正态分布假设，而gCNV对此类边界型缺失更敏感。在实际临床项目中，对于边界信号，手动复核和家系验证永远是不可省略的步骤。

5. 常见问题与排查技巧实录

5.1 假阳性SV断点大量聚集在重复区域

我刚开始跑LUMPY的时候，碰到过一种让人非常头疼的情况：几乎每个样本都会在特定染色体位置（例如染色体9的着丝粒两侧）报告几十个SV断点。最初我以为是样本污染，后来才发现是参考基因组中的segmental duplication（节段重复）捣的鬼。

排查思路：

• 首先比对每个断点区域与UCSC Genome Browser的segdup track是否重叠。
• 对断点做聚类可视化，看看是否呈成对聚集。
• 在注释阶段对segdup区域制定“低置信度”策略，要么直接过滤，要么改成人工复核数据单独标注。

节段重复区域的read比对是多重匹配的，任何依赖唯一比对信号的工具都会在断点检测上失真。这不是软件bug，而是比对原理的必然。后期我习惯直接用RepeatMasker和UCSC的简单重复track给变异注释增加一个repeat_mask列，临床分析人员自行决定是否关注这些位点。

5.2 WES样本的CNV结果出现“波浪形”深度

WES数据做CNV检测时，最经典的问题是深度信号呈现系统性波浪——长范围GC偏好造成的。捕获探针的GC含量分布越不均匀，这种波浪越严重。表现在最终结果上，就是CNV分段图到处是小的拷贝数增减，重复率很高，完全没法直接解读。

解决这个问题有几个层次：

• 第一层：在CollectReadCounts阶段就采用GATK的--interval-merging-rule OVERLAPPING_ONLY清洗探针重叠区域。
• 第二层：使用PoN（panel of normals）建立基线。用至少20个同批次或跨批次但来源相似的正常样本构建PoN，再对每个样本做标准化。这一步的效果远大于任何事后归一化。
• 第三层：用GATK gCNV的denoise模块重跑一遍，它会自动学习GC协变量和测序平台噪声特征。

在我实操的项目中，有了PoN之后，ExomeDepth和cn.MOPS的假阳性率至少下降了60%。所以强烈建议在pipeline中预留构建PoN的正常样本库。注意：PoN的样本越多样越好，但不能包含目标变异对应的异常样本，否则目标信号会被“学习”成正常背景。

5.3 不同工具结果如何判断谁更可信

遇到工具间冲突的结果时，很多新手会陷入工具崇拜，认为某个软件名气大，就默认它是正确的。我的建议是建立一套多信号证据打分：

• 如果一个SV同时证据链完整（方向异常、分裂read、覆盖度改变），给高分。
• 如果一个CNV边界与已知连锁不平衡模式或基因编码边界高度重合，给高分。
• 如果一个变异只在单一工具中出现，且该工具输出里还伴有大量低质量断点，则降低信任度，直接进入人工复核阶段。

一个很实用的小技巧：用BEDTools把变异区域与BAM的split read和discordant read信号做交集。实际动手比对时，可以写一行命令：

bedtools intersect \ -a candidate_svs.bed \ -b sample.split_reads.bed \ -wa -c > sv_split_support.tsv

当分裂read支持数大于等于2时，这个SV基本可信；如果是deletion，还应该看到覆盖深度在区域内的下降。如果一个SV有方向异常和覆盖深度改变，但没有分裂read支持，大概率是中等大小（200bp~10kb）的缺失，因为大片段缺失的同源末端重组通常不产生分裂read。用这种思维做证据链推理，比单纯依赖软件打分踏实得多。

5.4 性别表写错了导致的X染色体假缺失

CNV检测中一个很低级但很常见的坑是样本性别标注错误。男性样本被标成女性后，X染色体上半数的深度缺失会自动被解释为“X染色体杂合缺失”。同样，女性样本被标成男性后，X染色体的深度在PoN标准化后可能出现高拷贝数调高的假信号。

排查方式非常简单：在Stage 1阶段就利用BAM文件中的X/Y染色体的read深度估算遗传性别，然后与样本元数据表交叉比对。计算方法可以简单粗暴一点：

# 计算X和Y染色体的平均覆盖度 samtools view -b sample.bam chrX | samtools depth - | awk '{sum+=$3; n++} END {print "chrX mean depth:", sum/n}' samtools view -b sample.bam chrY | samtools depth - | awk '{sum+=$3; n++} END {print "chrY mean depth:", sum/n}'

对于正常男性样本，chrY的平均深度大约是常染色体的50%~70%（男性样本Y染色体的无假常染色体区比X和常染色体低），chrX大约是常染色体的50%。对于女性样本，chrY几乎为0，chrX与常染色体相当。一旦发现实际信号与样本表性别不符，立即返回数据源头确认，不要带着错误的性别标记进入下游分析。

6. 一个容易被忽视的元问题：pipeline版本与重复性管理

SV/CNV检测pipeline和SNV流程相比更脆弱，因为工具版本更新频繁、参考基因组版本和数据库版本不断变化。即使在同一个团队内部，几个月前跑出来的结果和新跑的结果也可能对不上。

我强烈建议把pipeline 1.2.3做成完全可复现的形式。具体可以做到这样几个层面：

• 使用容器封装每个工具：所有工具的版本、运行时环境固定在Docker或Singularity镜像里，镜像名称中带版本号，比如delly:1.1.6、manta:1.6.0。
• 使用工作流管理器：Snakemake或Nextflow流程中明确记录每个步骤的输入输出、参数和软件版本。这样即使将来重跑，也能在流程报告中回溯哪一步是用什么参数跑的。
• 对中间文件做哈希检查：在数据量允许时对每个中间文件的MD5进行记录，防止磁盘静默损坏影响结果。

实际执行中我还会在流程最后生成一个汇总网页或文本报告，包含样本ID、性别、BAM文件名、比对率、重复率、检测出的变异数目、变异类型直方图、各个工具的输出数量。这个报告既给合作医生看，也给我自己留档。有了这个意识之后，返工的概率大幅下降。

7. 写在最后的几个实操心得

pipeline 1.2.3这个三段式结构，我前前后后调了将近两年，中间踩过的坑远比我在这篇博客里写的多。最深刻的一条体会是：没有万能的参数组合，只有不断积累适合自己样本类型的基线。同一个工具和同一套参数，在不同测序平台、不同捕获试剂盒、不同人群标本上，效果可能差异巨大。所以我在每个新合作项目启动时，都会先跑一小批“阳性对照”样本（比如已知携带特定SV/CNV的细胞系或者家系样本），用它们校准阈值，再扩展到全队列。

另一个很实际的体会是：SV/CNV检测报告一定要留出“人工复核”这一层，而且这个环节不能省。很多自动化工具能帮你把范围缩小到十几个候选位点，但最终决定是否上报的，还是人对IGV图的判断、对reads分布的理解、对临床病史的把握。技术永远是辅助，讨论清楚生物学证据链才是终点。

最后分享一个小技巧：在Stage 2里我一般会故意放宽弱信号的阈值，宁可多留一些低置信度候选，也不会在早期就过滤得太死。原因很简单——后期过滤成本低，重跑成本高。宁可后期多花五分钟看一个IGV截图，也别因为前期参数太严格把一个真实变异漏掉。这个思路在临床和科研项目中都帮我避免过不少麻烦，建议你试试看。