邵世濱，閔自信，郭源旭，王權(quán)成，孫夢瑤，韓 燕，孫 健

（1．西安交通大學(xué)醫(yī)學(xué)部生物化學(xué)與分子生物學(xué)系，陜西西安 710061；2．山東醫(yī)學(xué)高等?？茖W(xué)校生物化學(xué)教研室，山東臨沂 276002）

非編碼 RNA（non-coding RNA，ncRNA）是基因組轉(zhuǎn)錄出的一類不同于mRNA的遺傳信息分子，包括小分子RNA（small RNA，sRNA）和長鏈非編碼 RNA（long non-coding RNA，lncRNA），廣泛地存在于各種生物中，而且隨著生物復(fù)雜度的升高，基因組中的非編碼序列的比例也相應(yīng)增大，提示ncRNA在生物進(jìn)化過程中的重要意義［1］。對真核細(xì)胞中ncRNA及其基因的發(fā)掘和功能研究，有可能揭示一個(gè)由ncRNA介導(dǎo)的遺傳信息傳遞方式和表達(dá)調(diào)控網(wǎng)絡(luò)，從不同于蛋白質(zhì)編碼基因的角度注釋和闡明基因組的結(jié)構(gòu)與功能，深入闡明生命活動(dòng)的本質(zhì)和規(guī)律［2］。

sRNA一般長度小于 200nt（nucleotide），由RNA聚合酶Ⅱ或RNA聚合酶Ⅲ所合成，種類眾多，通常與蛋白質(zhì)組成復(fù)合物，在細(xì)胞的生命活動(dòng)中起重要的作用［3］。按照細(xì)胞內(nèi)空間定位，sRNA分為細(xì)胞核小RNA（small nuclear RNA）和細(xì)胞質(zhì)小scRNA（small cytoplasmic RNA）兩類。其中，研究最多是microRNA（miRNA，miR），長約20～24nt一類負(fù)調(diào)控因子，起源于RNA聚合酶Ⅱ轉(zhuǎn)錄的內(nèi)源性單鏈sRNA［4］。其發(fā)現(xiàn)和深入研究進(jìn)一步揭示了基因表達(dá)調(diào)控的復(fù)雜性，功能幾乎包括了整個(gè)生物發(fā)育生長的過程，并不斷發(fā)現(xiàn)新的功能［5］。

由于ncRNA沒有編碼蛋白質(zhì)的讀碼框架，使得發(fā)現(xiàn)和鑒定sRNA有一定難度，常用方法有直接克隆測序技術(shù)、生物信息學(xué)預(yù)測、正向遺傳學(xué)方法。第二代高通量測序技術(shù)可以一次性產(chǎn)出百萬級的高質(zhì)量sRNA序列，極大地彌補(bǔ)原來的技術(shù)缺陷，成為目前最有效、最常用的sRNA鑒定方法，尤其對于特定時(shí)空的細(xì)胞表達(dá)出的全部ncRNA譜的研究，是闡明其參與的調(diào)控網(wǎng)絡(luò)的有效途徑［6］。

針對sRNA的高通量測序的Illumina的Solexa平臺，利用單分子陣列測試基因分型［7］，通過大量的平行測序，深度鑒定并定量出全基因組水平的sRNA譜，而且可減少因被檢測的核苷酸自身二級結(jié)構(gòu)，導(dǎo)致的一段區(qū)域缺失。但是，高通量測序數(shù)據(jù)會(huì)出現(xiàn)大量假陽性、假陰性及其他的雜質(zhì)（noise）。必須運(yùn)用生物信息學(xué)和統(tǒng)計(jì)學(xué)結(jié)合的方法，對這些數(shù)據(jù)進(jìn)行去規(guī)范化和標(biāo)準(zhǔn)化。

