貉外周血單核細(xì)胞轉(zhuǎn)錄組RNA-Seq 數(shù)據(jù)的de novo 拼接和信息比對研究

易 立，仝明薇，程悅寧，程世鵬

(中國農(nóng)業(yè)科學(xué)院特產(chǎn)研究所 特種經(jīng)濟(jì)動(dòng)物分子生物學(xué)國家重點(diǎn)實(shí)驗(yàn)室，吉林 長春 130112)

烏蘇里貉是貉屬的一種珍貴毛皮動(dòng)物，現(xiàn)階段我國廣泛飼養(yǎng)的烏蘇里貉亞種都是地產(chǎn)品種[1]。目前，我國貉的免疫相關(guān)分子信息挖掘亟待開展，如Toll 樣受體、干擾素、細(xì)胞因子及受體等貉免疫分子信息大部分尚未解析，采用傳統(tǒng)的依據(jù)同源基因(相似動(dòng)物)設(shè)計(jì)引物，克隆測序顯得較為繁瑣。采用2 代測序通過一次測序就能獲得大量序列信息，結(jié)合生物信息學(xué)比對，獲得全基因組序列信息及注釋信息[2-3]。

轉(zhuǎn)錄組是一個(gè)細(xì)胞中的所有轉(zhuǎn)錄本信息，包括轉(zhuǎn)錄本的數(shù)量、表達(dá)動(dòng)態(tài)、轉(zhuǎn)錄本序列信息。由于轉(zhuǎn)錄組(mRNA 測序)幾乎所有序列均為外顯子，去除了內(nèi)含子等無效信息，所以更為精簡有用。烏蘇里貉尚無全基因組序列，因此，建立貉外周血單核細(xì)胞(PBMCs)轉(zhuǎn)錄組測序及分析，對于研究其免疫系統(tǒng)特征及免疫機(jī)制具有重要意義。

本研究利用Illumina 平臺(tái)首次對貉PBMCs 轉(zhuǎn)錄組進(jìn)行序列分析，并對RNA-seq 數(shù)據(jù)進(jìn)行拼接和注釋，為后續(xù)開展貉進(jìn)化研究及免疫機(jī)制研究提供參考，為貉全基因組測序項(xiàng)目提供較好的注釋基礎(chǔ)。

1 材料和方法

1.1 主要實(shí)驗(yàn)材料 外周血樣品采集自10 只配種后公貉，每只約20 mL，共計(jì)200 mL；紅細(xì)胞裂解液購自康為世紀(jì)公司；TRIzol Reagent 購自Invitrogen公司；PolyATtract?mRNA Isolation System III with Magnetic Stand 試劑盒購自Promega 公司；TransScript?One-Step gDNA Removal and cDNA Synthesis SuperMix試劑盒購自TransGen Biotech 公司；QIAquick PCR Purification 試劑盒購自QIAGEN 公司。

1.2 貉PBMCs 總RNA 提取 按照紅細(xì)胞裂解液說明書方法，處理貉外周血并收獲PBMCs，參照TRIzol Reagent 說明書提取總RNA。利用1.2 %瓊脂糖凝膠電泳及分光光度計(jì)檢測確定總RNA 完整性和純度質(zhì)量。

1.3 cDNA 文庫的構(gòu)建和Illumina 測序 提取樣品總RNA，采用帶有Oligo(dT)的磁珠富集mRNA；在制備的mRNA 中加入fragmentation buffer 使其片段化為短片段，以其為模板，利用六堿基隨機(jī)引物合成cDNA 第一鏈，并加入緩沖液、dNTPs、RNase H和DNA polymerase I 合成cDNA 第二鏈，經(jīng)過QIAquick PCR Purification 試劑盒純化并加EB 緩沖液洗脫經(jīng)末端修復(fù)、加堿基A，加測序接頭，回收目的片段后進(jìn)行PCR 擴(kuò)增，從而完成整個(gè)文庫制備工作。構(gòu)建的文庫利用Illumina HiSeqTM2500 進(jìn)行測序，測序策略為PE100(Paired-end 技術(shù)，雙向測序100 bp)測序[4]。

