cell-free tumor DNA序列数据的生物信息学分析

Bioinformatics Analysis for Cell-Free Tumor DNA Sequencing Data

DOI: 10.1007/978-1-4939-7717-8_5

摘要

ctDNA (cell-free tumor DNA)是液体活检的主要生物标志物，可从血液、尿液或其他循环液体中提取，能够提供全面的肿瘤遗传信息，较组织活检更好地克服肿瘤异质性问题。下一代测序技术是近年来发展起来的一种广泛应用于ctDNA分析的技术。尽管处理ctDNA样本的技术已经成熟，但从噪声测序数据中检测低突变等位基因频率(MAF)变化的任务仍然具有挑战性。在本章中，作者将首先解释ctDNA测序数据分析的难点，回顾相关技术，然后提出一些新的生物信息学方法，以更好地分析ctDNA NGS数据。

Key words Liquid biopsy, Circulating tumor DNA, ctDNA, Gene fusion, CNV, Mutation visualization, OpenGene

---

介绍

ctDNA及其应用

ctDNA如何测序

ctDNA NGS数据分析的难点

无细胞游离肿瘤DNA只是cfDNA的一小部分，特别是对于早期癌症患者的样本。这一事实使检测肿瘤特异性突变变得困难。此外，PCR和测序错误、DNA氧化损伤和软件引入的伪象会产生高水平的噪声，并引入许多假阳性突变。

肿瘤特异性DNA的含量可以从低于0.01%到超过90%不等。ctDNA丰度的可变性与肿瘤负荷、分期、血管密度、细胞周转和对治疗的反应有关。理论上，任何部分的DNA改变都可以通过足够数量的分子深度测序检测到。然而，在异质混合物聚合酶链反应中，扩增偏置会导致群体倾斜，聚合酶错误会导致错误的碱基合并和重排。此外，在测序过程中产生的错误可能导致约0.1-1%的错误碱基调用[6]，这称为测序错误。表1显示了不同主要NGS平台的错误率。

文库准备也可能引入重大错误。例如，鸟嘌呤氧化是人工突变的一个重要来源，因为8-oxoG倾向于与腺嘌呤配对，而不是胞嘧啶。长时间的热培养，这在许多DNA提取和混合捕获协议中是常见的，可以显著增加G-->T替换。最近一项研究表明，DNA修复过程可以分别消除77%和82%的G-->T和C-->A错误。这项研究表明，DNA损伤会导致大量的错误。

除了样品制备和测序过程中引入的误差，软件和分析工具也会引入误差。特别是在同源序列和重复序列的参考基因组区域可以称为假阳性变异。

细胞游离DNA片段通常很短，在167 bp附近有一个紧凑的峰值。这一事实增加了两个不同的原始cfDNA片段共享相同序列的可能性，因此增加了删除这些重复的难度，因为重复数据删除算法将无法区分这些由放大引起的相同和重复读取。

总之，从嘈杂的ctDNA测序数据中检测低频突变具有挑战性。传统的工具不能很好地处理ctDNA分析任务，因此需要更专门的工具。

ctDNA测序数据分析管道

为了分析ctDNA测序数据，需要使用一系列的软件工具。例如，来自Illumina测序仪的原始测序数据以基调用(BCL)格式获得。这个BCL文件需要根据样本条形码进行解复用以分离FASTQ文件。然后用质量控制工具测量FASTQ文件，以确保它们满足质量要求，并过滤掉低质量和错误表示的读取。接下来，过滤后的FASTQ文件将使用对齐器与参考基因组进行对齐，输出应该是SAM/BAM文件。然后需要对BAM文件进行排序并删除重复项。然后，变量调用者需要处理BAM文件，并生成带有原始变量记录的VCF。接下来，这个VCF文件应该用dbSNP和COSMIC这样的数据库进行注释。基线技术将被应用于标记一些假阳性突变，然后支持每个突变的唯一读数将被计数，从而形成一个完整的VCF。然后，这个VCF文件将被过滤，生成一个干净的文件，并使用交互式分析工具进行可视化。最后对目标突变进行解释和报告。这些工具可以被安排成一个管道。图3展示了作者经常使用的ctDNA测序数据分析管道。

对于Illumina平台，使用bcl2fastq工具将BCL格式文件转换为FASTQ文件。Illumina平台通过对不同的样本使用不同的条形码来支持multiplexing，因此在转换的同时也进行了de-multiplexing。

一些附加的工具可以对FASTQ文件执行质量控制和数据过滤，例如FastQC和Trimmomatic。作者建议使用AfterQC，它是高度优化的ctDNA测序数据处理。AfterQC将在下一节中介绍。

许多校准器可以用于将DNA测序reads映射到参考基因组，如bowtie2和BWA。根据我们的实践，BWA在对齐质量和速度方面都有较好的性能。BWA是一个软件包，用于根据大型参考基因组绘制低发散序列。该算法由BWA-backtrack、BWA-SW和BWA-MEM三种算法组成。BWA-MEM通常推荐用于高质量的查询，因为它更快、更准确。但要注意，BWA和任何其他对齐器仍然可能引入不对中，特别是在有重复或同源序列的参考基因组区域。