為了研究sRNA對骨骼發(fā)育與骨形成的調(diào)控作用，我們選擇了軟骨生成不同階段的重要轉(zhuǎn)折點(diǎn)：出生后d0（新生）、d21（斷乳）、d42（性成熟）3個(gè)時(shí)間點(diǎn)的雌性SD大鼠股骨頭軟骨作為研究對象［8-9］，使用Solexa測序平臺鑒定了大于18nt的高質(zhì)量sRNA序列表達(dá)，依次為：14 813 676、17 054 400、17 025 148條，并作生物信息學(xué)分析。

1 材料與方法

1．1 主要試劑、儀器及分析平臺

1．1．1 主要試劑和儀器 Trizol?試劑（Invitrogen公司）；DEPC（Amresco公司）；Agarose（Solarbio公司）；Tris base（Solarbio公司）；RNA Marker RL6000（TaKaRa 公司）；RevertAidTM First Strand cDNA Synthesis Kit（Fermentas，MBI公司）；SYBR?Premix Ex TaqTMⅡ（TaKaRa 公司）；Primers（Invitrogen公司）。Solexa Genome AnalyzerⅡ測序儀（Illumina公司）。

1．1．2 主要軟件 ①SOAP（short oligonucleotide alignment program）［10］；②MIREAP（https：／／sourceforge．net／projects／mireap／）。

1．1．3 數(shù)據(jù)庫 ①SD大鼠基因組［11］；②miRBase20．0 database；③Rfam9．1database。

1．2 實(shí)驗(yàn)動(dòng)物及樣本 實(shí)驗(yàn)選用的SPF級SD大鼠購于第四軍醫(yī)大學(xué)實(shí)驗(yàn)動(dòng)物中心。動(dòng)物在SPF級動(dòng)物房繁育，環(huán)境維持室溫（25±3）℃，相對濕度40%～60%，每日人工照明12h。動(dòng)物墊料經(jīng)高壓蒸汽滅菌，動(dòng)物飲滅菌涼開水，喂經(jīng)鈷60照射滅菌的顆粒飼料，自由進(jìn)食及飲水。分別取d0，d21，d42齡遠(yuǎn)交系SD雌性大鼠各5只，斷頸處死后，立即截取右后肢股骨，迅速剝下股骨頭軟骨組織，用冰凍的DEPC-生理鹽水沖洗后，立即用Trizol?法提取總RNA，定量后將同一時(shí)間點(diǎn)的5個(gè)樣本等量合并，每份樣本都含有總RNA 20μg。

1．3 小分子RNA文庫的構(gòu)建及測序 ①使用聚丙烯酰胺凝膠電泳（PAGE）從20μg總RNA樣本分離純化出小于30nt的sRNA分子。②將分離出的sRNA分子與Solexa接頭分子連接作為合成cDNA的模板，與接頭分子互補(bǔ)的引物，進(jìn)行12個(gè)PCR循環(huán)，獲得長度約90bp的cDNA片段。③瓊脂糖凝膠電泳中分離純化cDNA擴(kuò)增產(chǎn)物。④使用Illumina Genome Analyzer（Illumina，San Diego，CA，USA）對純化的cDNA擴(kuò)增產(chǎn)物進(jìn)行序列分析，得到的所有數(shù)據(jù)直接被處理為高質(zhì)量的圖像文件。⑤去除接頭序列和污染序列后，讀取計(jì)算分析處理清潔序列［12］。

1．4 In silico分析策略 不同的樣本由于測序深度等原因，其總量并不一樣，不能直接對不同文庫的測序拷貝數(shù)的結(jié)果進(jìn)行比較。因此，須先進(jìn)行“TPM標(biāo)準(zhǔn)化（normalize）”，即某種特定轉(zhuǎn)錄本的表達(dá)量用TPM（transcript per million，TPM）來表示（指每一百萬的總轉(zhuǎn)錄本中，該轉(zhuǎn)錄本的拷貝數(shù)）。再進(jìn)行：①重復(fù)序列的合并：為了減少計(jì)算量，分別將3個(gè)文庫中的全部清潔序列sRNA完全相同序列合并，并計(jì)算拷貝數(shù)，獲得unique sRNA序列。②TPM（transcript per million，TPM）標(biāo)準(zhǔn)化處理3個(gè)文庫中每個(gè)unique sRNA的表達(dá)量。③SOAP：使用SOAP將unique sRNA序列與SD大鼠基因組比對，保留與基因組完全匹配（perfect match）的序列作進(jìn)一步分析。④同源比對：將與基因組完全匹配的unique sRNA序列與miRBase20．0數(shù)據(jù)庫比對，確定3個(gè)文庫中表達(dá)出的已知miRNA。⑤剔除已知non-coding sRNA：將剩余數(shù)據(jù)與RNA數(shù)據(jù)庫Rfam9．1比對，Rfam database包括除miRNA外的全部已知tRNA、rRNA、snRNA、snoRNA 和 其 他 non-coding RNA序列信息。剔除與Rfam database序列信息完全匹配的unique sRNA序列。

