一种基于聚类算法的高性能k-mer频次计数方法及系统

摘要

本发明公开了一种基于聚类算法的高性能k‑mer频次计数方法及系统，涉及生物信息学技术领域，该方法利用CPU读取基因组格式文件，得到长度为L的基因序列reads；通过CUDA流将文件传输至GPU，采用字符串拆分方法分裂基因序列，并进行k‑mer序列碱基坐标值转换；然后通过Kmeans聚类求取坐标数组质心，将k‑mer序列及其质心返回CPU；最后CPU对k‑mer出现频次进行计数，输出k‑mer频次分布结果。本发明大幅提高了k‑mer频次计数的速度，减小了对计算资源的占用，有助于生信分析人员以更快的速度、更短的时间获得准确的分析结果。

说明书

技术领域

本发明属于生物信息学领域，涉及一种基因组K-mer频次计数方法，尤其涉及一种基于聚类算法的k-mer频次计数方法。

背景技术

自2005年罗氏推出第一款二代测序仪罗氏454后，生命科学正式开始进入高通量测序时代。Illumina系列测序平台的推出，极大地降低了二代测序的价格，使得高通量测序在生命科学各个研究领域得到了广泛的普及。至今，第二代短读长测序技术在全球测序市场上仍然占有绝对的优势地位。

在生物信息学领域，计算每个k-mer(长度为k的子串)在一个长串中的出现次数是一个核心的子问题，包括基因组组装、测序读取的纠错、快速多序列比对和重复检测，以及评估基因组的大小、杂合程度和物种样品污染程度等。第二代测序技术虽然大大提高了测序的通量，但其获得的单条序列长度很短，往往只在50～300bp，并且因为测序的覆盖范围更深，导致基因组项目中需要处理的序列数量迅速增加，为后续基因组分析提供了巨大挑战。

目前关于k-mer计数提出了多种算法，但均有着无法避免的缺陷。如获取k-mer的过程(即分裂基因序列的过程)几乎都采用了遍历的方式，包括前向遍历与逆向遍历，其时间复杂度高；如对k-mer频次进行统计时，为了在查找时间上更高效，几乎全部采用了Hash表算法。即使是流行的k-mer计数软件jellyfish，也依然存在两个致命问题：首先，基于Hash表的计数算法中的Hash表均存在哈希冲突，使其算法效率降低，极端情况下可能出现直到耗尽最大探测次数仍然找不到合适位置的情况；其次，Hash表需要驻留在内存中以进行随机访问，jellyfish虽然实现了并行，提高了效率，但是对内存的占用也极为显著。

发明内容

本发明的目的在于克服现有技术的不足，提供一种基于聚类算法的高性能k-mer频次计数方法及系统，方法基于坐标偏移算法、无监督机器学习的聚类算法，使用唯一质心代替hash值表达k-mer序列，解决了Hash算法的低效率问题；采用CUDA数据流，配合异步操作，解决了GPU浪费及内存开销大的技术问题，提高计算速度；采用CPU+CUDA异构编程，保证CPU、GPU间高速通讯和协同计算，实现高准确度、高效率的k-mer计数。

本发明的目的是通过以下技术方案来实现的：

一种基于聚类算法的高性能k-mer频次计数方法，包括以下步骤：

CPU读取基因组格式文件，得到长度为L的基因序列reads；

CPU创建CUDA流，所述CUDA流包括以下GPU核函数及其执行的操作：

S1：按照预设子序列长度K将reads拆分为L-K+1段k-mer序列；

S2：以k-mer序列中每个碱基位置为坐标点，将拆分获得的多段k-mer序列分别转换为对应的坐标数组并添加坐标偏移；

S3：利用聚类算法计算出每个坐标数组的质心特征值，并将所有k-mer序列与k-mer序列对应的质心特征值返回至CPU；

CPU根据k-mer序列及对应的质心特征值，对所有k-mer出现的频次进行计数；

CPU输出k-mer频次分布结果。

具体的，所述基因组格式文件包括fasta或fastq格式文件。

