基于转录组测序数据分析及高通量GO注释理论的研究

刘粉香 杨文国 孙勤红

摘要 随着二代测序技术的快速发展，转录组测序在越来越多的动植物中完成，人们获得了大批量的转录组数据序列。如何从这些海量的序列数据中挖掘具有生物意义的信息已成为很多研究的关键所在，对未知基因的功能进行预测和注释就是其中一个重要的问题。转录组序列的功能注释是功能基因組学研究的一项重要内容，基因本体论（gene ontology，GO）注释目前是一种最重要的功能注释方式。介绍了利用生物信息学软件进行转录组测序数据分析过程，包括数据质量控制和过滤、从头拼接（De novo assembly）、同源比对以及大规模GO注释，为从事转录组测序特别是非模式植物转录组测序研究者在数据分析方面提供参考。

关键词 二代测序；转录组；从头拼接；GO注释

中图分类号 Q-3文献标识码 A文章编号 0517-6611（2018）31-0088-04

Abstract With the development of sequencing technology， the transcriptome sequencing has been completed in more and more plants.A large number of transcriptome sequence data were obtained.How to mine biologically meaningful information from these massive serial data has become the key point of many researches.Predicting and annotating the function of unknown genes is an important issue.Functional annotation of transcriptome sequences is an important part of functional genomics. Gene Ontology （GO） annotation is currently one of the most important functional annotation methods.We introduced the analysis of transcriptome sequencing data using bioinformatics software， including data quality control and filtering， De novo assembly， homology comparison and largescale GO annotation，which provided a reference for researchers engaged in transcriptome sequencing， especially nonmodel plant transcriptome sequencing in data analysis.

Key words Nextgeneration sequencing；Transcriptome；De novo assembly；GO annotation

广义上的转录组是指生物体细胞或组织在特定状态下所转录出来的所有RNA的总和，包括RNA（即mRNA）编码蛋白质和RNA（ncRNA，如rRNA、tRNA、microRNA等）非编码蛋白质；狭义上的转录组通常指所有mRNA的总和[1]。转录基因组学研究被转录的基因，是挖掘转录基因的功能基因极其重要的途径，功能基因组学研究在基因进化、遗传育种等研究中具有非常重要的意义[2]。转录组研究的技术手段大体上有EST序列构建、芯片技术和二代测序技术等。随着二代测序（next generation sequencing）技术的发展和应用，许多物种已经完成了转录组测序。早在2008年，Nagalakshmi 等[3]利用 RNA-Seq 技术进行了酵母转录组测序。近年来，越来越多的无参考基因组物种先后完成了转录组测序。2012年，Zhang 等[4]对不同发育阶段的6个麻竹花器官的转录组进行测序，并分析基因的差异表达，最后预测了81个转录因子家族在麻竹花组织发育过程中的差异表达。Mudalkar等[5]于2014年对亚麻转录组进行测序，并且在拼接得到的53 854个转录本序列数据中发现了19 379个SSR标记位点。同年，Upadhyay等[6]通过比较天冬根组织和叶组织转录组拼接结果，发现在根组织中特异表达的基因，从而推测其在体甾皂苷元合成中表达的基因。从目前公布的这些无参考基因组的物种转录组测序数据的研究成果[4-7]来看，转录组测序生物信息学分析的主要内容有：①功能注释、分类及代谢途径分析；②预测编码序列框（CDS）；③样品间基因差异表达（2个及2个以上样品）；④分子标记（SNPs、SSR）的研究进展。同时，这些研究也反映出转录组测序技术的几个突出优点：①任何物种都可以进行完整的转录组分析（无需了解物种的基因或基因组的信息，可以直接在任何物种中进行最全面的转录组分析）；②更准确的基因注释；③不仅可以检测已知的转录本，还可以识别新的基因、鉴定变异体。转录组测序作为一种更为精确的测定方法，在转录组学的应用中具有革命性的意义，开辟了转录组学研究的新纪元[8]。