1.4 質(zhì)量控制和序列的拼接 測序得到的某些原始序列含有帶接頭的、低質(zhì)量的序列，為保證信息分析質(zhì)量，對原始序列進(jìn)行過濾，得到clean reads，后續(xù)分析均基于clean reads。數(shù)據(jù)處理步驟包括：1)去除含N 比例大于10 %的read；2)去除超過50 %的堿基質(zhì)量值低于5 的read；3)去除接頭污染。對質(zhì)量控制前后整體序列的質(zhì)量進(jìn)行統(tǒng)計(jì)，質(zhì)量控制前后reads 的數(shù)量和長度通過perl 腳本進(jìn)行統(tǒng)計(jì)。

序列的de novo 拼接采用基于構(gòu)建De Bruijn 圖方法的Trintiy[5]，Trinity 能夠高效的構(gòu)建de novo 轉(zhuǎn)錄組。Trinity 利用de Bruijn 圖論的原理，針對轉(zhuǎn)錄本具有可變剪接的特點(diǎn)，組裝出全長轉(zhuǎn)錄本。

1.5 轉(zhuǎn)錄物功能注釋及分類 利用Trinotate 對開放閱讀框(ORF)和contigs 進(jìn)行功能注釋，利用Uniprot database、RNAMMER、eggNOG、GO、KEGG 對預(yù)測的序列進(jìn)行注釋。由于目前無貉基因組信息，因此本研究中的序列拼接參考BLASTX、COG(Clusters of orthologous groups of proteins)、GO(Gene ontology)數(shù)據(jù)庫的比對結(jié)果，利用BLASTX 將組裝序列與COG 直系同源簇?cái)?shù)據(jù)庫進(jìn)行比對分析，獲得貉表達(dá)基因的COG 功能注釋及其COG 功能分類。最后分別使用Blast2GO 軟件和WEGO 工具對轉(zhuǎn)錄物序列進(jìn)行GO 功能注釋及功能分類統(tǒng)計(jì)。

1.6 ORF 預(yù)測 對組裝的轉(zhuǎn)錄本利用TransDecoder鑒定編碼區(qū)域，并對長度大于300 bp 的序列進(jìn)行best ORF 查找。

1.7 基因表達(dá)分析 使用bowtie 軟件，并按照每個(gè)樣品各自特異的測序接頭信息把對應(yīng)的reads mapping 到已拼接的轉(zhuǎn)錄物上，以轉(zhuǎn)錄物序列中的reads 數(shù)表示基因的表達(dá)豐度；分析轉(zhuǎn)錄組整體表達(dá)情況，利用組裝的轉(zhuǎn)錄本，采用RSEM(RNA-Seq by Expectation-Maximization)軟件估計(jì)每個(gè)樣品表達(dá)量。

2 結(jié)果

2.1 總RNA 質(zhì)量檢測 按常規(guī)方法提取貉PBMCs樣品總RNA 并對其濃度進(jìn)行檢測，其OD260nm/OD280nm=2.05，濃度為212 ng/μL，表明提取的RNA質(zhì)量較好，可以滿足正常的實(shí)驗(yàn)要求。

2.2 測序數(shù)據(jù)初步分析 Illumina 測序共得到3.1 GB 的數(shù)據(jù)，共32 245 804 條未加工數(shù)據(jù)(Raw reads)，經(jīng)過質(zhì)量預(yù)處理后共得到有效數(shù)據(jù)(Valid reads)28 797 350 條。經(jīng)過Trinity 拼接去重復(fù)后，共得到118 868 條長度大于200 bp 的contigs，平均長度為525.53 nt，N50 值為674，最長的contig 為8 982 bp，最短的為305 bp，contigs 的長度統(tǒng)計(jì)分析見表1。

