CN108388772A - 一种利用文本比对分析高通量测序基因表达水平的方法 - Google Patents

Abstract

本发明属于生物信息学领域，提供了一种分析高通量测序序列基因表达量的方法，首先对测序序列进行编码、打散、随机组合，选取其中10万条序列作为查询序列分别与100万条序列进行比对，从每个查询序列随机选取9组20mer比对100万条序列去重后获得该序列的转录数量。利用查询序列首尾20mer从匹配的比对重叠群进行拼接。合并所有查询序列组比对表达数量即得到了该拼接序列的表达量，与用互补链进行比对得到的负链表达量相当。本方法可以有效用于高通量测序基因表达量及进行序列从头组装等分析。

Description

一种利用文本比对分析高通量测序基因表达水平的方法

技术领域

本发明属于生物信息学领域，涉及一种利用开放源代码操作系统命令行进行文本匹配，对高通量测序获得的短核苷酸序列进行相似性比对，并拼接匹配的重叠序列群，用于分析生物个体组织内基因表达水平的分析方法。

背景技术

高通量测序技术对数百万个DNA分子进行同时测序，使得对一个物种或者样品中的转录组和基因组进行细致全貌的分析成为可能。目前，常用的高通量测序技术主要Roche/454、ABI/SOLID测序技术、Illumina/Solexa测序技术、单分子测序技术及IonTorrrent等测序技术。RNA-Seq高通量测序，又称为转录组测序，是从2008年之后迅速开展的高通量测序方法，首先提取RNA并进行片段化后进行逆转录得到其互补DNA(cDNA)，双链cDNA加上接头并进行PCR扩增，之后进行高通量测序，每次可以获得百万条以上短序列(通常50～300个核苷酸)，从获得的cDNA序列可以了解各个基因的转录表达情况。因此，确切获得基因组或转录组信息是高通量测序的首要任务，进行DNA/RNA序列比对分析工具很多，多数需要复杂的算法来完成。其中，BLAST(Basic Local Alignment Search Tool)是由NCBI(National Center for Biotechnology Information)开发的序列相似搜索程序，是最常用的序列比对及基因注释软件。BLAST在序列数据库中快速查找与给定序列具有最优局部比对匹配的序列的算法，需要借助打分矩阵获得比对打分来确定比对序列之间(包括搜索目标双链序列)的相似程度。

目前，已有的通过高通量测序测定基因的表达水平的方法有CLC,Trinity,SOAP,Oases,ABySS,NextGENe,TopHAT,RSEM,eXpress,Sailfish,kallisto,NURD等，这些方法仍在不断改进，每种方法各有特点，算法原理各异，不同方法得到的结果显然是有差异的(同一算法设置参数不同其结果也相差很大)，因此，开发适合分析高通量测序基因表达水平的方法仍然是必要的。分析工具应具有快速方便，算法简单，序列分析过程可追溯及后续深入分析，本发明基于从海量序列数据中选取合适数据大小，利用序列文本查找目标数据集进行匹配获得相似序列重叠群，然后进行拼接延伸得到较长转录序列，重叠群里面的序列仍然可以进行序列差异分析(包括选择性剪接等)。

发明内容

本发明的目的在于提供一种简便高效的分析高通量测序基因表达水平的方法。

本发明的目的还在于提供一种可以同时分析DNA双链高通量测序序列比对及拼接的方法。

本发明的目的还在于提供一种系统兼容的分析高通量测序序列比对及拼接的方法和系统。

本发明运用开放源代码计算机操作系统的文本处理命令进行组合实现相似序列之间的比对，并进一步拼接获得较长的基因序列，合并各个比对相似序列数量即获得该基因序列的表达水平。高通量测序过程进行了核酸片段化，目前，基本上均对获得的cDNA进行双端测序，如果是同一cDNA分子，其中的片段被重复测序、不同部位被测序或者只是测互补链序列均是可能的，理论上双链被测序的概率几乎相同。假如用某一序列去查找高通量测序序列中的相似序列，显然直接用完全匹配进行比对更直接。本发明即是基于这一前提，利用开放源代码的文本处理命令直接进行序列匹配查找获得序列之间的比对相似性，结果自然也可以区分分别以双链互补序列查找获得其对应的匹配比对。本方法不需要借助参考基因组序列即可以获得匹配重叠序列的拼接(从头组装)，进行联网BLAST比对可以获得该拼接序列(基因)的注释，进一步合并重叠群的所有转录数从而获得组织内基因转录表达水平。