基因注释是基于“同源基因，功能相似”假设的基础[9-10]，利用生物信息学方法来搜索未知基因序列与公共数据库中序列的相似性，并通过与数据库中已注释的基因的的同源性来预测未知基因的功能。核酸数据库主要有GenBank（NCBI）、EMBL和DDBJ，蛋白质数据库主要有UniProt和PDB等，搜索比对软件主要有Blast系列软件等。目前基因功能分类主要有2种方法：KEGG功能分类和Gene Ontology（简称GO）分类。GO是国际标准的基因功能分类体系，它提供了一套动态更新的标准词汇表（controlled vocabulary）来全面描述生物体基因和基因产物的性质[11]。GO共有3个本体（ontology），分别描述的是分子功能（molecular function）、细胞组分（cellular component）和生物过程（biological process）[12]。GO的基本单位是term[13]（节点），每个term都对应一个属性。GO功能分析，一方面给出了基因GO功能的分类注释，另一方面给出了基因GO功能的显著性富集分析。GO功能分类注释给出了具有某个GO功能的基因数目统计量的基因列表。GO功能显著性富集分析给出了与基因组背景相比显著富集基因的GO功能条目，因而给出了显著相关的基因的生物学功能。该分析首先将所有基因映射到Gene Ontology数据库的各个term，计算每个term的基因数，然后使用超几何测试来识别GO条目，与整个基因组背景相比，显著富集的GO条目。转录组测序技术的应用和发展，将大大推动功能基因组学的发展。

尽管转录组测序已成为获得大量植物功能基因组数据的重要技术，但是非模式植物转录组研究也面临许多挑战。首先，从转录组测序中获得大量的短序列，数据分析时对计算机运算速度和内存有较高的要求。其次，由于缺乏参考基因组信息，非模式植物转录组的构建和量化必须依靠从头拼接（De novo assembly），错误拼接、不完整拼接、拼接得到的冗余数据都将影响下游分析的质量。另外，非模式植物转录组分析过程包括使用多个在线或本地化数据库、安装和使用Linux平台应用程序，以及选择和评估大规模计算参数等。所有这些都将给研究者带来不少困难。笔者以单端测序数据为例，详细介绍非模式植物转录组测序数据的分析过程，包括原始测序数据质量控制和从头开始拼接序列获得转录本序列（transcripts）、Blast同源比对、Blast2go进行大规模GO注释和基因功能预测等。这套非模式植物转录组分析流程为研究者在相关软件安装、使用方法以及注意事项等方面提供参考。

1 转录组测序数据分析

1.1 测序数据质量控制

笔者以鹰嘴豆（chickpea）的根及芽组织转录组测序数据为例介绍转录组测序数据分析过程、软件使用和结果说明。该数据包含31 028 774条长度为51 bp的原始序列，可根据数据号SRR063784直接从NCBI网站的SRA数据库下载[14]。

从SRA上下载的鹰嘴豆转录组数据为sra格式文件，这种文件不能直接使用软件进行分析，需要转化为fasta或fastq[15]格式文件才可以使用。所以，首先使用sratoolkit（http：//www.ncbi.nlm.nih.gov/Traces/sra）中的一个可执行程序fastq-dump，将下载的sra格式的序列文件（SRR063784.sra）转化为fastq格式的文件（SRR063784.fastq）。