2.3 轉(zhuǎn)錄物功能注釋及分類 為從整體上了解轉(zhuǎn)錄物序列信息，首先利用Uniprot database 數(shù)據(jù)庫對上述轉(zhuǎn)錄物數(shù)據(jù)進(jìn)行BLASTX 比對。隨后利用COGs、GO 及功能注釋分類體系對具有同源對比信息的轉(zhuǎn)錄物進(jìn)行比對和功能注釋(E-value<1e-05)。在COG功能分類體系中，共獲得8 261 個(gè)COG 功能注釋，涉及24 個(gè)COG 功能類別(圖1)。其中，一般功能基因(General function prediction only)的轉(zhuǎn)錄物比例最大；其次為蛋白質(zhì)翻譯后修飾與轉(zhuǎn)運(yùn)及分子伴侶相關(guān)基因(Posttranslational modification，protein turnover and chaperones)。在該轉(zhuǎn)錄組中，涉及貉生長發(fā)育相關(guān)的功能定義主要包括：氨基酸轉(zhuǎn)運(yùn)與代謝(E)、碳水化合物轉(zhuǎn)運(yùn)與代謝(G)、核酸轉(zhuǎn)運(yùn)與代謝(E)和脂質(zhì)轉(zhuǎn)運(yùn)與代謝(I)等物質(zhì)代謝過程、翻譯，核糖體結(jié)構(gòu)與生物合成(J)、蛋白質(zhì)翻譯后修飾與轉(zhuǎn)運(yùn)，分子伴侶(O)、信號(hào)轉(zhuǎn)導(dǎo)機(jī)制(T)、無機(jī)離子轉(zhuǎn)運(yùn)與代謝(P)、輔酶轉(zhuǎn)運(yùn)與代謝(G)、次級代謝產(chǎn)物生物合成、運(yùn)輸與分解代謝(Q)及防御機(jī)理(V)等多個(gè)生理生化過程。其中信號(hào)轉(zhuǎn)導(dǎo)機(jī)制(T)及防御機(jī)理(V)功能定義可能與免疫信號(hào)相關(guān)。S 為功能未知基因，從中可以發(fā)掘新基因。

圖1 貉PBMCs 轉(zhuǎn)錄組24 種COG 功能分類Fig.1 Twenty-four types of COG function classification of the PBMCs transcriptome of raccoon dog

對所有轉(zhuǎn)錄本進(jìn)行GO 功能分類(圖2)。轉(zhuǎn)錄物的GO 注釋歸為生物學(xué)過程、細(xì)胞組分和分子功能3 大功能，分別包含了23、14 和17 個(gè)功能亞類。在生物學(xué)過程功能類型中，主要的生命過程分別是細(xì)胞過程(Cellular process)和代謝過程(Metabolic process)，細(xì)胞殺傷(Cell killing)及節(jié)律過程(Rhythmic process)比例最低。在細(xì)胞組分功能類型中，細(xì)胞(Cell)和細(xì)胞部分(Cell part)所含比例最高；而胞外區(qū)部分(Extracellular region part)、合胞體(Symplast)、神經(jīng)元(Synapase)、神經(jīng)元組成(Synapase part)、病毒(Virion)、病毒組成(Virion part)均未檢測到。在分子功能類型中，蛋白結(jié)合(Binding)和轉(zhuǎn)錄調(diào)節(jié)活性(Transcription regulator)所含比例最高；抗氧化活性(Antioxidant)、輔助運(yùn)輸?shù)鞍谆钚?Auxiliary transport protein)、蛋白標(biāo)簽(Protein tag)、調(diào)節(jié)酶活性(Proteasome regulator)等均未檢測到。

圖2 貉PBMCs 轉(zhuǎn)錄組的54 種GO 功能分類Fig.2 54 types of GO classification of the PBMCs transcriptome of raccoon dog

通過KEGG 數(shù)據(jù)庫，本轉(zhuǎn)錄組共有16 650 條能夠在KEGG 數(shù)據(jù)庫中得到注釋，并分別定位到259條相應(yīng)的通路上。其中，T-細(xì)胞受體信號(hào)通路、B-細(xì)胞受體信號(hào)通路、補(bǔ)體-血凝聯(lián)級反應(yīng)通路、吞噬體通路、自然殺傷細(xì)胞介導(dǎo)的細(xì)胞毒性通路、產(chǎn)生IgA 的腸道免疫網(wǎng)絡(luò)通路、趨化因子信號(hào)通路、TLR 受體信號(hào)通路、NLR 受體信號(hào)通路、RIR 受體信號(hào)通路、溶酶體通路等通路與貉免疫與抗病相關(guān)。其中Toll 樣受體信號(hào)途徑(Toll-Like receptor signaling pathway)見圖3，在該信號(hào)途徑中，絕大部分的分子均能夠表達(dá)。

