基于轉(zhuǎn)錄組測序數(shù)據(jù)分析及高通量GO注釋理論的研究

劉粉香 楊文國 孫勤紅

摘要 隨著二代測序技術(shù)的快速發(fā)展，轉(zhuǎn)錄組測序在越來越多的動植物中完成，人們獲得了大批量的轉(zhuǎn)錄組數(shù)據(jù)序列。如何從這些海量的序列數(shù)據(jù)中挖掘具有生物意義的信息已成為很多研究的關(guān)鍵所在，對未知基因的功能進(jìn)行預(yù)測和注釋就是其中一個重要的問題。轉(zhuǎn)錄組序列的功能注釋是功能基因組學(xué)研究的一項重要內(nèi)容，基因本體論（gene ontology，GO）注釋目前是一種最重要的功能注釋方式。介紹了利用生物信息學(xué)軟件進(jìn)行轉(zhuǎn)錄組測序數(shù)據(jù)分析過程，包括數(shù)據(jù)質(zhì)量控制和過濾、從頭拼接（De novo assembly）、同源比對以及大規(guī)模GO注釋，為從事轉(zhuǎn)錄組測序特別是非模式植物轉(zhuǎn)錄組測序研究者在數(shù)據(jù)分析方面提供參考。

關(guān)鍵詞 二代測序；轉(zhuǎn)錄組；從頭拼接；GO注釋

中圖分類號 Q-3文獻(xiàn)標(biāo)識碼 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

廣義上的轉(zhuǎn)錄組是指生物體細(xì)胞或組織在特定狀態(tài)下所轉(zhuǎn)錄出來的所有RNA的總和，包括RNA（即mRNA）編碼蛋白質(zhì)和RNA（ncRNA，如rRNA、tRNA、microRNA等）非編碼蛋白質(zhì)；狹義上的轉(zhuǎn)錄組通常指所有mRNA的總和[1]。轉(zhuǎn)錄基因組學(xué)研究被轉(zhuǎn)錄的基因，是挖掘轉(zhuǎn)錄基因的功能基因極其重要的途徑，功能基因組學(xué)研究在基因進(jìn)化、遺傳育種等研究中具有非常重要的意義[2]。轉(zhuǎn)錄組研究的技術(shù)手段大體上有EST序列構(gòu)建、芯片技術(shù)和二代測序技術(shù)等。隨著二代測序（next generation sequencing）技術(shù)的發(fā)展和應(yīng)用，許多物種已經(jīng)完成了轉(zhuǎn)錄組測序。早在2008年，Nagalakshmi 等[3]利用 RNA-Seq 技術(shù)進(jìn)行了酵母轉(zhuǎn)錄組測序。近年來，越來越多的無參考基因組物種先后完成了轉(zhuǎn)錄組測序。2012年，Zhang 等[4]對不同發(fā)育階段的6個麻竹花器官的轉(zhuǎn)錄組進(jìn)行測序，并分析基因的差異表達(dá)，最后預(yù)測了81個轉(zhuǎn)錄因子家族在麻竹花組織發(fā)育過程中的差異表達(dá)。Mudalkar等[5]于2014年對亞麻轉(zhuǎn)錄組進(jìn)行測序，并且在拼接得到的53 854個轉(zhuǎn)錄本序列數(shù)據(jù)中發(fā)現(xiàn)了19 379個SSR標(biāo)記位點(diǎn)。同年，Upadhyay等[6]通過比較天冬根組織和葉組織轉(zhuǎn)錄組拼接結(jié)果，發(fā)現(xiàn)在根組織中特異表達(dá)的基因，從而推測其在體甾皂苷元合成中表達(dá)的基因。從目前公布的這些無參考基因組的物種轉(zhuǎn)錄組測序數(shù)據(jù)的研究成果[4-7]來看，轉(zhuǎn)錄組測序生物信息學(xué)分析的主要內(nèi)容有：①功能注釋、分類及代謝途徑分析；②預(yù)測編碼序列框（CDS）；③樣品間基因差異表達(dá)（2個及2個以上樣品）；④分子標(biāo)記（SNPs、SSR）的研究進(jìn)展。同時，這些研究也反映出轉(zhuǎn)錄組測序技術(shù)的幾個突出優(yōu)點(diǎn)：①任何物種都可以進(jìn)行完整的轉(zhuǎn)錄組分析（無需了解物種的基因或基因組的信息，可以直接在任何物種中進(jìn)行最全面的轉(zhuǎn)錄組分析）；②更準(zhǔn)確的基因注釋；③不僅可以檢測已知的轉(zhuǎn)錄本，還可以識別新的基因、鑒定變異體。轉(zhuǎn)錄組測序作為一種更為精確的測定方法，在轉(zhuǎn)錄組學(xué)的應(yīng)用中具有革命性的意義，開辟了轉(zhuǎn)錄組學(xué)研究的新紀(jì)元[8]。