获得原始数据后，需进行序列的从头拼接，这是后续研究的基础。原始数据中具有大量的测序接头序列、低质量碱基盒未检测碱基（用N表示）将严重影响后续组装的质量。所以，首先需要对测序数据做一些预处理，经过质控后得到的数据即为有效数据，也称为clean data。一般使用FastQC（http：//www.bioinformatics.babraham.ac.uk/projects/fastqc/）查看raw data的质量，为此可执行如下命令：./fastqc -o./ -f fastq SRR063784.fastq，其中，-o指定文件輸出路径，-f给出输入序列文件格式。FastQC输出的结果为一个压缩文件，解压后，打开文件夹中html格式文件可看到序列文件一些统计信息。统计信息包括每个碱基位点的平均质量值（per base sequence quality）、每条序列平均质量值的分布（per sequence quality scores）、序列GC含量（per sequence GC content）、序列是否含有接头（adapter content）等12项内容。通过结果报告概要（summary）就可以对数据的情况有一个初步的了解，每一项统计分析前都有一个标志，这种标志共有3种颜色：绿色、黄色和红色。绿色代表“通过”（pass），黄色代表“警告”（warn），红色代表“不合格”（fail），FastQC以此向用户指出需要注意序列数据哪些方面。

了解数据大致情况后，使用工具包NGS QC Toolkit （http：//59.163.192.90：8080 /ngsqctoolkit/）中的IlluQC.pl对raw data进行进一步过滤，为此可执行如下命令：perl IlluQC.pl -se SRR063784.fastq N A -s 20 -l 70 -o./SRR063784_NGS/，其中，-se给出输入的single-end的序列文件，N表示不过滤接头接头文库（FastQC结果显示reads不包含接头），A表示自动识别fastq文件的版本（不同版本采用不同的质量标示方案），-s设置Phred值，-l设置大于设定Phred值的read length占该序列长度的比例，-o指定输出文件路径。在执行上述命令时，当raw data中的reads的Phred值≥20（即base calling正确率要大于等于99%）的碱基数≥reads长度的70%时，reads被保留，否则被过滤掉。

程序运行结束后，所有输出的结果文件都保存在文件夹SRR063784_NGS中。其中，output_SRR063784.html中记录了raw data质量和数据过滤记录，SRR063784.fastq_filter是过滤后的序列（clean data）文件。过滤后，31 028 774条raw reads中有24 735 426条（79.72%）高质量reads保留下来，保留下来的clean data将用于从头拼接。

1.2 从头拼接

从头拼接是将De novo测序得到的序列拼接组装成连续较长的序列[16]。将这些拼接后得到的较长序列与公共数据库中公布的基因或蛋白质序列进行同源比对分析（Blast），最终可以确定基因序列。从头组装是进行无参考序列及短序列组装、快速获得表达基因的一种有效的方法。近年来，研究者们设计了各种适用于De novo assembly的软件。目前，常用的拼接软件有Trans-ABySS（http：//www.bcgsc.ca/platform/bioinfo/software/trans-abyss）、SOAPdenovo（http：//soap.genomics.org.cn/soapdenovo.html）、Trinity（http：//trinityrnaseq.sourceforge.net）、Velvet（http：//www.ebi.ac.uk/～zerbino/velvet）、Velvet/Oases（http：//www.ebi.ac.uk/～zerbino/oases）。

该研究使用Velvet结合Oases进行转录组序列的De novo 拼接。由于Velvet默认的K-mer值上限为3 若要使用的K-mer值大于3 则需要重新编译软件。例如，若将K-mer值上限设置为57，则可执行如下编译命令：make ‘MAXKMERLENGTH=57，另外，部分Velvet算法支持多核计算，对OPENMP选项进行编译后，这部分程序即可使用多核运行。如需编译OPENMP选项，可执行如下编译命令：make ‘OPENMP=1，编译好软件后，首先选择5个不同的K-mer值（27、31、37、41、47）进行单端测序序列的拼接，并执行如下velveth命令：./velveth chickpea 27，47，10 -short -fastq SRR063784_hq.fastq、./velveth chickpea 3 4 10 -short -fastq SRR063784_hq.fastq，其中，chickpea为输出文件名称；27，47，10表示输入多个K-mer值，27≤K≤47（K为奇数），10为K值步长（步长为偶数）；-fastq指出输入文件格式为fastq；-short指出输入数据类型。结果将产生5个文件夹，分别为chickpea_27、chickpea_31、chickpea_37、chickpea_41、chickpea_47。每个文件夹里包含2个文件，分别是Roadmap以及Sequences。