2.4 ORF 預(yù)測及基因表達(dá)分析 為獲得具有完整ORF 的重要功能基因序列，本研究選用TransDecoder對轉(zhuǎn)錄物進(jìn)行分析鑒定編碼區(qū)域，共獲得118 868條ORF 信息。依據(jù)序列中包含的reads 個(gè)數(shù)(表達(dá)量)，對轉(zhuǎn)錄物表達(dá)豐度進(jìn)行分析。轉(zhuǎn)錄物reads 個(gè)數(shù)分布范圍較廣，F(xiàn)PKM+1 取log10 為表達(dá)量，制作表達(dá)量密度分布圖(圖4)，表達(dá)量最高的序列在logFPKM 值為0.6～0.8，表明此次測序的重復(fù)率較低。

圖3 貉PBMCs 轉(zhuǎn)錄組Toll 樣受體信號(hào)途徑分子表達(dá)Fig.3 Toll-Like receptor signaling pathway of the PBMCs transcriptome of raccoon dog

圖4 貉PBMCs 轉(zhuǎn)錄組表達(dá)量密度圖Fig.4 Expression density of the PBMCs transcriptome of raccoon dog

2.5 本地BLAST 數(shù)據(jù)庫的建立及初步應(yīng)用 下載NCBI 本地BLAST 軟件，將拼接的118 868 條ORF信息通過Makeblastdb 命令進(jìn)行數(shù)據(jù)的初始化，建立貉PBMCs 轉(zhuǎn)錄組數(shù)據(jù)庫。利用本地BLAST 命令，對其進(jìn)行檢測，以貉Toll-1 基因?yàn)槔?，進(jìn)行本地BLAST，共搜尋到7 條匹對的組裝后序列，相似率均為100 %。表明本研究建立的數(shù)據(jù)庫可信。

3 討論

本研究在文庫構(gòu)建中，采用RNA 樣本片段化后再反轉(zhuǎn)處理的方法，充分利用了RNA 對二價(jià)陽離子的敏感性及穩(wěn)定性好、減少RNA 二級結(jié)構(gòu)等優(yōu)點(diǎn)，獲得了更均勻的覆蓋率和更全面的轉(zhuǎn)錄物信息[6]。采用雙端測序(Paired-end)方法進(jìn)行高通量測序，不但增加了測序的深度，而且提高了de novo 拼接的效率和準(zhǔn)確性。

本研究利用Uniprot database 數(shù)據(jù)庫對轉(zhuǎn)錄物數(shù)據(jù)進(jìn)行BLASTX 比對，結(jié)果顯示，有少數(shù)序列與其它物種蛋白序列并無匹配，分析原因如下：1)序列片段太短，無法獲得同源性比對結(jié)果。隨著序列長度的增加，能比對上的序列比例會(huì)明顯增加。芝麻[7]、水稻[8]等物種的轉(zhuǎn)錄組序列分析中均出現(xiàn)過此現(xiàn)象。2)注釋信息不完全。目前貉基因組學(xué)及轉(zhuǎn)錄組學(xué)研究較為空白，其生物信息數(shù)據(jù)庫仍在不斷的更新完善中，基因功能注釋信息不全面造成了部分序列暫時(shí)無法獲得功能注釋信息。3)貉物種或其組織轉(zhuǎn)錄組的特異性問題。貉基因組學(xué)及功能基因組學(xué)研究相對落后，目前常選用犬基因組作為參考序列，但由于犬與貉在遺傳進(jìn)化上距離相對較遠(yuǎn)，在貉轉(zhuǎn)錄本中可能存在某些屬特異性新基因，導(dǎo)致其同源序列難以被發(fā)現(xiàn)。

