基于高通量测序的结直肠癌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