5 系统的测试与扩展

5.1 MRUnit测试类编写

MRUnit是一个MapReduce的测试库，它能将已知的输入通过函数的形式直接传递给Map()函数，然后该函数直接被调用运行，或者检查Reduce()函数的输出和测试的指定输出匹配，看是否符合预期输出。MRUnit与标准的测试执行框架（如：JUnit）一起使用^[1]，因此可以在正常的开发环境中运行MapReduce作业的测试。

5.1.1 Map任务的测试类编写与调度

这里通过对Mapper函数传入一个样本名称，如SRR3226035，测试其程序是否调用流程分析脚本，并输出正常结果。

public class wgsMapperTest {
    /*
        测试Mapper
     */
    @Test
    public void wgs() throws IOException {
        Text value = new Text(“SRR3226035”);
        new MapDriver<LongWritable, Text, Text, Text>()
                .withMapper(new wgsMapper())
                .withInput(new LongWritable(0), value)
                .withOutput(new Text(“1”), value)
                .runTest();
}
}

5.1.2 Reduce任务的测试类编写与调度

这里通过对Reducer函数传入多个样本名称，如SRR3226035，测试其程序是否调用shell脚本，整合多个gvcf变异文件，并整合为一个VCF文件。

public class wgsReducerTest {
    /*
        测试Reducer
     */
    @Test
    public void merge() throws IOException {
        Text value1 = new Text("SRR3226035");
        Text value2 = new Text("SRR3226039");
        Text value3 = new Text("SRR3226042");
        new ReduceDriver<Text, Text, Text, Text>()
                .withReducer(new wgsReducer())
                .withInputKey(new Text("1"))
                .withInputValues(Arrays.asList(value1,value2,value3))
                .withOutput(new Text("1"), value3)
                .runTest();
    }
}

5.2 Hadoop测序平台的测试

人类的参考基因组数据文件很大，其序列长度大概在3GB，包含30亿个碱基信息。而一个人的高深度测序所需要准备的样本数据往往是这个数字的30倍多，大约达到100GB。如果直接用这样的数据来完成基因测序分析，会下载大量的数据以及花费大量的处理时间，而且对机器性能也要有很高的要求。

因此在保证WGS全基因组测序流程的的完整性和正确性的情况下，使用一种大肠杆菌E.coli K12的全基因组数据作为代替。这是美国国家生物技术信息中心的一个研究课题，主要是针对一种大肠杆菌的全基因组测序的实验项目。该E.coli K12生物的基因组比较单一，全基因组中仅包含NC_000913.3这一种染色体，数据量相比于人类来说算很小的，它的全基因组的文本数据总大小只有4.6Mb，非常适合用来测试基于Hadoop的全基因组测序分析平台。以后对人类基因组数据作分析时再作替换即可。

5.2.1 测试环境与测试数据准备

在测试环境的准备中，在上一章中着重介绍了Hadoop伪分布式的搭建和Hadoop分布式的搭建流程，由于机器的限制，就本次的测试选择Hadoop伪分布环境下进行，我可以在我自己学习用的笔记本上直接搭建好该环境，接着对测试样本数据进行下载和预处理。

1）处理E.coli K12的参考基因组序列
物种的基因组参考序列可以在NCBI上获取，下面是E.coli K12参考序列的ftp地址。为了接下来表达上的清晰和操作上的方便，通过bgzip将这个序列文件进行解压并重命名为E.coli_K12_MG1655.fa
$ gzip -dc GCF_000005845.2_ASM584v2_genomic.fna.gz > E.coli_K12_MG1655.fa
接着，使用Samtools软件为参考文件建立索引，方便其他数据分析工具（比如GATK）能够快速地获取fasta上的任何位点的序列信息。
$ /home/elon/biosoft/samtools/1.0/bin/samtools faidx E.coli_K12_MG1655.fa

2）下载E.coli K12的测序数据
基因组参考序列准备好之后，我们还需要准备该物种的测序样本数据。这里使用一组科研的E.coli K12实验数据E.coli K12测试数据 ，其作为一种供研究使用的模式生物，已经在NCBI上有许多的测序数据了。我这里选取某次科研组织的样本数据共9组，样本号序列分别为：SRR3226034、SRR3226035 … SRR3226042。

