CN104375807B - 基于众核协处理器的三级流水序列比对方法 - Google Patents

Abstract

本发明公开了一种基于众核协处理器的三级流水序列比对方法，目的是提高序列比对软件的比对速度。技术方案是：采用MIC众核协处理器使用多线程并行进行序列比对，并将MIC上序列比对过程中从主存读取序列到MIC、比对序列和将比对结果写回主存这三个串行步骤采用三级流水方式，即在进行每一轮序列比对的同时，读取下一组比对所需要的序列，同时将上一组比对的结果写回主存，使得读写操作和比对操作并行进行。本发明实现了序列读取，序列比对和比对结果写回的三个主要过程的并行执行，提高了比对效率，降低了比对时间，与两路八核CPU相比，本发明可加速比对过程2.3倍以上，且避免拷贝大量内存空间，提高了程序的时空效率。

Description

基于众核协处理器的三级流水序列比对方法

技术领域

本发明涉及生物信息领域序列比对的方法，尤其指一种基于众核协处理器的序列比对方法。

背景技术

分子生物学是从分子水平上研究生命现象物质基础的学科，通过研究生物分子的结构、功能和合成等方面的原理，从而使生物体的功能和性状在前所未有的分子细节上得到详尽的分析和理解，进而更加科学严谨地阐明生命现象的本质。

在分子生物学研究中，DNA的序列分析是进一步研究和改造目的基因的基础。DNA(脱氧核糖核酸)是一种生物大分子，一共分为四种碱基，记为A、T、C、G，这些大分子的排列顺序决定了某种遗传指令，这些遗传指令是建构细胞内其他的化合物，如蛋白质与核糖核酸的需要。带有蛋白质编码的DNA片段称为基因，即遗传物质，是DNA分子上具有遗传信息的特定核苷酸序列。基因经过转录、翻译，最终产生结构和功能各异的、表现生物体性状的蛋白质。

DNA序列分析的基础是对DNA分子进行测序，即确定DNA分子中A、T、C、G四种碱基的排列顺序。当前的DNA测序技术，一次实验最多只能直接测得不大于5000个碱基的排列顺序，形成多个DNA短序列(称为read)。而一般生物的基因组碱基数目都十分巨大，如人类基因组总长约为30亿个碱基对。这样，绝大多数生物的基因组都不能通过实验手段一次性获得，而必须借助于计算机技术进行后续拼接得到完整的基因组。

序列比对是目前广泛使用的DNA序列分析方法，它是将测序得到的read短序列直接与拼接完成的参考基因组进行比对，确定read在参考基因组中是否出现以及出现的具体位置。通过序列比对进行DNA序列分析，避免了对目标基因组进行组装，可以很大程度上节省序列分析的时间和工作量，提高序列分析的效率。

由于比对时read数量都较大，无法一次性全部存放到计算机主存中。所以目前常用的DNA序列比对方法均按照以下步骤进行：

步骤1：根据计算机主存可用空间大小，将read平均分为若干组，每一组所占空间大小不超过计算机主存容量；

步骤2：从磁盘上读取一组read到主存中；

步骤3：对读取到主存中的read逐个进行比对；

步骤4：将read比对结果写回磁盘；

步骤5：检查磁盘中是否还存在未比对的序列，如果存在，返回步骤2；如果不存在，结束比对过程。

目前比对主要使用的运算器件为计算机中央处理器CPU或者图形处理器GPU。

虽然与进行序列组装相比，序列比对可以节省大量时间，但是目前广泛使用的基于CPU的串行序列比对方法的速度仍比较慢，如在配备两路八核Intel 2.4GHz CPU的常用服务器上，采用李恒在论文《Fast and accurate short read alignment with Burrows-Wheeler Transform》中公布的基于BW(Burrows-Wheeler)变换的序列比对方法，对长度为100个碱基的8千万条序列进行比对，需要花费一天以上时间，很难满足后序的序列分析对于时间的要求，更是无法满足时效性要求较高的临床需求。

基于CPU的并行序列比对方法使得多个线程能够并行进行序列比对，有效地提高了序列比对的速度。但是目前绝大部分研究机构使用的是单节点服务器，CPU计算能力十分有限。而随着测序技术的发展，特别是新一代高通量测序技术的出现，单位时间内产生的read数量翻了几翻，基于CPU的并行序列比对软件也已经很难有效处理如此大量的read。

