CN107609350B - 一种二代测序数据分析平台的数据处理方法 - Google Patents

Abstract

本发明一种二代测序数据分析平台的数据处理方法，其中二代测序数据分析平台IMP将整个二代测序NGS处理流程实现为从输入FASTQ文件格式的短读长序列到输出标准VCF文件格式的变异检测的单个步骤，同时，还提供以标准SAM或BAM格式输出序列比对中间结果的选项，通过大量的内存访问、而不是使用缓慢的I/O来交换数据，可避免较慢的硬盘和SSD的I/O访问所需的数据搜索和加载时间，使哈希表写入或读出、删除重复比对记录，以及变异检测都更为迅速，在不影响分析质量的前提下，能实现快速的二代测序NGS数据分析，与现有方案相比速度提升达20倍。

Description

一种二代测序数据分析平台的数据处理方法

技术领域

本发明涉及一种二代测序数据分析平台的数据处理方法。

背景技术

随着人类基因组计划的顺利实施和测序技术的快速发展，测序的成本显著降低，而测序速度得到了显著提高，人类全基因组测序的测序成本已经降至$1000以内，DNA序列的数据量呈指数增长。如何快速的利用、表达这些数据，进而分析与解释基因序列里的潜在问题，从海量数据里发掘出对人类有利的信息，成为一个迫切需要解决的问题。应用越来越广泛的人类全基因组测序（WGS）产生的序列数据、以及对海量序列数据进行快速分析处理的持续需求，使数据分析形成了一个新的技术瓶颈，对二代测序技术的临床应用成为制约。

同时，为了推动精准医疗，二代测序技术的临床应用对数据分析工具有如下的要求。第一，对程序运行时间上的要求，数据分析方法速度要快。由于二代测序技术产出数据的通量越来越高，检测数据分析方法的速度需要与之相匹配，才能够达到快速确认，快速应对的目的。第二，对数据的私密性的要求，基因数据的隐秘性和安全性需要得到保证。第三，分析精度上的要求。

目前在国际上生物信息学领域有许多二代测序数据分析工具可用，其中最广泛使用的短读长序列比对的工具包括SOAP3-dp、BWA-aln、BWA-mem和Arioc等，最常用的变异检测工具包括GATK HaplotypeCaller、Samtools-mpileup和freebayes等。在二代测序数据，尤其是人类全基因组测序的分析流程中被广泛采用的做法是BWA-GATK流程，该流程在变异检测的准确性方面实现了高性能，然而，将整个流程应用于人类全基因组测序WGS中是非常耗时的。

GATK流程由若干个独立的模块组成，分别完成序列比对、排序、去除重复序列、以及最后的变异检测各项任务。其中：

步骤1、序列比对是最基本、最重要的操作，序列比对时，将输入的短读长序列匹配到参考序列上，并生成SAM格式的比对文件；

步骤2、排序，是对SAM格式的比对文件中的所有序列比对记录，按照其在参考序列上的比对位置重新排序，并产生新的BAM文件；

步骤3、去除重复序列，是为了去除PCR扩增过程中产生的重复序列。在制备文库的过程中，由于PCR扩增过程中会产生一些偏差，有的序列被过量扩增，这些扩增出来的完全相同的序列会被比对到基因组的相同位置，从而影响到变异检测的精确度。因此，这个步骤会对这些由PCR扩增过程中产生的重复序列进行标记或者去除后产生新的BAM文件并输出，该输出文件为步骤4变异检测的输入。

该数据分析流程模块化强、步骤清晰，但是在实际应用过程中，由于二代测序尤其是人类全基因组测序的数据量大，在每一个模块之间从硬盘读写文件的IO十分耗时，使得整个流程的工作时间很长。例如30倍人类全基因组数据的分析处理，通常需要超过20个小时。全基因组数据分析也可以采用超级计算机中心的超级计算机完成，但是超级计算机的租金很高，资源也很有限。

除了单机解决方案，全基因组数据分析也可以采用集群计算机方案，利用分布式计算资源，将计算任务拆分由数台计算机同时计算，再将结果汇总。集群方案又包括公有云和私有云。公有云方案是指采用基于云计算的网络服务平台，租用云计算存储和计算资源。其优点是不需要自己维护硬件，缺点是需要对海量基因数据进行网络传输和存储，同时在如何保护隐秘性和安全性的情况下开放基因数据也成为云平台方案的一大挑战。私有云方案例如搭建小型的服务器工作站，但是需要有专门的技术人员维护和管理，硬件成本和维护成本高。

发明内容

本发明的目的在于提供一种二代测序数据分析平台IMP的数据处理方法，所有的数据处理都是基于内存数据的存储和计算，从而避免了在多个处理步骤之间使用基于文件的中间结果的导入导出，减少了I/O开销，提高了运行效率。

本发明一种二代测序数据分析平台IMP的数据处理方法，包括如下步骤：

步骤1、二代测序数据分析平台IMP输入短读长序列文件和经过索引的参考序列；