本发明的目的通过下述技术方案实现：

一方面，本发明提供了一种分析高通量测序基因表达水平的方法，其包含以下步骤：

S1.获取待分析样本的高通量测序信息；

作为一种示范性的实施方式，所述的待分析样本为山茶花植物花期叶片及花蕾里面的花瓣组织总RNA提取物，并进一步纯化得到mRNA后进行高通量测序。当然，本发明的方法的待分析样本也可以包括动物、植物、微生物的DNA或RNA提取物，或者是环境样品如水、大气、土壤中的微生物的DNA或RNA提取物。

其中，所述的高通量测序信息可以通过Roche的454测序平台、Illumina的Hiseq/Miseq测序平台和ABI的SOLiD测序平台中的一种获得。

作为本发明一种示范性的实施方式，所述的待测样本为RNA提取物，并采用Illumina HiSeq的测序平台获取待分析样本的高通量测序信息。Illumina HiSeq高通量测序每次可获得双端测序序列近5000万条、单端约2500万条的150bp短序列。Illumina HiSeq测序前进行了RNA片段化，要从测序获得的短序列(本发明示范性的实施例中用IlluminaHiSeq 4000平台测序获得的150bp短序列)得到组织内基因转录表达水平就需要将重复测序的片段进行比对，同时拼接同一序列片段化后的多段序列，获得较长的基因序列，然后可以获知测序序列群体内基因表达水平。许多工具对此问题进行了探讨，而其算法较为复杂，操作系统程序兼容性不好，计算过程无法追溯。

S2.序列加编号，打散，随机组合。

具体地，序列加编号为提取步骤S1中获取的高通量测序序列，仅保留序列，每条序列加上编号，例如，有2500万条序列(单边测序序列)，则从第1条序列到第2500万条序列编号为00000001-25000000。当然，所述的编号也可以采用其他的方式进行。

编号完成后，则对序列进行打散和随机组合。作为一种优选的实施方式：本发明采用逐级按照每10万→100万→5万条切割序列文档后，随机排序再合并所有序列为一个文档。其中，按照100万条切割的每个100万条序列文档分成若干个目录，单独按照进行每1万条的随机排序后随机合并得到100万条序列，所有目录得到的文档再随机组合合并所有序列(25000000条)。然后按照每100万条切割随机排序后的文档，用于进行以下步骤的分析。

步骤S2的目的在于将测序的序列进行充分的随机分布，此后的分析过程则采用部分抽样进行分析，大大降低分析样本数量，减少分析工作量。

S3.选取上述步骤S2中的10万条序列作为查询序列，将10万条查询序列与100万条序列进行比对后统计查询序列表达量。

具体地，随机选取上述步骤S2中打散的100万条序列，再把这100万条序列按照每10万条分割，随机选取一个10万条序列作为查询序列。

在上述的10万条序列中，对于每一条序列，均进行以下操作：

1.每隔5个核苷酸，取20个连续的核苷酸序列(20mer)，将每条查询序列分为多个20mer的短片段。

2.在每条序列中的多个20mer的短片段中，随机选取至少9个20mer的短序列。作为优选的实施方式，随机选择9个20mer的短序列；选择的序列越少，则比对工作量越少。经本发明研究发现，随机选取9个20mer的短序列，已经完全满足分析需求。注：因每隔5个核苷酸，选取20个连续的核苷酸序列，所以连续的20mer的短序列之间会有序列的重叠。

3.采用随机选取的至少9个20mer的短序列，与100万条序列匹配比对；同时，获得9个20mer序列的互补链亦进行与100万条序列匹配比对，统计匹配序列编号去除重复编号，获得每条序列与100万条匹配后的数量。

4.分别计算两次(正链与互补链)比对的统计数，即可初步获得每个查询序列与100万条序列比对的匹配序列，即匹配该条查询序列的正链与互补链的表达量(每100万条序列转录物数量)。

作为示范性的实施方式，如采用Illumina的Hiseq/Miseq测序平台，其测得的每一条序列是150mer，则如果每隔5个核苷酸，取20个连续的核苷酸序列(20mer)，每条序列则分为27个20mer的短序列，在这27个20mer的短序列中，再随机选取9个20mer的短序列，与100万条序列匹配比对，统计匹配序列编号去除重复编号，获得每条序列与100万条匹配后的数量。

S4.将步骤S3中，10万条查询序列与100万条序列比对得到的序列进行重叠群拼接。

具体包括以下步骤：