其次运行velvetg，由于这里使用Velvet结合Oases进行转录组测序序列组装，所以运行velvetg时只设置1个参数。具体执行如下命令：./velvetg chickpea_27 -read_trkg yes，该命令中的-read_trkg参数要求结果给出更细致的拼接描述（yes表示打开该选项）。当程序运行结束时，屏幕上会显示nodes数n50的值、最长contig的长度（bp）以及总的组装序列的大小。同时，文件夹chickpea_27中将产生8个文件，分别是contigs.fa、LastGraph、Pregraph、Sequences、Graph2、Log、Roadmaps和stats.txt。contigs.fa即为拼接得到的contigs文件，Log文件记录Velvet运行情况（包括开始时间、软件版本、执行命令、运行结果），stats.txt文件则记录对拼接得到的每一条contig的描述。对velveth产生的其他4个文件夹进行同样的操作（分别运行velvetg），最终产生5个组装结果。

比较这5个拼接结果的n50长度、contigs的數目（nodes）和contigs的平均长度这3个参数，选择最好的拼接结果。如图1所示，当K-mer为37时，拼接得到的n50长度最长（620 bp）、最大的contig长度最长（7 339 bp）、contigs的平均长度较长（202 bp），所以最终选择K-mer值为37时的拼接结果进行后续分析。

最后运行oases对Velvet拼接得到的contigs进行进一步的拼接，最终获得转录本（transcripts）。运行oases的前提是安装并运行了Velvet，并且需要将Velvet所在的文件夹命名为“velvet”或者指明Velvet的路径，为此可执行如下命令：

make ‘VELVET_DIR=～/software/velvet，值得注意的是oases默认的K-mer值上限为3 若使用的K-mer值大于3 则在使用软件前需重新编译K-mer的值。若将K-mer值上限设置为75，可执行如下命令：make ‘MAXKMERLENGTH=75，运行oases时执行如下命令：oases chickpea_37，运行结束后文件夹chickpea_37中产生2个文件，分别是transcripts.fa和contig-ordering.txt。transcripts.fa为包含组装得到的transcripts文件，而contig-ordering.txt记录了每一个transcripts中contigs的组成情况（图1）。

将拼接得到的contig或scaffold从大到小排序，累加其长度，当累加长度达总contig或scaffold长度50%的时候，最后一个contig或scaffold的长度即为n50的值。

1.3 基因注释与功能分类

基因注释是通过比对已知数据库中已被注释的同源基因的信息推断未知基因的功能。Blast+（ftp：//ftp.ncbi.nlm.nih.gov/ blast/executables/blast+/LATEST/）中Blastx的功能是将输入核苷酸序列翻译成蛋白，并将其与蛋白质数据库比对，最后输出几个相似度高的结果。该研究使用Blastx将拼接得到的transcripts比对到nr数据库（NCBI非冗余蛋白质数据库）。

要进行本地Blast搜索，首先需要从NCBI的 ftp站点下载并格式化数据库nr.gz。将下载的nr.gz放在目录ncbi-blast-2.2.25+/bin/中，解压后，利用文件夹bin/中的可执行文件makeblastdb格式化数据库，为此可执行如下命令：makeblastdb –in nr –dbtype prot -parse_seqids -out nrdb，其中，-in（nr）输入待格式化的文件（nr），-dbtype（prot）给出数据库类型（蛋白质数据库），-parse_seqids启动序列ID解析，-out（nrdb）指定输出文件名。

格式化数据库后，即可运行Blastx将拼接得到的transcipts比对到本地nr数据库，为此执行如下命令：./blast+/bin/blastx -query transcripts.fa -out transcripts.xml -db ～/software/blast+/bin /nrdb -outfmt 5 -evalue 1.0E-6 -max_target _seqs 10 -num_threads 20，上述命令中，-query给出输入待比对数据文件路径及数据文件名（transcripts.fa），-out指定输出文件名（transcripts.xml），-db 指定用于比对的数据库名称（nrdb），-outfmt指定 输入数据格式（xml格式），-evalue设置输出结果的E-value值，-num_threads：使用多线程运算。