基因注釋是基于“同源基因，功能相似”假設(shè)的基礎(chǔ)[9-10]，利用生物信息學(xué)方法來搜索未知基因序列與公共數(shù)據(jù)庫中序列的相似性，并通過與數(shù)據(jù)庫中已注釋的基因的的同源性來預(yù)測未知基因的功能。核酸數(shù)據(jù)庫主要有GenBank（NCBI）、EMBL和DDBJ，蛋白質(zhì)數(shù)據(jù)庫主要有UniProt和PDB等，搜索比對軟件主要有Blast系列軟件等。目前基因功能分類主要有2種方法：KEGG功能分類和Gene Ontology（簡稱GO）分類。GO是國際標(biāo)準(zhǔn)的基因功能分類體系，它提供了一套動態(tài)更新的標(biāo)準(zhǔn)詞匯表（controlled vocabulary）來全面描述生物體基因和基因產(chǎn)物的性質(zhì)[11]。GO共有3個本體（ontology），分別描述的是分子功能（molecular function）、細(xì)胞組分（cellular component）和生物過程（biological process）[12]。GO的基本單位是term[13]（節(jié)點(diǎn)），每個term都對應(yīng)一個屬性。GO功能分析，一方面給出了基因GO功能的分類注釋，另一方面給出了基因GO功能的顯著性富集分析。GO功能分類注釋給出了具有某個GO功能的基因數(shù)目統(tǒng)計量的基因列表。GO功能顯著性富集分析給出了與基因組背景相比顯著富集基因的GO功能條目，因而給出了顯著相關(guān)的基因的生物學(xué)功能。該分析首先將所有基因映射到Gene Ontology數(shù)據(jù)庫的各個term，計算每個term的基因數(shù)，然后使用超幾何測試來識別GO條目，與整個基因組背景相比，顯著富集的GO條目。轉(zhuǎn)錄組測序技術(shù)的應(yīng)用和發(fā)展，將大大推動功能基因組學(xué)的發(fā)展。

盡管轉(zhuǎn)錄組測序已成為獲得大量植物功能基因組數(shù)據(jù)的重要技術(shù)，但是非模式植物轉(zhuǎn)錄組研究也面臨許多挑戰(zhàn)。首先，從轉(zhuǎn)錄組測序中獲得大量的短序列，數(shù)據(jù)分析時對計算機(jī)運(yùn)算速度和內(nèi)存有較高的要求。其次，由于缺乏參考基因組信息，非模式植物轉(zhuǎn)錄組的構(gòu)建和量化必須依靠從頭拼接（De novo assembly），錯誤拼接、不完整拼接、拼接得到的冗余數(shù)據(jù)都將影響下游分析的質(zhì)量。另外，非模式植物轉(zhuǎn)錄組分析過程包括使用多個在線或本地化數(shù)據(jù)庫、安裝和使用Linux平臺應(yīng)用程序，以及選擇和評估大規(guī)模計算參數(shù)等。所有這些都將給研究者帶來不少困難。筆者以單端測序數(shù)據(jù)為例，詳細(xì)介紹非模式植物轉(zhuǎn)錄組測序數(shù)據(jù)的分析過程，包括原始測序數(shù)據(jù)質(zhì)量控制和從頭開始拼接序列獲得轉(zhuǎn)錄本序列（transcripts）、Blast同源比對、Blast2go進(jìn)行大規(guī)模GO注釋和基因功能預(yù)測等。這套非模式植物轉(zhuǎn)錄組分析流程為研究者在相關(guān)軟件安裝、使用方法以及注意事項等方面提供參考。