1.利用查询序列的5'和3'两端20mer序列进行匹配，切除原始查询序列后找出剩下最长的序列进行组装拼接；

2.将重叠群(接近10万条，部分序列可能没有获得重叠群)相互比对，去除重复的重叠群；

3.将获得的去除重复后的所有重叠群序列再次进行至少9x20mer随机片段比对；

4.比对结束后，匹配的序列组，先在序列组内按照序列号的大小进行排序，然后在不同的序列组间，按组内第一个序列号相同的序列组合并再比对的重叠群；

5.分别计算单独的重叠群和2个及以上的重叠群，找出比对的表达量进行合并得到重叠群的表达量。

S5.重叠群序列的比对以及注释提取；

具体地，将组装拼接得到重叠群序列与美国国立生物技术信息中心(NCBI，http://www.ncbi.nlm.nih.gov/)远程BLAST服务器核酸数据库进行比对获得该序列的比对及注释信息。

优选地，每次上传500条进行比对，保存比对结果进行注释提取。

通常，每条重叠群序列进行远程BLAST比对不一定返回有效比对结果也就得不到其注释，获得注释的序列根据前面步骤分析获得的该序列对应的重叠群的表达量即可获得所分析组织内该基因的转录表达水平。

作为优选的实施方式，步骤S1-S5，均采用开放源代码操作系统进行操作。本发明采用类Unix系统中的Unix命令行完成上述的操作。

另一方面，本发明还提供了一种分析高通量测序基因表达水平的系统，该系统包括：

1)高通量测序系统：所述的高通量测序系统用于获取待分析样本的高通量测序信息；

2)数据处理系统：所述的数据处理系统用于实现上述分析高通量测序基因表达水平的方法的步骤；

作为一种优选的实施方式，所述的数据处理系统为开放源代码操作系统。

具体地，所述的开放源代码操作系统可以选自类Unix系统。作为优选的实施方式，本发明采用类Unix系统中的Unix命令行完成上述的操作。

类Unix系统(Unix-like System):指各种传统的Unix系统(比如FreeBSD、OpenBSD、SUN公司的Solaris)以及各种与传统Unix类似的系统(例如Minix、Linux等)。它们虽然有的是自由软件，有的是商业软件，但都相当程度地继承了原始UNIX的特性，有许多相似处，并且都在一定程度上遵守POSIX规范。作为一种示范性的实施方式，本发明的开放源代码操作系统选自FreeBSD(版本11.0)。

另一方面，本发明还提供了一种计算机可读介质，其上存储有计算机指令，这些指令被处理执行时实现上述分析高通量测序基因表达水平的方法的步骤。

本发明相对于现有技术具有如下的优点及效果：

本发明的主要优点在于简化高通量测序序列比对及拼接过程的复杂性，高效快速获得组织高通量测序的基因表达情况，为进一步开展基因表达调控分析打下重要基础。具体有以下优势：

1.序列随机排序充分：序列多次进行打散随机排序，从获得的重叠群其表达量可以看出随机排序效果显著，从100万条序列进行比对和拼接就能够获得基因表达水平，大大减少了分析的工作量；

2.序列进行文本比对、拼接，方法简单、快捷、透明、高效：上述提及的比对方法直接明了，利用每条查询序列每隔5个碱基选取20mer的短序列，之后再从中随机选用几个短序列(例如9x20mer)与100万目标序列比对，比对过程简单快速，每一步骤透明可控。用查询序列首尾20mer进行匹配重叠群的拼接也简单快捷，分多个窗口可显著提高拼接速度，拼接长度可以随意控制，从实例看到经过两次拼接可以获得超过1000mer的序列。

3.可以同时分析cDNA双链表达水平：用上述的至少9x20mer的互补链进行比对可以获得查询序列的互补链表达量，从实例可以看出两者的表达量相当，与基因表达实际情况相一致，也说明通常不需要进行互补链的表达量分析；

4.本方法利用开源系统无程序兼容性问题：本方法基于开放源代码操作系统，因而本发明所使用的文本处理命令组合在这些系统上不存在兼容性问题。

附图说明

图1是本发明的序列分析流程图。

图2品种1叶片10万条序列与1个100万条序列(包含自身10万条)部分序列比。

图3品种1叶片10万条序列与1个100万条序列比对结果排序去重复后匹配序列。

图4.叶片10万条序列与100万条序列比对匹配序列第一次拼接部分序列。

图5.叶片10万条序列与100万条序列比对匹配序列第二次拼接部分序列。

