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。如果最终结果出现假阳性我可以快速定位是哪个信号出了问题。比如BSRbin-wise scaled read depth信号波动太大问题出在Stage 1的样本QC而如果是某一个工具单独报出的变异过多问题在Stage 2的参数设置。整条链路的可解释性大幅提升。二是便于类似样本的批量扩展。精准医学项目通常不是做一个样本而是几十上百个。分阶段设计配合批量脚本之后每加一个新样本只需要跑一遍Stage 1Stage 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、ExomeDepthWES专用和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清洗探针重叠区域。第二层使用PoNpanel 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几乎为0chrX与常染色体相当。一旦发现实际信号与样本表性别不符立即返回数据源头确认不要带着错误的性别标记进入下游分析。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截图也别因为前期参数太严格把一个真实变异漏掉。这个思路在临床和科研项目中都帮我避免过不少麻烦建议你试试看。 
长河君 · 2026/10/13 11:03:47 · 阅读 1.2k SV/CNV检测pipeline 1.2.3：从信号提取到临床注释的完整指南
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。如果最终结果出现假阳性我可以快速定位是哪个信号出了问题。比如BSRbin-wise scaled read depth信号波动太大问题出在Stage 1的样本QC而如果是某一个工具单独报出的变异过多问题在Stage 2的参数设置。整条链路的可解释性大幅提升。二是便于类似样本的批量扩展。精准医学项目通常不是做一个样本而是几十上百个。分阶段设计配合批量脚本之后每加一个新样本只需要跑一遍Stage 1Stage 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、ExomeDepthWES专用和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清洗探针重叠区域。第二层使用PoNpanel 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几乎为0chrX与常染色体相当。一旦发现实际信号与样本表性别不符立即返回数据源头确认不要带着错误的性别标记进入下游分析。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截图也别因为前期参数太严格把一个真实变异漏掉。这个思路在临床和科研项目中都帮我避免过不少麻烦建议你试试看。