具体的，所述步骤S2具体包括以下子步骤：

S201，对拆分后得到的k-mer序列，在ATCG笛卡尔坐标系中，以A、T为垂直坐标的正负轴，G、C为水平坐标的正负轴，以k-mer序列中每个碱基位置为坐标点，将k-mer序列中的每一个碱基都表达为坐标系中的坐标点，从而将每段k-mer序列转换为ATCG笛卡尔坐标系中对应的坐标数组；

S202，基于火星坐标系的坐标偏移算法，在每个坐标数组各个碱基的坐标值中添加坐标偏移，将坐标点转换为浮点型。

具体的，所述步骤S3具体包括以下子步骤：

S301，以碱基坐标点之间欧几里得距离的均方差作为准则函数，对每个坐标数组进行Kmeans聚类计算分析，获得每个坐标数组的质心特征值，并将其作为每个坐标数组对应的k-mer序列质心表达；

S302，将基因组格式文件中reads序列的所有k-mer序列都转换为质心表达后，将所有k-mer序列与k-mer序列对应的质心特征值异步返回至CPU。

具体的，所述聚类算法包括K均值聚类、均值漂移聚类和基于高斯混合模型的最大期望聚类算法。其中，K均值聚类(Kmeans)是优选算法。

具体的，所述CUDA流还包括内存复制指令，内存复制指令用于CPU与GPU之间的数据传输，并降低GPU内存或cache。

具体的，所述内存复制指令包括同步内存复制指令和异步内存复制指令。其中，同步内存复制指令用于从CPU内存中同步向GPU传输基因组格式文件，提高计算速度；异步内存复制指令用于将GPU处理的k-mer序列及对应的质心特征值异步返回至CPU，减少GPU浪费和内存开销。

一种基于聚类算法的高性能k-mer频次计数系统，采用上述的基于Kmeans的基于聚类算法的高性能k-mer频次计数方法实现，系统包括k-mer预处理模块、坐标转换模块、Kmeans计算模块和k-mer频次统计模块。其中，

K-mer预处理模块用于通过CPU从磁盘读取测序文件，通过CUDA流将其传输至GPU，并利用GPU按照预设子序列长度K将测序文件拆分为多段k-mer序列；

坐标转换模块用于以k-mer序列中每个碱基位置为坐标点，将拆分获得的多段k-mer序列分别转换为对应的坐标数组并添加坐标偏移；

Kmeans计算模块用于利用Kmeans聚类计算出每个坐标数组的质心特征值，并将所有k-mer序列与k-mer序列对应的质心特征值返回至CPU；

k-mer频次统计模块用于通过CPU根据k-mer序列的对应的质心特征值进行计数，统计出测序文件的所有k-mer序列频次。

具体的，坐标转换模块将拆分获得的多段k-mer序列分别转换为对应的坐标数组并添加坐标偏移处理过程具体为：坐标转换模块在ATCG笛卡尔坐标系中，以A、T为垂直坐标的正负轴，G、C为水平坐标的正负轴，并以k-mer序列中每个碱基位置为坐标点，将对拆分后得到的k-mer序列中的每一个碱基都表达为坐标系中的坐标点，从而将每段k-mer序列转换为ATCG笛卡尔坐标系中对应的坐标数组；最后坐标转换模块基于火星坐标系的坐标偏移算法，在每个坐标数组各个碱基的坐标值中添加坐标偏移，将坐标点转换为浮点型。

具体的，Kmeans计算模块利用Kmeans聚类计算出每个坐标数组的质心特征值，并将所有k-mer序列与k-mer序列对应的质心特征值返回至CPU的处理过程具体为：Kmeans计算模块先以碱基坐标点之间欧几里得距离的均方差作为准则函数，对每个坐标数组进行Kmeans聚类计算分析，获得每个坐标数组的质心特征值，并将其作为每个坐标数组对应的k-mer序列质心表达；然后Kmeans计算模块将测序文件中reads序列的所有k-mer序列都转换为质心表达后，将所有k-mer序列与k-mer序列对应的质心特征值返回至CPU。