對轉(zhuǎn)錄本進(jìn)行功能注釋(E-value<1e-05)。在COG 及KEGG 功能分類體系中，有多種注釋涉及貉免疫功能，因此，后續(xù)的研究重點(diǎn)將放在進(jìn)一步分析每種注釋中匹配的基因，發(fā)掘參與貉免疫應(yīng)答和防治疾病過程中的關(guān)鍵基因。

根據(jù)研究物種是否有參考基因組信息，RNASeq 序列拼接的策略也分為基于參考基因組的序列拼接和de novo 序列拼接[9-10]，或當(dāng)基因組信息不完整時(shí)，將二者結(jié)合起來進(jìn)行序列拼接。De novo 策略使用的軟件有Trinity 和Oases 等，其中Trinity 可以發(fā)現(xiàn)更多的全長轉(zhuǎn)錄本，其靈敏度甚至能夠近乎于基于基因組信息的拼接。目前已經(jīng)發(fā)表的利用de novo 拼接方法進(jìn)行的非模式生物轉(zhuǎn)錄組研究涵蓋了昆蟲、植物等許多物種[11-12]。

本研究首次對貉PBMCs 進(jìn)行了較全面的轉(zhuǎn)錄組測序，PBMCs 作為重要的免疫細(xì)胞，其轉(zhuǎn)錄本中可能包含更豐富的免疫相關(guān)基因信息，為后續(xù)開展貉進(jìn)化研究、免疫機(jī)制研究提供參考并為貉全基因組測序項(xiàng)目提供較好的注釋基礎(chǔ)。

[1]樸厚坤.特種經(jīng)濟(jì)動(dòng)物養(yǎng)殖業(yè)的現(xiàn)狀與展望[J].吉林農(nóng)業(yè)大學(xué)學(xué)報(bào)，1995，(1)：100-105.

[2]杜玲，劉剛，陸健，等.高通量測序技術(shù)的發(fā)展及其在生命科學(xué)中的應(yīng)用[J].中國畜牧獸醫(yī)，2014，12：109-116.

[3]岳桂東，高強(qiáng)，羅龍海，等.高通量測序技術(shù)在動(dòng)植物研究領(lǐng)域中的應(yīng)用[J].中國科學(xué)：生命科學(xué)，2012，02：107-124.

[4]Marioni J C,Mason C E,Mane S M,et al.RNA-seq:An assessment of technical reproducibility and comparison with gene expression arrays[J].Genome Res,2008,18:1509-1517.

[5]Haas B J,Papanicolaou A,Yassour M,et al.De novo transcript sequence reconstruction from RNA-seq using the Trinity platform for reference generation and analysis[J].Nature Protocols,2013,8(8):1494-1512.

[6]Wang Z,Gerstein M,Snyder M.RNA-Seq:A revolutionary tool for transcriptomes[J].Nat Rev Genet,2009,10(1):57-63.

[7]Wen Liang-wei,Xiao Qiong-qi,Lin Hai-wang,et al.Characterization of the sesame(Sesamum indicum L.)global transcriptome using Illumina paired-end sequencing and development of ESTSSR markers[J].BMC Genomics,2011,12:451.

[8]Ting Ting-lu,Guo Jun-lu,Dan Lin-fan,et al.Function annotation of the rice transcriptome at single-nucleotide resolution by RNA-seq[J].Genome Res,2010,20:1238-1249.

[9]Birzele F,Schaub J,Rust W,et al.Into the unknown:Expression profiling without genome sequence information in CHO by next generation sequencing[J].Nucleic Acids Res,2010,38(12):3999-4010.

[10]Sultan M,Schulz M H,Richard H,et al.A global view of gene activity and alternative splicing by deep sequencing of the human transcriptome[J].Science,2008,5891:956-960.

[11]Bin Chen,Yu Juan-zhang,Zheng Bo-he,et al.De novo transcriptome sequencing and sequence analysis of the malaria vector Anopheles sinensis(Diptera:Culicidae)[J].Parasit Vectors,2014,7(1):314.

[12]Rastogi S,Meena S,Bhattacharya A,et al.De novo sequencing and comparative analysis of holy and sweet basil transcriptomes[J].BMC Genomics,2014,15:588.