图6第二次拼接得到的单一序列表达量大于等于100的正负链表达量柱状图。

图7第二次拼接得到含有2个以上匹配序列组合并后最高表达的50个序列的正负链表达量柱状图。

具体实施方式

下面结合实施例及附图对本发明作进一步详细的描述，但本发明的实施方式不限于此。

重叠群(contig)：通过重叠部分将相邻reads组装形成的单元称为contigs。转录组测序的原始测序包含了很多的reads(一次测序中仪器读取的核苷酸长度)，通过序列的拼接，具有重叠区的reads会被组装成更大的片段，称之为contig；；

单独的重叠群：重叠群片段拼接后再次进行比对没有再找到匹配片段的重叠群；

两个以上的重叠群：重叠群片段拼接后再次进行比对可以再找到两个或以上匹配片段的重叠群

实施例1 山茶花高通量测序

山茶花开花期花枝成熟全展叶片和未打开花苞的花瓣，由公司提供测序服务，包括总RNA提取及建库,双端测序(Paired-End,Illumina HiSeq 4000)。序列格式fastq，提交6G高质量数据(clean data)及7G未处理原始测序数据(raw data)，每条序列长度150mer，合并双端得到每个样品测序序列约5000万条。

实施例2 测序序列加编号，打散，随机组合

提取实施1中获取的高通量测序序列，仅保留序列，每条序列加上编号本次测序有5000万条(双端测序各分别产生约2500万条序列)，一端的测序从第1条序列到第2500万条序列编号为00000001-25000000。然后采用逐级按照每10万→100万→5万条进行随机排序并按随机方式合并序列切割序列文档后随机排序再合并所有序列为一个文档。其中，按照100万条切割的每个100万条序列文档分成若干个目录，单独按照进行每1万条的随机排序后随机合并得到100万条序列，所有目录得到的文档再随机组合合并所有序列(25000000条)。然后按照每100万条切割随机排序后的文档，用于进行以下分析过程。其中，每个100万分成若干个目录单独进行每1万条的随机排序，再随机组合合并2500万条序列。

实施例3 在100万条序列中，选取10万条作为查询序列进行比对，得出每条查询序列的表达量

随机选取上述步骤实施例2中的打散的100万条序列，在把这100万条序列按照每10万条分割，随机选取一个10万条序列作为查询序列。

在上述的10万条查询序列中，对于每一条序列，进行以下操作：

1.每隔5个核苷酸，取20个连续的核苷酸序列(20mer)，每条查询序列可分为27x20mer和短序列；

2.在每条查询序列中，随机选取9个20mer的短序列；

3.采用随机选取的至少9个20mer的短序列，与100万条序列匹配比对，同时，至少9个20mer短序列的互补链亦进行与100万条序列匹配比对，统计匹配序列编号去除重复编号，获得每条序列与100万条匹配后的数量；

4.统计与正链与互补链查询序列比对的匹配序列数即获得每个150mer序列与100万条序列比对的表达量；

作为示范，以下表1随机示出了叶片的24条查询序列与100万条序列比对的匹配序列数量。

表1.叶片24条查询序列的20mer与1个叶片100万条序列比对的匹配序列数量

注：A:27个20mer；B～D:27个20mer序列从左至右的9个20mer；E～I:27个20mer分5次随机选取9个20mer。

从表1可以看出，用9个随机的20mer短序列与27个所有的20mer序列比对得到的比对匹配数量很一致(A,E～I)，而不进行随机选取的3个9个20mer的匹配数量与27个20mer的差异较明显(A,B～D)。因此，为了加快比对进程，后面的比对结果全部选取了150mer序列中随机选取9个20mer进行比对。

表2.24条序列9个随机20mer与多个100万条及一个300万条序列比对的匹配序列数量

注：RNA-Seq序列，A～C:分别为不同100万条品种1叶片的序列；D～300万条品种1叶片的序列；E：品种2叶片100万条的序列；F、G：品种1花瓣100万条(F-双端读1/G-读2)。

表2用同样的叶片(品种1)测序的24条查询序列与其叶片测序的3个100万条(A/B/C)、1个300万条(D)和品种2的1个叶片100万条(E),及2个品种1花瓣的100万条(F-双端读1/G-读2)进行比对，每个100万条是随机选取的，300万条序列只选了一个。从表2可以看到，24条查询序列与三个100万品种1叶片序列(A、B、C)的匹配数量一致很高，其中，较高表达的24条序列与300万条的匹配数量也成倍增加，低水平表达的则不然，也说明基因表达并不是简单的线性增加。初步也可以看到两个品种的叶片表达数量是相似的，有些序列与花瓣则有明显的差异。从表1和表2结果可以看出，用100万条序列进行比对可反映高通量测序基因表达水平。