本发明的有益效果：

1.本发明创造性地提出了一种基于Kmeans的高性能CUDA异构异步算法，完成k-mer序列的频次统计。相较于以往的统计算法，无监督机器学习与异构编程结合的方法在保证准确率的基础上大大地提高了运算速度。相较于传统的滑窗算法获取k-mer序列，本方法直接拆分基因序列，便于GPU计算；同时用坐标点转换、Kmeans求取质心来取代Hash值表达k-mer，避免了传统算法中可能出现的哈希冲突；采取CUDA流形式传输数据，降低内存开销，避免了jellyfish等算法内存不足的问题；使用CUDA异构编程，替代jellyfish算法中使用CAS操作提高并行的方法，采取GPU加速的策略，协调CPU、GPU共同运算，提高并行速度。总而言之，本方法是一种集高精准度、低算力需求、高效率于一体的k-mer计数方案，大幅提高了涉及k-mer分析的生物信息学分析流程的速度，有助于生信分析人员以更快的速度、更短的时间获得准确的分析结果。

2.依托于基因测序技术以及生物信息学等科学领域的快速发展，本方法为低成本、快速高效的基因组k-mer分析提供了有力保障，在真正意义上实现了以计算机科学技术为工具，对生物信息进行储存、检索与分析。同时，生物信息学是21世纪自然科学的核心领域之一，本方法作为生物信息学中DNA-seq分析的关键技术，在基因组学中具有极高技术优势和市场价值，适合于基因分析的实际应用与技术上的进一步推广。

附图说明

图1是本发明的方法步骤流程图；

图2是系统整体框架图；

图3是k-mer唯一化表达流程图；

图4是系统程序工作流程图。

具体实施方式

为了对本发明的技术特征、目的和有益效果有更加清楚的理解，现对本发明的技术方案精选以下详细说明。显然，所描述的实施案例是本发明一部分实施例，而不是全部实施例，不能理解为对本发明可实施范围的限定。基于本发明的实施例，本领域普通技术人员在没有作出创造性劳动前提下所获得的其他所有实施例，都属于本发明的保护范围。

实施例一：

本实施例中，如图1所示，一种基于聚类算法的高性能k-mer频次计数方法，包括以下步骤：

CPU读取基因组格式文件，得到长度为L的基因序列reads；

CPU创建CUDA流，所述CUDA流包括以下GPU核函数及其执行的操作：

S1：按照预设子序列长度K将reads拆分为L-K+1段k-mer序列；

S2：以k-mer序列中每个碱基位置为坐标点，将拆分获得的多段k-mer序列分别转换为对应的坐标数组并添加坐标偏移；

S3：利用聚类算法计算出每个坐标数组的质心特征值，并将所有k-mer序列与k-mer序列对应的质心特征值返回至CPU；

CPU根据k-mer序列及对应的质心特征值，对所有k-mer出现的频次进行计数；

CPU输出k-mer频次分布结果。

由于受目前测序水平的限制，基因组测序时需要先将基因组打断成DNA片段，然后再建库测序。reads(读长)指的是测序仪单次测序所得到的碱基序列，也就是一连串的ATCGGGTA之类的，它不是基因组中的组成。不同的测序仪器,reads长度不一样。对整个基因组进行测序，就会产生成百上千万的reads。测序得到的原始图像数据经base calling转化为序列数据，结果以fastq文件格式存储，fastq文件为用户得到的最原始文件，面存储reads的序列以及reads的测序质量。在进行高通量测序时，在芯片上的每个反应，会读出一条序列，是比较短的，叫read，它们是原始数据。

均值漂移聚类是基于滑动窗口的算法，来找到数据点的密集区域。这是一个基于质心的算法，通过将中心点的候选点更新为滑动窗口内点的均值来完成，来定位每个组/类的中心点。然后对这些候选窗口进行相似窗口进行去除，最终形成中心点集及相应的分组。该算法处理流程大致如下：

1.确定滑动窗口半径r，以随机选取的中心点C半径为r的圆形滑动窗口开始滑动。均值漂移类似一种爬山算法，在每一次迭代中向密度更高的区域移动，直到收敛。