2 結(jié) 果

2．1 各樣本小分子RNA文庫表達(dá)的基本情況 使用高通量的Solexa平臺測序技術(shù)，從d0、d21和d42三個(gè)大鼠軟骨生成不同階段的股骨頭軟骨細(xì)胞sRNA文庫中，共讀取到長度18nt至30nt的sRNA拷貝數(shù)，分別為：16 641 832條、18 560 892條和19 298 789條。為了簡化測序數(shù)據(jù)，分別對3個(gè)文庫讀取到的所有sRNA進(jìn)行分組和歸類，序列完全相同的sRNA標(biāo)簽為同一序列，并統(tǒng)計(jì)其拷貝數(shù)。在剔除了低質(zhì)量測序結(jié)果和接頭序列，并SD大鼠基因組比對后，文庫中包含了各自所表達(dá)的全部清潔序列sRNA［13］，其中，d0包含有11 169 201條sRNA，d21包含有13 814 641條sRNA，d42包含有15 792 490條（表1）。

表1 sRNA測序的基本結(jié)果Tab．1 The summary of data produced by small RNA sequencing

統(tǒng)計(jì)d0、d21、d42的3個(gè)文庫中，所有檢測到的全部清潔序列sRNA長度，20～24nt的sRNA占比分別為71．94%、72．85%、86．39%；21～23nt的sRNA 占比分別為58．67%、60．17%、74．69%；22nt的sRNA占比分別為：29．71、28．19、37．76%，峰值均為22nt（圖1A）。

2．2 In silico分析結(jié)果 ①各文庫small RNA的重復(fù)序列分析：3個(gè)文庫中的全部清潔序列sRNA中完全相同序列合并，并計(jì)算拷貝數(shù)，獲得unique sRNA序列見表2。②各文庫sRNA的基因組定位分析見表3。③各文庫sRNA表達(dá)的miRNA：3個(gè)文庫中，sRNA的長度分布峰值均為22nt，清潔序列的占比分別高達(dá)49．9%、47．2%，55．5%（圖1B）；成熟 miRNA的種類在unique sRNA種類中，占比分別為0．69%、0．78%和0．63%（表2）。④已知的ncRNA保留的與大鼠基因組完全匹配的sRNA序列，按照其生物起源和注釋，分為rRNA、tRNA、snRNA、snoRNA、exon、intron等（表2）。⑤未知的ncRNA：剔除已知miRNA、miRNA＊以及其它已知non-coding RNA后，剩余的無法與miRBase20．0和Rfam9．1匹配的未知sRNA，命名為“unann”。在d0、d21、d42文庫中，unann占unique sRNA序列標(biāo)簽的比例，分別為54．24%、57．04%、54．09%，序列拷貝數(shù)分別是2 432 498、3 322 033、4 179 434（表2）。

圖1 3個(gè)關(guān)節(jié)軟骨文庫中mRNA的長度特征Fig．1 The length characteristics of rat articular cartilage miRNAs at different development stages

表2 大鼠軟骨細(xì)胞中表達(dá)的sRNA數(shù)據(jù)統(tǒng)計(jì)Tab．2 The distribution of small RNAs in rat cartilage cells

表3 sRNA基因組定位分析Tab．3 Mapping statistics of sRNA

