基于高通量測序的結直腸癌miRNA表達譜分析*

梁高峰，李俊生，何向峰，陳寶安#

1)東南大學醫學院附屬中大醫院血液科 南京210009 2)河南科技大學醫學技術與工程學院 洛陽471003 3)東南大學醫學院附屬中大醫院普外科 南京210009

微小RNA(microRNA，簡稱miRNA)是真核生物內源性的小分子非編碼RNA，能夠通過轉錄后水平調節基因的表達。人miRNA 總數可能占總基因數的1%，人類全部基因的1/3 可能受miRNA 的調控［1］。研究證實，miRNA 與多種疾病有關，如糖尿病［2］、心臟病［3］、癌癥［4］、神經系統疾病［5］等。該研究采用SOLiD 高通量測序技術分析了10 對癌組織和癌旁正常組織miRNA 表達譜，篩選出兩者間的差異miRNA，利用數學模型來綜合分析這些差異表達的miRNA 對其靶mRNA 的調控作用及結直腸癌組織中受miRNA 調控的KEGG 通路和GO 生物過程，為進一步研究結直腸癌發生發展的機制提供理論依據。

1 材料與方法

1．1 材料

1．1．1 標本來源 收集東南大學附屬中大醫院組織庫2010年至2011年存檔的冰凍組織標本20例，包括結直腸癌組織10例、癌旁正常組織10例，均參照文獻［6］標準進行診斷和分類。

1．1．2 主要試劑 Trizol(Qiagen，德國)，mirVanaTMmiRNA 提取分離試劑盒(Ambion，美國)，SOLiDTM小RNA 表達檢測試劑盒(Ambion，美國)。

1．2 總RNA 提取及質量鑒定 將每個冰凍組織標本切5～6 片，片厚8 μm，用Trizol 提取總RNA，然后用mirVanaTMmiRNA 提取分離試劑盒分離小RNA，總RNA 完整性采用10 g/L 瓊脂糖凝膠電泳鑒定。

1．3 SOLiD 測序及分析 使用SOLiDTM小RNA 表達檢測試劑盒將樣品中的小RNA 轉換成雙鏈的cDNA。為避免小RNA 樣品提取過程中產生的miRNA前體斷裂片段的影響，將這些比對上的片段再與miRbase 數據庫中的成熟miRNA 序列進行最后的比對，從中篩選出miRNA 及其SOLiD 測序片段數。

