本发明涉及基因测序数据的一种压缩方式，特别涉及基因测序数据的压缩方法。
背景技术：
高通量测序技术(hts)的发展和飞速进步使得将基因组信息用作多个领域的日常实践成为可能。例如，随着最新高通量测序仪器的发布，人类全基因组测序(wgs)的成本已经下降到仅1,000美元。很快测序成本就会下降到大约100美元。这些在降低测序成本方面取得的成就打开了个性化医疗的大门，使得对患者的基因组信息进行测序和分析可以像如今的标准血液检测一样频繁。然而，由于存储和处理数据相关的it成本，测序数据的不断增长已然严重阻碍了公共卫生测序更广泛地传播。
时至今日，一个单一的测序系统每年可以传送超过18,000个人类全基因组，按照目前的测序能力，每年产生近乎5pb的数据(1pb＝2^50b)，因此，随着基因测序技术的普及不断增加，基因组数据的高效存储和传输将变得至关重要。然而，缺乏适当的表示方法和有效的压缩技术严重限制了基因组数据在科学和公共卫生目的潜力的。
在信息论里面，数据压缩是指在不丢失有用信息的前提下，通过按照一定的算法对数据进行重新组织以缩减其数据量，从而减少存储空间，提高其传输、存储和处理效率的技术。比较常用的压缩算法包括广泛用于文本或通用数据压缩的lempel-ziv(lz)压缩方法及其变种，用于图像压缩的gif,jpeg算法，以及用于音视频压缩的mpegaudio/video压缩算法等。然而，由于基因测序数据的特性和普通的文本或音视频数据较大差异，直接使用这些通用的压缩算法对基因测序数据进行直接压缩的效果并不理想，故而需要根据基因测试数据的格式和统计特性，研究更有效的压缩技术。
目前，基因测序数据一般用文本格式进行存储。常用的格式包含fastq格式和sam格式。其中，fastq用于存储未经定序(alignment)的原始基因测试数据，sam用于存储经过定序(alignment)的基因测序数据，下面分别对fastq于sam格式进行详细介绍。
fastq：
fastq是基于文本的，保存生物序列(通常是核酸序列)和其测序质量信息的标准格式。其序列以及质量信息都是使用一个ascii字符标示，最初由sanger开发，目的是将fasta序列与质量数据放到一起，目前已经成为高通量测序结果的事实标准。
fastq文件中每个序列通常有如下四行：
序列标识以及相关的描述信息(readname)，以‘@’开头；
第二行是序列片段的序列信息(seq)；
第三行以‘+’开头，后面是序列标示符、描述信息，或者什么也不加；
第四行，是质量信息(qual)，与第二行的序列相对应，每一个序列都有一个质量评分，根据评分体系的不同，每个字符的含义表示的数字也不相同。
fastq文件序列例子：
@@hwusi-eas100r:6:73:941:1973#0/1
gatttggggttcaaagcagtatcgatcaaatagtaaatccatttgttcaactcacagttt
+
！”*((((***+))％％％++)(％％％％).1***-+*”))**55ccf>>>>>>ccccccc65
其中，序列标识的格式一般由测序仪生产厂商制定。比如上述例子的序列标识根据illumina格式的解释如下：
序列标识：@hwusi-eas100r:6:73:941:1973#0/1，其中：
质量编码：指的是一个碱基的错误概率的对数值。其关系如下：
q＝-10log10p
其中，p是碱基的错误概率，q是对应的质量编码值，其最初在phred拼接软件中定义与使用，其后在许多软件中得到使用。
对于每个碱基的质量编码标示，不同的测序厂商采用不同的方案，目前有5种方案：
sanger，质量编码值的范围从0到93，对应的ascii码从33到126，但是对于测序数据(rawreaddata)质量得分通常小于60，序列拼接或者mapping可能用到更大的分数。
solexa/illumina1.0,质量编码值的范围从-5到63，对应的ascii码从59到126，对于测序数据，得分一般在-5到40之间；
illumina1.3+，质量编码值的范围从0到62，对应的ascii码从64到126，(原始测序数据一般得分在0到40之间)；
illumina1.5+，其中0到2作为另外的标示。
illumina1.8+基本返回sanger格式。ascii码＝质量编码值+33。
sam格式：
sam与fastq文件格式类似，sam文件格式也使用文本文件来表示基因测序的数据。一般情况下，sam格式用于表示经过aligner定序软件定序处理后的基因测序数据。因此，同fastq文件相比，sam格式增加了每个序列的定位信息。sam输出的结果中每一行都包括十二项通过tab分隔(\t)，从左到右分别是：
qname：序列标识，一般和fastq文件相同；
flag：概括出一个合适的标记，各个数字分别代表；
rname：比对上的参考序列的名字(染色体)；
pos：比对上序列在参考序列上的位置(染色体上的位置)；
mapq：比对质量，越高则比对结果约可靠；
cigar：代表比对结果的cigar字符串，如37m1d2m1i，这段字符的意思是37个匹配，1个参考序列上的删除，2个匹配，1个参考序列上的插入，m代表的是alignmentmatch(可以是错配)；
rnext：mate序列所在参考序列的名称；下一个片段比对上的参考序列的编号，没有另外的片段，这里是’*‘，同一个片段，用’＝‘；
pnext：mate序列在参考序列上的位置；下一个片段比对上的位置，如果不可用，此处为0；
tlen：估计出的片段的长度，当mate序列位于本序列上游时该值为负值。template的长度，最左边得为正，最右边的为负，中间的不用定义正负，不分区段(single-segment)的比对上，或者不可用时，此处为0；
seq：序列片段的序列信息，如果不存储此类信息，此处为’*‘，注意cigar中m/i/s/＝/x对应数字的和要等于序列长度；一般情况下该信息和fastq文件相同。
qual：ascii码格式的序列质量；序列的质量信息，格式同fastq一样。
可选的字段(field)。
sam格式中第1、10和11项的信息和fastq对应的信息格式相同。一遍情况下，aligner输出的sam文件第1、10项的信息同对应的fastq输入原始数据文件相同，qual信息可能相同，也可能会通过basequalityrecalibration算法修改而更准确地翻译序列的质量。其他项的数据为比对信息，由aligner通过对比算法增加。sam文件序列数据的一个示例如下：
hwi-st170:265:5:44:14178:183344#01451624213763m1i35m18568439490
cctgtatacatagtaatcaaagtgtaccactggtcggtgtttgtgttcaggcccctgttgggtaatgtgcatgtgaagacctcaggtggtatagttttg
cee？@f@be@ggegfbhhedeeedeebedhhbghggfhhdfhhhggggfffeeehfhfgfhhhhhfhhhfhhhhghghehhhhhhhhhfhhhhhhhhhhrg:z:du23m01_durocxt:a:unm:i:4sm:i:37am:i:37x0:i:1x1:i:0xm:i:3xo:i:1xg:i:1md:z:20t22c1a52。
目前，基因测序数据的压缩还是一个研究课题，比较常用的基因测序数据压缩格式包括bam和cram。
bam文件格式使用通用文件压缩方式gzip(gnuzip)的一种变种gbzf(blockedgnuzipformat)进行压缩。为了实现随机访问，bam文件通过bgzf压缩方式，将基因测序数据分为许多数据块，然后通过gzip算法进行压缩，压缩后的数据块大小限制为64kb。由于gzip至少一种通用的文本压缩方式，使用该方式进行基因测序数据压缩并不高效。同时，gbzf的压缩方式进一步影响了压缩率。
相比于bam，cram通过一些办法改善了基因测序数据的压缩效率。首先，cram采用一种tokenizer的做法，把基因测序数据的不同字段通过tokenizer的方式拆分为不同的子数据流，然后，对不同的子数据流，cram会首先尝试不同的压缩方式，比如run-lengthcoding，ransgolombrice，哈夫曼编码，然后选择最好的压缩方式对其进行压缩。比如上述例子的序列标识：
@hwusi-eas100r:6:73:941:1973#0/1，其中：
cram会将其通过tokenizer的办法拆分为”@hwusi-eas100r”,“:”,“6”,“:”,“73”,“:”,“941”,“:”,“1973”,“#”,“0”,“/”,”1”等多个子数据流分别进行压缩编码。
通过对基因测序文件进行分析后可以发现，在对系列标识进行tokenizer操作后，大部分的字段信息对所有的系列都是相同的。在这种情况下，编码变得非常简单，只需要简单的run-len编码就可以很高效的进行编码。然而，有些字段为比较随机的数字序列，需要更好的办法进行压缩。又由于这些数字序列的值的变化范围比较大，直接对数字序列进行熵编码比较困难，前述的cram或其他基因数据压缩算法通常通过把这些数字序列都被转化为8比特长的字符序列，再使用通用的字符压缩算法比如进行编码的方式解决此问题。然而，这样处理的缺点是转化后的字符序列没法保持原来数字序列的统计特性，导致这些方法的压缩效率往往并不理想。
对于质量数据，最直接的办法是把它当作字符串，利用通用字符串压缩算法进行编码。然而，这样做没有充分利用质量数据的上下文相关性对质量信息进行更有效的压缩。
有鉴于此，本发明人特别研制出一种优化的基因测序数据的压缩方法，本案由此产生。
技术实现要素：
本发明的目的在于提供一种基因测序数据的压缩方法，以提高压缩效率。
为了实现上述目的，本发明的技术方案如下：
基因测序数据的压缩方法，其中序列标识压缩包括的步骤如下：
将质量数据序列改写成{q,r}的形式，其中q为质量数据，r为该质量数据的重复次数；
分别用不同的统计属性对质量数据q和重复次数r进行熵编码。
进一步，采用基于上下文建模编码的方式对质量数据q和重复次数r进行编码，具体包括：
设待编码的序列为s＝{s0,s1,s2,…,sn,…}，对每个待编码的符号sn,首先选择与它的概率分布函数相关的信息作为它的上下文模型c；可选的上下文模型选择包含：
1)在sn之前出现的符号{si|n-l<i<n}；
2)符号sn的位置n；
3)和序列s分布相关其他序列；
在选择了上下文模型后，通过预扫描的办法或自适应的办法，建立不同上下文模型下sn的概率分布函数p(sn|c)；
根据p(sn|c)对sn进行熵编码。
优选的，上述方法还包括：
对系列标识进行tokenize操作；
对tokenize操作后得到的数据{xi|i＝1,…,n}进行移位处理，得到高位数据{hi＝xi>>l|i＝1,…,n}和低位数据{li&(1<<l-1}|i＝1,…,n},其中，l是一个大于或等于0的整数，代表移位量，>>和<<分别代表右移和左移操作；
对高位数据{hi|i＝1,…,n}进行熵编码；
对低位数据{li|i＝1,…,n}不编码进行直接传输。
进一步，对高位数据{hi|i＝1,…,n}进行熵编码具体包括：
使用预扫描的办法先建立起经验概率函数p(h)，
p(h)＝|{hi|hi＝h；i＝1,…,n}|/n；其中，操作符|·|表示计算集合里面元素的个数；
再使用得到的经验概率函数p(h)利用算术编码对高位数据{hi|i＝1,…,n}进行编码。
进一步，利用经验概率函数p(h)对高位数据{hi|i＝1,…,n}中的h1进行编码之后，在每次编码hi,i＝2,…,n之前，都利用之前编码过的hi-1的值对p(h)进行更新，更新方式为：
若h＝hi-1，则p’(h)＝(kp(h)+1)/(k+1),；
否则，p’(h)＝kp(h)/(k+1)；
其中，k为一个大于0的常数。
本发明具有以下优点：
一、对于质量数据，通过改写质量数据序列得到质量数据q和重复次数r，再针对q和r采用不同的统计属性进行编码，提高了压缩效率。
二、对于序列标识，通过对okenize操作后得到的数据进行位移处理，得到高位数据与低位数据，针对于高位数据进行熵编码，而低位数据(字符序列)不编码，从而在熵编码的精度限制下，充分利用了原来数据的统计特性提高压缩效率。
附图说明
图1是本发明基因测序数据的压缩方法的流程图；
图2是本发明序列标识的编码方式示意图；
图3是本发明质量数据序列编码方式示意图。
具体实施方式
如图1所示，本实施例揭示的一种基因测序数据的压缩方法，该方法可以用在fastq或sam文件中系列标识和质量数据的压缩，包括：
其中序列标识压缩包括的步骤如下：
s101、对系列标识进行tokenize操作，以将不同字段拆分为不同的子数据流；
s102、对tokenize操作后得到的数据{xi|i＝1,…,n}进行移位处理，得到高位数据{hi＝xi>>l|i＝1,…,n}和低位数据{li&(1<<l-1}|i＝1,…,n},其中，l是一个大于或等于0的整数，代表移位量，>>和<<分别代表右移和左移操作；比如对16进制数字0x7b5a,l＝16，进行移位操作后高位数据为0x7b,低位数据为0x5a；
s103、对高位数据{hi|i＝1,…,n}进行熵编码(entropycoding)；
s104、对低位数据{li|i＝1,…,n}不编码进行直接传输；
对高位数据{hi＝xi>>l|i＝1,…,n}的熵编码可以灵活地使用不同的熵编码方案进行编码，比如rans，算术编码算法等。使用熵编码最关键的问题是需要知道数据{hi|i＝1,…,n}的概率密度函数p(h)，如图2所示，一种可行的做法是使用预扫描(2-pass)的办法先建立起经验概率函数p(h)；
p(h)＝|{hi|hi＝h；i＝1,…,n}|/n，其中，操作符|·|表示计算集合里面元素的个数；
然后再使用得到的经验概率函数p(h)利用算术编码对{hi|i＝1,…,n}进行编码。在使用这个办法的时候需要同时把得到的经验概率函数p(h)放入编码后的码流，并传输给解码器，以便解码器对hi进行解码。
可使用自适应的熵编码器对hi进行编码，以进一步提高编码准确性；自适应的方式首先利用经验概率函数p(h)对高位数据{hi|i＝1,…,n}中的h1进行编码之后，在每次编码hi,i＝2,…,n之前，都利用之前编码过的hi-1的值对p(h)进行更新，更新方式为：
若h＝hi-1，则p’(h)＝(kp(h)+1)/(k+1),；
否则，p’(h)＝kp(h)/(k+1)；
其中，k为一个大于0的常数。
如图3所示，质量数据压缩的步骤如下：
s105、充分利用质量信息经常重复的特点，将质量数据序列改写成{q,r}的形式，其中q为质量数据，r为该质量数据的重复次数(如果该质量数据仅出现一次，则r＝0，以此类推)；比如以下质量数据，
！”*((((***+))％％％％％％％ccccccc65
可以改写为：
{！,0}{',1}{*,0}{(,3){*,2}{+,0}{),1}{％,6}{c,6}{6,0}{5,0}
s106、分别用不同的统计属性对质量数据q和重复次数r进行熵编码。
为了充分利用质量信息之间的相关性达到更好的编码效果，本实施例优先选择基于上下文建模编码(contextmodellingcoding)的方式对质量数据q和重复次数r进行编码，具体步骤为：
设待编码的序列为s＝{s0,s1,s2,…,sn,…}，对每个待编码的符号sn,首先选择与它的概率分布函数相关的信息作为它的上下文模型c；可选的上下文模型选择包含：
1)在sn之前出现的符号{si|n-l<i<n}；
2)符号sn的位置n；
3)和序列s分布相关其他序列；比如，如果选择sn-1，sn-2和n作为编码符号sn的context,则c为一个三元数组(sn-1,sn-2，n)；
在选择了上下文模型后，通过预扫描的办法或自适应的办法，建立不同上下文模型下sn的概率分布函数p(sn|c)；
根据p(sn|c)对sn进行熵编码。比如哈夫曼算法，算术编码算法等。
据质量数据和重复次数的分布特点和关联性，本例选择(n,qn-1,rn-1)作为编码qn的上下文，(n,qn)作为编码rn的上下文，同时，为了避免使用太多数量的上下文，可以通过量化的办法对上下文的数量进行压缩。比如对位置信息n，我们可以把每l个位置作为一个context。在这种情况下context变为(n‘,qn-1,rn-1)，其中n’＝rounding(n/l)，rounding()为取整操作。对context里面的质量信息也可以做同样的处理。
以上仅为本发明的具体实施例，并非对本发明的保护范围的限定。凡依本案的设计思路所做的等同变化，均落入本案的保护范围。