步骤2、序列对比时每次读取一定长度的短读长序列放入缓存，采取多线程工作模式逐一比对，将输入的短读长序列匹配到参考序列上，对每一条成功比对上的短读长序列产生一条或多条序列比对记录：

步骤3、同时，序列比对记录经过数据压缩、通过哈希表去除重复序列和排序处理后写入内存，供后续的变异检测模块使用，序列比对记录在内存中通过哈希表寻址，每一条序列比对记录以全局比对位置作为哈希表键值，用来计算该序列比对记录的哈希值，具有相同哈希表键值的序列比对记录形成链表，根据全局比对位置进行排序；

步骤4、将经过排序和去除重复序列的序列比对记录从内存中输出成SAMDedup文件或者BAM文件；

步骤5、变异检测模块将参考序列分段，采用多CPU并行处理，最终输出VCF文件。

进一步的，步骤2中所述的一条序列比对记录，包括以CORE数据结构表示的必须的字段，和以EXT数据结构表示的可选的字段，所述CORE数据结构仅包括变异检测所需的字段，而EXT数据结构包括如果指定为输出BAM文件，则要写入BAM文件的其他字段；

A.对于单端和双端短读长序列，所述的CORE数据结构包含字段SIZE、OFFSET、POSITION、TLEN、FLAG、MAPQ、CIGAR、SUBREAD、QUAL，其中：

SIZE：当前序列比对记录在内存中所需的总字节数；

OFFSET：当前序列比对记录在当前内存块中的相对地址；

POSITION：该短读长序列在参考序列中的比对位置，是从0到（N-1）的全局位置，其中N是参考序列的长度

TLEN：该值仅用于双端短读长序列的序列比对，表示该序列的比对位置与相应的mate序列的比对位置之间的距离；

FLAG：位标志，每一个bit代表一种比对情况，与标准SAM文件里的FLAG一致；

MAPQ：由比对算法给出的比对质量，MAPQ的值范围为0到60，并使用单个字节表示；

CIGAR：简要比对信息表达式，以参考序列为基础，使用数字加字母表示比对结果；

SUBREAD：比对序列的子序列，指完全重建该短读长序列本身所需的原始序列的子集，子序列中的每个碱基对使用3位进行无损编码；

QUAL：序列的质量信息，使用Rice编码进行无损压缩；

B.对于单端和双端短读长序列，所述单端短读长序列的EXT数据包含字段RNAME、MD、QNAME、AS、XS、NM和RNEXT，其中：

RNAME：当参考序列包括多条染色体时，表示该单端短读长序列在参考序列中对应的染色体名称；

MD：该字符串用于表示从该单端短读长序列完全重构相应位置的参考序列的子序列；

QNAME：表示该单端短读长序列的名称；

AS：表示序列比对的分数；

XS：序列比对给出多个序列比对记录时，将该单端短读长序列映射到参考序列的不同位置，XS表达第二位的序列比对记录的分数；

NM：从该单端短读长序列到参考序列的编辑距离，即从该单端短读长序列变换到参考序列对应位置的子序列所需的编辑次数；

所述双端短读长序列还包括附加字段RNEXT、PNEXT。其中：

RNEXT：是该双端短读长序列的mate在参考序列中对应的染色体名称；

PNEXT：是该双端短读长序列的mate在参考序列中的比对位置。

进一步的，采用多CPU并行处理时使用分块存储的内存管理办法：

首先，对基因比对记录的处理是按基因区域分段进行的，各段的操作相对独立，将参考基因数据和序列比对记录按照其在参考基因所在的位置进行分块存储，每一区域块对应一段固定长度的参考基因；

然后，每个区域块里的数据再按照其在块中的相对位置建立本地哈希表，分块后的数据连同本地哈希表保存到操作系统下的共享内存缓冲区中，以方便后续测序数据分析流程以多进程方式访问这些数据；

二代测序数据分析平台IMP对共享内存缓冲区执行写入或读出序列比对记录时，首先计算出序列比对记录所属的区域块，然后进入该区域块，根据该区域块的本地哈希表的键值，查找到该键值对应的缓冲区地址查找到相应的比对记录，保证每个区域块的数据总量不超过一个节点处理器的内存资源，以避免IMP进程在对区域块进行操作时进行QPI内存访问；

二代测序数据分析平台IMP同时运行多个进程，每个进程分别对不同的基因区段进行计算，各个进程的输出结果再按顺序拼接成最终输出文件，由于每个单独的进程仅对分配到的区域块进行操作，运行时，每个进程和对应的区域块分配到同一节点处理器上进行计算。

进一步的，步骤3中序列比对记录进行数据压缩，包括如下步骤：

（1）短读长序列基于参考序列的压缩

为每个短读长序列提取一个子序列，该子序列通过使用全局比对位置得到的该位置处的相应参考基因组序列和通过扩展的CIGAR字符串来完全重建整个短读长序列；

CIGAR操作符包括：

S、软剪切，用于表示序列头和尾被截取的部分

M、用于表示当前位匹配或不匹配

I、用于表示当前位相对于参考序列有插入碱基

D、用于表示当前位相对于参考序列有碱基被删除