基于GPU的序列比对软件，利用GPU具有大量计算核心的特性，使用其对序列比对进行加速，有效地增强了服务器的计算能力，与基于CPU的并行序列比对方法相比，进一步提高了速度。GPU指令集设计与CPU相比较为简单，但是能够快速处理简单的浮点和整型计算。而目前广泛使用的李恒在论文《Fast and accurate short read alignment withBurrows-Wheeler Transform》中公布的基于BW(Burrows-Wheeler)变换的序列比对方法运算过程复杂，程序分支多，当GPU中的一个核心遇到分支时，与其同组的其他核心均要等待该分支处理完毕才能继续并行执行，很大程度上降低了序列比对的效率。

MIC(Many Integrated Core)是Intel公司开发的一款众核协处理器，具有与传统x86CPU兼容的指令集，能够快速处理分支判断等复杂指令。而且，每个MIC协处理器配备了50个以上的计算核心，每个计算核心可以启动4个硬件线程，并行规模可达200个线程以上，其核心主频约为1.1GHz，并且包含512位宽度的向量处理单元，配备6GB以上片上存储空间，其单卡双精度浮点计算峰值超过1TFlops。MIC与GPU相比，在峰值计算能力相当的情况下，能够有效地处理复杂指令，提高复杂程序的运行效率。因此，本发明采用MIC众核协处理器加速序列比对。

但是MIC仍然存在一定局限，其计算核心只能访问它自身的片上存储空间，无法直接访问计算机主板上的主存空间，只能由CPU来负责主机与MIC之间的数据传输。一种解决方案是在运行前将程序和所需的数据事先拷贝到MIC上，然后登录上MIC上的操作系统，启动程序，程序的输出写在MIC上，等运行完毕后再将结果手动拷贝到主存。以人类基因组为例，序列比对过程的输入、输出，以及中间文件所占空间超过10GB，但是目前MIC卡上的存储空间无法满足需求，所以该方案困难非常大。另一种解决方案是从CPU端启动程序，在程序运程过程中由CPU将序列比对的输入数据拷贝到MIC上，MIC完成比对后，再由CPU将比对结果拷贝到主存，并将下一批序列拷贝到MIC上，如此循环，直到所有序列比对完毕。但是CPU与MIC之间频繁的数据传输占用了大量的程序运行时间，很大程度上影响了程序的运行效率。如何提高计算效率，降低数据传输开销是利用MIC加速序列比对的难点。目前还没有利用MIC进行序列比对的技术方案的公开报导。

发明内容

本发明要解决的技术问题是提出一种基于MIC众核协处理器的三级流水序列比对方法，提高序列比对软件的比对速度。

本发明的技术方案是：采用MIC众核协处理器使用多线程并行进行序列比对，并将MIC上序列比对过程中从主存读取序列到MIC、比对序列和将比对结果写回主存这三个串行步骤采用三级流水方式，即在进行每一轮序列比对的同时，读取下一组比对所需要的序列，同时将上一组比对的结果写回主存，使得读写操作和比对操作并行进行。

具体技术方案如下：

主要变量定义：

M_CPU：计算机主存可用空间大小。

M_DNA：DNA短序列所占空间大小。

M_MIC：MIC上可用空间大小。

M_REF：参考基因组大小。

M_SEQ：MIC上每块缓存空间大小。

步骤1：CPU根据计算机主存可用空间大小M_CPU，以及DNA短序列所占空间大小M_DNA，将DNA短序列即read平均分为L组，L为正整数， 表示对“x”向上取整；

步骤2：在MIC上声明三个指针变量：Seqs_ptr、Read_ptr、和Write_ptr，并根据MIC上可用空间大小M_MIC以及参考基因组大小M_REF，在MIC上分别为三个指针分配同样大小的缓存空间，MIC上每块缓存空间大小M_SEQ＝(M_MIC-M_REF)/3，其中Seqs_ptr指向的空间存储当前一组正在比对的序列，Read_ptr指向的空间存储下一组将要比对的序列，Write_ptr指向的空间存储上一组序列比对的结果；