拼接的结果中有42 203条transcripts参与比对，其中38 622条（91.5%）transcripts获得相似性搜索结果（基因注释）。此次比对获得的hits在大豆中的分布最多（47 520），其次是鹰嘴豆（33 898）。这样的结果表明，一方面参与比对的序列与豆科植物基因表现出显著的相似性，另一方面表明公共数据库中可获得的鹰嘴豆的基因组资源依然较少[13]。

Blast+只是一种预测新基因功能的基本工具，仅通过Blast的结果无法得到新基因的GO注释信息。可以将Blast搜索结果文件（xml文件）作为Blast2Go[17]的输入数据，使用Blast2Go软件进行GO注释，最终得到与输入序列相关的GO注释信息，并将GO注释信息分为molecular function、cellular component和biological process 3类及其子类。

2 高通量GO注释工具Blast2Go

目前，能进行基因产物功能注释的生物信息学软件或生物信息学方法有很多[18]，但是对非模式物种测序序列进行大规模功能注释的软件不多。在获得Blast结果后，如果再到基因本体论网站查询相关的GO注释信息，将会浪费大量的时间[19]。Blast2Go是一款用于大规模GO注释的工具，Blast2Go是一套在植物基因组研究中对未知基因功能分析的综合软件，其主要特点是：①综合多种注释策略，输出格式多样，支持多种注释数据库，包括GO、Enzyme Codes、InterPro以及KEGG；②直观的图形化界面，可输出多种结果统计图；③综合处理数据，除对序列做GO注释，还可以进行KEGG Pathway分析等，并能根据用户的设置进行分析；④可进行大规模数据的本地自动化注释，可一次性处理20 000条序列的分析。Blast2Go的注释进程包括3个步骤：Blast、Mapping和Annotation。

2.1 启动Blast2Go

进入Blast2Go主页（http：//www.blast2go.com/），下载适合计算机内存容量的版本，下载后得到图形化界面程序blast2go*.jnlp。运行Blast2Go有3个必要条件：①网络连接；②JAVA运行环境（JRE）；③配置本地数据库（本地数据库包含了执行Mapping步骤的必要信息）。若使用Blast2Go Pro（Blast2Go的付費版本），则可以使用Blast2Go提供的在线数据库，无需再配置本地数据库。在lunix下打开Blast2Go运行界面，可执行命令：Javaws -Xnosplash blast2go*.jnlp，打开Blast2Go运行界面后，在运行Blast2Go之前，需设置数据库。可供选择的数据库有3类：①公共数据库；②本地数据库（事先本地化的数据库）；③Pro Server（Blast2Go Pro用户可选）。

2.2 Blast步骤

启动Blast2Go后，可直接输入Blast的结果文件（xml格式），也可以直接输入拼接后的结果文件进行Blast比对。用户可选择的Blast方式有3种：①在NCBI运行Blast；②使用本地Blast（Blast+，需本地化数据库）；③使用CloudBLAS进行Blast。

使用NCBI的Blast+进行本地Blast比对时，可以选择的Blast程序有4种，即Blastx、Blastp、Blastn和tBlastx[13]。用户可以根据自己的需要设置E-value值，同时Blast2Go提供数目众多的数据库供用户选择，如nr、nt、swissprot、refseq_ protein、est等。选择合适的Blast程序及比对数据库，设置E-value值和最大hits数后，点击“start”便开始Blast比对步骤：①直接输入Blast结果（xml文件）；②输入序列文件；③选择Blast运行方式；④Blast设置，包括选择Blast运行程序、选择比对数据库、设置e值、设置Blast hits数以及输出文件格式；⑤查看Blast结果统计图。

2.3 Mapping步骤

