检测多组学高通量测序数据中的个体间污染与样本错配

摘要

本方案描述了一种质量控制框架的实施方法，通过验证个体内部样本对的遗传一致性，来检测下一代测序数据中的个体间污染和样本错配问题。

摘要

通过下一代测序技术对患者生物样本进行高通量处理，并将分子数据与患者层面和样本层面的临床数据进行比较，需要在整个生物样本保管链中精确追踪和匹配样本标识符，这对于实现生物标志物临床试验结果的可靠解读至关重要。除了追踪样本和数据处理流程中的各个步骤外，生物信息学解决方案还可用于确认样本是否来源于同一患者。本文展示了一种利用生物信息学流程识别来自同一受试者的匹配样本的方法。该分析流程适用于任意两组或更多组需进行比较并验证患者样本来源的NGS数据集。通过基于全基因组范围样本比对的评分算法，用户可判断两个样本是否来源于同一个体。具体而言，该方法利用选定连锁不平衡区域内的单核苷酸多态性（SNPs）进行样本识别与比对。研究确定了宽松和严格条件下用于匹配与不匹配样本筛选的阈值组合。该方案的有效性已在超过2,000例患者的临床肿瘤组织和血液样本的质量控制与验证中得到验证，涵盖多种组学检测模式。

引言

大规模收集和分析临床样本需要在样本保管链的各个环节进行精确追踪，因为正确匹配来自相同或不同检测方式的分子数据以及临床患者层面和样本层面的数据，对于准确解读和科学决策至关重要。尽管在良好临床实践规范下已努力优化样本处理流程，但在从活检/样本提取、制备与处理步骤直至数据分析的各个阶段，仍可能发生样本混淆或标签错误（图1）。随着样本数量和处理步骤的增加，样本混淆和交叉污染的可能性也随之上升。这可能导致对样本与患者关系错误的数据进行分析，从而影响后续分析和结论，因此在临床基因组学研究中是一个必须重视的问题。在临床研究中，样本的错误识别会显著影响整体研究结果，尤其对于样本量较小的研究更是如此¹。个体之间的污染会导致在比较同一患者多个样本时，检测差异的能力下降，并产生假阳性结果。样本混淆会影响遗传关联分析的检出效能，并可能导致在全基因组关联分析中低估复杂性状的遗传力²。

癌症研究是开展大规模基因组和转录组分析的领域之一³，特别是用于监测患者间和患者内部的基因组及表型异质性。癌症研究的一个特点是，来自同一患者的样本可能携带不同的突变和拷贝数变异，因而表现出独立的变异等位基因频率⁴。尤其是在解读多种组学数据时，正确整合来自同一个个体内多模态数据集至关重要，因此需要监控个体间的交叉污染^5,6,7,8。针对癌症基因组图谱计划（TCGA）和肺部基因组研究联盟（LGRC）数据集的研究发现，样本错误识别率平均为3%，在某些研究中甚至高达约20%^2,9,10,11。这些实例凸显了监控样本错换和交叉污染现象的重要性¹²。除了在每个实验步骤中进行常规监控和质量控制外，对测序结果的比较分析可作为最终的质量核查手段，确保在进入数据分析与解读之前实现准确的样本匹配。

已有多种生物信息学方法被建立，用于鉴定样本是否来源于同一个体^{1,4,13,14,15,16}。早期方法利用短串联重复序列来验证样本的一致性¹⁷。如今，基于RNA和DNA水平的高通量测序数据，可通过单核苷酸多态性对成对样本进行比对分析¹⁸。这些方法在适用性方面有所不同，例如适用于RNA测序¹⁹或全外显子组测序数据⁵，在实现方式上也存在差异，例如可用于跨测序通道的比对²⁰，同时在使用便捷性上也各有特点。尽管基于20–45个单核苷酸多态性即可识别来自同一受试者的样本，但在癌症研究中常用的低至中等测序深度策略通常需要整合大量SNP位点¹。