2.每一次滑动到新的区域，计算滑动窗口内的均值来作为中心点，滑动窗口内的点的数量为窗口内的密度。在每一次移动中，窗口会想密度更高的区域移动。

3.移动窗口，计算窗口内的中心点以及窗口内的密度，知道没有方向在窗口内可以容纳更多的点，即一直移动到圆内密度不再增加为止。

4.步骤一到三会产生很多个滑动窗口，当多个滑动窗口重叠时，保留包含最多点的窗口，然后根据数据点所在的滑动窗口进行聚类。

本实施例采用均值漂移聚类算法进行聚类分析，方法的实现主要分为以下三个方面进行：

第一方面：提供了一种无监督机器学习唯一化表达k-mer序列的方法。

这一方面又分为了三个层次：

第一层次为序列拆分层。针对测序下机预处理后fasta、fastq等测试文件，摒弃传统的滑窗算法分裂基因序列，采取字符串拆分方式，若一条read的长度为L，指定子序列长度为K，则从read序列的起始位置开始，将read拆分为L-K+1段长度为K的k-mer子序列。例如一段预处理后的clean read为TGCTTACTCTGCGTTGATACCACTGCTTAGATCGGAAGAGCACACGTCTGAA，长度为52bp，指定K值为5，则可将其拆分为TGCTT、GCTTA、CTTAC……CTGAA等48段子序列。

第二层次为坐标转换层。在GPU进行数据处理的过程中，使用GPU计算的程序必须具有以下特点：需要处理的数据量比较大，数据以数组或矩阵形式有序存储，并且对这些数据要进行的处理方式基本相同，各个数据之间的依赖性或者说耦合很小，需要复杂数据结构的计算如树，图等，则不适用于使用GPU进行计算。

因此，本实施例将拆分后得到的k-mer序列，将其中的每一个碱基转换为ATCG笛卡尔坐标系中的坐标数组。设计A、T为垂直坐标的正负轴，G、C为水平坐标的正负轴。考虑碱基的位置，将k-mer序列的每一个碱基表达为坐标系中的点，使得碱基的坐标点散布在笛卡尔坐标系的坐标轴上。例如一个k-mer序列：TGCTT，T为第一个碱基，其坐标为(0，-1)，G为第二个碱基，其坐标为(2，0)，以此类推，可将第一段k-mer序列转换为{(0，-1)，(2，0)，(-3，0)，(0，-4)，(0，-5)}的坐标数组。另外，考虑到类似CTATC和TCACT、ACTTC和ATCCT、CTTCA和TCCTA这类对称或近似对称的序列会产生相同质心，进一步参考我国火星坐标系的坐标偏移算法，将各碱基的坐标转换为浮点型，从而使不同的k-mer序列都有唯一化的质心，据此进行区分。

第三层次为质心计算层。在将k-mer序列转换为坐标点后，即可运用均值漂移聚类进行质心计算，先设定滑动窗口半径r，随机在k-mer序列中选择一个碱基坐标点为中心点，开始以半径为r的圆形滑动窗口进行滑动，计算出这组坐标点的质心，其可以作为该段k-mer序列的质心特征值，唯一化表达一段k-mer序列。至此，可以将测序文件中reads的所有k-mer序列转换为质心表达，用于后续k-mer出现频次的统计。

第二方面：提供了一种CUDA流配合异步操作进行Kmeans计算的方法。

为了提高k-mer计算的速度，采用CUDA流加速，并配合异步操作，添加核函数启动、内存复制等指令，使得内存向GPU传输测序文件的同时，GPU能够分别进行基因序列的拆分、k-mer序列的预处理转换和Kmeans聚类计算，提高运算性能。

CUDA流表示一个GPU操作队列，并且该队列中的操作将以指定的顺序执行。我们可以在流中添加一些操作，如核函数启动，内存复制等。将这些操作添加到流的顺序也就是他们的执行顺序。然而，在硬件中没有流的概念，而是包含一个或多个引擎来执行内存复制操作，以及一个引擎来执行核函数。