为了后续变异检测软件的需要，扩展CIGAR操作符时，将“M”分为“X”和“=”，其中“X”表示不匹配，“=”表示匹配；在子序列中记录该序列里相应位置的碱基，根据参考序列、扩展后的CIGAR操作符、以及子序列，可以完全重构出当前序列的全部碱基；对于匹配“=”和删除“D”，则不需要记录任何碱基；

（2）短读长序列的质量信息采用差异Rice编码进行无损压缩

短读长序列的质量信息字符串的第一个字符按原始数据编码，其余的质量信息字符按与之前一个字符的差异进行编码；

在对整个短读长序列的质量信息进行Rice编码后，如果发现Rice编码的压缩率达不到阈值要求，则编码器自动切换成原始编码模式；将编码后的比特流的第一位作为标志位，用于区分质量信息是采用原始编码还是Rice编码。

进一步的，步骤3中通过哈希表去除重复序列，包含如下步骤：

将序列比对记录所对应的DNA片段通过哈希表写入内存，在内存中记录DNA片段的起始位置和长度、以及对应于该DNA片段的最佳短读长序列在内存中的检索位置，该DNA片段的起始位置和长度可以唯一确定一个DNA片段，对于每一个DNA片段，只保留一个或一对最佳短读长序列，来自于同一DNA片段的短读长序列都被认为是重复序列，通过序列的质量信息来确定是否为最佳短读长序列，短读长序列的每一个碱基都有一个质量值，将短读长序列所有碱基的质量值相加，所得的质量值越大，则认为该短读长序列越佳；

PCR重复序列的处理是通过以下两个布尔型变量来控制的：

标记重复序列MarkDuplicate：缺省为true，表示来自同一DNA片段的所有短读长序列中，除了质量最佳的那一个或一对，其他的都会被标记为PCR duplicate，并显示在程序输出的SAM或BAM文件里；

去除重复序列RemoveDuplicate：缺省为false，表示被标记为PCR duplicate的那些短读长序列仍然会被保留在输出的SAM或BAM文件里，如果设为true，则会直接从SAM/BAM文件里删除；

在去除重复序列的过程中，DNA片段的哈希表操作和比对序列的哈希表操作交替进行的；

需要保存一个新的比对序列时，先根据该比对序列的比对位置和软剪切计算出相应的DNA片段，然后在DNA片段的哈希表里查找该DNA片段记录是否已经存在，如果已经存在，则比较当前比对序列与已有比对序列的质量信息，将质量信息较低的标记为重复序列，如果已有的比对序列为重复序列，则需要根据DNA片段上保存的地址信息找到该比对序列，并修改其标记，如果当前正在保存的新的比对序列被判断为重复序列，则在写入内存前直接修改其标记；如果不存在,则将该DNA片段写入内存，不需要标记任何比对序列为重复序列；标记完重复序列后，再将标记过的当前序列的比对记录写入内存；

从内存中导出比对记录时，仅导出未被标记为重复序列的比对记录，从而完全地去除重复序列。

进一步的，步骤3中通过哈希表进行排序，包括如下具体步骤：

所述的序列比对记录以全局比对位置作为哈希表键值，将一条新的序列比对记录插入到内存中时，将其与已经使用相同键值插入哈希表的序列比对记录进行比较，当新的序列比对记录的比对位置与现有的序列比对记录的比对位置相等时，对于双端测序数据和单端测序数据处理步骤分别为：

（1）对于双端测序数据，将其mate的比对位置与现有的序列比对记录的mate的比对位置进行比较，按照mate的比对位置排序，如果不相等，则认为两条记录不是重复的，并将新记录写入内存。比对位置和mate的比对位置均相同的情况下，按照写入内存的顺序排列，如果两条记录的mate的比对位置相等，则比较两条记录的比对质量MAPQ，如果新记录的比对质量MAPQ等于或低于现有记录的比对质量MAPQ，则将其丢弃，否则将写入内存；

（2）单端测序数据只有一条短读长序列而没有mate，所以直接比较两条记录的比对质量MAPQ，如果新记录的比对质量MAPQ比现有记录的高，才按照比对位置排序写入内存，比对位置相同的情况下按照写入内存的顺序排列。

本发明一种二代测序数据分析平台的数据处理方法，其中二代测序数据分析平台IMP将整个二代测序NGS处理流程实现为从输入FASTQ文件格式的短读长序列到输出标准VCF文件格式的变异检测的单个步骤，同时，还提供以标准SAM或BAM格式输出序列比对中间结果的选项，通过大量的内存访问、而不是使用缓慢的I/O来交换数据，可避免较慢的硬盘和SSD的I/O访问所需的数据搜索和加载时间，使哈希表写入或读出、删除重复比对记录，以及变异检测都更为迅速，在不影响分析质量的前提下，能实现快速的二代测序NGS数据分析，与现有方案相比速度提升达20倍。

附图说明

图1为本发明二代测序数据分析平台IMP的数据处理分析流程图；

图2为本发明通过哈希表将序列比对记录写入内存的流程图；

图3为多个CPU的多路服务器通常使用NUMA非统一内存访问的内存管理架构示意图；