通过点击下载数据样本，并把最后SRR_number替换成指定的样本号。这些样本数据来自Illumina MiSeq测序平台，read长度是300bp，测序类型Pair-End，说明样 本有R1和R2双端测序数据。

从NCBI上下载下来的测序数据是一种NCBI的特殊格式SRA，它具较高的压缩率。我们需要对下载下来的SRA数据进行转换，得到我们需要的FASTQ格式数据，需要用到NCBI的官方工具包sratoolkit^[2]，解压之后就可以使用了。

sratoolkit是一个工具包，所有的执行程序都在它解压后的bin文件夹下，我们要把SRA转换为fastq，只需直接通过工具包中的fastq-dump即可完成。
$ /home/elon/biosoft/sratoolkit/2.9.0/bin/fastq-dump --split-files SRR_number.sra

然后我们就会得到这个E.coli K12数据的read1和read2了：SRR_number_1.fastq和SRR_number_2.fastq。在得到R1和R2序列对之后，可以用bgzip命令将其压缩为.gz文件，这样可以节省空间，并且不会对接下来的数据分析产生影响。

$ bgzip -f SRR_number_1.fastq
$ bgzip -f SRR_number_2.fastq

在得到gz格式的FASTQ样本数据之后，接下来就可以进行下面具体的测序分析流程了。

5.2.2 系统平台的测试

在本机的伪分布式的Hadoop环境中，用IDEA开发平台进行测试。首先通过以下指令运行Hadoop集群：
$ start-dfs.sh && start-yarn.sh

编辑Mapper阶段的输入文本sample.txt，每行指定一个样本数据：

SRR3226034
SRR3226035
SRR3226036
SRR3226037
SRR3226038
SRR3226039
SRR3226040
SRR3226041
SRR3226042

接下来，需要通过以下命令将样本数据上传到HDFS上的指定文件夹：
$ hadoop fs -put sample.txt /wgsv2-input/

最后，调用wgsDriver驱动程序，这是整个项目的入口程序。从这里开始整个项目的MapReduce分析流程，并在指定文件夹中打印日志和脚本信息，最终得到VCF变异检测文件，并已上传到HDFS的指定位置。

对gVCF文件和VCF文件的解读在第二章的相关文件格式介绍中已经提及，其中的样本数据就是来自这一组测试数据得到的最终变异集合。

5.3 测序平台的分析与优化

5.3.1 测序平台与传统测序流程的比较

与传统的处理流程相比，对于这九组样本数据，在传统的分析流程中，需要手动编写九组测序的脚本，这个重复性的工作在FreeMarker的引入之后可以简化了。只需编写一套测序分析模板脚本，通过模版生成引擎，得到九组不同的测序分析脚本。

Map阶段测序分析运行时间比对如图5-1所示，横坐标表示的是用于测序的样本个数，纵坐标表示的是测序流程总的运行时间。在传统的测序分析是进行串行处理的方式，在处理单个样本对数据时，耗时大概12分钟左右，因此对于样本个数逐个递增的情况其总耗时时间呈线性递增；而对于分布式环境的Map任务执行，各个map任务在集群中是同时提交并行执行的，在每一种测序样本数处理情况中，以处理用时最长那一组的测序时间作为总运行时间。
因此由图5-1可知，随着样本数的增多，传统的测序方法的运行时间呈线性增加，而基于Hadoop的分布式测序所需要的时间基本保持不变。

Reduce阶段测序分析运行时间比对如图5-2所示，在Reducer规约阶段需要将上一步骤中各个样本产生的单个变异检测样本合并，而这里都是经由一个脚本程序处理，因此传统基因测序时间和分布式基因测序时间相近。

5.3.2 增加测序流程处理的时间戳标记

在每个分析流程执行的时候，在各个Shell命令中加入time命令，可以在日志中打印该流程所花费的时间。比如：程序执行的开始时间和程序执行的结束时间，一共花费了多少时间等，这样可以对分析流程中各个程序的执行时间和执行流程顺序更加了解。

5.3.3 引入Log4j日志框架

在整个项目中引入Log4j日志框架，将不同的日志类型的信息导向不同的文件中存储，便于项目的开发和调试。下面是我的log4j的配置文件：

### Log4j建议只使用四个级别，优先级从高到低分别是ERROR、WARN、INFO、DEBUG。###
log4j.rootLogger = stdout,D,E,I