品种1叶片10万条查询序列与包含自身序列的叶片100万条进行比对获得的部分(9个)比对结果如图2，比对过程大约一周(每天约1.5万条，可同时运行6个相同任务计算机运行其他程序不受影响)，对比对结果进行排序去重复(剩下92025条)，以及去除没有匹配的(9701条)，得到90298条匹配(排序见图3)，最高匹配序列6615个，大于1000个的2371个，小于1000的87927个。

实施例4 重叠群拼接

将实施例3中，10万条查询序列与100万条序列比对得到的序列进行重叠群拼接；具体地，将150mer查询序列首尾20mer序列作为匹配序列与找到的每一行内各自的匹配序列进行拼接，连接后得到加长的拼接序列，得到90298条拼接序列，拼接过程运行时间较长约10天(按照每个窗口1万条分别进行拼接时间只需要3天)。为了进一步合并拼接序列，去除了含有10个连续polyA/C/G/T的拼接序列(测序序列预处理的时候通常会去除的mRNA的尾部序列，如果多个序列出现超出20个A/T/G/C也会拼接的，去除是为了不影响拼接)，其中没有得到重叠群的序列有7912条，序列超过150mer的有77954条(重叠群序列见图4)，拼接长度最长410mer。

将获得的重叠群序列(77954条)进行再次相互比对，即所有序列进行相互间比对，每个序列仍然进行9x20mer短序列与所有拼接序列比对，没有得到匹配的序列有16478条，得到超过1条以上重叠群的有37256条，合并第一个序列号后为16754条，将这些序列选取第一个序列头尾20mer进行拼接得到新的重叠群，长度最长为1174mer(部分序列截图见图5)，长度超过500mer的有9729条。因为是按照第一个序列号组合的，新的16754条重叠群仍然有匹配的序列群，所以需要进一步拼接，按照第一次拼接的方法用获得的拼接序列进行20mer分割随机选取9个去找到匹配的重叠群，再次拼接后得到10629个不再出现匹配，6125个仍具有一个以上的匹配群序列，而且最多匹配序列号只有218个。

上面的拼接中，得到的单一拼接序列包括再次拼接得到的16478条及将含有2个以上的16754条合并第一个序列号再比对匹配拼接后得到的10629条，仍有6125条序列含有2条及以上的匹配序列，如果进一步合并可以去除重复的匹配序列。因此，从10万条序列与100万条序列的比对匹配后拼接可以大致获得3万条左右的拼接序列。利用这些拼接序列的匹配序列按照拼接过程合并其表达数量，就得到它们的表达量计数。由于测序不排除正负链，所以计算表达量应该将正负链的表达进行计数。下面介绍得到的拼接序列的表达量。

第二次拼接得到的16478条单一拼接序列正负链的表达量相当，最高3369条，大于100的有53条，而大于1000的只有16条(即千分之一)，小于等于10的最多有14762条。限于篇幅，图6列出了大于100的所有序列正负链表达量。

第二次拼接得到的含有2条以上匹配序列的16754条序列作为拼接序列合并其表达量，正负链也是相当的，而表达量普遍较高，大于1000的有547条，大于等于10000有20条，大于100的有3040，小于等于100的有13739条，小于等于10的4204条。图7也列出了前50条序列的表达量。表达量较高是因为合并了多个匹配序列组，如图7的第一个序列(正链表达量为25520)拼接合并了86组首序列号相同的序列组，其中单组表达序列数大于500的就有58组(最高为758个)。

实施例5 重叠群序列的比对以及注释提取

单一拼接序列与含有2个以上拼接序列组的拼接序列其是否真实反映了某个基因，需要进一步利用该序列与NCBI核酸数据库进行blast比对，从返回比对序列及其注释可以初步推知。

每次用500条序列(fasta格式)与NCBI的blastn服务器进行比对，保存text格式的比对结果，提取得到了所有拼接序列的blast比对及其注释，完成所有3万条拼接序列的远程blalstn比对只需要半天。blastn程序只是返回比对打分大于70％以上的比对，前面的两组拼接序列中再次拼接得到的16478条单一序列blastn比对返回8889条相似打分(Identities)高于70％的结果，而再次拼接大于多于两条拼接序列的16754条序列进行blastn比对返回12255条，说明较长的序列可以得到更多匹配的blastn比对。