1 轉(zhuǎn)錄組測序數(shù)據(jù)分析

1.1 測序數(shù)據(jù)質(zhì)量控制

筆者以鷹嘴豆（chickpea）的根及芽組織轉(zhuǎn)錄組測序數(shù)據(jù)為例介紹轉(zhuǎn)錄組測序數(shù)據(jù)分析過程、軟件使用和結(jié)果說明。該數(shù)據(jù)包含31 028 774條長度為51 bp的原始序列，可根據(jù)數(shù)據(jù)號SRR063784直接從NCBI網(wǎng)站的SRA數(shù)據(jù)庫下載[14]。

從SRA上下載的鷹嘴豆轉(zhuǎn)錄組數(shù)據(jù)為sra格式文件，這種文件不能直接使用軟件進(jìn)行分析，需要轉(zhuǎn)化為fasta或fastq[15]格式文件才可以使用。所以，首先使用sratoolkit（http：//www.ncbi.nlm.nih.gov/Traces/sra）中的一個可執(zhí)行程序fastq-dump，將下載的sra格式的序列文件（SRR063784.sra）轉(zhuǎn)化為fastq格式的文件（SRR063784.fastq）。

獲得原始數(shù)據(jù)后，需進(jìn)行序列的從頭拼接，這是后續(xù)研究的基礎(chǔ)。原始數(shù)據(jù)中具有大量的測序接頭序列、低質(zhì)量堿基盒未檢測堿基（用N表示）將嚴(yán)重影響后續(xù)組裝的質(zhì)量。所以，首先需要對測序數(shù)據(jù)做一些預(yù)處理，經(jīng)過質(zhì)控后得到的數(shù)據(jù)即為有效數(shù)據(jù)，也稱為clean data。一般使用FastQC（http：//www.bioinformatics.babraham.ac.uk/projects/fastqc/）查看raw data的質(zhì)量，為此可執(zhí)行如下命令：./fastqc -o./ -f fastq SRR063784.fastq，其中，-o指定文件輸出路徑，-f給出輸入序列文件格式。FastQC輸出的結(jié)果為一個壓縮文件，解壓后，打開文件夾中html格式文件可看到序列文件一些統(tǒng)計信息。統(tǒng)計信息包括每個堿基位點(diǎn)的平均質(zhì)量值（per base sequence quality）、每條序列平均質(zhì)量值的分布（per sequence quality scores）、序列GC含量（per sequence GC content）、序列是否含有接頭（adapter content）等12項內(nèi)容。通過結(jié)果報告概要（summary）就可以對數(shù)據(jù)的情況有一個初步的了解，每一項統(tǒng)計分析前都有一個標(biāo)志，這種標(biāo)志共有3種顏色：綠色、黃色和紅色。綠色代表“通過”（pass），黃色代表“警告”（warn），紅色代表“不合格”（fail），F(xiàn)astQC以此向用戶指出需要注意序列數(shù)據(jù)哪些方面。