图4为多线程应用程序在多路服务器上运行的示意图；

图5为本发明基于共享内存和分块哈希表的内存管理示意图；

图6为本发明基于区域块的变异检测计算示意图。

以下结合附图和具体实施例对本发明作进一步详述。

具体实施方式

本发明一种二代测序数据分析平台IMP的数据处理方法，使用哈希表寻址，通过哈希表支持数据排序、去除重复序列的功能，通过数据的无损压缩避免过多的使用内存，所有的数据处理都是基于内存数据的存储和计算，各模块内部以及不同模块之间的多线程并行处理。

如图 1所示，本发明一种二代测序数据分析平台IMP的数据处理方法，具体包括如下步骤：

步骤1、二代测序数据分析平台IMP输入短读长序列文件和经过索引的参考序列；

步骤2、序列对比时每次读取一定长度的短读长序列放入缓存，采取GPU或CPU上的多线程工作模式逐一比对，将输入的短读长序列匹配到参考序列上，对每一条成功比对上的短读长序列产生一条或多条序列比对记录，将序列比对记录直接输出成原始SAM文件：

所述的一条序列比对记录，包括以CORE数据结构表示的必须的字段，和以EXT数据结构表示的可选的字段，所述CORE数据结构仅包括变异检测所需的字段，而EXT数据结构包括如果指定为输出BAM文件，则要写入BAM文件的其他字段；

A.对于单端和双端短读长序列，所述的CORE数据结构包含字段SIZE、OFFSET、POSITION、TLEN、FLAG、MAPQ、CIGAR、SUBREAD、QUAL，其中：

SIZE：当前序列比对记录在内存中所需的总字节数；

OFFSET：当前序列比对记录在当前内存块中的相对地址；

POSITION：该短读长序列在参考序列中的比对位置，是从0到（N-1）的全局位置，其中N是参考序列的长度

TLEN：该值仅用于双端短读长序列的序列比对，表示该序列的比对位置与相应的mate序列的比对位置之间的距离；

FLAG：位标志，每一个bit代表一种比对情况，与标准SAM文件里的FLAG一致；

MAPQ：由比对算法给出的比对质量，MAPQ的值范围为0到60，并使用单个字节表示；

CIGAR：简要比对信息表达式，以参考序列为基础，使用数字加字母表示比对结果；

SUBREAD：比对序列的子序列，指完全重建该短读长序列本身所需的原始序列的子集，子序列中的每个碱基对使用3位进行无损编码；

QUAL：序列的质量信息，使用Rice编码进行无损压缩；

B.对于单端和双端短读长序列，所述单端短读长序列的EXT数据包含字段RNAME、MD、QNAME、AS、XS、NM和RNEXT，其中：

RNAME：当参考序列包括多条染色体时，表示该单端短读长序列在参考序列中对应的染色体名称；

MD：该字符串用于表示从该单端短读长序列完全重构相应位置的参考序列的子序列；

QNAME：表示该单端短读长序列的名称；

AS：表示序列比对的分数；

XS：序列比对给出多个序列比对记录时，将该单端短读长序列映射到参考序列的不同位置，XS表达第二位的序列比对记录的分数；

NM：从该单端短读长序列到参考序列的编辑距离，即从该单端短读长序列变换到参考序列对应位置的子序列所需的编辑次数；

所述双端短读长序列还包括附加字段RNEXT、PNEXT。其中：

RNEXT：是该双端短读长序列的mate在参考序列中对应的染色体名称；

PNEXT：是该双端短读长序列的mate在参考序列中的比对位置。

对于具有可变长度的数据字段，该字段的长度也被编码为单独的变量。总的来说，对比记录的大小（其指示存储器中当前对比记录的总字节数）被添加到对比记录的开头；

步骤3、如图2所示，同时，序列比对记录经过数据压缩，通过哈希表进行标记重复序列和排序等处理后写入内存，供后续的变异检测模块使用，序列比对记录在内存中通过哈希表寻址，每一条序列比对记录以全局比对位置作为哈希表键值，用来计算该序列比对记录的哈希值，具有相同哈希表键值的序列比对记录形成链表，根据全局比对位置进行排序：

所述的哈希表是根据键值而直接访问在内存存储位置的数据结构，通过哈希函数将输入数据对应的键值转换为哈希表键值，然后通过哈希表键值将所需查询的数据映射到表中一个位置来访问记录，加快查找速度；

步骤3.1、压缩比对数据

由于原始的序列比对数据的数据量太大（对于30x WGS数据（包括序列本身和质量数据）的SAM文件约300~400 G），即使对于今天的高性能计算机也难以直接放入计算机内存中。为了避免过多的使用内存，压缩比对数据具体包括如下内容：

（1）短读长序列基于参考序列的压缩

为每个短读长序列提取一个子序列，该子序列通过使用全局比对位置得到的该位置处的相应参考基因组序列和经过扩展的CIGAR字符串来完全重建整个短读长序列；

CIGAR操作符包括：

S、软件切，用于表示序列头和尾被截取的部分

M、用于表示当前位匹配或不匹配

