RNA-seq技術(shù)在野生動物研究中的應(yīng)用

肖 慧 吳 盡 趙敏蝶 丁 新 劉學(xué)東

(東北林業(yè)大學(xué)，哈爾濱，150040)

1 RNA-seq技術(shù)的發(fā)展歷史及其在生物學(xué)研究中的作用

1964年Holley第一次成功獲得一個完整基因的核苷酸序列［6］，此后核酸測序方法就不斷地快速發(fā)展。隨著高通量測序時代的到來，大規(guī)模并行測序(massive parallel sequencing，MPS)平臺如Roche公司(454 GS-FLX)、Illumina公司(Genome Analyzer II)和 ABI公司(ABSOLiD)徹底變革了測序技術(shù)，也改變了轉(zhuǎn)錄組的研究方法，產(chǎn)生了RNA測序技術(shù)(RNA-seq)。

轉(zhuǎn)錄組研究能夠從整體水平研究基因功能以及基因結(jié)構(gòu)，揭示特定生物學(xué)過程以及疾病發(fā)生過程中的分子機(jī)理。RNA-seq作為一種新的高效、快捷的研究手段正在改變著人們對轉(zhuǎn)錄組的認(rèn)識。在生物學(xué)研究中，RNA-seq可以進(jìn)行轉(zhuǎn)錄本結(jié)構(gòu)研究(基因邊界鑒定、可變剪切研究等)、轉(zhuǎn)錄本結(jié)構(gòu)變異研究(如基因融合、編碼區(qū) SNP研究)、非編碼區(qū)域功能研究(Non-coding RNA、microRNA前體研究等)、基因表達(dá)水平研究以及全新轉(zhuǎn)錄本的發(fā)現(xiàn)［7］。這些研究為人類認(rèn)識和了解生物機(jī)體和各種功能提供了重要的方法和途徑。同樣，RNA-seq在野生動物的研究中也得到推廣和使用。

2 RNA-seq在野生動物研究中的方法

近年來，隨著轉(zhuǎn)錄組測序技術(shù)的成熟與發(fā)展，許多研究人員利用該技術(shù)對各種野生動物的轉(zhuǎn)錄本進(jìn)行研究，建立了較為全面的相關(guān)野生動物轉(zhuǎn)錄組數(shù)據(jù)庫［8－10］，為進(jìn)一步研究野生動物生長過程中的基因結(jié)構(gòu)和功能的變化及代謝通路的調(diào)控等提供實驗依據(jù)。研究的具體方法通常分為以下幾部分。

2.1 樣品分離和文庫制備

首先，研究者應(yīng)根據(jù)研究對象及研究目的，選擇性地進(jìn)行樣品采集和處理?？紤]到動物福利，盡量做到低損傷取樣甚至無損傷取樣，楊曉光［8］為了建立馬鹿(Cervus elaphus)鹿茸増生區(qū)茸皮和軟骨轉(zhuǎn)錄組數(shù)據(jù)庫，分別采集3頭3～4歲雄性東北馬鹿的生長期鹿茸樣本，單獨提取每個樣本中的總RNA后再進(jìn)行混合，以求盡可能囊括馬鹿增生區(qū)所有基因。楊秀峰［10］為了獲取狼(Canis lupus)和家犬(Canis lupus familiaris)的轉(zhuǎn)錄組數(shù)據(jù)，選取狼和家犬的血液樣本進(jìn)行轉(zhuǎn)錄組測序。以上樣品相對于動物其他組織器官來說更容易獲取，對動物損傷較低且可恢復(fù)。在樣本采集之后，利用生物公司提供的試劑盒進(jìn)行總RNA提取，隨后將總RNA片段化，連接到特定的接頭序列并反轉(zhuǎn)錄，得到的cDNA片段進(jìn)行克隆性擴(kuò)增，制備文庫以備高通量測序使用。

2.2 高通量測序