3 討 論

高通量測序技術(shù)進(jìn)行sRNA測序，能夠從基因組的角度研究sRNA的分布以及其在基因組上的染色質(zhì)修飾情況，是目前發(fā)現(xiàn)特定發(fā)育階段的sRNA最有效的方法，也是分析比較野生型和突變體中sRNA合成途徑重要基因表達(dá)差異與功能的有效手段，還能夠發(fā)現(xiàn)nat-siRNA和nat-miRNA等特殊作用方式的sRNA［14］。同時(shí)，避免了一代測序所固有的細(xì)菌克隆步驟，導(dǎo)致很小片段序列和極低豐度序列在克隆過程中的丟失現(xiàn)象，并能夠精確的讀取表達(dá)豐度極低的sRNA的拷貝數(shù)［15］。而且極大的提高了sRNA的發(fā)現(xiàn)效率，但是也產(chǎn)生了大量未知的候選sRNA的鑒定工作［15］。

大鼠從出生到性成熟，有著非常簡短而迅速的成長過程。僅僅6周時(shí)間，大鼠就歷經(jīng)了嬰幼兒期、兒童期和青春期這3個(gè)出生后重要的發(fā)育階段。這一過程正是骨骼發(fā)育與骨形成、軟骨細(xì)胞增殖與分化的重要過程［16］。因此，本研究中，我們以大鼠骨骼發(fā)育與骨形成的d0、d21、d42三個(gè)關(guān)鍵時(shí)點(diǎn)的股骨頭軟骨組織來構(gòu)建sRNA文庫。

但是，單只成年大鼠關(guān)節(jié)軟骨組織非常有限，軟骨組織中細(xì)胞比例也較低，RNA含量也較低，從軟骨中提取高質(zhì)量總RNA的難度相應(yīng)也較高，sRNA占總RNA比例又低。因此，從一只大鼠股骨頭軟骨能夠提取的總RNA量，不能達(dá)到Solexa測序構(gòu)建sRNA文庫所要求的單個(gè)樣本總RNA量20μg的低限。所以，我們只能將同一時(shí)間點(diǎn)的5個(gè)雌性SD大鼠樣本等量合并。

從3個(gè)文庫序列分布情況來看，基本為20～24nt的序列，大部分為21～23nt，峰值為22nt的基本特征，完全符合高質(zhì)量sRNA文庫的特征。超過總拷貝數(shù)62%的清潔序列來自于成熟miRNA序列。一般認(rèn)為，21nt和22nt miRNA傾向于與Ago蛋白構(gòu)成沉默復(fù)合體，前者使靶基因轉(zhuǎn)錄本裂解，后者使靶基因轉(zhuǎn)錄本翻譯阻遏；24nt miRNA的功能傾向于甲基化組蛋白，使靶基因轉(zhuǎn)錄抑制。三個(gè)sRNA文庫的miRNA這種分布模式，可能暗示著有功能不同或者來源各異的miRNA，參與了對骨骼發(fā)育和骨形成關(guān)鍵階段軟骨細(xì)胞增殖與分化的調(diào)控［1，5］。

但是，總體上看，3個(gè)文庫中成熟miRNA的種類在的占比卻僅僅只有0．69%、0．78%和0．63%，只是整個(gè)sRNA種類的極小部分；有一半以上的unique sRNA序列無法與miRBase20．0和Rfam9．1匹配，在基因組上找不到匹配位點(diǎn)，說明這些未知的unann sRNA可能是外源的，或者是通過復(fù)雜的轉(zhuǎn)錄后加工、修飾而形成的［2］。這些都提示：sRNA的組成、來源以及功能是極其復(fù)雜的，我們對于sRNA的認(rèn)識才剛剛開始。本研究為sRNA調(diào)控軟骨細(xì)胞分化的分子機(jī)制提供了大量線索和證據(jù)。