上述实施例为本发明较佳的实施方式，但本发明的实施方式并不受上述实施例的限制，其他的任何未背离本发明的精神实质与原理下所作的改变、修饰、替代、组合、简化，均应为等效的置换方式，都包含在本发明的保护范围之内。

Claims (10)

1.一种分析高通量测序基因表达水平的方法，其特征在于，包含以下步骤：

S1.获取待分析样本的高通量测序信息；

S2.序列加编号，打散，随机组合；

S3.选取上述步骤S2中的10万条序列作为查询序列，将10万条查询序列与100万条序列进行比对后统计查询序列表达量；

S4.将步骤S3中，10万条查询序列与100万条序列比对得到的序列进行重叠群拼接；

S5.重叠群序列的比对以及注释提取。

2.根据权利要求1所述的方法，其特征在于，步骤S1中，所述的待分析样本为动物、植物、微生物的DNA或RNA提取物；或者，所述的待分析样本为水、大气、土壤中的微生物的DNA或RNA提取物。

3.根据权利要求1所述的方法，其特征在于，步骤S1中，所述的高通量测序信息采用Roche的454测序平台、Illumina的Hiseq/Miseq测序平台和ABI的SOLiD测序平台中的一种获得；更优选Illumina的Hiseq测序平台。

4.根据权利要求1所述的方法，其特征在于，步骤S2中，打散，随机组合的方法为：采用逐级按照每10万→100万→5万条进行随机排序并按随机方式合并序列，其中，每个100万分成若干个目录单独进行每1万条的随机排序，再随机组合合并所有序列。

5.根据权利要求1所述的方法，其特征在于，步骤S3中，随机选取上述步骤S2中打散的100万条序列，在100万条序列中随机选取10万条序列作为查询序列，针对每一条查询序列进行以下操作：

1.每隔5个核苷酸，取20个连续的核苷酸序列，将每条查询序列分为多个20mer的短片段；

2.在每条查询序列的多个20mer的短片段中，随机选取至少9个20mer的短序列；优选地，随机选择9个20mer的短序列；

3.采用随机选取的至少9个20mer短序列，与100万条序列匹配比对，同时，至少9个20mer短序列的互补链亦进行与100万条序列匹配比对，统计匹配序列编号去除重复编号，获得每条序列与100万条匹配后的数量；

4.分别计算正链与互补链比对的统计数即获得每个查询序列与100万条序列比对的匹配序列，即匹配该条查询序列的正链与互补链的表达量。

6.根据权利要求1所述的方法，其特征在于，步骤S4中，其具体步骤如下：

1.利用查询序列的5'和3'两端20mer序列进行匹配，切除原始序列后找出剩下最长的序列进行组装拼接从而得到重叠群；

2.将重叠群相互比对，去除重复的重叠群；

3.再次进行至少9x20mer随机片段比对；

4.对匹配的序列组，先在序列组内按照序列号的大小进行排序，然后在不同的序列组间，按组内第一个序列号相同的序列组合并重叠群；

5.分别计算单独的重叠群和2个及以上的重叠群，找出比对的表达量进行合并得到重叠群的表达量。

7.根据权利要求1所述的方法，其特征在于，步骤S5中，将步骤S4组装拼接得到的重叠群序列与NCBI的BLAST服务器核酸数据库进行比对获得序列的比对及注释信息。

8.根据权利要求1-7任一所述的方法，其特征在于，采用采用开放源代码操作系统进行操作；

优选地，所述的开放源代码操作系统选自类Unix系统；更优选地，所述的开放源代码操作系统选自FreeBSD、OpenBSD、Solaris、Minix或Linux；更优选地，所述的开放源代码操作系统选自FreeBSD。

9.一种分析高通量测序基因表达水平的系统，其特征在于，所述的系统包括：

1)高通量测序系统：用于获取待分析样本的高通量测序信息；

2)数据处理系统：所述的数据处理系统用于实现权利要求1-7所述的方法的步骤；

优选地，所述的数据处理系统为开放源代码操作系统；

更优选地，所述的开放源代码操作系统选自类Unix系统；

还更优选地，所述的开放源代码操作系统选自FreeBSD、OpenBSD、Solaris、Minix或Linux尤其优选地，所述的开放源代码操作系统选自FreeBSD。

10.一种计算机可读介质，其上存储有计算机指令，其特征在于，该指令被处理执行时实现权利要求1-7所述的方法的步骤。