在上述得到的cDNA片段兩端加上接頭，使用二代高通量測序儀測序得到足夠的多的reads序列。目前各測序平臺廣泛采用的一種測序方法是雙末端測序法，該方法對文庫中的每一個cDNA分子的兩端均進(jìn)行序列測定，這樣每個cDNA分子就可獲得兩個經(jīng)過測序的序列片段(reads)，相對于單端測序增加了物理覆蓋度［11－12］，通過對這兩個序列片段進(jìn)行標(biāo)記就可以在測序得到的數(shù)據(jù)中識別一對雙末端配對序列，顯著增強(qiáng)了對數(shù)據(jù)分析的能力。通常情況下，測序深度在10×～15×以上時覆蓋度和測序錯誤率控制均得以保證。在轉(zhuǎn)錄組測序中，雙末端測序不僅解決了測序reads不夠長的問題，還能發(fā)現(xiàn)新的結(jié)構(gòu)變異，區(qū)別不同的剪接體，鑒定由染色體重組造成的融合基因［11－14］等問題。

2.3 轉(zhuǎn)錄組裝配

接下來需要進(jìn)行轉(zhuǎn)錄組的裝配。模式生物常具有參考基因組DNA序列信息，可以對測序結(jié)果進(jìn)行基于參考基因組的比對和拼接，再通過進(jìn)行基因組定位和注釋來獲得基因組尺度的轉(zhuǎn)錄圖譜。這種方法雖然方便，但其明顯缺陷在于新轉(zhuǎn)錄序列的丟失。而非模式生物缺乏參考基因組，只能自體組裝(De novo assembly)，這時可以使用一些軟件如Velvet與ABySS自體組裝表達(dá)序列標(biāo)簽(EST)，或在近緣生物數(shù)據(jù)的幫助下引導(dǎo)裝配［15－16］。為提高裝配質(zhì)量，策略之一是盡量增加reads的覆蓋度，以及混合使用不同類型的reads。

2.4 生物信息學(xué)分析

轉(zhuǎn)錄組拼接后得到轉(zhuǎn)錄本，為了獲取轉(zhuǎn)錄本中的信息，接下來需要對轉(zhuǎn)錄本進(jìn)行功能注釋和表達(dá)定量分析。

2.4.1 轉(zhuǎn)錄組功能注釋及代謝通路分析

②依法劃定飲用水水源性坑塘保護(hù)區(qū)，禁止一切直接或間接污染水體的行為。制定防止污染和人為破壞的管理辦法，從源頭上保證水源的可持續(xù)利用和農(nóng)村經(jīng)濟(jì)可持續(xù)發(fā)展。

研究人員一般會用Blast程序?qū)D(zhuǎn)錄組功能進(jìn)行注釋［17］。常用的參考序列數(shù)據(jù)庫包括NCBI中的非冗余核酸數(shù)據(jù)庫(non-redundant protein database，nr database)、Swiss-Prot等。除了對每條轉(zhuǎn)錄本進(jìn)行功能注釋之外，對轉(zhuǎn)錄組的注釋還有COG(cluster of orthologous group)注釋［18］、GO(gene ontology) 注釋以及PATHWAY注釋。COG注釋可以根據(jù)基因功能和進(jìn)化關(guān)系對轉(zhuǎn)錄組中的序列進(jìn)行分類，進(jìn)而可以宏觀地認(rèn)識和比較物種的轉(zhuǎn)錄組構(gòu)成。GO注釋是從分子功能、細(xì)胞組成以及生物學(xué)過程3個方面對基因進(jìn)行注釋。Pathway指代謝通路，Pathway注釋主要指構(gòu)建轉(zhuǎn)錄組包含的生物代謝通路和調(diào)控關(guān)系網(wǎng)絡(luò)。目前Pathway分析主要有 KEGG［19］和 BioCyc［20］。以 KEGG 為例，KEGG是系統(tǒng)分析基因產(chǎn)物在細(xì)胞中的代謝途徑以及這些基因產(chǎn)物的功能的數(shù)據(jù)庫，利用KEGG可以進(jìn)一步研究基因網(wǎng)絡(luò)在生物學(xué)上的復(fù)雜功能。