了解數(shù)據(jù)大致情況后，使用工具包NGS QC Toolkit （http：//59.163.192.90：8080 /ngsqctoolkit/）中的IlluQC.pl對raw data進(jìn)行進(jìn)一步過濾，為此可執(zhí)行如下命令：perl IlluQC.pl -se SRR063784.fastq N A -s 20 -l 70 -o./SRR063784_NGS/，其中，-se給出輸入的single-end的序列文件，N表示不過濾接頭接頭文庫（FastQC結(jié)果顯示reads不包含接頭），A表示自動識別fastq文件的版本（不同版本采用不同的質(zhì)量標(biāo)示方案），-s設(shè)置Phred值，-l設(shè)置大于設(shè)定Phred值的read length占該序列長度的比例，-o指定輸出文件路徑。在執(zhí)行上述命令時，當(dāng)raw data中的reads的Phred值≥20（即base calling正確率要大于等于99%）的堿基數(shù)≥reads長度的70%時，reads被保留，否則被過濾掉。

程序運(yùn)行結(jié)束后，所有輸出的結(jié)果文件都保存在文件夾SRR063784_NGS中。其中，output_SRR063784.html中記錄了raw data質(zhì)量和數(shù)據(jù)過濾記錄，SRR063784.fastq_filter是過濾后的序列（clean data）文件。過濾后，31 028 774條raw reads中有24 735 426條（79.72%）高質(zhì)量reads保留下來，保留下來的clean data將用于從頭拼接。

1.2 從頭拼接

從頭拼接是將De novo測序得到的序列拼接組裝成連續(xù)較長的序列[16]。將這些拼接后得到的較長序列與公共數(shù)據(jù)庫中公布的基因或蛋白質(zhì)序列進(jìn)行同源比對分析（Blast），最終可以確定基因序列。從頭組裝是進(jìn)行無參考序列及短序列組裝、快速獲得表達(dá)基因的一種有效的方法。近年來，研究者們設(shè)計了各種適用于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結(jié)合Oases進(jìn)行轉(zhuǎn)錄組序列的De novo 拼接。由于Velvet默認(rèn)的K-mer值上限為3 若要使用的K-mer值大于3 則需要重新編譯軟件。例如，若將K-mer值上限設(shè)置為57，則可執(zhí)行如下編譯命令：make ‘MAXKMERLENGTH=57，另外，部分Velvet算法支持多核計算，對OPENMP選項進(jìn)行編譯后，這部分程序即可使用多核運(yùn)行。如需編譯OPENMP選項，可執(zhí)行如下編譯命令：make ‘OPENMP=1，編譯好軟件后，首先選擇5個不同的K-mer值（27、31、37、41、47）進(jìn)行單端測序序列的拼接，并執(zhí)行如下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為奇數(shù)），10為K值步長（步長為偶數(shù)）；-fastq指出輸入文件格式為fastq；-short指出輸入數(shù)據(jù)類型。結(jié)果將產(chǎn)生5個文件夾，分別為chickpea_27、chickpea_31、chickpea_37、chickpea_41、chickpea_47。每個文件夾里包含2個文件，分別是Roadmap以及Sequences。

其次運(yùn)行velvetg，由于這里使用Velvet結(jié)合Oases進(jìn)行轉(zhuǎn)錄組測序序列組裝，所以運(yùn)行velvetg時只設(shè)置1個參數(shù)。具體執(zhí)行如下命令：./velvetg chickpea_27 -read_trkg yes，該命令中的-read_trkg參數(shù)要求結(jié)果給出更細(xì)致的拼接描述（yes表示打開該選項）。當(dāng)程序運(yùn)行結(jié)束時，屏幕上會顯示nodes數(shù)n50的值、最長contig的長度（bp）以及總的組裝序列的大小。同時，文件夾chickpea_27中將產(chǎn)生8個文件，分別是contigs.fa、LastGraph、Pregraph、Sequences、Graph2、Log、Roadmaps和stats.txt。contigs.fa即為拼接得到的contigs文件，Log文件記錄Velvet運(yùn)行情況（包括開始時間、軟件版本、執(zhí)行命令、運(yùn)行結(jié)果），stats.txt文件則記錄對拼接得到的每一條contig的描述。對velveth產(chǎn)生的其他4個文件夾進(jìn)行同樣的操作（分別運(yùn)行velvetg），最終產(chǎn)生5個組裝結(jié)果。