本实施例中，CUDA流采用的是异步传输，即cudaMemcpyAsync/cudaMemcpy2DAsync/cudaMemcpy3DAsync，没有内存同步复制，当第一个文件传输至GPU后立即开始核函数运算，同时启动第二个文件的传输，即文件的传输是异步进行的，相互独立。

第三方面：提供了一种CUDA异构编程进行k-mer计数的方法。

采用前两方面所述的方法，可以使用GPU多核并行实现k-mer序列的唯一化表达。为了提高运算速度，充分发挥GPU、CPU的性能，采取CUDA异构编程，进行CPU、GPU协同计算。建立其CPU与GPU之间的多次通讯，CPU负责测序文件的读取传入，依托数据流形式将其传输给GPU，GPU采取多核并行，负责调取核函数完成基因序列分裂、坐标转换偏移、Kmeans计算等任务。计算完成后的质心点及对应的k-mer则返回到内存，依托CPU根据质心进行计数，从而完成k-mer频次的统计。

实施例二：

本实施例在实施例一的方法基础上，使用聚类算法中的Kmeans聚类分析算法替换均值漂移聚类算法进行质心计算，具体质心计算流程如下：

在将k-mer序列转换为坐标点后，即可运用Kmeans聚类，以样本点(坐标点)间欧几里得距离的均方差为准则函数，计算出这组坐标点的质心，其可以作为该段k-mer序列的质心特征值，唯一化表达一段k-mer序列。至此，可以将测序文件中reads的所有k-mer序列转换为质心表达，用于后续k-mer出现频次的统计。

其中，关于Kmeans聚类分析算法，其相关概念包括：(1)K值，即要得到的簇的个数；(2)质心：每个簇的均值向量，即向量各维取平均即可；(3)距离量度：常用欧几里得距离和余弦相似度(先标准化)。该算法的大致流程如下：

1、首先确定一个k值，即我们希望将数据集经过聚类得到k个集合。

2、从数据集中随机选择k个数据点作为质心。

3、对数据集中每一个点，计算其与每一个质心的距离(如欧式距离)，离哪个质心近，就划分到那个质心所属的集合。

4、把所有数据归好集合后，一共有k个集合。然后重新计算每个集合的质心。

5、如果新计算出来的质心和原来的质心之间的距离小于某一个设置的阈值(表示重新计算的质心的位置变化不大，趋于稳定，或者说收敛)，我们可以认为聚类已经达到期望的结果，算法终止。

6、如果新质心和原质心距离变化很大，需要迭代3～5步骤。

本实施例中，就质心特征值的聚类分析效果，将本实施例采用的K均值聚类(Kmeans聚类分析算法)同均值漂移聚类、基于高斯混合模型的最大期望聚类算法进行对比，对比结果如下表所示：

表1质心特征值的聚类分析效果比对表

由上述对比结果可以看出，均值漂移聚类基于滑动窗口，虽然也能返回中心值但是GPU运算效率低，而高斯混合模型需要返回高斯分布参数代替KMeans值，执行效率稍差。Kmeans聚类分析算法收敛速度快，当结果簇是密集的，而簇与簇之间区别明显时,它的效果较好。此外，该算法主要需要调参的参数仅仅是簇数k。可见，根据本发明的应用场景而言，采用Kmeans聚类分析算法进行质心计算的方法是本发明的优选实施例。

实施例三：

在实施例一提供的方法基础上，本发明还进一步提供了一种基于聚类算法的高性能k-mer频次计数系统，系统包括k-mer预处理模块、坐标转换模块、Kmeans计算模块和k-mer频次统计模块。其中，

K-mer预处理模块用于通过CPU从磁盘读取测序文件，通过CUDA流将其传输至GPU，并利用GPU按照预设子序列长度K将测序文件拆分为多段k-mer序列。

坐标转换模块用于以k-mer序列中每个碱基位置为坐标点，将拆分获得的多段k-mer序列分别转换为对应的坐标数组并添加坐标偏移。