2.4.2 表達(dá)定量與表達(dá)差異分析

以注釋和read映射為基礎(chǔ)，基于瀏覽器對于數(shù)據(jù)質(zhì)量可視化和特定事件進(jìn)行解釋非常重要。不過，它們只提供了有關(guān)研究的質(zhì)量畫面，大量數(shù)據(jù)及其相關(guān)細(xì)節(jié)并不能簡單地通過這種方式表現(xiàn)出來。因此，大部分RNA-seq第二階段的內(nèi)容是關(guān)于全基因組轉(zhuǎn)錄事件的自動定量研究。其研究內(nèi)容包括定量已知元件(即，已注釋的基因或外顯子)，與檢測尚未在數(shù)據(jù)庫中注釋為外顯子的新轉(zhuǎn)錄區(qū)。通常，定量步驟是進(jìn)行任何差異表達(dá)研究的基礎(chǔ)。RNA-seq的基因表達(dá)分析技術(shù)是基于對reads的計數(shù)，對低表達(dá)的基因也能夠檢測，具有靈敏度高、分辨率高、無飽和區(qū)等優(yōu)點［3，21－22］。考慮到樣本大小的影響(如 reads 數(shù)量隨基因長度不同而不同;測序深度不同，測序獲得的reads總數(shù)也不同)，Mortazavi提出了RPKM(reads per kilobases per million reads)和FPKM(fragments per kilobase per million reads)作為標(biāo)準(zhǔn)化定量指數(shù)來消除這兩種系統(tǒng)差異。經(jīng)過這種歸一化處理，不同長度、不同測序深度下的基因表達(dá)量具有可比性［23］。另外可以在KEGG通路上進(jìn)行富集分析，之后再進(jìn)行表達(dá)差異分析。

基因表達(dá)差異分析是指找出不同時間點、不同組織或者不同處理條件下具有差異表達(dá)的基因。而轉(zhuǎn)錄組研究的最終目的就是要定量多個樣本的表達(dá)來獲得差異基因表達(dá)，確定樣本特異性可變剪接異構(gòu)體和它們差異性豐度。為了分析野生動物基因表達(dá)與作用因素之間的關(guān)系，可以用統(tǒng)計學(xué)的方法對高通量測序得到的基因表達(dá)量進(jìn)行分析比對，找出樣本間差異性顯著的基因，再進(jìn)行進(jìn)一步的分析研究。研究人員通常會采用假設(shè)檢驗算法(P值)對樣本之間的差異基因進(jìn)行篩選，然后再通過錯誤發(fā)現(xiàn)率(false discovery rate，F(xiàn)DR)方法來校準(zhǔn)P值。符合以下標(biāo)準(zhǔn)的基因被認(rèn)為存在基因表達(dá)差異:當(dāng)比對某基因在兩樣本間的FDR值＜0.001，且log2Ratio值＞1時可認(rèn)為基因存在表達(dá)差異。

3 RNA-seq在野生動物研究中的應(yīng)用

3.1 構(gòu)建野生動物基因藍(lán)圖

RNA-seq最大的優(yōu)點是不局限于檢測已知基因組序列的轉(zhuǎn)錄組，還可以用于全基因組序列未知的非模式生物的轉(zhuǎn)錄組測序分析，這一點極大地拓寬了轉(zhuǎn)錄組學(xué)的研究對象范圍，為研究人員提供更多的研究思路。2008年Vera等［24］使用454 GS-20測序技術(shù)對非模式生物——慶網(wǎng)蛺蝶(Melitaea cinxia)進(jìn)行了轉(zhuǎn)錄組測序研究，這是第一個利用自體組裝進(jìn)行轉(zhuǎn)錄組研究的范例。隨后包括野生動物在內(nèi)的大量生物轉(zhuǎn)錄組測序研究不斷出現(xiàn)，RNA-seq技術(shù)極大地推動野生動物轉(zhuǎn)錄組研究，這其中包括國外研究的鱘魚(Acipenser fulvescens)［25］、虹鱒 (Oncorhynchus mykiss)［26］、響尾蛇(Crotalus adamanteus)［27］、食蟹猴(Macaca fascicu-laris)［28］和國內(nèi)研究的斑海豹(Phoca largha)［29］、絨山羊［30］、馬鹿［8］和狼［10］等。