比較這5個拼接結(jié)果的n50長度、contigs的數(shù)目（nodes）和contigs的平均長度這3個參數(shù)，選擇最好的拼接結(jié)果。如圖1所示，當(dāng)K-mer為37時，拼接得到的n50長度最長（620 bp）、最大的contig長度最長（7 339 bp）、contigs的平均長度較長（202 bp），所以最終選擇K-mer值為37時的拼接結(jié)果進(jìn)行后續(xù)分析。

最后運(yùn)行oases對Velvet拼接得到的contigs進(jìn)行進(jìn)一步的拼接，最終獲得轉(zhuǎn)錄本（transcripts）。運(yùn)行oases的前提是安裝并運(yùn)行了Velvet，并且需要將Velvet所在的文件夾命名為“velvet”或者指明Velvet的路徑，為此可執(zhí)行如下命令：

make ‘VELVET_DIR=～/software/velvet，值得注意的是oases默認(rèn)的K-mer值上限為3 若使用的K-mer值大于3 則在使用軟件前需重新編譯K-mer的值。若將K-mer值上限設(shè)置為75，可執(zhí)行如下命令：make ‘MAXKMERLENGTH=75，運(yùn)行oases時執(zhí)行如下命令：oases chickpea_37，運(yùn)行結(jié)束后文件夾chickpea_37中產(chǎn)生2個文件，分別是transcripts.fa和contig-ordering.txt。transcripts.fa為包含組裝得到的transcripts文件，而contig-ordering.txt記錄了每一個transcripts中contigs的組成情況（圖1）。

將拼接得到的contig或scaffold從大到小排序，累加其長度，當(dāng)累加長度達(dá)總contig或scaffold長度50%的時候，最后一個contig或scaffold的長度即為n50的值。

1.3 基因注釋與功能分類

基因注釋是通過比對已知數(shù)據(jù)庫中已被注釋的同源基因的信息推斷未知基因的功能。Blast+（ftp：//ftp.ncbi.nlm.nih.gov/ blast/executables/blast+/LATEST/）中Blastx的功能是將輸入核苷酸序列翻譯成蛋白，并將其與蛋白質(zhì)數(shù)據(jù)庫比對，最后輸出幾個相似度高的結(jié)果。該研究使用Blastx將拼接得到的transcripts比對到nr數(shù)據(jù)庫（NCBI非冗余蛋白質(zhì)數(shù)據(jù)庫）。

要進(jìn)行本地Blast搜索，首先需要從NCBI的 ftp站點(diǎn)下載并格式化數(shù)據(jù)庫nr.gz。將下載的nr.gz放在目錄ncbi-blast-2.2.25+/bin/中，解壓后，利用文件夾bin/中的可執(zhí)行文件makeblastdb格式化數(shù)據(jù)庫，為此可執(zhí)行如下命令：makeblastdb –in nr –dbtype prot -parse_seqids -out nrdb，其中，-in（nr）輸入待格式化的文件（nr），-dbtype（prot）給出數(shù)據(jù)庫類型（蛋白質(zhì)數(shù)據(jù)庫），-parse_seqids啟動序列ID解析，-out（nrdb）指定輸出文件名。

格式化數(shù)據(jù)庫后，即可運(yùn)行Blastx將拼接得到的transcipts比對到本地nr數(shù)據(jù)庫，為此執(zhí)行如下命令：./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給出輸入待比對數(shù)據(jù)文件路徑及數(shù)據(jù)文件名（transcripts.fa），-out指定輸出文件名（transcripts.xml），-db 指定用于比對的數(shù)據(jù)庫名稱（nrdb），-outfmt指定 輸入數(shù)據(jù)格式（xml格式），-evalue設(shè)置輸出結(jié)果的E-value值，-num_threads：使用多線程運(yùn)算。