本文描述了一种利用单核苷酸多态性（SNP）连锁不平衡区段进行匹配样本质量控制的方法的实施过程及相应调整¹⁵。该方法已被证明具有较低的误报率和错误匹配率，其工作流程支持不同检测模式之间的比较，例如全外显子组测序与RNA测序样本之间的比对，同时也适用于多种数据格式。为了在临床试验样本的大规模数据集中推广应用该方法，生物信息学分析流程采用通用工作流语言（Common Workflow Language, CWL）实现^21,22。由于CWL具有良好的可读性和类似YAML的语法结构，编程经验有限的科研人员也能轻松理解工作流的整体架构和分析结果。CWL的另一关键特性是其分散/聚合（scatter/gather）功能，可实现流程的并行化处理，从而充分利用分配的计算资源。用户可设定特定步骤的执行条件，进一步提升分析的灵活性。CWL还可与完整的工作流管理系统中的其他组件（如数据库存储、图形用户界面和作业调度器）集成，构建一个强大平台，用于创建、运行和维护可重复的科学分析流程。因此，该实现方式有助于用户更便捷地访问工作流，并在既定的工作流管理系统框架下实现数据集的高通量处理。

此外，研究了在匹配与不匹配样本之间调整选择参数阈值的影响，并确定了宽松和严格筛选错配样本的阈值。展示了修改这些参数对样本对选择的影响，以及其在不同组学模式内部和跨模式应用的适用性。有效优化这些参数将使用户能够调整其分析结果解释的严格程度。该工作流程已应用于一组包含数千个样本的大规模临床数据集。

方案

伦理声明：本项关于个体间交叉污染的分析是根据罗氏公司负责任的数据再利用流程，对已完成的I期和II期临床研究中的个体患者层面数据进行回顾性分析而完成的，并符合各研究的主知情同意书要求。每项研究在开展前均已获得伦理委员会/机构审查委员会的批准。所有参与者均已签署知情同意书，同意参加这些研究。

生物信息学工作流程
注意：生物信息学工作流程的实施始于来自下一代测序数据的原始 fastq 文件，例如全基因组、全外显子组或全转录组测序。此处描述的各个步骤已整合到 CWL 工作流程中。

1. 所需参考材料

1. 进行生物信息学分析流程时，需提供以下内容：
   1. 人类参考基因组（GRCh38），格式为.fasta。
   2. 相应的索引文件，即.fasta.fai。
   3. 相应字典文件，格式为.dict。
   4. 与目标基因组区域匹配的单倍型图谱，用于SNP及连锁不平衡区段（例如，Picard build_fingerprint_maps，SCR_006525）。
2. 确保单倍型图谱的头部信息与参考基因组一致。

2. 比对参考基因组、排序与索引

1. 运行 CWL 工作流程，以生成下述步骤所需的输出结果（工具列表见材料表）（图2）。
2. 通过提供成对的 fastq 文件目录或已比对的 bam 文件目录，启动 CWL 工作流程。
3. 提供模式（例如 R1/R2），作为 CWL 命令中的输入正则表达式模式。
4. 此外，还需提供工作流程的运行位置、基因组参考序列以及单倍型映射文件。
   注意：对于全对全比较，将使用完整的输入目录。若仅需比较部分文件，请提供一个包含待比较文件名的逗号分隔文件。高级选项允许选择运行该流程所需的随机存取内存（RAM）和中央处理器（CPU）数量。
5. 使用比对算法将 fastq 文件比对至人类参考基因组。
   注意：根据测序方式的不同，DNA 测序结果使用 BWA-MEM²³，RNA 测序结果则使用支持剪接比对的 STAR 工具²⁴。以下为 STAR 比对的示例代码：
   STAR \
   --readFilesCommand zcat \
   --runThreadN 8 \
   --outSAMmapqUnique 60 \
   --outSAMattributes All \
   --outReadsUnmapped Fastx \
   --outTmpDir /tmp/STARtmp/ \
   --runDirPerm All_RWX \
   --outSAMtype BAM Unsorted \
   --outFileNamePrefix /FILENAME_ \
   --outSAMattrRGline ID:FILEID.L001 SM:SAMPLE \
   --genomeDir /REF/GENOME/DIR \
   --readFilesIn /path/to/FILENAME.R1.fastq.gz /path/to/FILENAME.R2.fastq.gz
6. 使用 SAMtools sort 命令（samtools sort -o FILENAME_OUT.bam FILENAME_IN.bam）按读段坐标对生成的比对二进制比对矩阵（BAM）文件进行排序。
7. 使用 SAMtools index 命令（samtools index FILENAME_OUT.bam）对排序后的 BAM 文件建立索引。
8. 使用 Picard MarkDuplicates 标记并去除重复序列。
9. 使用 SAMtools 对 BAM 文件重新建立索引，并包含读组（RG）信息。
   注意：工作流程正常运行要求 BAM 文件必须包含 RG 标签。

3. 提取指纹