Hale等通過RNA-seq技術(shù)成功獲得了鱘魚的轉(zhuǎn)錄組序列，確定超過5000個表達(dá)序列標(biāo)簽并鑒定877個候選SNP，大約每460個堿基就有一個SNP。楊秀峰應(yīng)用RNA-seq技術(shù)對2只家犬和3只狼的血液轉(zhuǎn)錄組進(jìn)行測序，該研究通過Illumina HiSeqTM2000平臺進(jìn)行測序，組裝后獲得了26212個基因，其中新基因1989個;總共鑒定出33229個轉(zhuǎn)錄本，其中新轉(zhuǎn)錄本1993個。這些研究為構(gòu)建野生動物基因的藍(lán)圖增添了基礎(chǔ)數(shù)據(jù)。

3.2 差異表達(dá)分析

RNA-seq的另一個優(yōu)勢是它可以捕捉不同組織或狀態(tài)下的轉(zhuǎn)錄組動態(tài)變化，通過分析不同因素作用下的基因在RNA水平表達(dá)差異性，可以將那些顯著差異表達(dá)的基因與某些生物學(xué)功能聯(lián)系起來。例如，通過RNA-seq技術(shù)對狼和家犬的血液轉(zhuǎn)錄組進(jìn)行分析，GO富集分析發(fā)現(xiàn)6個在狼中高表達(dá)的基因富集到了狼的先天性免疫系統(tǒng)上，推測這可能與狼在某些病毒的抵抗能力上大于家犬相關(guān)［10］。另外，我們研究組對馬鹿鹿茸角快速生長期鹿茸生長點茸皮和軟骨進(jìn)行高通量測序，發(fā)現(xiàn)茸皮和軟骨各自特異性表達(dá)6961和2776條基因，通過GO功能富集分析及KEGG通路富集分析，這些差異表達(dá)基因主要參與細(xì)胞結(jié)構(gòu)、細(xì)胞代謝、蛋白質(zhì)相互作用、催化活性等生物學(xué)進(jìn)程。其中涉及注釋了5328個基因的236條通路，這些基因和通路對茸皮和軟骨組織特異性生長發(fā)育具有重要的調(diào)控作用［8］。這些差異基因的發(fā)現(xiàn)為進(jìn)一步研究野生動物體內(nèi)的分子生物學(xué)機(jī)制奠定基礎(chǔ)。

3.3 非編碼區(qū)域功能研究

非編碼RNA(ncRNA)指的是未被翻譯成蛋白質(zhì)的RNA分子，重要的非編碼RNA有轉(zhuǎn)運RNA(tRNAs)、核糖體RNA(rRNAs)以及小RNA如microRNAs(miRNA)、siRNAs等，在一系列與細(xì)胞存活有關(guān)的活動中發(fā)揮重要功能［31］。Berezikov［32］等使用大規(guī)模平行測序來比較人類和黑猩猩(Pan troglodytes)大腦的microRNA含量，發(fā)現(xiàn)了447個新的miRNA基因，其中許多新的miRNA在靈長類之間并不保守，以此為基礎(chǔ)探討miRNA的進(jìn)化過程以及人類與黑猩猩的大腦進(jìn)化和功能的差異。陳艷霞等通過Solexa高通量測序?qū)β谷总浌呛腿灼そM織的miRNA進(jìn)行了全面的鑒定與特征分析，首次通過同源比對鑒定了鹿茸軟骨和茸皮miRNA共684個，其中611個哺乳動物保守miRNA和73個鹿茸新的候選miRNA。通過對這些miRNA靶基因功能注釋，發(fā)現(xiàn)在快速生長期鹿茸軟骨和茸皮中的很多基因參與細(xì)胞或細(xì)胞器構(gòu)成、核酸與蛋白質(zhì)生物合成、催化活性、代謝過程、細(xì)胞増殖、配體/受體相互作用及多種信號通路。這些基因和通路在鹿茸快速生長發(fā)育過程中起到重要調(diào)控作用［9］。

