基因组分析工具包：Apache Spark
基因组分析工具包：Apache Spark
自2000年人类基因组计划（Human Genome Project）产生人类基因组首份草案序列以来，测序成本从几乎每个基因组的1亿美元左右急剧下降到今天的约1,000美元。在同一时期，我们看到Apache Hadoop等大数据技术的存储和处理能力大幅增长。因此，使用Hadoop生态系统中的工具进行基因组学分析就水到渠成，Cloudera与Broad Institute及其他行业合作伙伴就借着这股东风，发布了他们运行在Apache Spark上的第4版基因组学分析工具套装（Genome Analysis Toolkit，GATK）的alpha版本。
基因组数据处理就像任何行业的数据处理一样，都是关于运行数据管线的。图1显示了DNA测序的管线。
图1. DNA测序流程
流水线从DNA样本开始；由一台机器排序，排序后导出一个包含DNA序列片段（由字母A，C，G和T组成）的文件。原始序列数据并不是非常有用，因为此时并未包含序列片段在基因组中的位置信息。因此，需要使用一款称为对齐器的软件将待测序列与参考基因组序列进行比对（该参考基因组序列是人类基因组计划的产物）。在这个阶段，个体的DNA序列是已知的，但随之而来的问题是：“这个待测序列与其他个体有什么不同？”流水线上的这一步等同于寻找待测序列中存在的与参考序列不同的变体：在概念上，您可以将其视为在待测序列和参考序列上运行Unix diff命令。输出是针对个体的一组变体识别。 图1中的流水线在此处停止，但实际上，变体识别数据是研究人员下游分析的原材料。
基因组分析工具包（GATK）涵盖了流水线的变体发掘部分。变体发掘本身包含许多步骤，而GATK提供了用于运行这些步骤的工具（这些步骤在GATK最佳实践文档中进行了描述）。GATK版本4中的新增功能提供了这些工具的Spark实现，这将允许在Spark集群上按比例运行流水线上的变体发现部分。为了理解所涉及的内容，让我们更详细地看看GATK。
更详细的了解GATK
GATK具有一种称为Walker的概念，本质上它是一个回调迭代器API。例如，ReadWalker遍历基因组中的所有读取，并在每个读取序列的回调接口上调用一个方法（一个读取序列是已与参考基因组进行比对的简短DNA片段。）。Walker API最吸引人的地方在于它使用户可以轻松编写自己的工具。最简单的程序之一是计算文件中的读取序列的数量，如下所示：
publicfinalclassCountReadsextendsReadWalker{
privatelongcount=0;
@Override
publicvoidapply(GATKRead read,ReferenceContext referenceContext,FeatureContext featureContext){
++count;
@Override
publicObjectonTraversalDone(){
returncount;apply()方法简单地增加一个计数器; 那么一旦达到文件结尾onTraversalDone()就被调用来报告计数。CountReads是一个串行程序，因此只有一个进程用于对文件（甚至是多个文件）中的读取序列进行计数。它的优点就是简单，但是缺点也很明显，就是速度慢，在处理TB量级的数据时这一点尤为明显。
Spark等价物（GATK）的结构是不同的，因为代替回调函数，工具实现它们自己的输入处理。下面是CountReads的一个分布式版本（为了简单起见，在GATK源代码中已经从原始代码中删除了一些样板代码）。
publicfinalclassCountReadsSparkextendsGATKSparkTool{
@Override
protectedvoidrunTool(finalJavaSparkContext ctx){
finalJavaRDD reads=getReads();
finallongcount=reads.count();
System.out.println(count);在这种情况下，读取序列将作为Spark RDD（JavaRDD<GATKRead>count()）进行访问，然后调用内置函数count()触发Spark作业来计算RDD中的条目数。输入被分成多个部分（默认情况下，每个部分的大小均为128MB），并且Spark作业为每个并行分割运行一个任务。使用这种方法，之前使用Walker版本运行需要花费数小时的作业仅需要几分钟内就可完成，即便是具有少量节点的适度Spark群集。
计数读取是非常枯燥琐碎的，其并行处理也非常鸡肋。GATK中的大多数工具都比较复杂，不能轻易移植到Spark上。让我们更详细地看看其中的一个：Mark Duplicates算法。
在Spark中标记重复序列
测序过程本身就是一个嘈杂的过程，而且经常发生相同的DNA片段多次测序，产生重复读取序列。所以需要删除这些重复项目以减少不必要的额外工作。Mark Duplicates算法的作用就是查找并标记这些相同的序列。
如何判断两个（或更多）读取序列是否重复？读取序列由几个区域组成对齐信息：如果这些区域对于一组读取序列是相同的，则它们被认为是重复的。重复读取序列是不相同的，因此算法根据读取的其他方面（如质量测量）对每个重复读取序列进行评分，具有最高分数的读取序列保持不变，其他的片段被标记为重复。
Mark Duplicates算法归结为通过对齐信息字段对所有读取序列进行排序。使用Walker API来实现这一点很困难，因为这些读取序列并不全都适合内存。该实现使用自定义的磁盘分类算法，但即使这样也只能利用一台机器。
借助Spark，我们可以访问高度调整的分布式排序实现，并且这隐藏在一个漂亮的Java API之中。让我们看看Mark Duplicates实现的核心部分。我们从由读取分组和名字分组的读取序列开始（文件通常已经按照这种方式排序，但如果没有，则需要进行初始排序）。
JavaPairRDD<String,Iterable>keyedReads=...;接下来，我们将每次读取的对齐信息字段提取到一个字符串中，并为该值构建PairedEnds对象。读数通常是成对的，一对中的每个成员来自DNA片段的任一末端进行测序。一个PairedEnds对象只是一对读取的包装。
JavaPairRDD<String,PairedEnds>keyPairs=
keyedReads.flatMapToPair(keyedRead->{...});然后我们通过操作来执行一个组，以将具有相同对齐信息的所有读取（在PairedEnds对象中）放在一起。这些组就是重复的集合。
JavaPairRDD<String,Iterable>keyedPairs=
keyPairs.groupByKey(numReducers);最后，对于每个组，可以进行评分和标记过程：
JavaRDD markedDups=markPairedEnds(keyedPairs);用Spark创作工具
Mark Duplicates工具将读取的最终RDD写入输出文件，以供下游其他GATK工具进一步处理。由于输出是RDD，因此另一种选择是在单个Spark作业中组合工具，以便中间步骤不需要在文件系统上实现。
在最新的GATK4 alpha版本中，并非所有工具都已移植到Spark中，因此还无法将整个测序流水线作为单个Spark作业运行。然而，这将在未来发生变化，目标是在Spark集群上运行对齐器（如BWA-MEM）和GATK变体发掘管线。
由于许多项目采用Spark进行大规模计算，因此将工具集成到基因组管线中的可能性正在增加。ADAM是第一个将Spark作为基因组学平台的项目，该项目还使用Apache Parquet为基因组数据定义文件格式。作为选项，GATK4可以读取和写入ADAM Parquet格式化数据。还有一个有趣的可能性是混合和匹配管道中的工具。例如，您可以使用来自ADAM的变体识别，或来自Hammer Lab（纽约西奈山医院的Jeff Hammerbacher实验室）的变体识别。
展望未来，另一个目标是将由管道生成的变体调用数据集加载到Hive（如Parquet格式）或Apache Kudu（孵化）等Hadoop本机存储引擎中。这将使用户能够利用Apache Impala（孵化），Ibis或Spark等工具进行分析，甚至可以为这些框架之上的科学家构建工具。
参与开源基因组学
GATK4还很年轻，但其接受程度已经非常令人鼓舞。作为其潜力的一个例子，Broad研究所的一个研究小组正在研究拷贝数变异，这是基因组中复杂的变体，其中一段DNA被重复。他们在GATK3上编写一个Spark工具，由于其计算复杂性，它在GATK3上没有尝试过，按照他们的估计，其运行速度比它运行在GATK3上快一到两个数量级。
针对新的Spark实现进行算法优化和运行时调优的小组数量令人印象深刻：它们包括来自Cloudera，Cray，Google，IBM和Intel的开发人员。由于它的工作负载很大，基因组学对于尖端的硬件和软件来说是一个很好的测试平台，所以它已经被用来发现错误（例如在Java 8上运行时的Kryo中的这个错误）也许并不奇怪。
Spark正在履行其作为普通分布式计算结构的承诺，该结构既可以在云端也可以在本地运行。我们在Cloudera希望其他开发者能够参与像Spark这样基于GATK的项目。要了解更多关于GATK，请访问GitHub的仓库https://github.com/broadinstitute/ GATK。
腾讯云开发者