Blast步骤完成后，接着可以进行Mapping步骤。Mapping是一个检索与Blast得到的hits相关的GO terms的进程。Blast2Go进行3种不同的Mapping方式：①Blast结果中的基因序列号（accession number）用来检索基因名称，检索会用到2个由NCBI提供的Mapping文件（gene-infor、gene2accession）；②Blast结果中的GI identifiers用于重新检索在UniProt ID号，检索使用来自PIR（the protein information resource，蛋白质信息资源数据库）[20]非冗余参考蛋白质数据库的Mapping文件，这个非冗余参考蛋白质数据库搜罗了来自PSD、UniProt、Swiss-Prot、TrEMBL、RefSeq、GenPept以及PDB数据库的蛋白质信息；③Blast结果中的基因序列号（accession number）直接在GO数据库中的DBXRef Table中进行搜索。

2.4 Annotation步骤

Mapping步骤结束后，进入Annotation注释步骤。通过Annotation步骤，将Mapping步骤中获得的GO terms分配到各个输入序列，得到与输入序列相关的GO注释信息，并将GO注释信息分为molecular function，cellular component和biological process这3类及其子类。利用大量的序列数目和GO terms的结果数目，通过GO slim（GO联合会提供的简化本体论术语）将得到的GO terms归类到更高层次的terms，从而可以在更高的层次上研究基因的功能。

2.5 利用Blast2Go在GO注释结果中挖掘信息

利用Blast2Go还可以进行KEGG Pathway分析。KEGG（kyoto encyclopedia of genes and genomes）是系统分析基因功能、基因组信息数据库，KEGG可以查询整合代谢途径（pathway），这样有利于研究者将基因及表达信息作为一个整体网络进行研究。在Blast2Go注释的过程中，会给出相关unigene的EC（enzyme code）号。在代谢通路中，EC号是节点（酶）的识别符，即通过EC号，可以找到unigene参与的生物学通路（pathway），因此能推断出对应的unigene如何参与生命活动及其在生命活动中发挥的作用（图2）。

3 讨论

目前，绝大多数已报道的转录组研究资料仅介绍了某个物种的转录组研究成果，很少有资料介绍转录组分析中使用的软件及软件的详细使用方法。该研究以NCBI网站SRA数据库下载的Illumina测序平台产生的数据（sra文件）为例，使用工具包NGS QC Toolkit中的IlluQC.pl对raw data（31 028 774条raw reads）进行过滤得到clean data（24 735 426条clean reads）。随后使用Velvet/Oases进行转录组拼接，最后进行基因注释和功能分类。最终，拼接得到42 203条transcripts中，有38 622条（91.5%）transcripts获得相似性搜索结果，这表明转录组测序技术是功能基因组学研究的有利手段。

该研究详细介绍了转录组测序数据（singleend）分析的流程，但研究者在具体的数据分析过程中，可能还会遇到各种各样的问题。如测序中出现的错误会影响到从头拼接的质量，所以在质量控制时，会根据数据质量情况对reads末端碱基进行适当的剪切（trimming）。其次，该研究使用的是Singleend reads，所以在进行拼接时，可以直接运行velvet。

在组装Pairedend reads时，由于velvet软件只能采用两端序列混合在一起的fasta或fastq文件，因此需先使shuffleSe quences_fastq.pl或shuffleSequen ces_fasta.pl将paired-end数据结合在一起。大多数拼接軟件使用的算法最初都是为基因组测序设计的，但由于可变剪切的存在，一个基因通常都会编码多个转录本，这给真核生物转录组拼接带来巨大的挑战[16]。

另外，由于一般实验室计算机内存限制无法一次性完成所有数据的GO注释，可以将拼接后得到的转录本大文件（transcript.fa）分成几个大小合适的fasta文件进行基因注释及GO分类，在查看annotation结果图（Statistics -> Annotation Statistics）时可分别将注释结果以txt格式输出（save->export as text），最终将结果汇总即可。