另外，由于生物信息學(xué)預(yù)測標(biāo)準(zhǔn)模糊、算法不同產(chǎn)生了大量假陽性數(shù)據(jù)和pre-miRNA長度的可變性使得預(yù)測結(jié)果必然存在遺漏，而且生物信息學(xué)無法預(yù)測沒有基因組或轉(zhuǎn)錄組數(shù)據(jù)庫信息的物種和物種特異的sRNA，所以需要生物信息學(xué)分析與實(shí)驗(yàn)結(jié)合對結(jié)果進(jìn)行驗(yàn)證［17-18］。

［1］NGREITZ JM，SIROKMAN K，MCDONEL P，et al．RNARNA interactions enable specific targeting of noncoding RNAs to nascent Pre-mRNAs and chromatin sites［J］．Cell，2014，159：188-199．

［2］HWARTZ S，BERNSTEIN DA，MUMBACH MR，et al．Transcriptome-wide mapping reveals widespread dynamic-regulated pseudouridylation of ncRNA and mRNA［J］．Cell，2014，159：148-162．

［3］RANOV P，OZSOLAK F，KIM SW，et al．New class of genetermini-associated human RNAs suggests a novel RNA copying mechanism［J］．Nature，2010，466：642-646．

［4］HUTVAGNER G，MCLACHLAN J，PASQUINELLI AE，et al．A cellular function for the RNA-interference enzyme Dicer in the maturation of the let-7 small temporal RNA［J］．Science，2001，293：834-838．

［5］LAURESSERGUES D，COUZIGOU JM，CLEMENTE HS，et al．Primary transcripts of microRNAs encode regulatory peptides［J］．Nature，2015，520（7545）：90-93．

［6］COLLINS LJ．The RNA infrastructure：an introduction to ncRNA networks［J］．Adv Exp Med Biol，2011，722：1-19．

［7］COKUS SJ，F(xiàn)ENG S，ZHANG X，et al．Shotgun bisulphite sequencing of the Arabidopsis genome reveals DNA methylation patterning［J］．Nature，2008，452：215-219．

［8］SUN J，ZHONG N，LI Q，et al．MicroRNAs of rat articular cartilage at different developmental stages identified by Solexa sequencing［J］．Osteoarthr Cartilage，2011，19：1237-1245．

［9］ZHONG N，SUN J，MIN Z，et al．MicroRNA-337is associated with chondrogenesis through regulating TGFBR2 expression［J］．Osteoarthr Cartilage，2012，20：593-602．

［10］LI R，LI Y，KRISTIANSEN K，et al．SOAP：short oligonucleotide alignment program［J］．Bioinformatics，2008，24：713-714．

［11］GIBBS RA，WEINSTOCK GM，METZKER ML，et al．Genome sequence of the Brown Norway rat yields insights into mammalian evolution［J］．Nature，2004，428：493-521．

［12］HAFNER M，LANDGRAF P，LUDWIG J，et al．Identification of microRNAs and other small regulatory RNAs using cDNA library sequencing［J］．Methods，2008，44：3-12．

［13］GLAZOV EA，COTTEE PA，BARRIS WC，et al．A microRNA catalog of the developing chicken embryo identified by a deep sequencing approach［J］．Genome Res，2008，18：957-964．

［14］LU C，JEONG DH，KULKARNI K，et al．Genome-wide analysis for discovery of rice microRNAs reveals natural antisense microRNAs（nat-miRNAs）［J］．PNAS，2008，105：4951-4956．

［15］CHEN X，LI QB，WANG J，et al．Identification and characterization of novel amphioxus microRNAs by Solexa sequencing［J］．Genome Biol，2009，10：R78．2-R78．13．

［16］FOX JG，ANDERSON LC，LOEW FM．Laboratory Animal Medicine［M］．San Diego：Academic Press，2002：1265-1325．

［17］CZECH B，MALONE CD，ZHOU R，et al．An endogenous small interfering RNA pathway in Drosophila［J］．Nature，2008，453：798-802．

［18］ADAI A，JOHNSON C，MLOTSHWA S，et al．Computational prediction of miRNAs in Arabidopsis thaliana［J］．Genome Res，2005，15：78-91．