1. 使用已排序和索引的 BAM 文件，通过 Picard ExtractFingerprints 工具识别 SNP 指纹。
2. 仅出于样本间比较的目的，临时存储生成的变异调用格式（VCF）文件。

4. 相似性评分的计算

1. 使用 Picard CrossCheckFingerprints 计算基于连锁不平衡区段的相似性对数似然比分数（LOD），结果以 crosscheck_metrics 文件格式输出，文件名为 crosscheck_metrics.txt。
   注意：该工作流程支持对所有可能的样本对组合进行交叉比较，或根据预定义列表对特定样本进行比较。
2. 删除中间 VCF 文件。
3. crosscheck_metrics 文件为每对测试样本提供四种比较结果。LOD 分数的解读如下：
   LOD 分数 > 0：样本很可能来自同一个体（样本匹配）
   LOD 分数 ≤ 0：样本很可能来自不同个体
4. 通过直方图可视化 LOD 分数的分布，以验证样本匹配的阈值 cutoff 值，例如使用 R 或 Python 软件（图 3）。
5. 建议由包含经验丰富的生物信息学专家的协作团队共同评估证据是否足以明确鉴定样本身份，并考虑是否需要调整所用阈值，例如纳入全部三个 LOD 分数，或与已知来自不同个体的样本的 LOD 分数分布进行比较。
   注意：样本混淆可能导致多个意外匹配（同一供体样本的 LOD 分数 < 0）和错配（不同供体样本的 LOD 分数 > 0）。复杂情况的解读需要跨职能团队与生物信息学专家密切沟通。

5. 代码可用性

该计算工作流程将在 Github 上提供：https://github.com/Roche/sample-matching-workflow。

结果

CWL 工作流的实现
本研究实现了一种基于先前已建立方法的工作流，用于识别样本匹配情况。该方法利用单核苷酸多态性（SNP）的连锁不平衡区域来识别样本错配问题¹⁵。已有研究显示，该方法的分类准确率可达 0% 的假匹配率（FMR）和 0.01% 的假非匹配率（FFR）。与其他方法相比，该方法在高覆盖度和中等覆盖度下与 NGSCheckmate 表现相当，在低覆盖度以及基因组区域重叠较少的情况下优于 NGSCheckmate。而在与 Conpair 和 BAMixChecker 的比较中，结果尚不明确^13,15,25。本文中，该工作流被实现为 CWL 格式，对 LOD 阈值进行了探究与优化，并应用于组织来源样本中 RNA 测序对与 DNA 测序对之间的比较，以及组织样本与外周血样本之间的跨模态比较（图 3、表 1）。该工作流的实现支持对所有可能的样本对组合进行交叉比较，或基于预定义列表在特定样本间进行选择性比较。

该工作流程的输入使用一组选定的单倍型。这些单倍型用于计算连锁不平衡区域中的单核苷酸多态性（SNP）。通过计算不同样本对之间这些SNP区域的对数优势比（LOD）得分，可区分匹配与不匹配的样本。先前研究表明，LOD得分在LOD < -5 和 LOD > 5 范围内能够正确分类匹配的样本对¹⁵。此外，还计算了额外的LOD得分（LOD_SCORE_TUMOR_NORMAL、LOD_SCORE_NORMAL_TUMOR），其考虑了任一样本的肿瘤样本中可能存在的杂合性缺失情况（即在一个样本中为杂合的区域在另一样本中被检测为纯合）。

输入参数和阈值对匹配/不匹配率的影响
本工作流程的评估突显了影响其性能和准确性的三个关键方面。首先，单倍型图谱所覆盖的基因组区域的选择被证明是至关重要的一步。这些区域的选择直接影响匹配过程的区分能力。其次，读段比对策略的组合以及用于指纹提取的特定单倍型图谱显著影响最终的分析结果。上游处理步骤中的差异可能引入细微偏差，并传递至匹配评分中（图4A–B）。第三，确定样本匹配的阈值需经过仔细评估和选择，这一点至关重要。最优阈值可能因具体的数据类型（例如全外显子组测序与全转录组测序）以及所评估的基因组区域而有显著差异。不同的阈值可调节方法的严格程度（高假阳性率与高假阴性率之间权衡）（图4C）。为应对大规模临床样本队列的需求，该方法被改进为可根据所使用的LOD评分组合及比对对象，定义宽松和严格的样本匹配评分标准。第一种阈值设定方法（I）将三个LOD评分（LOD_SCORE、LOD_SCORE_TUMOR_NORMAL、LOD_SCORE_NORMAL_TUMOR）中任一为正值的情况视为匹配。通过纳入TUMOR_NORMAL和NORMAL_TUMOR评分信息，可减轻由于癌症样本拷贝数变异导致的杂合性缺失所带来的影响。相反，为减少不匹配样本中的假阳性，采用了更严格的第二类阈值（II），其包含两种替代性过滤标准：（a）仅当LOD_SCORE为正值时才归类为匹配；（b）对于某一给定样本，即使其LOD_SCORE本身为负值，只要该值大于该样本与其他已知非同源样本（基于记录的患者来源）之间所有成对LOD_SCORE的最大值，仍可归类为匹配。