拼接的結(jié)果中有42 203條transcripts參與比對，其中38 622條（91.5%）transcripts獲得相似性搜索結(jié)果（基因注釋）。此次比對獲得的hits在大豆中的分布最多（47 520），其次是鷹嘴豆（33 898）。這樣的結(jié)果表明，一方面參與比對的序列與豆科植物基因表現(xiàn)出顯著的相似性，另一方面表明公共數(shù)據(jù)庫中可獲得的鷹嘴豆的基因組資源依然較少[13]。

Blast+只是一種預(yù)測新基因功能的基本工具，僅通過Blast的結(jié)果無法得到新基因的GO注釋信息?？梢詫last搜索結(jié)果文件（xml文件）作為Blast2Go[17]的輸入數(shù)據(jù)，使用Blast2Go軟件進(jìn)行GO注釋，最終得到與輸入序列相關(guān)的GO注釋信息，并將GO注釋信息分為molecular function、cellular component和biological process 3類及其子類。

2 高通量GO注釋工具Blast2Go

目前，能進(jìn)行基因產(chǎn)物功能注釋的生物信息學(xué)軟件或生物信息學(xué)方法有很多[18]，但是對非模式物種測序序列進(jìn)行大規(guī)模功能注釋的軟件不多。在獲得Blast結(jié)果后，如果再到基因本體論網(wǎng)站查詢相關(guān)的GO注釋信息，將會浪費(fèi)大量的時間[19]。Blast2Go是一款用于大規(guī)模GO注釋的工具，Blast2Go是一套在植物基因組研究中對未知基因功能分析的綜合軟件，其主要特點(diǎn)是：①綜合多種注釋策略，輸出格式多樣，支持多種注釋數(shù)據(jù)庫，包括GO、Enzyme Codes、InterPro以及KEGG；②直觀的圖形化界面，可輸出多種結(jié)果統(tǒng)計圖；③綜合處理數(shù)據(jù)，除對序列做GO注釋，還可以進(jìn)行KEGG Pathway分析等，并能根據(jù)用戶的設(shè)置進(jìn)行分析；④可進(jìn)行大規(guī)模數(shù)據(jù)的本地自動化注釋，可一次性處理20 000條序列的分析。Blast2Go的注釋進(jìn)程包括3個步驟：Blast、Mapping和Annotation。

2.1 啟動Blast2Go

進(jìn)入Blast2Go主頁（http：//www.blast2go.com/），下載適合計算機(jī)內(nèi)存容量的版本，下載后得到圖形化界面程序blast2go*.jnlp。運(yùn)行Blast2Go有3個必要條件：①網(wǎng)絡(luò)連接；②JAVA運(yùn)行環(huán)境（JRE）；③配置本地數(shù)據(jù)庫（本地數(shù)據(jù)庫包含了執(zhí)行Mapping步驟的必要信息）。若使用Blast2Go Pro（Blast2Go的付費(fèi)版本），則可以使用Blast2Go提供的在線數(shù)據(jù)庫，無需再配置本地數(shù)據(jù)庫。在lunix下打開Blast2Go運(yùn)行界面，可執(zhí)行命令：Javaws -Xnosplash blast2go*.jnlp，打開Blast2Go運(yùn)行界面后，在運(yùn)行Blast2Go之前，需設(shè)置數(shù)據(jù)庫?？晒┻x擇的數(shù)據(jù)庫有3類：①公共數(shù)據(jù)庫；②本地數(shù)據(jù)庫（事先本地化的數(shù)據(jù)庫）；③Pro Server（Blast2Go Pro用戶可選）。

2.2 Blast步驟

啟動Blast2Go后，可直接輸入Blast的結(jié)果文件（xml格式），也可以直接輸入拼接后的結(jié)果文件進(jìn)行Blast比對。用戶可選擇的Blast方式有3種：①在NCBI運(yùn)行Blast；②使用本地Blast（Blast+，需本地化數(shù)據(jù)庫）；③使用CloudBLAS進(jìn)行Blast。