4 展望

自RNA-seq技術(shù)問世以來，已經(jīng)成為生物研究中不可或缺的技術(shù)手段，并且隨著高通量測序的快速發(fā)展，測序成本的不斷降低，轉(zhuǎn)錄學(xué)研究逐漸成為研究熱點。RNA-seq技術(shù)在野生動物研究中的應(yīng)用已有數(shù)年，NGS測序技術(shù)的快速發(fā)展又為野生動物轉(zhuǎn)錄組學(xué)的研究提供了必要的技術(shù)手段。尤其對于那些沒有基因組序列信息的野生動物，利用RNA-seq技術(shù)可以通過比較基因組學(xué)對其進(jìn)行基因組數(shù)據(jù)注釋，這極大拓寬了RNA-seq技術(shù)在野生動物領(lǐng)域的研究，在后續(xù)的野生動物研究中RNA-seq將發(fā)揮越來越大的作用。

［1］ Chu Yongjun，Corey D R．RNA sequencing:platform selection，experimental design，and data interpretation ［J］．Nucleic Acid Therapeutics，2012，22(4):271 －274．

［2］ Maher C A，Kumar-Sinha C，Cao Xuhong A，et al．Transcriptome sequencing to detect gene fusions in cancer［J］．Nature，2009，458(7234):97－101．

［3］ Wang Zhong，Gerstein M，Snyder M．RNA-seq:a revolutionary tool for transcriptomics ［J］．Nature Reviews Genetics，2009，10(1):57－63．

［4］ 劉紅亮，鄭麗明，劉青青，等．非模式生物轉(zhuǎn)錄組研究 ［J］．遺傳，2013，35(8):955－970．

［5］ Ingolia N T，Brar G A，Rouskin S A，et al．The ribosome profiling strategy for monitoring translation in vivo by deep sequencing of ribosome-protected mRNA fragments ［J］．Nature Protocols，2012，7(8):1534－1550．

［6］ Holley R W．Alanine transfer RNA，in Nobel lectures in Adokctdar biology 1933－1975［S］．Elsevier North Holland:New York，NY，USA．1977:285－300．

［7］ 祁云霞，劉永斌，榮威恒．轉(zhuǎn)錄組研究新技術(shù):RNA-seq及其應(yīng)用 ［J］．遺傳，2011，33(11):1191－1202．

［8］ 楊曉光．馬鹿(Cervus elaphus)鹿茸角頂端茸皮與軟骨組織轉(zhuǎn)錄組研究［D］．哈爾濱:東北林業(yè)大學(xué)，2015．

［9］ 陳艷霞．馬鹿(Cervus elaphus)鹿茸快速生長期生長點軟骨和茸皮組織microRNA表達(dá)譜研究 ［D］．哈爾濱:東北林業(yè)大學(xué)，2015．

［10］ 楊秀峰．基于高通量測序的狼和家犬血液轉(zhuǎn)錄組研究［D］．曲阜:曲阜師范大學(xué)，2016．

［11］ Ozsolak F，Milos P M．RNA sequencing:advances，challenges and opportunities［J］．Nature Reviews Genetics，2011，12(2):87－98．

［12］ Maher C A，Palanisamy N，Brenner J C，et al．Chimeric transcript discovery by paired-end transcriptome sequencing［J］．Proceedings of the National Academy of Sciences of the United States of America，2009，106(30):12353－12358．

［13］ Au K F，Jiang Hui，Lin Lan，et al．Detection of splice junctions from paired-end RNA-seq data by SpliceMap［J］．Nucleic Acids Research，2010，38(14):4570－4578．