### 输出信息到控制台 ###
log4j.appender.stdout = org.apache.log4j.ConsoleAppender
log4j.appender.stdout.Target = System.out
log4j.appender.stdout.layout = org.apache.log4j.PatternLayout
log4j.appender.stdout.layout.ConversionPattern = [%-5p] %d{yyyy-MM-dd HH:mm:ss,SSS} method:%l%n%m%n

### 输出DEBUG 级别以上的日志到 $PROJECT_HOME/wgsv2/wgs-logs/debug.log ###
log4j.appender.D = org.apache.log4j.DailyRollingFileAppender
log4j.appender.D.File = ./wgs-logs/debug.log
log4j.appender.D.Append = true
log4j.appender.D.Threshold = DEBUG
log4j.appender.D.layout = org.apache.log4j.PatternLayout
log4j.appender.D.layout.ConversionPattern = %-d{yyyy-MM-dd HH:mm:ss}  [ %t:%r ] - [ %p ]  %m%n

### 输出ERROR 级别以上的日志到 $PROJECT_HOME/wgsv2/wgs-logs/error.log ###
log4j.appender.E = org.apache.log4j.DailyRollingFileAppender
log4j.appender.E.File =./wgs-logs/error.log
log4j.appender.E.Append = true
log4j.appender.E.Threshold = ERROR
log4j.appender.E.layout = org.apache.log4j.PatternLayout
log4j.appender.E.layout.ConversionPattern = %-d{yyyy-MM-dd HH:mm:ss}  [ %t:%r ] - [ %p ]  %m%n

### 输出INFO 级别以上的日志到=$PROJECT_HOME/wgsv2/wgs-logs/TemplateEngine-info.log ###
log4j.appender.I = org.apache.log4j.DailyRollingFileAppender
log4j.appender.I.File =./wgs-logs/info.log
log4j.appender.I.Append = true
log4j.appender.I.Threshold = INFO
log4j.appender.I.layout = org.apache.log4j.PatternLayout
log4j.appender.I.layout.ConversionPattern = %-d{yyyy-MM-dd HH:mm:ss}  [ %t:%r ] - [ %p ]  %m%n

从配置文件中可以看到，日志优先级共四种，优先级从高到低分别是ERROR、WARN、INFO、DEBUG，对于不同日志等级的日志文件进行分文件存储。这样对于后期的系统调试至关重要。

5.4 基于Hadoop基因测序平台的扩展

在对不同物种进行基于Hadoop平台的全基因组测序流程开发时，只需对WGS处理模板进行重构，适配特定的物种即可，在本章的测试样本中，我使用的是大肠杆菌的K12生物来进行测序流程构建，而将其换作是其他生物的全基因组数据完全可以。下图是测序平台的扩展性示意图，由于每个样本的测序分析流程都是封装在Map任务和Reduce任务中的，因此可以同时并行化的开展对不同物种的测序流程，Hadoop大数据测序平台扩展如图5-3所示。

在该Hadoop基因测序平台中，对于其他物种如人类，只需将人类的全基因组测序样本上传至HDFS分布式存储系统中，将人类的参考基因组数据在集群的各个主机上进行分布，再设计人类基因测序流程的模板文件即可适配该平台。因为各个Map任务和Reduce任务是独立存在，并隔离在不同的container中执行的，因此各个子任务之间互不干扰，达到完全并行化计算处理。
##5.5 本章小结
本章主要对搭建的基于Hadoop的WGS分析流程平台的一个测试，在保证WGS全基因组测序流程的的完整性和正确性的情况下，用 E.coli K12生物的数据作为代替，测试平台的可用性和健壮性等。对传统WGS分析流程与Hadoop大数据框架整合，从测序源到产生最终的变异结果gVCF和VCF文件的过程进行测试。
在这次测试中，通过观察运行过程的日志信息，发现尝试分析最耗时间的步骤基本聚集在比对和变异检测这两步中，最后对于测序平台的扩展性进行论述，说明其如何适配不同物种的基因测序流程，并进行并行化计算等。

---

1. Apache. Apache MRUNIT™[EB/OL]. http://mrunit.apache.org/, 2018. ↩

2. The sratoolkit Toolkit [EB/OL]. https://github.com/ncbi/sra-tools/wiki/Downloads, 2018. ↩