对齐过程将生成一个SAM文件，其中包含可以立即转换为BAM的对齐信息，BAM是SAM的二进制标识。这个BAM文件通常是无序的，应该对其进行排序并建立索引。最常用的BAM文件排序和索引工具是Samtools，还有其他一些工具可以更快地对BAM进行排序。例如，Sambamba是一个使用SAM/BAM数据的高性能工具。Sambamba是用D语言编写的，其源代码可在以下网站获得: https://github.com/lomereiter/sambamba。

在BAM文件被排序和索引之后，一个可选的过程是应用重新排列来改进插入和删除(INDELs)的检测。一些工具，如ABRA[16]可以执行基于程序集的重新排列以输出更干净的indel，但这些工具通常很慢。现在可以应用BAM文件的质量控制来评估数据的对齐质量并检测不必要的偏差。这个过程可以用Qualimap这样的工具来完成。

后续处理就是重复数据删除。Samtools rmdup和Picard markduplicate (http://picard.sourceforge.net))通常用于根据读的映射坐标和质量分数识别和分解读重复。由于cfDNA片段较短，长度分布紧接近167 bp，不同原始DNA片段的许多reads可能共享相同的映射坐标，不应视为重复。因此，我们不建议使用Samtools rmdup或Picard markduplates进行重复数据删除，我们将在下一节中讨论新的方法和策略。

变异调用是BAM操作(排序、重新排列、去重)之后的关键过程。众所周知，癌症基因组具有广泛的突变，包括单核苷酸变异(SNVs)、多核苷酸变异(MNVs)、小插入和删除(INDELs)以及复杂变异(如拷贝数变异(CNVs)和基因融合)。许多不同的调用者，如GATK HaplotypeCaller， FreeBayes (https://github.com/ekg/FreeBayes)， MuTect2和VarScan2，可以用来调用SNV, MNV和小型INDELs。根据我们的经验，GATK HaplotypeCaller和FreeBayes并不擅长从超深测序数据中调用ctDNA的低频体细胞突变，因为它们最初是设计用于基因分型和发现遗传多态性的。MuTect2在调用体细胞突变方面要好得多，特别是在肿瘤-正常配对数据时。然而，它仅能很好地处理组织测序数据，但对检测ctDNA测序数据中的低频突变不够敏感。VarScan2在检测低频突变方面非常敏感，但可能报告大量假阳性突变。因此，我们无法找到一个完美的变异调用器来检测ctDNA测序等超深NGS数据中的低频突变。目前我们建议使用VarScan2，结合严格的变量过滤。请注意，一些可变调用程序，如GATK HaplotypeCaller，不能很好地随深度扩展，通常是下采样(随机删除数据的部分)来提高计算性能。然而，下采样会显著降低检测低等位基因频率突变的敏感性，不建议用于ctDNA测序数据分析。

变量调用过程完成后，得到原始VCF文件。该VCF文件可以使用ANNOVAR等注释工具进行注释，获取编码序列和蛋白质变化，并与dbSNP、ClinVar、COSMIC等数据库进行比较。

突变基线将用于注释每个变体，说明该变体在过去的数据中被记录了多少次。该信息可用于过滤由软件工件和其他常规系统错误引起的假阳性突变。基线技术将在下一节中介绍。

为了更准确地计算每个突变的支持读数，我们可以将具有相同映射坐标的读视为一个唯一的读。一个名为MrBam (https://github.com/OpenGene/MrBam)的工具用于计算每个突变的惟一引用支持和惟一替代支持。

在完成唯一读取计数之后，我们获得一个完整的VCF文件。这个VCF文件中的记录可以添加到突变基线中。这个VCF文件可以根据不同的条件进行过滤，以尽可能多地去除假阳性突变。在筛选过程中，通常使用由重要临床靶点(即癌症可药物突变靶点)组成的白名单，以避免重要靶点突变被意外过滤掉。

另一方面，可以使用像MutScan (https://github.com/OpenGene/MutScan))这样的工具可视化被调用的变体，以生成用于交互分析的突变可视化。对癌症诊断和治疗很重要的突变将被人工解释。

除了snv和INDELs，另外两种重要的癌症诊断变异是基因融合和拷贝数变异(CNV)。这些工具中的大多数只能使用已排序的BAM文件。例如，DELLY和Factera可用于检测基因融合，CNVkit (https://github.com/etal/CNVkit) 可用于检测靶向DNA测序中的基因扩增。一个例外是，作者开发的工具FusionDirect可以直接使用FASTQ文件来检测目标融合。

作者创建了一个开源项目来演示这个管道，该项目可在GitHub (https://github.com/OpenGene/ctdna-pipeline)上获得。通过研究它，读者可以了解如何安装工具，准备所需的数据库和参考数据，并使用FASTQ文件尝试管道进行测试。

在上面介绍的管道中，超过一半的工具是常用的软件(例如，BWA、Samtools和VarScan2)，而其余的工具是由作者开发的(例如，MutScan、AfterQC和MrBam)。这些新开发的工具是高度优化的ctDNA测序数据分析。这些工具中的大多数都是GitHub组织OpenGene (https://github.com/OpenGene))下的开源项目。我们将在下一节中介绍其中一些。

---

新方法

由于肿瘤特异性DNA只是cfDNA的一小部分，ctDNA中体细胞突变的突变等位基因频率(MAF)通常很低。为了检测具有如此低MAF的突变，我们应该应用目标捕获和超深测序(例如，10000×或更深)。然而，这种超深测序中的测序错误和实验错误（即PCR错误）可能会导致高水平的背景噪声，并且很难从ctDNA NGS数据中检测到高灵敏度和特异性的突变。此外，基因融合的检测也很困难，因为cfDNA片段通常很短，而肿瘤特异性DNA片段太少。由于肿瘤细胞中拷贝数的变化只会导致cfDNA的拷贝数总数的微小差异，因此检测拷贝数变异(CNV)比检测融合更加具有挑战性。

在本节中，我们将介绍一些新方法来部分解决上面列出的问题。其中一些是由作者开发的，并已在我们的常规管道中使用。

更好的数据预处理

数据预处理是为下游分析获取更清晰数据的重要步骤。对于NGS原始数据（FASTQ格式），有必要丢弃低质量读取、去除接头并应用其他过滤器。此外，还需要质量控制（QC）方法来确保数据符合质量要求。

一些好的工具可以执行质量控制，例如具有每碱基和每序列质量分析功能的FastQC和具有FASTA/FASTQ统计功能的PRINSEQ，而其他一些工具可以执行读取修剪，例如Trimmomatic和SolexaQA。由于数据过滤的方式取决于QC结果，并且过滤的数据也需要后过滤QC，因此仍然需要一个具有丰富QC和过滤功能的工具。

由于cfDNA片段通常较短（~167 bp），2×150对端序列将导致读对重叠。基于这个事实，我们可以对配对末端测序数据进行重叠分析。当DNA模板长度小于测序长度的两倍时，这对读取将重叠。请注意，重叠区域中的每个碱基实际上被测序了两次，因此这些碱基对的不一致可能反映了测序错误。

AfterQC是作者开发的一种工具，用于解决许多实际测序数据质量控制和过滤问题。除了常规的质量控制功能，如每个循环的碱基数和质量统计，AfterQC还提供了许多新功能，如自动修剪和重叠分析。例如，我们发现一些测序器（如Illumina NextSeq系列）可能会输出大量具有高质量分数的polyX读取。

AfterQC可以使用其polyX过滤器去除它们，而普通质量的过滤器不能。我们还发现，如果扩增或测序过程存在严重的链偏差，序列读取将显示K-MER计数偏差（即ATCGATCG及其反向补体CGATCGAT的计数显著不同）。基于这一发现，AfterQC提供了基于K-MER计数的链偏差分析。该工具的另一个主要贡献是对成对末端测序数据进行重叠分析，可用于分析测序错误率，并将其用于错误基础校正或消除。对于单个或一对FASTQ文件的每个输入，AfterQC输出一个HTML报告，其中包含质量控制和数据筛选摘要，以及一个交互式图形列表。表2 显示了AfterQC和其他NGS质量控制或过滤工具的特性比较。

AfterQC旨在批量处理FASTQ文件。它通过一个文件夹，其中包含所有FASTQ文件（可以是单端或成对端输出），这些文件通常是不同样本的测序运行数据，并将每个FASTQ或成对文件传递到QC和过滤管道。首先，AfterQC将运行bubble检测，以查找测序过程中出现的bubble；其次，将进行预过滤QC，以使用每个周期的基本内容和质量曲线对数据进行剖面分析；第三，AfterQC将根据数据质量评测执行自动读取修剪；第四，每次读取都将通过气泡过滤器、polyX过滤器、质量过滤器和重叠分析过滤器进行过滤，未通过这些过滤器的将被视为错误读取而丢弃；第五，基于重叠分析的误差校正将应用于配对末端测序数据；最后，AfterQC将存储正确的读取，执行过滤后QC分析，并生成HTML报告。

（AfterQC检测bubble的原理是什么？）

AfterQC可以处理FASTQ数据的自动修剪。有两种修剪策略，本地策略和全局策略。有些工具，如Trimmomatic，应用本地策略，逐reads进行裁剪。然而，局部修剪策略有一些缺点。第一个缺点是局部修剪仅使用质量信息进行修剪，而不能利用全局统计信息来发现异常cycle。第二个缺点是局部修剪会导致未对齐的修剪，这意味着重复读取可能会被不同地修剪，从而导致某些重复数据消除工具（如Picard）失败。大多数重复数据消除工具仅通过对具有相同映射位置的读取进行群集来检测重复数据。相反，AfterQC执行全局修整策略，即以相同的方式修整所有读取。使用一种算法来确定要在前部和尾部修剪多少个循环，该算法基于每个循环基本内容曲线和基本质量曲线的分段。

AfterQC的一个主要优点是重叠分析。让T表示测序DNA模板的长度，S表示成对双端测序长度的长度，那么如果T≤S、 重叠长度将是2S-T、 如果Sed（R1o，R2o）

基于重叠分析，AfterQC可以检测不匹配。如果不匹配对的质量分数不平衡，这意味着一个基数具有高质量分数（即>Q30），而另一个基数的质量分数非常低（即

重叠分析可用于自动适配器切割。在重叠分析过程中，我们得到了每对最佳局部对齐的最佳偏移量O。该对的重叠长度可以使用偏移量O直接计算。如果O为负值，重叠区域外的碱基将被视为适配器序列的一部分，然后自动切割。

AfterQC是一个开源工具：https://github.com/OpenGene/AfterQC。它是用Python和C++实现的，启用了PyPy支持。AfterQC为每个输入生成一个独立的HTML报告，图按Plotly绘制。示例报告可在以下网址找到：http://opengene.org/AfterQC/report.html.

分子条形码测序及其数据分析

PCR和测序方法引入的系统性错误阻碍了ctDNA NGS深度测序的潜力。分子索引与深度测序相结合有望打破PCR和测序错误的限制，并能够检测罕见和极罕见的突变。

自2007年以来，人们提出并报告了用分子条形码标记单个模板。分子条形码或分子索引有多种名称，如唯一标识符（UID）、唯一分子标识符（UMI）、引物ID、双链条形码等。它们通常设计为一串完全随机的核苷酸（如NNNNNNNN）、部分退化的核苷酸（例如NNNRNYNN）或定义的核苷酸（当模板分子有限时）。UID或UMI可以通过连接或通过PCR或逆转录过程中的引物引入靶模板。

用UID或双链条形码标记DNA片段已被证明可以减少错误并提高测序准确性，因为根据共享相同UID的一致读取结果，可以将真正的突变与PCR错误或测序错误区分开来。目前，经典的基于标签的方法有SafeSeq、CircleSeq和duplex测序。SafeSeq是一种基于“条形码”的单链标记方法。基于剪切点的单链标签的替代方法是circle sequencing，该方法利用Phi29 DNA聚合酶的链置换活性，在扩增前串联产生多个循环DNA分子副本。然而，这两种方法都无法区分真变异体和在最初几轮PCR扩增过程中引入的人工变异体。相比之下，双链测序通过标记dsDNA的两条链来解决这些类型的错误，利用DNA作为双链实体自然存在的事实，一个分子相互编码其补体的序列信息。表3比较了SafeSeq、CircleSeq和双工测序的声称错误率。

分子条形码测序数据的分析可分为三个步骤。

第一步是提取UID。需要注意的是，连接到原始DNA模板上的条形码通常是由DNA合成技术制成的，通常错误率很高。例如，如果设计了8-nt条码，由于合成错误，我们仍然有机会获得7-nt或9-bt条码。为了解决这个问题，通常使用由几个碱基（通常是三到五个碱基）组成的固定序列来表示UID和原始DNA序列的边界。分割算法应在设计位置附近寻找该标志，通常算法应允许一个碱基错配，以实现DNA合成或测序容错。通过使用特殊适配器，一些分子条形码方法将UID放置在多路复用索引位置（Illumina TrueSeq的I7或I5索引）。在这种情况下，UID提取要容易得多，因为它可以直接从样本索引中提取。这个过程是用FASTQ数据完成的。

第二步是对来自相同原始DNA的读取进行聚类。这些读取应该共享非常相似的UID和映射协调。但由于PCR和测序错误的存在，它们不需要完全相同。通常可以容忍一个基替换不匹配，松散的聚类方法可以允许INDEL或多个替换不匹配。此过程通常使用排序的BAM文件完成，但也可以使用基于序列聚类算法的FASTQ文件完成。

最后一步是为每个读取集群生成一致读取。首先，同一集群中的读取应该对齐在一起。这一过程可以通过像Clustal这样的多序列比对工具来完成。完整的多序列比对通常很耗时，如果我们限制不匹配替换和INDEL的数量，一些简单的方法可以运行得更快。校准完成后，可以通过从前到后扫描来生成一致性读数。对于每个职位，该职位中的所有基地将根据其质量分数投票给共识基地。对于基础完全相同的职位，该共识基础的质量分数可以调高一点，反之亦然，对于没有共识的职位，结果基础的质量得分可以调低一点。如果只有两个读取被聚集，如果相同位置的两个基数不同，但都有高质量分数，则可以用N或零质量分数掩盖此位置。

baseline methods

NGS数据有不同类型的错误。有些错误，如测序错误和PCR错误，是随机的，可以发生在任何基因组位置的任何核苷酸上，尽管有一些偏差。有些错误更为常见，比如由基因组高度重复区域的错位引起的错误。这些常规错误可以通过基线技术消除。

基线技术是从尽可能多的样本中组合和存储所有相关检测到的突变和其他相关信息，然后对这些数据进行统计，并为查询和更新提供接口。基线数据通常存储在数据库中，因此可以使用标准SQL语言进行插入、更新、删除和查询。可以使用两种不同类型的数据库：面向行的数据库和面向列的数据库。面向行的数据库是关系数据库的主流形式，如MySQL和PostgreSQL，而面向列的数据库则鲜为人知，如Infobright和MonetDB。面向行的数据库可以支持在线事务处理（OLTP），并针对关系查询进行了高度优化，而面向列的数据库可以提供更高的数据压缩率。

基线应该存储每个突变及其染色体、位置、参考和替代碱基，以及突变读取数和总深度。利用这个基线，我们可以计算检测到多少次具有特定变化的特定位置的突变，其平均MAF是多少，以及突变的读取数是多少。

由于在许多不同类型的癌症中可以检测到一些突变，一个更好的解决方案是用健康人的数据序列建立一个特定的基线。然后，该基线可用于筛选假阳性突变。当一个变异被调用时，它的基线重复数将被评估。若基线重复数太高，则该突变可被视为假阳性，需要仔细评估。

基线的另一个用途是检测热点突变，包括体细胞突变和种系突变。通过从与肿瘤个体建立的基线中挖掘热点突变，我们可以找到有潜力成为生物标记物的目标突变。

直接扫描FASTQ数据的目标变异检测

NGS数据的常规突变检测管道通常涉及多个工具的不同步骤。由于应用了不同的过滤器，这些工具可能会导致信息丢失，并可能最终导致漏检真突变，尤其是MAF低的突变。这种由数据分析引起的假阴性在临床应用中是不可接受的，因为它会使患者错过更好的治疗机会。

相反，也应避免这些关键突变的假阳性检测，因为它可能导致昂贵但无效的治疗，甚至可能导致严重的不良反应。常规的NGS管道可以检测到大量替换和INDEL，并不可避免地产生误报。特别是，由于比对物的参考基因组定位不准确，在基因组的高重复区域中调用的INDEL中有很大一部分是假阳性。

作者开发了一些工具，可以通过扫描原始FASTQ数据检测目标突变，而无需进行任何校准和变量调用。一个工具是MutScan，它基于容错字符串搜索算法，并通过滚动哈希和bloom过滤器对速度进行了高度优化。MutScan可以在无参考模式下运行，以检测程序中预定义的目标突变。通过提供VCF文件及其相应的参考FastA文件，MutScan可以扫描VCF中的所有变体，并通过为每个变体创建HTML文件来可视化它们。

MutScan超灵敏、超快速。它可以在只支持一个变异读取的情况下捕获变异。它能跑50×，如果只扫描预定义的癌症药物靶点，则比常规管道（AfterQC+BWA+Samtools+VarScan2）更快。此外，MutScan生成的交互式HTML报告有助于可视化和验证目标突变。图5显示了MutScan的交互式突变连锁反应。

MutScan位于：https://github.com/OpenGene/MutScan。它是用C++编写的，支持多线程。它支持单端和配对端数据，对于后者，它将尝试通过质量调整和纠错来合并每对数据。

作者开发的另一个工具是FusionDirect，它可以直接从原始FASTQ数据中检测基因融合。此工具还可以直接处理FASTQ文件，不需要对齐。它可以输出融合位点（基因和位置），以及支持融合的读数。图6给出了FusionDirect的输出示例。

FusionDirect需要一个包含四列（染色体、起始位置、结束位置、基因名称）的BED文件。如果未提供此文件，FusionDirect将使用内置BED文件，其中包含大多数具有高度临床重要性的融合基因。

FusionDirect的网址为：https://github.com/OpenGene/FusionDirect.jl。它是用Julia编写的，这是一种允许高性能技术计算的新语言。FusionDirect基于OpenGene Julia库构建(https://github.com/OpenGene/OpenGene.jl），提供常规NGS相关文件格式（即FASTQ/FastA/VCF）的基本序列和变量表示以及I/O功能。

去重和唯一支持读取计数

在确定被调用变量的可信度时，最重要的证据是其支持读取的数量和质量。为了计算支持读取的数量，我们需要识别并折叠重复的读取。

有一些工具可以消除PCR重复。Picard MarkDuplicates比较SAM/BAM文件中读取和读取对的五个主要位置中的序列。标记重复读取后，此工具使用一种算法来区分主要读取和重复读取，该算法通过读取的基本质量分数的总和对读取进行排序。然而，当该工具与某些野生型读取共享映射协调时，可能会导致不必要的肿瘤衍生变异读取的删除。

CAPP-seq引入了另一种方法。它会折叠那些具有完全相同序列的读取，但具有超低质量分数的读取除外。与Picard MarkDuplicates相比，此方法删除的读取更少，因此损失更少。然而，它通常受到排序错误的影响，因此处理的数据的重复级别仍然很高。

上面介绍的分子条形码测序是一种新的方法，似乎可以有效地消除PCR重复。由于UID连接是在任何扩增发生之前进行的，因此来自相同原始DNA的读取将共享相同的UID。基于UID和读取序列的聚类，可以检测到PCR重复，一致读取生成过程将删除重复的读取。表4比较了现有的重复数据消除工具。

上述方法在调用变体之前检测重复。另一种策略是在变量调用完成后检测重复，它将具有相同映射位置（开始和结束）的读取折叠为唯一读取，并为每个变异提供支持引用和替代基的读取数。这种独特的读取计数方法可以提供更精确的支持读取计算。应用此策略后，我们可以应用损失较小的重复数据消除方法，如CAPP-seq方法，为变量调用保留更多信息。如果变量调用方能够处理重复数据，我们甚至可以在变量调用之前跳过重复数据消除。

MrBam是为这种独特的读取计数任务而设计的工具。它区分了共享相同映射协调的单个读取或多个读取生成的结果读取。对于成对末端测序数据，它区分突变位于读对重叠或非重叠区域的情况。MrBam将给出以下条件组合的唯一读取数：支持引用或替代、由单个或多个读取聚集，以及位于重叠或非重叠区域。

MrBam的结果可用于筛选从ctDNA测序数据中调用的变体。根据我们的经验，要报告一个突变，我们需要至少两个支持它的唯一读取对，每一对应该在其重叠区域中有这个突变，或者是由多个对生成的一致对。由于高比率测序错误和ctDNA测序数据的极端深度，仅由非重叠区域的少数单次读取支持的突变通常为假阳性。

---

Cell-Free DNA的甲基化分析

甲基化变化在不同的癌症类型中很常见，通常发生在癌症发展的早期，通常抑制肿瘤抑制基因的表达。与突变相比，异常的DNA甲基化可能为血液中的肿瘤DNA提供一个更为一致的、因而广泛适用的标记。

有大量公开信息描述肿瘤组织中的DNA甲基化模式及其对患者预后的影响。当肿瘤DNA流入血流时，这些模式也可在血浆和血清中检测到。

肿瘤特异性ctDNA甲基化可用于量化肿瘤DNA，提供有关肿瘤负担水平的信息，并揭示肿瘤中的甲基化模式。基于DNA甲基化的生物标记物可以纳入患者护理和管理，但对临床实践的影响很小，例如甲基化ctDNA最近在确定癌症预后和在手术后或化疗期间的疾病监测中的应用。甲基化ctDNA检测也被开发出来，以满足癌症筛查所需的严格标准。

下一代测序平台允许以单个碱基分辨率构建DNA甲基化的基因组图。用亚硫酸氢钠脱氨基处理基因组DNA，使非甲基化胞嘧啶（C）转化为尿嘧啶（U），而甲基化C残基保持不变。在随后的聚合酶链反应（PCR）中，U最终转化为胸腺嘧啶（T）。全基因组亚硫酸氢盐测序（WGBS）和减少代表性亚硫酸氢酯测序（RRBS）是全基因组甲基化研究的两种经典方法。

WGBS（BS seq；MethylC seq）理论上涵盖了所有C信息。在这种方法中，基因组DNA被纯化并剪切成片段。碎片DNA进行端部修复；腺嘌呤碱基被添加到DNA片段的3' 端（A尾），甲基化适配器被连接到DNA片段。在亚硫酸氢钠处理和PCR扩增之前，对DNA片段进行大小选择，并对产生的文库进行测序。WGBS的主要优点是它能够评估几乎每个CpG位点的甲基化状态，包括低CpG密度区域，如基因间“基因沙漠”、部分甲基化域和远端调控元件。它还可以确定绝对DNA甲基化水平，并揭示甲基化序列的背景。

RRBS的开发成本低于WGBS，后者集成了Msp1限制酶消化、亚硫酸氢盐转化和下一代测序，用于分析特定片段的甲基化模式。研究发现，40到220 bps之间的MspI消化片段的大小选择覆盖了85%的CGI，大部分位于启动子中，启动子仅占哺乳动物基因组的1–3%，因此显著减少了测序数量。基于RRBS的协议比WGBS更具成本效益，因为这些方法侧重于富集靠近限制酶识别序列的富含CpG的区域。然而，这些方案可能在相对较少研究的基因间和远端调控元件方面缺乏覆盖。

最近，基于靶捕获的亚硫酸氢盐测序方法也已经开发出来，一些试剂盒，如NimbleGen SeqCap Epi已经商业化，以提供靶向甲基化分析。由于cfDNA中肿瘤DNA的比例较低，通常需要进行超深测序，因此进行目标捕获亚硫酸氢盐测序的能力对于分析ctDNA样本的甲基化信息非常重要。

ctDNA甲基化分析的主要应用之一是检测早期癌症。血浆中的循环甲基化SEPT9 DNA被开发为结直肠癌的生物标记物，SHP-1启动子2（SHP1P2）的甲基化被报道为非小细胞肺癌（NSCLC）的生物标记。这些生物标记物通常比蛋白质生物标记物（即癌胚抗原，CEA）更敏感，有潜力应用于癌症筛查或早期癌症检测。

ctDNA甲基化分析的另一个主要应用是确定未知原发癌（CUP）的组织起源。这种应用基于这样一个事实，即不同的人类组织和细胞具有不同的DNA甲基化模式。最近，一种鉴定甲基化单倍型区的方法被开发出来，用于从血浆DNA进行肿瘤组织起源定位。

分析亚硫酸氢盐测序（BS-seq）数据的生物信息学管道与分析正常测序数据不同。分析BS序列数据的关键步骤是质量控制、绘图、甲基化评分、差异甲基化评估等。

BS-seq数据的QA过程与正常测序数据的相同过程，包括质量分析、适配器修剪和低质量读取过滤。然而，要注意亚硫酸氢盐处理会导致T的过度表达和C的不足表达，这可能被传统的QC工具认为是有偏见的。因此，传统的QC工具，如FastQC，不是处理BS序列数据质量控制的好选择。BseQC和MethyQA是更好的选择，因为它们专门用于BS序列数据。

将BS-seq读数映射到参考基因组具有挑战性，因为序列与参考基因组不完全匹配，并且由于亚硫酸氢盐处理，文库复杂性降低。此外，每个给定的T都可能是真正的基因组T或转化的非甲基化C。由于这些原因，传统的比对工具如BWA和Bowtie不适合将BS序列读数映射到参考。一些BS-seq专用对齐器已经开发出来，通常可以分为两个通配符对齐器和三个字母对齐器。通配符如BSMAP通过将C替换为Y（胞嘧啶或胸腺嘧啶的IUPAC代码）来操作，而三字母比对符如Bismark在测序读取和参考中都将C转换为T。

一旦比对完成，就可以计算胞嘧啶或基因组区域的甲基化分数，以找到差异甲基化胞嘧啶（DMC）和差异甲基化区域（DMR）。胞嘧啶甲基化评分是通过聚集重叠读数并计算C或T的比例来计算的，这称为β评分。这个过程可以通过Bismark和GBSA等工具实现。像Methylkit这样的软件提供了一种将基因组划分为小仓的策略，平均β-得分作为bin得分。然后，可以应用Fisher精确检验（FET）等统计检验来评估样本之间DMC/DMR的统计相关性。这部分工作也可以用Methylkit完成，这是一个用于分析DNA甲基化的全面R包(https://code.google.com/p/methylkit).)。

最近发表了一些新的BS-seq数据甲基化分析方法。例如，Gao等人提出了一种搜索具有高度协调甲基化的基因组区域的方法。这种方法基于紧密耦合的CpG位点块，称为甲基化单倍型块（MHB）。然后可以在块水平（MHL）上进行甲基化分析，基于MHL分析的结果比基于分析单个CpG位点的结果要好得多，这意味着该方法可以用于鉴定来源组织。

亚硫酸氢盐测序作为分析DNA甲基化的黄金方法，已被研究多年，并开发了许多方法和工具。由于迫切需要建立用于癌症筛查和组织起源鉴定的甲基化分析，BS-seq数据分析将引起研究人员的更多关注。我们不能在本章中讨论BS seq的所有方面。可以在OMIC在线工具中找到BS seq数据分析工具和管道的集合(https://omictools.com/bsseq-category.)). （此链接无法访问）

---

机器学习方法

机器学习（ML）技术在许多领域中创建数据模型非常流行，它也可以应用于ctDNA数据分析。最适用的方法是监督学习方法，它基于标记数据的训练来构建分类器。在本小节中，我们将展示如何使用ML技术构建带有ctDNA测序数据的分类器。

一个ML应用是对cfDNA数据和非cfDNA数据进行分类。CfDNA具有一定的片段模式，可以带来测序数据开始周期的非随机碱基含量曲线。Chandrananda等人于2014年首次报道了cfDNA片段模式，其一个核苷酸分辨率。他们在cfDNA片段的两侧发现了一些高频率的10个核苷酸基序，而在切割位点的cfDNA的前两个碱基可以决定其他八个碱基中的大部分。他在2015年的进一步研究表明，这些断裂模式与染色体上的非随机生物分裂有关。DNA切割位点两侧的十个位置显示出一致的模式，其中特定核苷酸优先于核小体核心和连接区。图7显示了血浆cfDNA测序数据的碎片模式。

由于cfDNA的这种片段模式是稳定和独特的，它可以用于区分cfDNA的数据和其他种类的样品的数据。作者开发了一个名为CfdnaPattern的开源工具，用于训练分类器，如SVM、KNN或随机森林，以预测FASTQ是否从cfDNA中测序。使用0.632 + bootstrapping与3000多个FASTQ文件进行交叉验证，结果平均准确率为99.8%，使用随机森林、线性SVM或KNN分类器获得。这个工具是用Python编写的，使用了广泛使用的Python机器学习包scikit-learn。这个工具可从https://github.com/OpenGene/CfdnaPattern访问。

另一个ML应用是预测突变是体细胞还是种系。通常，肿瘤和正常样本都被测序，正常样本可以用作参考，以确定肿瘤样本中称为种系或体细胞突变的突变。但在某些情况下，我们可能没有与肿瘤样本匹配的正常样本，然后我们可以应用ML方法根据支持参考文献和突变的读数对突变进行分类。

DeepSomatic是提供此类功能的工具。它可以用深度神经网络对体细胞和种系突变进行分类。如果读取数大于256，则提取覆盖突变的所有读取并采样到256个读取。然后将突变位点周围的这些读取碱基编码为2D图像，每个像素包含以下通道：读取碱基及其质量分数、参考碱基以及插入或删除的长度。然后，用五个常规层构建深度卷积神经网络（CNN）。利用肿瘤正常配对数据对模型进行训练和验证，然后交叉验证评估表明，该模型的平均准确率高于99.9%。DeepSomatic也是一个开源工具，可在https://github.com/OpenGene/DeepSomatic访问。

---

数据模拟

调整生物信息学管道和训练软件参数需要具有已知地面真相的测序数据，而这实际上很难从真实测序数据中获得。特别是，对于ctDNA测序应用，其目的是从超深度测序数据中检测低频变异，很难区分所谓的变异是真阳性还是假阳性，这是由测序或其他过程的错误引起的。在这些情况下，具有配置变化的模拟数据可用于诊断和验证生物信息学程序。

尽管已经开发了许多下一代测序模拟器，但大多数模拟器缺乏模拟一些实际特征的能力，例如目标捕获测序、拷贝数变化、基因融合、扩增偏差和测序错误。作者开发了SeqMaker，这是一种现代NGS模拟器，能够模拟不同类型的变化，并集成了放大偏差和测序误差。目标捕获测序只需使用捕获面板描述文件即可支持，测序错误率、平均重复水平、DNA模板长度分布和质量分布等其他特性可以通过简单的JSON格式配置文件轻松配置。通过整合测序误差和扩增偏差，SeqMaker能够模拟更真实的下一代测序数据。可配置的变体和捕获区域使SeqMaker非常有助于生成数据，以训练生物信息学管道，用于体细胞突变调用等应用。表5比较了SeqMaker和其他NGS模拟器的特性。

SeqMaker是一种工具，可生成SNV、INDEL、CNV和基因融合的测序读数，并集成测序误差和PCR偏差。该工具使用JSON格式的配置文件来描述测序模拟设置，使用BED格式的TSV文件来配置捕获的目标区域。首先，模拟器从整个基因组或面板文件配置的目标区域中提取DNA片段，并在此过程中模拟CNV。其次，根据配置文件中配置的变异列表，DNA片段将被改变以模拟SNV、INDEL和基因融合。第三，将在每个DNA片段上模拟测序过程以生成NGS读数，测序误差和扩增偏差也将在该过程中模拟。最后，生成的读取被写入FASTQ文件。

SeqMaker是用Julia编写的，源代码可从GitHub获得：https://github.com/OpenGene/SeqMaker.jl/。目前，它仅支持Illumina平台。需要更多的努力来构建其他平台的模拟器，尤其是PacBio和Nanopore平台等新一代测序仪。

---

讨论

作为癌症领域的一种创新方法，液体活检在癌症诊断、监测和筛查中具有当前或潜在的应用。无细胞肿瘤DNA作为液体活检的主要成分，已广泛用于肿瘤患者的个性化药物指导。对于那些不适合通过手术或针头穿刺获取组织样本的患者，ctDNA测序为他们诊断肿瘤提供了新的机会。

由于ctDNA应该进行非常深入的测序，因此通常会考虑到成本，使用小基因板进行靶捕获。然而，小面板有一些缺点。小面板不允许检测目标区域外的突变，难以检测大规模拷贝数变化，并且难以计算总突变负担（TMB），这通常需要大面板或整个外显子组测序。随着测序成本的下降，不难推测整个外显子组甚至全基因组深度测序将变得负担得起，并被更广泛地用于ctDNA测序。然后将获得非常大的测序数据，对这些数据的数据处理和分析将非常具有挑战性。

---

结论

在本章中，我们介绍了ctDNA的概念和应用，解释了分析ctDNA NGS数据的困难，回顾了一些相关工具，并提出了一些新的方法或工具。人们应该认识到，cfDNA中的体细胞突变通常具有非常低的MAF，因为肿瘤特异性DNA片段通常是整个cfDNA的一小部分。人们还应该意识到，在实验和测序步骤中可能会发生错误，软件也可能会引入错位或假阳性变量调用等伪影。

---

未来的工作

尽管我们已经讨论了用于ctDNA NGS数据分析的生物信息学的许多方面，但仍存在上文未讨论的主题。

数据压缩是我们在本章中没有讨论的一个关键主题。由于ctDNA通常需要超深度测序，它通常产生非常大的数据。想象一下，如果10000？如果应用WES，我们将为单个样本获得超过500 Gb的数据，从而得到大于1 TB的未压缩原始文件。存储或传输这样大的文件将非常具有挑战性，并且迫切需要提供高压缩比的方法。从信号处理的角度来看，ctDNA测序数据是高度冗余的，因为它非常深，并且有可能被高比率压缩。然而，由于三个原因，压缩此类数据仍然不容易：排序错误导致的不一致读取、不同的质量分数以及无损压缩的要求。与通用压缩器（如gzip和bzip2）相比，当前的方法（如DSRC）表现出了更好的性能，但压缩比的改善仍不令人满意。一些新型压缩机，如gtz(https://github.com/Genetalks/gtz）已经开发出来，但它们仍然没有针对深度测序数据进行优化。我们认为，完美的深度测序数据压缩器应该实现局部从头组装或应用基于参考的策略来实现更高的压缩比。

另一个有待讨论的话题是CNV检测。由于肿瘤特异性DNA只是cfDNA的一小部分，肿瘤细胞中的拷贝数变化只会导致ctDNA测序数据中的微小拷贝数差异。例如，如果肿瘤特异性DNA占整个cfDNA的1%，并且肿瘤细胞中的拷贝数倍数为5，那么整个cfDNA数据中的拷贝数量将为104%，略高于平均水平。目前的CNV检测器，如CNVkit，并不是用来处理ctDNA测序数据的，也不够灵敏，无法检测CNV中的这种细微变化。更好的CNV检测器仍有待开发，这将为深度和目标捕获的ctDNA测序数据提供更好的标准化。

最近，一些针对癌症免疫学的新方法正在吸引人们。一个主题是预测癌症免疫治疗的结果，特别是PD-1/PD-L1检查点抑制剂。肿瘤突变负荷（TMB）已被证明与癌症免疫疗法的反应有关。然而，TMB通常使用组织全外显子组测序数据计算，并且由于低MAF和高噪声水平，使用ctDNA计算TMB仍然具有挑战性。需要针对基于ctDNA的TMB计算进行优化的方法，这一主题可以在未来讨论。另一个与癌症免疫治疗相关的话题是新抗原的发现。2016年12月，帕克癌症免疫治疗研究所和其他机构宣布成立肿瘤新抗原选择联盟。这个联盟包括来自30个非营利机构的研究人员，旨在确定能够最好地预测患者肿瘤DNA新抗原的软件。目前，计算预测能够在患者中引发有效抗肿瘤反应的新抗原仍然是一个偶然事件。从患者的ctDNA进行同样的预测甚至更具挑战性。新抗原预测研究将是学术界和工业界的热门话题，其进展和结果将在未来讨论。

---

end~