步骤3：CPU初始化循环变量i为零；

步骤4：CPU将磁盘中L组read中的第i组读入主存中；

步骤5：CPU根据MIC上Seqs_ptr指向的空间的大小M_SEQ，将读入主存中的read平均分为M小组，M为正整数，

步骤6：CPU将循环变量m置为零；

步骤7：CPU将主存中M个小组read中的第m组读到MIC的Seqs_ptr指向的空间内；

步骤8：根据MIC上可用计算核心数Core_MIC，以及MIC上每个计算核心支持的最大硬件线程数Thread_MIC，在MIC上同时启动N+2(N>0，为整数，N＝(Core_MIC-1)*Thread_MIC，其中MIC卡需要保留一个核心处理主存与MIC卡之间数据调度)个线程，线程编号为0到N+1，N+2个线程并行执行以下步骤：

步骤8.1：第0到第N-1号线程并行比对Seqs_ptr对应空间中的read，并将比较结果写入Seqs_ptr对应空间中；比对方法采用李恒在论文《Fast and accurate short readalignment with Burrows-Wheeler Transform》中公布的基于BW(Burrows-Wheeler)变换的方法，这一步与背景技术步骤3的区别在于是在MIC上由N个线程并行地对进行比对，而不是在CPU上由单个线程完成，将Seqs_ptr对应空间中的所有read比对完毕后，转到步骤9；

步骤8.2：第N号线程将循环变量m加1，判断m是否等于M，如果m不等于M，执行步骤8.2.1，如果m等于M，结束第N号线程，转到步骤9；

步骤8.2.1：第N号线程将主存中M个小组read中的第m组读到MIC的Read_ptr对应空间内，读取完毕后，转到步骤9；

步骤8.3：第N+1号线程判断Write_ptr对应空间是否为空，如果Write_ptr对应空间不为空，执行步骤8.3.1，如果Write_ptr对应空间为空，结束第N+1号线程，转到步骤9；

步骤8.3.1：第N+1号线程将Write_ptr对应空间中的read比对结果写回主存，写完成后转到步骤9；

步骤9：同步第0到第N+1号线程，同步完成后，在MIC上的多线程部分执行完毕，以下步骤为单线程执行；

步骤10：MIC进行指针交换，在MIC上声明临时指针tmp_ptr，将Seqs_ptr值赋给tmp_ptr，将Read_ptr值赋给Seqs_ptr，将Write_ptr值赋给Read_ptr，将tmp_ptr值赋给Write_ptr，将tmp_ptr值置为空。

步骤11：MIC判断Seqs_ptr对应空间是否为空，如果不为空，转步骤8，如果为空，执行步骤12；

步骤12：MIC将Write_ptr对应空间中的read比对结果写回主存；

步骤13：CPU将内存中第i大组read比对的结果写回磁盘，并清空相应内存空间；

步骤14：CPU将循环变量i加1；

步骤15：CPU判断i是否等于L，如果i不等于L，转步骤4，如果i等于L，执行步骤16；

步骤16：MIC释放Seqs_ptr、Read_ptr、Write_ptr指向的空间；

步骤17：结束比对。

采用本发明可以达到以下技术效果：

本发明通过多线程并行技术，利用新型运算器件MIC众核协处理器进行序列比对，其中

步骤8在MIC上利用N+2(N>0)个线程并行进行序列比对，提高了比对速度，在多线程下可以取得接近线性的加速比，并实现了序列读取，序列比对和比对结果写回这三个主要过程的并行执行，提高了比对效率，降低了比对时间，与两路八核CPU相比，本发明可加速比对过程2.3倍以上。

步骤9在MIC上通过交换指针的方式实现数据交换，而不是直接对变量进行赋值，避免拷贝大量内存空间，提高了程序的时空效率。

附图说明

图1是本发明总体流程图；

图2是对图1中步骤8的分解示意。

具体实施方式

国防科大采用配备两路八核2.4GHz CPU和一块57核1.1GHz MIC卡的服务器作为环境，服务器硬盘大小为43TB，内存大小为132GB，MIC卡片上存储空间大小为6GB，输入数据为人类基因组，参考基因组占空间大小为3GB，DNA短序列占空间大小为240GB，包括8千万条序列，验证本发明的效果：