使用NCBI的Blast+進(jìn)行本地Blast比對時，可以選擇的Blast程序有4種，即Blastx、Blastp、Blastn和tBlastx[13]。用戶可以根據(jù)自己的需要設(shè)置E-value值，同時Blast2Go提供數(shù)目眾多的數(shù)據(jù)庫供用戶選擇，如nr、nt、swissprot、refseq_ protein、est等。選擇合適的Blast程序及比對數(shù)據(jù)庫，設(shè)置E-value值和最大hits數(shù)后，點(diǎn)擊“start”便開始Blast比對步驟：①直接輸入Blast結(jié)果（xml文件）；②輸入序列文件；③選擇Blast運(yùn)行方式；④Blast設(shè)置，包括選擇Blast運(yùn)行程序、選擇比對數(shù)據(jù)庫、設(shè)置e值、設(shè)置Blast hits數(shù)以及輸出文件格式；⑤查看Blast結(jié)果統(tǒng)計圖。

2.3 Mapping步驟

Blast步驟完成后，接著可以進(jìn)行Mapping步驟。Mapping是一個檢索與Blast得到的hits相關(guān)的GO terms的進(jìn)程。Blast2Go進(jìn)行3種不同的Mapping方式：①Blast結(jié)果中的基因序列號（accession number）用來檢索基因名稱，檢索會用到2個由NCBI提供的Mapping文件（gene-infor、gene2accession）；②Blast結(jié)果中的GI identifiers用于重新檢索在UniProt ID號，檢索使用來自PIR（the protein information resource，蛋白質(zhì)信息資源數(shù)據(jù)庫）[20]非冗余參考蛋白質(zhì)數(shù)據(jù)庫的Mapping文件，這個非冗余參考蛋白質(zhì)數(shù)據(jù)庫搜羅了來自PSD、UniProt、Swiss-Prot、TrEMBL、RefSeq、GenPept以及PDB數(shù)據(jù)庫的蛋白質(zhì)信息；③Blast結(jié)果中的基因序列號（accession number）直接在GO數(shù)據(jù)庫中的DBXRef Table中進(jìn)行搜索。

2.4 Annotation步驟

Mapping步驟結(jié)束后，進(jìn)入Annotation注釋步驟。通過Annotation步驟，將Mapping步驟中獲得的GO terms分配到各個輸入序列，得到與輸入序列相關(guān)的GO注釋信息，并將GO注釋信息分為molecular function，cellular component和biological process這3類及其子類。利用大量的序列數(shù)目和GO terms的結(jié)果數(shù)目，通過GO slim（GO聯(lián)合會提供的簡化本體論術(shù)語）將得到的GO terms歸類到更高層次的terms，從而可以在更高的層次上研究基因的功能。

2.5 利用Blast2Go在GO注釋結(jié)果中挖掘信息

利用Blast2Go還可以進(jìn)行KEGG Pathway分析。KEGG（kyoto encyclopedia of genes and genomes）是系統(tǒng)分析基因功能、基因組信息數(shù)據(jù)庫，KEGG可以查詢整合代謝途徑（pathway），這樣有利于研究者將基因及表達(dá)信息作為一個整體網(wǎng)絡(luò)進(jìn)行研究。在Blast2Go注釋的過程中，會給出相關(guān)unigene的EC（enzyme code）號。在代謝通路中，EC號是節(jié)點(diǎn)（酶）的識別符，即通過EC號，可以找到unigene參與的生物學(xué)通路（pathway），因此能推斷出對應(yīng)的unigene如何參與生命活動及其在生命活動中發(fā)揮的作用（圖2）。

3 討論