Kmeans计算模块用于利用Kmeans聚类计算出每个坐标数组的质心特征值，并将所有k-mer序列与k-mer序列对应的质心特征值返回至CPU。

k-mer频次统计模块用于通过CPU根据k-mer序列的对应的质心特征值进行计数，统计出测序文件的所有k-mer序列频次。

其中，坐标转换模块将拆分获得的多段k-mer序列分别转换为对应的坐标数组并添加坐标偏移处理过程具体为：坐标转换模块在ATCG笛卡尔坐标系中，以A、T为垂直坐标的正负轴，G、C为水平坐标的正负轴，并以k-mer序列中每个碱基位置为坐标点，将对拆分后得到的k-mer序列中的每一个碱基都表达为坐标系中的坐标点，从而将每段k-mer序列转换为ATCG笛卡尔坐标系中对应的坐标数组；最后坐标转换模块基于火星坐标系的坐标偏移算法，在每个坐标数组各个碱基的坐标值中添加坐标偏移，将坐标点转换为浮点型。

Kmeans计算模块利用Kmeans聚类计算出每个坐标数组的质心特征值，并将所有k-mer序列与k-mer序列对应的质心特征值返回至CPU的处理过程具体为：Kmeans计算模块先以碱基坐标点之间欧几里得距离的均方差作为准则函数，对每个坐标数组进行Kmeans聚类计算分析，获得每个坐标数组的质心特征值，并将其作为每个坐标数组对应的k-mer序列质心表达；然后Kmeans计算模块将测序文件中reads序列的所有k-mer序列都转换为质心表达后，将所有k-mer序列与k-mer序列对应的质心特征值返回至CPU。

本实施例中，系统整体框架如图2所示，整体可分为三个方面：k-mer预处理，Kmeans计算、k-mer频次统计。

整体系统在工作时，CPU首先从磁盘中读入fasta、fastq等测序文件，以数据流的方式传输到GPU中，GPU调用函数，将reads按照指定K值拆分为k-mer序列。然后将k-mer序列的每个碱基转换成坐标，并参照火星坐标系的偏移算法调整为浮点型，之后对每段k-mer的坐标进行Kmeans计算，得出表达该段序列的质心点，返回至CPU中，对k-mer序列及出现的频次进行计数，最后输出k-mer分布结果。

本实施例中，系统主要核心功能是实现k-mer唯一化表达。如图3所示，k-mer唯一化表达大体上可以分为三步，目的是要得到能够表达不同k-mer序列的唯一化质心。第一步使用字符串拆分方式将输入的fasta、fastq文件按照指定的K值，从reads的起始位置开始，分别将reads拆分为L-K+1段(L为read读长)。第二步将得到的k-mer序列按照位置关系转换成笛卡尔坐标轴上的点，并添加必要的坐标偏移算法，将其调整为浮点型。第三步则在上一步的基础上，以坐标点之间欧几里得距离的均方差作为准则函数，依托无监督机器学习的Kmeans聚类，将每段k-mer用质心表达出来。至此，完成了一段k-mer序列到质心特征值的转换。

本实施例中，系统程序工作流程如图4所示，整体系统运行可分以下步骤：

步骤1：CPU从磁盘读入fasta、fastq文件；

步骤:2：CUDA流将文件传输至GPU；

步骤3：字符串拆分方法分裂基因序列；

步骤4：k-mer序列碱基坐标值转换；

步骤5：Kmeans聚类求取坐标数组质心；

步骤6：k-mer序列及其质心返回CPU；

步骤7：CPU对k-mer出现频次进行计数；

步骤8：输出k-mer频次分布结果。

以上显示和描述了本发明的基本原理、主要特征和本发明的优点。本行业的技术人员应该了解，本发明不受上述实施例的限制，上述实施例和说明书中描述的只是说明本发明的原理，在不脱离本发明精神和范围的前提下，本发明还会有各种变化和改进，这些变化和改进都落入要求保护的本发明范围内。本发明要求保护的范围由所附的权利要求书及其等效物界定。