参考文献

[1] COSTA V，ANGELINI C，DE FIES I，et al.Uncovering the complexity of transcriptomes with RNASeq[J].Journal of biomedicine and biotechnolog，2010，2010：1-19.

[2] 刘红亮，郑丽明，刘青青，等.非模式生物转录组研究[J].遗传，2013，35（8）：955-970.

[3] NAGALAKSHMI U，WANG Z，WAERN K，et al.The transcriptional landscape of the yeast genome defined by RNA sequencing[J].Science，2008，320（5881）：1344-1349.

[4] ZHANG X M，ZHAO L，LARSONRABIN Z，et al.De novo sequencing and characterization of the floral transcriptome of Dendrocalamus latiflorus（Poaceae：Bambusoideae）[J].PLoS One，201 7（8）：1-15.

[5] MUDALKAR S，GOLLA R，GHATTY S，et al.De novo Transcriptome analysis of an imminent biofuel crop，Camelina sativa L.using Illumina GAIIX sequencing platform and identification of SSR markers[J].Plant Mol Biol，2014，84（1/2）：159-171.

[6] UPADHYAY S，PHUKAN U J，MISHRA S，et al.De novo leaf and root transcriptome analysis identified novel genes involved in Steroidal sapogenin biosynthesis in Asparagus racemosus [J].BMC Genomics，2014，15：1-13.

[7] LOGACHEVA M D，KASIANOV A S，VINOGRADOV D V，et al.De novo sequencing and characterization of floral transcriptome in two species of buckwheat（Fagopyrum）[J].BMC Genomics，201 12：1-17.

[8] 井赵斌，魏琳，俞靓，等.转录组测序及其在牧草基因资源发掘中的应用前景[J].草业科学，201 28（7）：1364-1369.

[9] 周华，张新，刘腾云，等.高通量转录组测序的数据分析与基因发掘[J].江西科学，201 30（5）：607-611.

[10] 黄子夏，柯才焕，陈军.大规模GO注释的生物信息学流程[J].厦门大学学报（自然科学版），201 51（1）：139-143.

[11] WANG Z Y，FANG B P，CHEN J Y，et al.De novo assembly and characterization of root trascriptome using Illumina paired-end sequencing and development of cSSR markers in sweetpotato（Ipomoea batatas）[J].BMC Genomics，2010，11（1）：726-739.

[12] 郝大程，马培，穆军，等.中药植物虎杖根的高通量转录组测序及转录组特性分析[J].中国科学，201 42（5）：398-412.

[13] HARRIS M A，CLARK J，IRELAND A，et al.The Gene Ontology（GO）database and informatics resource[J].Nucleic acids research，2004，32：258-261.

[14] GARG R，PATEL R K，TYAGI A K，et al.De novo assembly of chickpea transcriptome using short reads for gene discovery and marker identification[J].DNA Research，201 18（1）：53-63.

[15] COCK P J A，FEILDS C J，GOTO N，et al.The Sanger FASTQ file format for sequences with quality scores，and the Solexa/Illumina FASTQ variants[J].Nucleic acids research，2010，38（6）：1767-1771.

[16] CLARKE K，YANG Y，MARSH R，et al.Comparative analysis of de novo transcriptome assembly[J].Science China life science，2013，56（2）：156-162.

[17] CONESA A，GTZ S.Blast2Go：A comprehensive suite for functional analysis in plant genomics[J].International journal of plant genomics，2008，2008：1-12.

[18] KUMAR S，DUDLEY J.Bioinformatics software for biologist in the genomics era[J].Bioinformatics，2007，23（14）：1713-1717.

[19] 王成剛，莫志宏.整合BLAST搜索与GO注释的软件GoBlast[J].中国生物化学与分子生物学报，2006，22（12）：1003-1006.

[20] 胡绍军.蛋白质组学数据库信息资源的开发与利用[J].图书馆学研究，2006（7）：77-82.