1．4 基于測序數據的數學模型的構建 借助miRNA 靶基因預測數據庫PicTar(http://pictar．bio．nyu．edu)，獲得結直腸癌組織中所檢測到的miRNA的所有靶基因，使用公式［1］計算所有測序得到的miRNA 對mRNA 的抑制率。

其中，RmRNAk表示miRNA 對靶mRNAk的抑制率;CountmiRNAi表示所有測序miRNA 中第i個miRNA 的數量;i，調控mRNAk的miRNA 數量，i=1，2，3，……，l;nij，miRNAi的靶mRNA 的數量;weight，PicTar 得分;Ntotal，小RNA 樣品中所檢測到的總miRNA 數量。

結直腸癌組織中，所有miRNA 集合對靶mRNA抑制率差異的比較采用Z-test 方法［7-8］。P0假設使用如下公式計算:

當Z ＞1．96 或＜－1．96，拒絕P0假設，則結直腸癌組織中miRNA 對mRNA 的抑制率與癌旁正常組織比較，差異有統計學意義。

1．5 DAVID 分析 根據所有miRNA 集合對mRNA 的整體抑制情況，獲得表達水平受miRNA 影響的蛋白。采用DAVID 分析這些蛋白對KEGG 通路和GO 生物過程的影響。

2 結果

2．1 SOLiD 測序結果 結直腸癌組織和癌旁正常組織中小RNA 片段占1．13%和1．14%。比對結果表明，每個樣本所比對上的片段數都達107以上，長度為21～24 nt，尤以22 nt 的片段最多(圖1)。結直腸癌組織和癌旁正常組織所得測序片段的分布具有明顯的相似性，而成熟miRNA 的序列長度集中分布在18～24 nt。結直腸癌組織和癌旁正常組織中分別檢出627 和451個miRNA(表1)。

2．2 結直腸癌組織中差異表達的miRNA 結果見圖2。結直腸癌組織中表達上調的miRNA 中，46個miRNA 的表達高于2 倍;表達下調的miRNA 中，84個miRNA 的表達低于1/2。

Z-test 結果(圖3、表2)表明，結直腸癌組織中較多的miRNA 表達水平受到影響。超過1/2(194/313)的miRNA 表達差異明顯，其中80個miRNA 表達水平明顯上升，114個明顯下降。

基于Z-test 方法所篩選的miRNA 不僅包括了倍數分析所篩選的miRNA，如miR-101、miR-135b、miR-222、miR-451、miR-4315 等，而且還包括部分表達量變化不明顯但測序結果中讀長很大的miRNA，如let-7b、let-7i、miR-125b、miR-126 等。

圖1 與前體序列一致的SOLiD 測序序列的長度分布

表1 SOLiD 測序結果

圖2 結直腸癌組織中miRNA 表達差異的倍數分析

圖3 結直腸癌組織中miRNA 表達差異的Z-test 分析

2．3 miRNA 在染色體上的分布 結果見圖4。

2．4 miRNA 對其靶基因的調節 在結直腸癌組織中有7 273個mRNA 受到miRNA 的調控。在受調控的mRNA 中，865種顯著上調，占11．9%;973種顯著下調，占13．4%。見圖5。

2．5 結直腸癌組織中KEGG 通路和GO 生物過程的系統性調節 DAVID 分析結果表明:在結直腸癌組織中，眾多KEGG 通路受到調控，而只有部分受到顯著調控(圖6)。顯著調控的KEGG 通路主要涉及代謝、遺傳信息處理、環境信息處理和細胞過程四大分支，但其中一些KEGG 通路歸根于環境信息處理中的信號轉導，減弱或增強外界環境對細胞的刺激作用;在顯著調控的環境信息處理途徑中，WNT信號通路和MAPK 信號通路均受到顯著的影響，而這兩者又是結直腸癌發展過程中經常發生異常的信號通路。

而GO 分析［9］的結果表明，在這些變化的miRNA 中，顯著調控的GO 過程包括:細胞內的生物合成、細胞代謝、主動的細胞內加工、細胞周期調控、程序性細胞凋亡、細胞內信號轉導、細胞內的轉運等，可大致分為代謝過程、細胞增殖與細胞周期、細胞死亡、細胞通訊、物質運輸五大類(圖7)，大多數顯著調控的GO 過程歸根于調控過程;進一步結合KEGG 公共數據庫的結果可以發現，一些與結直腸腸癌相關的信號通路，例如Wnt 信號通路、MAPK 信號通路和TGF-β 信號通路，均受到miRNA 的調節。

表2 2種分析方法所得結直腸癌miRNA 表達譜的比較

圖4 癌旁正常組織(左)和結直腸癌組織(右)中miRNA 的染色體分布差異

圖5 結直腸癌組織中miRNA 對靶基因的抑制率

圖6 結直腸癌組織中miRNA 對KEGG 通路的調控

圖7 結直腸癌組織中miRNA 對GO 生物過程的調控

3 討論

結直腸癌的發生涉及多個癌基因、抑癌基因和信號轉導通路的改變［10］。作者采用SOLiD 技術對10 對結直腸癌和癌旁正常組織進行高通量測序，得到了一批顯著差異表達的miRNA。另外，作者比較了倍數分析法和Z-test 2種方法，發現miRNA 的倍數分析法存在著不容忽視的缺點:例如，倍數分析未能引入統計量，不能反映倍數變化的統計學水平。為了全面理解miRNA 在細胞中的調控作用，作者根據高通量測序數據將Z-test 引入到表達譜差異的分析中，借助數學模型綜合分析miRNA 的調控作用，該模型包含樣品中所有檢測到的miRNA。結果表明，miRNA 對部分mRNA 的抑制率達到顯著水平。miRNA 抑制的靶mRNA 中，少部分mRNA 的miRNA 抑制率發生顯著變化。該方法比常規倍數分析法更能真實反映miRNA 的表達變化。

應用不同的miRNA 表達譜篩選方法所得到的篩選結果亦存在區別。通過比較倍數分析法和Ztest 的結果，作者發現2種方法用于鑒別miRNA 表達水平增減趨勢時，所得結果一致。但在篩選表達量變化明顯的miRNA 時，基于Z-test 方法所篩選的miRNA 不僅包括了倍數分析所篩選的miRNA，而且還包括部分表達量變化不明顯但測序結果中拷貝數卻很大的miRNA。這說明在篩選miRNA 時，使用倍數分析法會遺漏部分倍數變化不大但覆蓋度較高的miRNA。這也說明了Z-test 用于miRNA 表達譜分析的優勢:不但考慮miRNA 倍數變化，更考慮其在miRNA 群體中的覆蓋度，最為重要的是具有統計學意義，避免了篩選miRNA 時主觀因素的影響。

結直腸癌發生發展過程中，miRNA 染色體分布的變化可反映結直腸癌發病過程中對miRNA 轉錄的差異調控。同樣，miRNA 的表達差異會對細胞內KEGG 通路和GO 生物過程產生重要影響，進而影響細胞內的信號轉導，如MAPK、WNT、TGF-β 信號通路等［11］。細胞調控miRNA 的差異表達，進而調節細胞過程。對細胞內GO 調控過程的調節可事半功倍地調控細胞效應，是miRNA 調控基因表達的經濟之體現，是對細胞信號通路的主動調節、自適應調節的過程。而受調控的KEGG 通路間的差異，反映了不同病理、生理狀態細胞間的差異。這說明miRNA 表達譜的異常變化對結直腸癌的發生起到重要的調節作用;同時，也說明這些異常表達的miRNA也可能參與其他相關癌癥的生物過程，也就是說，其中的一部分miRNA 是不同癌癥發生過程中都發生變化的［12］。

綜上，作者建立了基于miRNA 表達譜的結直腸癌基因調控系統分析體系，miRNA 表達譜變化是細胞對長期外部刺激的自適應調節過程，結合miRNA表達譜數據，借助生物信息學工具研究miRNA 對KEGG 通路和GO 生物過程的調控，并可據此系統地分析不同疾病過程中miRNA 表達譜的變化。

［1］Lim LP，Glasner ME，Yekta S，et al．Vertebrate microRNA genes［J］．Science，2003，299(5612):1540

［2］Locke JM，Da Silva Xavier G，Dawe HR，et al．Increased expression of miR-187 in human islets from individuals with type 2 diabetes is associated with reduced glucosestimulated insulin secretion［J］．Diabetologia，2014，57(1):122

［3］Kuster DW，Mulders J，Ten Cate FJ，et al．MicroRNA transcriptome profiling in cardiac tissue of hypertrophic cardiomyopathy patients with MYBPC3 mutations［J］．J Mol Cell Cardiol，2013，65:59

［4］Slattery ML，Herrick JS，Mullany LE，et al．An evaluation and replication of miRNAs with disease stage and colorectal cancer-specific mortality［J］．Int J Cancer，2015，137(2):428

［5］Maciotta S，Meregalli M，Torrente Y．The involvement of microRNAs in neurodegenerative diseases［J］．Front Cell Neurosci，2013，7:265

［6］Hamilton SR，Aaltonen LA．World Health Organization classification of tumor:pathology and genetics of tumours of the digestive system［M］．Lyon:IARC Press，2003:103

［7］Pisanic TR，Blackwell JD，Shubayev VI，et al．Nanotoxicity of iron oxide nanoparticle internalization in growing neurons［J］．Biomaterials，2007，28(16):2572

［8］Kal AJ，Van Zonneveld AJ，Benes V，et al．Dynamics of gene expression revealed by comparison of serial analysis of gene expression transcript profiles from yeast grown on two different carbon sources［J］．Mol Biol Cell，1999，10(6):1859

［9］Borchert GM，Lanier W，Davidson BL．RNA polymerase Ⅲtranscribes human microRNAs［J］．Nat Struct Mol Biol，2006，13(12):1097

［10］黃辰，劉利英，倪磊，等．腫瘤miRNAs 調控機制的研究進展與展望［J］．西安交通大學學報:醫學版，2015，36(3):221

［11］Tang M，Mao D，Xu L，et al．Integrated analysis of miRNA and mRNA expression profiles in response to Cd exposure in rice seedlings［J］．BMC Genomics，2014，15:835

［12］Kara M，Yumrutas O，Ozcan O，et al．Differential expressions of cancer-associated genes and their regulatory miRNAs in colorectal carcinoma［J］．Gene，2015，567(1):81