如图1所示，具体实施步骤如下：

步骤1：CPU根据计算机主存可用空间大小M_CPU＝45GB(操作系统及其他服务占用了一定内存，另外还需要预留一部分内存作程序运行时使用，所以可用内存大小为45GB，小于安装内存大小132GB)，以及DNA短序列所占空间大小M_DNA＝240GB，将DNA短序列即read平均分为L＝6组，L为正整数，

步骤2：在MIC上声明三个指针变量：Seqs_ptr、Read_ptr、和Write_ptr，并根据MIC上可用空间大小M_MIC＝4.5GB(MIC内存大小为6GB，其中1.5GB空间用于存储MIC片上操作系统以及运行时使用)以及参考基因组大小M_REF＝3GB，在MIC上分别为三个指针分配同样大小的缓存空间，每块空间大小M_SEQ＝(M_MIC-M_REF)/3＝0.5GB，其中Seqs_ptr指向的空间存储当前一组正在比对的序列，Read_ptr指向的空间存储下一组将要比对的序列，Write_ptr指向的空间存储上一组序列比对的结果；

步骤3：CPU初始化循环变量i为零；

步骤4：CPU将磁盘中L＝6组read中的第i组读入主存中；

步骤5：CPU根据MIC上Seqs_ptr指向的空间的大小M_SEQ＝0.5GB，将读入主存中的read平均分为M＝80小组，M为正整数，

步骤6：CPU将循环变量m置为零；

步骤7：CPU将主存中M＝80个小组read中的第m组读到MIC的Seqs_ptr指向的空间内；

步骤8如图2所示：根据MIC上可用计算核心数Core_MIC＝57，以及MIC上每个计算核心支持的最大硬件线程数Thread_MIC＝4，在MIC上同时启动N+2＝226(N>0，为整数，N＝(Core_MIC-1)*Thread_MIC＝(57-1)*4＝224，其中MIC卡需要保留一个核心处理主存与MIC卡之间数据调度)个线程，线程编号为0到N+1＝225，N+2＝226个线程并行执行以下步骤：

步骤8.1：第0到第N-1＝223号线程并行比对Seqs_ptr对应空间中的read，并将比较结果写入Seqs_ptr对应空间中；比对方法采用李恒在论文《Fast and accurate shortread alignment with Burrows-Wheeler Transform》中公布的基于BW(Burrows-Wheeler)变换的方法，这一步与背景技术步骤3的区别在于是在MIC上由N＝224个线程并行地对进行比对，而不是在CPU上由单个线程完成，Seqs_ptr对空间中的所有read比较完毕后，转到步骤9；

步骤8.2：第N＝224号线程将循环变量m加1，判断m是否等于M＝80，如果m不等于M＝80，执行步骤8.2.1，如果m等于M＝80，结束第N＝224号线程，转到步骤9；

步骤8.2.1：第N＝224号线程将主存中M＝80个小组read中的第m组读到MIC的Read_ptr对应空间内，读取完毕后，转到步骤9；

步骤8.3：第N+1＝225号线程判断Write_ptr对应空间是否为空，如果Write_ptr对应空间不为空，执行步骤8.3.1，如果Write_ptr对应空间为空，结束第N+1＝225号线程，转到步骤9；

步骤8.3.1：第N+1＝225号线程将Write_ptr对应空间中的read比对结果写回主存，写完成后转到步骤9；

步骤9：同步第0到第N+1＝225号线程，同步完成后，在MIC上的多线程部分执行完毕，以下步骤为单线程执行；

步骤10：MIC进行指针交换，在MIC上声明临时指针tmp_ptr，将Seqs_ptr值赋给tmp_ptr，将Read_ptr值赋给Seqs_ptr，将Write_ptr值赋给Read_ptr，将tmp_ptr值赋给Write_ptr，将tmp_ptr值置为空。

步骤11：MIC判断Seqs_ptr对应空间是否为空，如果不为空，转步骤8，如果为空，执行步骤12；

步骤12：MIC将Write_ptr对应空间中的read比对结果写回主存；

步骤13：CPU将内存中第i大组read比对的结果写回磁盘，并清空相应内存空间；

步骤14：CPU将循环变量i加1；