I、用于表示当前位相对于参考序列有插入碱基

D、用于表示当前位相对于参考序列有碱基被删除

为了后续变异分析软件的需要，扩展CIGAR操作符时，将“M”分为“X”和“=”，其中“X”表示不匹配，“=”表示匹配；同时，在扩展CIGAR字符串中具有相应长度的“S”，“X”和“I”时，说明序列中相应位置的碱基不存在于参考序列中，我们在子序列中记录该序列里相应位置的碱基，根据参考序列、扩展后的CIGAR操作符、以及子序列，可以完全重构出当前序列的全部碱基；对于匹配“=”和删除“D”，则不需要记录任何碱基；

下面给出一个示例，其中子序列由对应于位置3处的“1X”的碱基A、对应位置11处的“1I”的碱基T，和对应位置38处的“1X”的碱基C组成：

参考序列：GTGTTTAATACATTTAAATTTATATAGTTACTGATAAGTTAGATTC

短读长序列：GTATTTAATATCATTTAAATTTATATATTACTGATAACTTAGATTC

CIGAR序列：2M1X7M1I16M1D10M1X8M

子序列：ATC

子序列读取中的每个碱基对由三位编码，如下所示：

A，100

C，101

G，110

T，111

N，011

则，子序列ATC的编码为100,111,101（二进制），即0x13D（十六进制），用一个无符号长整型64位（UINT64）来表示21个碱基对并写入内存；这样，子序列的编码直接通过位操作完成，同时，UINT64的数量也被写入内存，以便对子序列进行正确解码；

（2）短读长序列的质量信息采用差异Rice编码进行无损压缩

短读长序列的质量信息字符串的第一个字符按原始数据编码，其余的质量字符按与之前一个字符的差异进行编码；

由于短读长序列的质量信息字符的范围有限，所以质量字符的值用6个比特位进行原始编码，压缩率为75％。在对整个序列的质量信息进行Rice编码后，如果发现Rice编码的压缩率达不到75%，则编码器将自动切换成原始编码模式。为了区分质量信息是采用原始编码还是Rice编码，将编码后的比特流的第一位做为标志位，1表示该质量信息采用Rice编码，0表示该质量信息采用原始编码。

步骤3.2、通过哈希表去除重复序列

将序列比对记录所对应的DNA片段通过哈希表写入内存，在内存中记录DNA片段的起始位置和长度、以及对应于该DNA片段的最佳短读长序列在内存中的检索位置，该DNA片段的起始位置和长度可以唯一确定一个DNA片段，对于每一个DNA片段，只保留一个或一对最佳短读长序列，来自于同一DNA片段的短读长序列都被认为是重复序列，通过序列的质量信息来确定是否为最佳短读长序列，短读长序列的每一个碱基都有一个质量值，将短读长序列所有碱基的质量值相加，所得的质量值越大，则认为该短读长序列越佳；

PCR重复序列的处理是通过以下两个布尔型变量来控制的：

标记重复序列MarkDuplicate：缺省为true，表示来自同一DNA片段的所有短读长序列中，除了质量最佳的那一个或一对，其他的都会被标记为PCR duplicate，并显示在程序输出的SAM或BAM文件里；

去除重复序列RemoveDuplicate：缺省为false，表示被标记为PCR duplicate的那些短读长序列仍然会被保留在输出的SAM或BAM文件里，如果设为true，则会直接从SAM/BAM文件里删除；

在当前的程序设定里，不论这些重复序列是否被保留在SAM/BAM文件，在做变异检测时，这些重复序列都没有被用到的。在去除重复序列的过程中，DNA片段的哈希表操作和比对序列的哈希表操作是交替进行的。

需要保存一个新的比对序列时，先根据该比对序列的比对位置和软剪切计算出相应的DNA片段，然后在DNA片段的哈希表里查找该DNA片段记录是否已经存在，若如果已经存在，则比较当前比对序列与已有比对序列的质量信息，将质量信息较低的标记为重复序列，如果已有的比对序列为重复序列，则需要根据DNA片段上保存的地址信息找到该比对序列，并修改其标记（FLAG），如果当前正在保存的新的比对序列被判断为重复序列，则在写入内存前直接修改其标记；如果不存在,则将该DNA片段写入内存，不需要标记任何比对序列为重复序列；标记完重复序列后，再将标记过的当前序列的比对记录写入内存；

从内存中导出比对记录时，仅导出未被标记为重复序列的比对记录，从而实现完全的去除重复序列；

内存中的序列比对记录将会被导出用于输出SAMDedup文件、输出BAM文件、变异检测；针对双端测序数据，当对同一组数据进行多次测试时，在比对位置和mate的比对位置均相同的情况下，比对序列的顺序可能有不同，但是数据的内容是一致的，多线程多CPU的数据处理模式不会导致多次测试时程序结果不同；

步骤3.3、通过哈希表进行排序

所述的序列比对记录以全局比对位置作为哈希表键值，将一条新的序列比对记录插入到内存中时，将其与已经使用相同键值插入哈希表的序列比对记录进行比较，当新的序列比对记录的比对位置与现有的序列比对记录的比对位置相等时，对于双端测序数据和单端测序数据处理步骤分别为：