在本应用中，为了建立宽松阈值，只要任意样本对符合上述任一标准（I、IIa、IIb）被认定为匹配，则该样本对即被视为匹配。这种方法可确保所有被识别为不匹配的结果具有高度可信度，但代价是部分真实的不匹配可能被误判为匹配（即假阴性）。将宽松阈值与更严格的阈值进行比较，结果显示被归类为不匹配的样本对比例发生了变化。在所有分析的研究中，不同方法之间的差异分别为：3.9%（三个LOD评分中任意一个为正值（I））、13.3%（LOD_SCORE必须为正值（IIa））、9.2%（将LOD_SCORE与样本的非匹配样本对进行比较（IIb）），以及3.6%（综合考虑上述任意条件来判定匹配）（图4C）。

基因组区域覆盖度的影响
当覆盖较广泛的基因组区域时（全基因组测序（WGS）或WGS样本与其他检测方法的比较），匹配与不匹配样本之间的LOD得分差异最大，从而有利于阈值的选择（图5A）。在全外显子组测序和RNA测序比较中，LOD得分更接近于零，且阈值设定方法会影响结果，这凸显了在基因组覆盖度较低的技术方法中评估阈值严格性的重要性。已知配对样本中具有阳性LOD得分的结果分布见图5B。Javed等（2020）已证明，当使用连锁不平衡区段时，仅需0.02%的基因组重叠即可有效区分匹配与不匹配的样本¹⁵。

验证
该方法在其他乳腺癌、结直肠癌和肺癌的全外显子组测序（WES）及RNA测序数据集上进行了验证，这些数据集中已知部分样本来源于同一受试者（图6A）。来自同一受试者的样本对显示出100%的匹配率（图6B），而与其他已知来源于不同受试者的样本进行比对时，均显示100%的不匹配率。在所有这些数据集中均未观察到假阳性或假阴性结果。

总之，实施质量控制流程可通过提供标准化且可重复的方法，促进对新一代测序样本对的个体间比较。所得到的LOD评分阈值在基因组区域重叠较大的样本中具有较低的假阳性率和假阴性率；对于测序深度较低或基因组重叠较少的样本，还可进一步应用阈值优化策略。

图1： 样本交换和错贴标签情况的示意图。 （A）两个个体之间各一个样本发生交换。（B）从活检取样到测序数据分析的样本处理步骤示意图。使用BioRender制作。Voith von Voithenberg, L. (2026) https://BioRender.com/xhbp178。 请点击此处查看此图的放大版本。

图2：在CWL中运行样本匹配质控工作流程所需的输入文件和参数的用户界面示意图。用于文件和参数输入的图形用户界面。请点击此处查看该图的放大版本。

图 3：样本匹配工作流程获得的结果。 展示了一组代表性 DNA 测序（全基因组测序和全外显子组测序）样本的 LOD 分数分布（左图），DNA 测序与 RNA 测序之间的比较（中图），以显示错配样本的行为与正确匹配样本分布的显著差异，以及一组代表性较大的 RNA 测序配对样本队列（右图），其中数据已知来源于不同个体（错配，浅红色）或同一个体（匹配，浅绿色）。缩写：LOD = 对数优势比。 请点击此处查看该图的放大版本。

图 4：观察到的 LOD 分数差异示例。（A）结合不同的序列比对方法和单倍型图谱对一组已知不匹配和匹配的样本进行分析时，以及（B）通过整合肿瘤与正常组织信息进行评分时的结果。（C）采用不同严格程度的阈值方法所定义的样本匹配与不匹配数量的分布情况。缩写：LOD = 对数优势比。请点击此处查看该图的放大版本。