步骤15：CPU判断i是否等于L＝6，如果i不等于L＝6，转步骤4，如果i等于L＝6，执行步骤16；

步骤16：MIC释放Seqs_ptr、Read_ptr、Write_ptr指向的空间；

步骤17：结束比对。

统计运行数据后发现，与仅使用两路八核CPU并行进行序列比对相比，本发明可以将序列比对速度提高2.3倍以上。

Claims (1)

1.一种基于众核协处理器的三级流水序列比对方法，其特征在于包括以下步骤：

步骤1：CPU根据计算机主存可用空间大小M_CPU，以及DNA短序列所占空间大小M_DNA，将DNA短序列即read平均分为L组，L为正整数， 表示对“x”向上取整；

步骤2：在MIC即众核协处理器上声明三个指针变量：Seqs_ptr、Read_ptr和Write_ptr，并根据MIC上可用空间大小M_MIC以及参考基因组大小M_REF，在MIC上分别为三个指针分配同样大小的缓存空间，MIC上每块缓存空间大小M_SEQ＝(M_MIC-M_REF)/3，其中Seqs_ptr指向的空间存储当前一组正在比对的序列，Read_ptr指向的空间存储下一组将要比对的序列，Write_ptr指向的空间存储上一组序列比对的结果；

步骤3：CPU初始化循环变量i为零；

步骤4：CPU将磁盘中L组read中的第i组读入主存中；

步骤5：CPU根据MIC上Seqs_ptr指向的空间的大小M_SEQ，将读入主存中的read平均分为M小组，M为正整数，

步骤6：CPU将循环变量m置为零；

步骤7：CPU将主存中M个小组read中的第m组读到MIC的Seqs_ptr指向的空间内；

步骤8：根据MIC上可用计算核心数Core_MIC，以及MIC上每个计算核心支持的最大硬件线程数Thread_MIC，在MIC上同时启动N+2个线程，线程编号为0到N+1，N>0，为整数，N＝(Core_MIC-1)*Thread_MIC，其中MIC卡需要保留一个核心处理主存与MIC卡之间数据调度，N+2个线程并行执行以下步骤：

步骤8.1：第0到第N-1号线程并行比对Seqs_ptr对应空间中的read，并将比较结果写入Seqs_ptr对应空间中；比对方法采用李恒在论文《Fast and accurate short readalignment with Burrows-Wheeler Transform》中公布的基于BW变换的方法，将Seqs_ptr对应空间中的所有read比对完毕后，转到步骤9；

步骤8.2：第N号线程将循环变量m加1，判断m是否等于M，如果m不等于M，执行步骤8.2.1，如果m等于M，结束第N号线程，转到步骤9；

步骤8.2.1：第N号线程将主存中M个小组read中的第m组读到MIC的Read_ptr对应空间内，读取完毕后，转到步骤9；

步骤8.3：第N+1号线程判断Write_ptr对应空间是否为空，如果Write_ptr对应空间不为空，执行步骤8.3.1，如果Write_ptr对应空间为空，结束第N+1号线程，转到步骤9；

步骤8.3.1：第N+1号线程将Write_ptr对应空间中的read比对结果写回主存，写完成后转到步骤9；

步骤9：同步第0到第N+1号线程，同步完成后，在MIC上的多线程部分执行完毕，以下步骤为单线程执行；

步骤10：MIC进行指针交换，在MIC上声明临时指针tmp_ptr，将Seqs_ptr值赋给tmp_ptr，将Read_ptr值赋给Seqs_ptr，将Write_ptr值赋给Read_ptr，将tmp_ptr值赋给Write_ptr，将tmp_ptr值置为空；

步骤11：MIC判断Seqs_ptr对应空间是否为空，如果不为空，转步骤8，如果为空，执行步骤12；

步骤12：MIC将Write_ptr对应空间中的read比对结果写回主存；

步骤13：CPU将内存中第i大组read比对的结果写回磁盘，并清空相应内存空间；

步骤14：CPU将循环变量i加1；

步骤15：CPU判断i是否等于L，如果i不等于L，转步骤4，如果i等于L，执行步骤16；

步骤16：MIC释放Seqs_ptr、Read_ptr、Write_ptr指向的空间；

步骤17：结束比对。