（1）对于双端测序数据，将其mate的比对位置与现有的序列比对记录的mate的比对位置进行比较，按照mate的比对位置排序，如果不相等，则认为两条记录不是重复的，并将新记录写入内存。比对位置和mate的比对位置均相同的情况下，按照写入内存的顺序排列，如果两条记录的mate的比对位置相等，则比较两条记录的比对质量MAPQ，如果新记录的比对质量MAPQ等于或低于现有记录的比对质量MAPQ，则将其丢弃，否则将写入内存；

（2）单端测序数据只有一条短读长序列而没有mate，所以直接比较两条记录的比对质量MAPQ，如果新记录的比对质量MAPQ比现有记录的高，才按照比对位置排序写入内存，比对位置相同的情况下按照写入内存的顺序排列；

通过上述步骤（1）或者（2），所有的序列比对记录都在内存中进行了排序，如图 2所示，该排序按三级实现，按当前序列比对记录的比对位置排序，在当前序列的比对位置相同的情况下，按mate序列的比对位置（仅对于双端测序序列）排序，在mate序列的比对位置均相同的情况下，按照写入内存的顺序排列；

步骤4、将经过排序和去除重复序列的序列比对记录从内存中输出成SAMDedup文件或者BAM文件；

步骤5、变异检测模块将参考序列分段，采用多CPU并行处理，最终输出VCF文件：

如图3所示，具有多个CPU的多路服务器通常使用NUMA非统一内存访问的内存管理架构，每个节点处理器使用本节点的内存控制器管理本地内存；当每个应用程序进程只需要使用一个节点处理器，且所使用的数据内存不超过一个节点所配置的内存时，该应用程序只会在一个节点处理器中运行，不会使用QPI 远程访问其他节点处理器的内存，这样会获得最佳的计算性能；

如图4所示，当多线程应用程序需要同时在多个节点处理器上运行，或者因为使用的数据量太大，需要访问不同节点的内存时，运行应用程序的节点处理器需要通过QPI快速通道互联对其他节点进行远程内存访问。这种对其他节点的内存操作，在NUMA 非统一内存访问的内存管理架构下非常低效和缓慢。操作系统会尝试把数据移动或复制到同一个节点。但如果需要使用的内存超过一个节点拥有的资源，是没法达到一个优化的分配方案的。而且这种情况下，程序所需的数据不是集中在一些物理地址空间，而是分散在整个多节点地内存中，在访问过程中需要触发大量远程内存访问，造成QPI出现瓶颈，导致程序运行效率严重下降。

在测序数据分析过程中，为了满足二代测序数据分析平台IMP测序数据分析算法对内存数量和计算性能的需求，也希望二代测序数据分析平台IMP能在多路多核CPU的服务器上以多线程方式运行。然而，由于上述原因，直接在多个CPU上运行的IMP线程需要通过QPI访问分散在整个内存的基因比对记录和参考数据内存的访问效率很低，严重影响整个系统的性能。所以在多CPU平台上使用全局的哈希表得不到良好的运算性能。

为了解决这个问题，如图5所示，本发明使用分块存储的内存管理办法：

首先，对基因比对记录的处理（比如变异检测）是按基因区域分段进行的，各段的操作相对独立，将参考基因数据和序列比对记录按照其在参考基因所在的位置进行分块存储，每一区域块对应一段固定长度的参考基因；

然后，每个区域块里的数据再按照其在块中的相对位置建立本地哈希表，分块后的数据连同本地哈希表保存到Linux系统下的共享内存缓冲区中（或Windows系统下的映射文件共享内存），以方便二代测序数据分析平台IMP以多进程方式访问这些数据；

二代测序数据分析平台IMP对共享内存缓冲区执行写入或读出序列比对记录时，首先计算出序列比对记录所属的区域块，然后进入该区域块，根据该区域块的本地哈希表的键值，查找到该键值对应的缓冲区地址查找到相应的比对记录，保证每个区域块的数据总量不超过一个节点的内存资源，以避免IMP进程在对区域块进行操作时进行QPI内存访问，从而提高内存访问效率。

下面以变异检测为例说明二代测序数据分析平台IMP如何使用上述内存管理方法来实现高通量多线程基因数据的分析。如图6所示，为了充分发挥多CPU平台的潜力，在进行变异检测时，二代测序数据分析平台IMP同时运行多个变异检测进程，每个进程分别对不同的基因区段进行计算，各个进程的输出结果再按顺序拼接成最终输出文件。由于每个单独的变异检测进程仅对分配到的区域块进行操作，运行时，每个进程和对应的区域块都可以分配到同一节点处理器上进行计算，从而获得最佳性能。对基因对比记录的其他操作，比如SAM/BAM输出也可以用类似的原理实现并行操作。在这种情况下，因为I/O输出瓶颈, 整个输出文件进程可以在一个CPU节点上运行，在进程内部实现同时多个解压缩解码输出线程就足够了。每个线程分别对不同的基因区段的哈希表进行输出，各个线程的输出比对记录最终按顺序拼接成输出的SAM/BAM文件。

