基于轉錄組測序數據分析及高通量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.