［14］ Edgren H，Murumagi A，Kangaspeska S，et al．Identification of fusion genes in breast cancer by paired-end RNA-sequencing ［J］．Genome Biology，2011，12(1):R6．

［15］ Birol I，Jackman S D，Nielsen C B，et al．De novo transcriptome assembly with ABySS ［J］．Bioinformatics，2009，25(21):2872－2877．

［16］ Zerbino D R，Birney E．Velvet:algorithms for de novo short read assembly using de bruijn graphs［J］．Genome Research，2008，18(5):821－829．

［17］ Altschul S F，Gish W，Miller W，et al．Basic local alignment search tool［J］．Journal of Molecular Biology，1990，215(3):403－410．

［18］ Tatusov R L，F(xiàn)edorova N D，Jackson J D，et al．The COG database:an updated version includes eukaryotes［J］．BMC Bioinformatics，2003，4(1):41．

［19］ Ogata H，Goto S，F(xiàn)ujibuchi W，et al．Computation with the KEGG pathway database［J］．Biosystems，1998，47(1/2):119 －128．

［20］ Karp P D，Ouzounis C A，Moore-Kochlacs C，et al．Expansion of the BioCyc collection of pathway/genome databases to 160 genomes［J］．Nucleic Acids Research，2005，33(19):6083－6089．

［21］ 't Hoen P A C，Ariyurek Y，Thygesen H H，et al．Deep sequencing-based expression analysis shows major advances in robustness，resolution and inter-lab portability over five microarray platforms［J］．Nucleic Acids Research，2008，36(21):e141．

［22］ Shendure J．The beginning of the end for microarrays? ［J］．Nature Methods，2008，5(7):585 －587．

［23］ Mortazavi A，Williams B A，Mccue K，et al．Mapping and quantifying mammalian transcriptomes by RNA-seq ［J］．Nature Methods，2008，5(7):621－628．

［24］ Vera J C，Wheat C W，F(xiàn)escemyer H W，et al．Rapid transcriptome characterization for a nonmodel organism using 454 pyrosequencing［J］．Molecular Ecology，2008，17(7):1636－1647．

［25］ Hale M C，McCormick C R，Jackson J R，et al．Next-generation pyrosequencing of gonad transcriptomes in the polyploid lake sturgeon(Acipenser fulvescens):the relative merits of normalization and rarefaction in gene discovery ［J］．BMC Genomics，2009，10(1):203．

［26］ Salem M，Rexroad C E，Wang Jiannan，et al．Characterization of the rainbow trout transcriptome using Sanger and 454-pyrosequencing approaches［J］．BMC Genomics，2010，11(1):564．

［27］ Rokyta D R，Wray K P，Lemmon A R，et al．A high-throughput venom-gland transcriptome for the eastern diamondback rattlesnake(Crotalus adamanteus)and evidence for pervasive positive selection across toxin classes［J］．Toxicon，2011，57(5):657 －671．

［28］ Huh J W，Kim Y H，Park S J，et al．Large-scale transcriptome sequencing and gene analyses in the crab-eating macaque(Macaca fascicularis)for biomedical research ［J］．BMC Genomics，2012，13(1):163．

［29］ Gao Xianggang，Han Jiabo，Lu Zhichuang，et al．Characterization of the spotted seal Phoca largha transcriptome using Illumina pairedend sequencing and development of SSR markers［J］．Comparative Biochemistry and Physiology Part D:Genomics and Proteomics，2012，7(3):277－284．

［30］ 劉紅亮．基于高通量測序的遼寧絨山羊轉(zhuǎn)錄組研究與應(yīng)用［D］．楊凌:西北農(nóng)林科技大學(xué)，2013．

［31］ Ponting C P，Oliver P L，Reik W．Evolution and functions of long noncoding RNAs［J］．Cell，2009，136(4):629 －641．

［32］ Berezikov E，Thuemmler F，Van Laake L W，et al．Diversity of microRNAs in human and chimpanzee brain ［J］．Nature Genetics，2006，38(12):1375－1377．