目前，絕大多數(shù)已報道的轉(zhuǎn)錄組研究資料僅介紹了某個物種的轉(zhuǎn)錄組研究成果，很少有資料介紹轉(zhuǎn)錄組分析中使用的軟件及軟件的詳細(xì)使用方法。該研究以NCBI網(wǎng)站SRA數(shù)據(jù)庫下載的Illumina測序平臺產(chǎn)生的數(shù)據(jù)（sra文件）為例，使用工具包NGS QC Toolkit中的IlluQC.pl對raw data（31 028 774條raw reads）進(jìn)行過濾得到clean data（24 735 426條clean reads）。隨后使用Velvet/Oases進(jìn)行轉(zhuǎn)錄組拼接，最后進(jìn)行基因注釋和功能分類。最終，拼接得到42 203條transcripts中，有38 622條（91.5%）transcripts獲得相似性搜索結(jié)果，這表明轉(zhuǎn)錄組測序技術(shù)是功能基因組學(xué)研究的有利手段。

該研究詳細(xì)介紹了轉(zhuǎn)錄組測序數(shù)據(jù)（singleend）分析的流程，但研究者在具體的數(shù)據(jù)分析過程中，可能還會遇到各種各樣的問題。如測序中出現(xiàn)的錯誤會影響到從頭拼接的質(zhì)量，所以在質(zhì)量控制時，會根據(jù)數(shù)據(jù)質(zhì)量情況對reads末端堿基進(jìn)行適當(dāng)?shù)募羟校╰rimming）。其次，該研究使用的是Singleend reads，所以在進(jìn)行拼接時，可以直接運(yùn)行velvet。

在組裝Pairedend reads時，由于velvet軟件只能采用兩端序列混合在一起的fasta或fastq文件，因此需先使shuffleSe quences_fastq.pl或shuffleSequen ces_fasta.pl將paired-end數(shù)據(jù)結(jié)合在一起。大多數(shù)拼接軟件使用的算法最初都是為基因組測序設(shè)計的，但由于可變剪切的存在，一個基因通常都會編碼多個轉(zhuǎn)錄本，這給真核生物轉(zhuǎn)錄組拼接帶來巨大的挑戰(zhàn)[16]。

另外，由于一般實(shí)驗(yàn)室計算機(jī)內(nèi)存限制無法一次性完成所有數(shù)據(jù)的GO注釋，可以將拼接后得到的轉(zhuǎn)錄本大文件（transcript.fa）分成幾個大小合適的fasta文件進(jìn)行基因注釋及GO分類，在查看annotation結(jié)果圖（Statistics -> Annotation Statistics）時可分別將注釋結(jié)果以txt格式輸出（save->export as text），最終將結(jié)果匯總即可。

參考文獻(xiàn)

[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] 劉紅亮，鄭麗明，劉青青，等.非模式生物轉(zhuǎn)錄組研究[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] 井趙斌，魏琳，俞靚，等.轉(zhuǎn)錄組測序及其在牧草基因資源發(fā)掘中的應(yīng)用前景[J].草業(yè)科學(xué)，201 28（7）：1364-1369.

[9] 周華，張新，劉騰云，等.高通量轉(zhuǎn)錄組測序的數(shù)據(jù)分析與基因發(fā)掘[J].江西科學(xué)，201 30（5）：607-611.

[10] 黃子夏，柯才煥，陳軍.大規(guī)模GO注釋的生物信息學(xué)流程[J].廈門大學(xué)學(xué)報（自然科學(xué)版），201 51（1）：139-143.

[11] WANG Z Y，F(xiàn)ANG 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] 郝大程，馬培，穆軍，等.中藥植物虎杖根的高通量轉(zhuǎn)錄組測序及轉(zhuǎn)錄組特性分析[J].中國科學(xué)，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，F(xiàn)EILDS 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].中國生物化學(xué)與分子生物學(xué)報，2006，22（12）：1003-1006.

[20] 胡紹軍.蛋白質(zhì)組學(xué)數(shù)據(jù)庫信息資源的開發(fā)與利用[J].圖書館學(xué)研究，2006（7）：77-82.