图 5：比较不同二代测序模式的LOD分数分布示例。 （A）血液和肿瘤组织DNA测序与肿瘤组织RNA测序之间，预期错配与匹配样本的LOD分数分布。（B）不同模式组合下匹配样本的LOD分数分布。缩写：LOD = 对数优势比。 请点击此处查看该图的放大版本。

图 6：乳腺癌全外显子组测序（WES）和RNA测序数据集分析中LOD分值的分布。 去标识化的乳腺癌数据集来源于Caris Life Sciences，源自全面的肿瘤分子谱分析。（A）肿瘤WES样本间（左）和RNA测序样本间（右）LOD分值的对数发生频率。（B）来自同一受试者的预期样本对的LOD分值分布（WES为上排，RNA测序为下排）。缩写：LOD = 对数优势比。请点击此处查看该图的放大版本。

左侧组值	右侧组值	结果	LOD分数	LOD分数_ 肿瘤_正常	LOD分数_ 正常_肿瘤
样本 1	样本 1	预期匹配	38.119266	29.649485	29.649485
样本 1	样本 2	预期不匹配	-2.552644	-4.57422	5.283698
样本 2	样本 1	预期不匹配	-2.552644	5.283698	-4.57422
样本 2	样本 2	预期匹配	12.328737	8.796457	8.796457

表1：运行Crosscheck指纹分析获得的示例结果。 该表格展示了一对样本通过样本匹配方法进行比对的示例结果。

讨论

目前已有多种用于识别样本匹配的方法^{1,4,13,14,15,16}。本文描述了一种基于单核苷酸多态性（SNP）连锁不平衡区段的方法实现，该方法适用于多种组学模态，且具有较低的假阳性率和假阴性率¹⁵。该方法采用通用工作流语言（CWL）实现，以便在标准化的工作流环境中对多个数据集进行高通量处理。工作流评估确定了应用该方法时需考虑的三个关键方面。其中关键步骤之一是选择单倍型图谱所覆盖的基因组区域。此外，结合不同的单倍型图谱进行读段比对与指纹提取，可能影响分析结果。同时，阈值的审慎评估与选择至关重要，其可能依赖于数据模态及所覆盖的基因组区域，并可能导致样本匹配判定更为宽松或更为严格。

为了评估大量临床样本，该方法经过调整，以定义宽松和严格样本匹配评分的阈值组合。通过考虑包含任意一种检出限评分（LOD_SCORE、LOD_SCORE_TUMOR_NORMAL、LOD_SCORE_NORMAL_TUMOR）或两种替代过滤标准的综合评分，对该方法进行了优化，从而实现对样本更为严格的筛选。该方法的应用仅限于那些在多个基因组区域具有单核苷酸多态性（SNP）信息的样本，例如新一代测序数据。此外，进行比较分析和样本匹配时，至少需要来自同一受试者的成对样本。其他临床信息，如通过靶向方法获得的突变状态或患者元数据（如性别），也可用于进一步验证高维分子数据与患者水平临床数据之间的匹配关系。

在工作流管理环境中实施该方法，并结合并行数据存储的可能性，可实现对样本进行高通量的质量控制分析，以检测个体间的污染。因此，该方法在不同数据集和样本间的可及性与可重复性得以提升。该方法中阈值判定标准具有可调节性，能够处理和分析肿瘤突变负荷低或高的癌症样本，以及存在拷贝数变异的样本；这些变异可能导致杂合性缺失，从而影响基因型似然性。

该方法可广泛适用于任何涉及来自人类个体的高通量测序数据且每个人拥有多个样本的项目。其应用范围可涵盖针对个体患者的个性化医疗方案，到针对不同疾病领域收集高维分子数据的大型临床试验。该方法可与质控工作流程相结合，用于检测跨物种污染，并可整合高维分子数据集与临床信息，从而融入任何高通量测序数据的质控分析流程中。

材料

本文使用的材料清单

姓名	公司	目录编号	评论
FastQC	v0.11.9	SCR_014583
MultiQC	v1.8	SCR_014982
BWA-MEM	v0.7.17	SCR_010910
STAR	v2.7.9a	SCR_004463
SAMtools	V1.12，v1.19.2	SCR_005227
-  faidx
-  分选
-  索引
- 添加替换项
皮卡德	V2.25.5，v3.0.0	SCR_006525
- 创建序列字典
- 标记重复序列
- 构建指纹图谱
- 提取指纹
- 交叉核对指纹
CWL	v1.2	SCR_015528
R	R v4.3.1	SCR_001905
dplyr v1.1.4