以上所述，并非对本发明的技术范围作任何限制，故凡是依据本发明的技术实质对以上实施例所作的任何细微修改、等同变化与修饰，均仍属于本发明技术方案的范围内。

Claims (5)

1.一种二代测序数据分析平台IMP的数据处理方法，其特征在于包括如下步骤：

步骤1、二代测序数据分析平台IMP输入短读长序列文件和经过索引的参考序列；

步骤2、序列对比时每次读取一定长度的短读长序列放入缓存，采取多线程工作模式逐一比对，将输入的短读长序列匹配到参考序列上，对每一条成功比对上的短读长序列产生一条或多条序列比对记录；

步骤3、同时，序列比对记录经过数据压缩、通过哈希表去除重复序列和排序处理后写入内存，供后续的变异检测模块使用，序列比对记录在内存中通过哈希表寻址，每一条序列比对记录以全局比对位置作为哈希表键值，用来计算该序列比对记录的哈希值，具有相同哈希表键值的序列比对记录形成链表，根据全局比对位置进行排序；

步骤4、将经过排序和去除重复序列的序列比对记录从内存中输出成SAMDedup文件或者BAM文件；

步骤5、变异检测模块将参考序列分段，采用多CPU并行处理，最终输出VCF文件；

采用多CPU并行处理时使用分块存储的内存管理办法：

首先，对基因比对记录的处理是按基因区域分段进行的，各段的操作相对独立，将参考基因数据和序列比对记录按照其在参考基因所在的位置进行分块存储，每一区域块对应一段固定长度的参考基因；

然后，每个区域块里的数据再按照其在块中的相对位置建立本地哈希表，分块后的数据连同本地哈希表保存到操作系统下的共享内存缓冲区中，以方便后续测序数据分析流程以多进程方式访问这些数据；

二代测序数据分析平台IMP对共享内存缓冲区执行写入或读出序列比对记录时，首先计算出序列比对记录所属的区域块，然后进入该区域块，根据该区域块的本地哈希表的键值，查找到该键值对应的缓冲区地址查找到相应的比对记录，保证每个区域块的数据总量不超过一个节点处理器的内存资源，以避免IMP进程在对区域块进行操作时进行QPI内存访问；

二代测序数据分析平台IMP同时运行多个进程，每个进程分别对不同的基因区段进行计算，各个进程的输出结果再按顺序拼接成最终输出文件，由于每个单独的进程仅对分配到的区域块进行操作，运行时，每个进程和对应的区域块分配到同一节点处理器上进行计算。

2.根据权利要求1所述的一种二代测序数据分析平台IMP的数据处理方法，其特征在于：步骤2中所述的一条序列比对记录，包括以CORE数据结构表示的必须的字段，和以EXT数据结构表示的可选的字段，所述CORE数据结构仅包括变异检测所需的字段，而EXT数据结构包括如果指定为输出BAM文件，则要写入BAM文件的其他字段；

A.对于单端和双端短读长序列，所述的CORE数据结构包含字段SIZE、OFFSET、POSITION、TLEN、FLAG、MAPQ、CIGAR、SUBREAD、QUAL，其中：

SIZE：当前序列比对记录在内存中所需的总字节数；

OFFSET：当前序列比对记录在当前内存块中的相对地址；

POSITION：该短读长序列在参考序列中的比对位置，是从0到（N-1）的全局位置，其中N是参考序列的长度

TLEN：该值仅用于双端短读长序列的序列比对，表示该序列的比对位置与相应的mate序列的比对位置之间的距离；

FLAG：位标志，每一个bit代表一种比对情况，与标准SAM文件里的FLAG一致；

MAPQ：由比对算法给出的比对质量，MAPQ的值范围为0到60，并使用单个字节表示；

CIGAR：简要比对信息表达式，以参考序列为基础，使用数字加字母表示比对结果；

SUBREAD：比对序列的子序列，指完全重建该短读长序列本身所需的原始序列的子集，子序列中的每个碱基对使用3位进行无损编码；

QUAL：序列的质量信息，使用Rice编码进行无损压缩；

B.对于单端和双端短读长序列，所述单端短读长序列的EXT数据包含字段RNAME、MD、QNAME、AS、XS、NM和RNEXT，其中：

RNAME：当参考序列包括多条染色体时，表示该单端短读长序列在参考序列中对应的染色体名称；

MD：该字符串用于表示从该单端短读长序列完全重构相应位置的参考序列的子序列；

QNAME：表示该单端短读长序列的名称；

AS：表示序列比对的分数；

XS：序列比对给出多个序列比对记录时，将该单端短读长序列映射到参考序列的不同位置，XS表达第二位的序列比对记录的分数；

NM：从该单端短读长序列到参考序列的编辑距离，即从该单端短读长序列变换到参考序列对应位置的子序列所需的编辑次数；

所述双端短读长序列还包括附加字段RNEXT、PNEXT；

其中：

RNEXT：是该双端短读长序列的mate在参考序列中对应的染色体名称；

PNEXT：是该双端短读长序列的mate在参考序列中的比对位置。

3.根据权利要求1所述的一种二代测序数据分析平台IMP的数据处理方法，其特征在于步骤3中序列比对记录进行数据压缩，包括如下步骤：

（1）短读长序列基于参考序列的压缩

为每个短读长序列提取一个子序列，该子序列通过使用全局比对位置得到的该位置处的相应参考基因组序列和通过扩展的CIGAR字符串来完全重建整个短读长序列；

CIGAR操作符包括：

S、软剪切，用于表示序列头和尾被截取的部分

M、用于表示当前位匹配或不匹配

I、用于表示当前位相对于参考序列有插入碱基

D、用于表示当前位相对于参考序列有碱基被删除

为了后续变异检测软件的需要，扩展CIGAR操作符时，将“M”分为“X”和“=”，其中“X”表示不匹配，“=”表示匹配；在子序列中记录该序列里相应位置的碱基，根据参考序列、扩展后的CIGAR操作符、以及子序列，可以完全重构出当前序列的全部碱基；对于匹配“=”和删除“D”，则不需要记录任何碱基；

（2）短读长序列的质量信息采用差异Rice编码进行无损压缩

短读长序列的质量信息字符串的第一个字符按原始数据编码，其余的质量信息字符按与之前一个字符的差异进行编码；

在对整个短读长序列的质量信息进行Rice编码后，如果发现Rice编码的压缩率达不到阈值要求，则编码器自动切换成原始编码模式；将编码后的比特流的第一位作为标志位，用于区分质量信息是采用原始编码还是Rice编码。

4.根据权利要求1所述的一种二代测序数据分析平台IMP的数据处理方法，其特征在于步骤3中通过哈希表去除重复序列，包含如下步骤：

将序列比对记录所对应的DNA片段通过哈希表写入内存，在内存中记录DNA片段的起始位置和长度、以及对应于该DNA片段的最佳短读长序列在内存中的检索位置，该DNA片段的起始位置和长度可以唯一确定一个DNA片段，对于每一个DNA片段，只保留一个或一对最佳短读长序列，来自于同一DNA片段的短读长序列都被认为是重复序列，通过序列的质量信息来确定是否为最佳短读长序列，短读长序列的每一个碱基都有一个质量值，将短读长序列所有碱基的质量值相加，所得的质量值越大，则认为该短读长序列越佳；

PCR重复序列的处理是通过以下两个布尔型变量来控制的：

标记重复序列MarkDuplicate：缺省为true，表示来自同一DNA片段的所有短读长序列中，除了质量最佳的那一个或一对，其他的都会被标记为PCR duplicate，并显示在程序输出的SAM或BAM文件里；

去除重复序列RemoveDuplicate：缺省为false，表示被标记为PCR duplicate的那些短读长序列仍然会被保留在输出的SAM或BAM文件里，如果设为true，则会直接从SAM/BAM文件里删除；

在去除重复序列的过程中，DNA片段的哈希表操作和比对序列的哈希表操作交替进行的；

需要保存一个新的比对序列时，先根据该比对序列的比对位置和软剪切计算出相应的DNA片段，然后在DNA片段的哈希表里查找该DNA片段记录是否已经存在，如果已经存在，则比较当前比对序列与已有比对序列的质量信息，将质量信息较低的标记为重复序列，如果已有的比对序列为重复序列，则需要根据DNA片段上保存的地址信息找到该比对序列，并修改其标记，如果当前正在保存的新的比对序列被判断为重复序列，则在写入内存前直接修改其标记；如果不存在,则将该DNA片段写入内存，不需要标记任何比对序列为重复序列；标记完重复序列后，再将标记过的当前序列的比对记录写入内存；

从内存中导出比对记录时，仅导出未被标记为重复序列的比对记录，从而完全地去除重复序列。

5.根据权利要求1所述的一种二代测序数据分析平台IMP的数据处理方法，其特征在于步骤3中通过哈希表进行排序，包括如下具体步骤：

所述的序列比对记录以全局比对位置作为哈希表键值，将一条新的序列比对记录插入到内存中时，将其与已经使用相同键值插入哈希表的序列比对记录进行比较，当新的序列比对记录的比对位置与现有的序列比对记录的比对位置相等时，对于双端测序数据和单端测序数据处理步骤分别为：

（1）对于双端测序数据，将其mate的比对位置与现有的序列比对记录的mate的比对位置进行比较，按照mate的比对位置排序，如果不相等，则认为两条记录不是重复的，并将新记录写入内存；

比对位置和mate的比对位置均相同的情况下，按照写入内存的顺序排列，如果两条记录的mate的比对位置相等，则比较两条记录的比对质量MAPQ，如果新记录的比对质量MAPQ等于或低于现有记录的比对质量MAPQ，则将其丢弃，否则将写入内存；

（2）单端测序数据只有一条短读长序列而没有mate，所以直接比较两条记录的比对质量MAPQ，如果新记录的比对质量MAPQ比现有记录的高，才按照比对位置排序写入内存，比对位置相同的情况下按照写入内存的顺序排列。

Images