基于高通量测序的结直肠癌组织中环状RNA 功能研究与预后分析

傅中懋，罗 再，戎泽印，章建明，李腾飞，余志龙，黄 陈

上海交通大学附属第一人民医院普外科，上海200080

结直肠癌是当今世界最常见的恶性消化道肿瘤之一，每年新增结直肠癌患者人数为185 万，新增死亡人数为88 万［1］。在美国，结直肠癌的发病率与病死率在男性和女性罹患的恶性肿瘤中均排第三［2］。据我国肿瘤中心登记的数据显示，2015 年结直肠癌在总新发肿瘤中发病率与病死率均排第五［3］。结直肠由可治愈的局部恶性病变发展到转移性结直肠癌需要经过一个漫长的时间窗，因此早期诊断的意义重大。目前，结直肠癌检查主要为直肠指诊和结肠镜检查。由于侵入式检查会给患者带来不适感，患者往往不愿意定期检查［4］。计算机断层扫描检查虽然为非侵入性，但由于其局限于形态成像，可能会将残余粪便误诊为肿瘤而导致假阳性结果［5］，因此血源性生物标志物被认为是一种侵入性最小的选择。近年来，癌胚抗原（carcinoembryonic antigen，CEA）、糖类抗原19-9（carbohydrate antigen 19-9，CA19-9）、CA242 等已被证实与结直肠癌的诊断和预后相关［6］。但是这些指标的特异性往往有所欠缺，例如在其他恶性肿瘤如卵巢癌、胰腺癌，甚至在良性疾病如炎症性肠病中，CEA 的水平也会升高［7］。因此为了改善患者预后，深入研究结直肠癌的发生机制，挑选新型生物标志物至关重要。

随着二代测序技术与生物信息技术的发展，人们对环状RNA（circular RNA，circRNA）的研究与认识越发深入。circRNA在哺乳动物中普遍存在，作为一种不具有5'端帽子结构和3'端多腺苷酸尾部结构的特殊环状闭合RNA，其结构稳定，不易被核酸外切酶降解［8］。由于circRNA的广泛性、组织发育阶段特异性、高度保守性等特点，其具有成为新型临床诊断标志物的潜力［9］。circRNA因其在真核生物中所具有的细胞或组织特异性以及发育阶段特异性，已被证实参与不同疾病的致病过程［10］。同时，circRNA 结构特异性决定了其特有的生物学功能，如调节蛋白翻译、基因转录等［11］。目前，越来越多的证据表明circRNA 作为竞争性内源性RNA 在细胞质参与miRNA 结合竞争，从而在致癌过程中发挥着重要作用，如circNHSL1 通过miR-1306-3p/SIX1/vimentin 途径促进癌症进程［12］。但是当前在结直肠癌中关于表达差异circRNA 的选择、circRNA 的生物学功能及其与患者预后关系的研究仍很少。

本研究通过对12 对结直肠癌患者的新鲜结直肠癌组织和癌旁组织进行高通量测序分析，筛选出表达差异的circRNA 并进行功能预测，借助实时荧光定量PCR（realtime fluorescent quantitative PCR，RT-qPCR）技术验证测序结果的准确性，构建circRNA 与微RNA（microRNA，miRNA）、mRNA 的互作网络，并初步研究差异表达的circRNA与患者预后的关系，以期为结直肠癌发生机制与新型临床标志物的筛选提供研究基础。

1 对象与方法

1.1 研究对象

收集2017 年1 月—2018 年12 月在上海交通大学附属第一人民医院普外科行结直肠癌根治术的12 例患者的癌组织和癌旁组织（距离癌组织>5 cm）用于高通量测序。收集2015年1月—2017年12月行结直肠癌根治术的44 例患者的癌组织和癌旁组织（距离癌组织>5 cm），用于RTqPCR 实验。纳入标准：①符合结直肠癌诊断标准。②患者病历资料完整，随访数据完整。③术前未经过新辅助放射治疗和化学治疗。④患者无其他系统严重疾病。排除标准：①合并其他部位原发恶性肿瘤患者。②患者临床、随访资料不完整。组织标本取出后即刻放入液氮盒或−80 ℃冰箱保存。研究获上海交通大学附属第一人民医院伦理委员会批准（批件编号2019SQ274）。所有患者或其家属均签署知情同意书。

1.2 circRNA高通量测序

RNA 高通量测序由上海欧易生物医学科技有限公司提供。使用mirVana miRNA （AM1561，中国Ambion）分离试剂盒提取样品总RNA，用ribo-zero 试剂盒（RS-122-2301，美国Illumina）捕获并去除核糖体RNA，用RNA 碎片化试剂（AM8740，美国ABI）将RNA 片段化，配置一链合成体系通过六碱基随机引物合成第一链cDNA，随后在此基础上配置二链合成体系，用上标Ⅱ反转录酶（18064014，美国Invitrogen）增加cDNA 产量，其中dUTP 代替dTTP 合成第二链cDNA，用UNG 酶法将富含dUTP 的第一链消化，将第二链cDNA 进行纯化，将纯化后的cDNA 末端补平，加A 尾修饰和加测序接头，qPCR 扩增构建测序文库。将构建好的文库使用安捷伦2100 生物分析仪（5067-1511，美国Aglient）检测RNA完整性，RNA 完整性参数（RIN）≥7 的样品进行后续分析，在Illumina 测序平台（HiSeq 2500，美国Illumina）对这些文库进行测序。最后根据筛选标准：基因表达值倍数变化（fold change，FC）>2 或FC<−2，P<0.05 进行差异显著circRNA筛选。

1.3 差异表达circRNA的功能预测分析

1.3.1 差异表达circRNA GO 富集分析 基因本体（gene ontology，GO）分析主要通过分析细胞组分、分子功能、生物过程3个方面来描述基因行使的主要功能。在得到差异表达的circRNA 后，根据circRNA 来源基因的GO 功能注释信息，计算每个GO条目中差异表达的circRNA数量，超几何分布检验差异circRNA 在GO 条目中富集的程度，以P<0.05 为阈值筛选显著富集的GO 条目，结合GO 注释结果，预测差异表达的circRNA所行使的主要生物学功能。

1.3.2 差异表达circRNA KEGG 信号通路分析 circRNA信号通路分析主要依赖于京都基因和基因组百科全书（Kyoto Encyclopedia of Genes and Genomes，KEGG）公共数据库，通过基因的KEGG 注释结果研究疾病信号通路和转导。在得到差异表达的circRNA 后，计算统计每个信号通路中差异circRNA 来源基因的富集程度，根据P<0.05进行KEGG富集通路筛选分析。

1.4 RT-qPCR验证分析

在筛选出的具有表达差异的circRNA 中各挑选2 个显著上调与显著下调的RNA，在44 对结直肠癌组织与癌旁组织中进行RT-qPCR 验证circRNA 含量；其中挑选的上调的circRNA 为表达差异倍数最大前5 中的hsa_circ_0023984 与hsa_circ_0008192，挑选的下调的circRNA 为表达差异倍数最大的hsa_circ_0020093 与hsa_circ_0069922。4 个circRNA 和内参甘油醛-3-磷酸脱氢酶（glyceraldehyde-3-phosphate dehydrogenase，GAPDH） 基因的引物序列如表1 所示。首先使用PrimeScriptTMRT 试剂盒（RR036A，日本TaKaRa） 反转录500 ng RNA 为cDNA，然后用SYBR®Premix Ex TaqTM（RR820A，日本TakaRa）扩增cDNA 进行RT-qPCR。将GAPDH 作为内部对照，用公式2-ΔΔCt计算circRNA的相对表达水平。

表1 RT-qPCR 引物序列Tab 1 Primer sequences of RT-qPCR

1.5 circRNA-miRNA靶向作用预测

挑选上述实验的差异表达circRNA， 通过CircInteractome 分析网站挑选具有潜在结合circRNA 能力的miRNA，根据结合位点个数及context+score percentile评分综合筛选出最优的的miRNA。接着通过MicroRNA、TargetScan 与TangetMiner 网站对筛 选 的miRNA 进行靶基因预测，取3 个预测结果的交集，得到miRNA 靶基因。最后根据度分析指标在Cytoscape 软件上构建circRNAmiRNA-mRNA网络。

1.6 患者预后分析

将表达上调的hsa_circ_0023984、hsa_circ_0008192与表达下调的hsa_circ_0020093 与hsa_circ_0069922 按照表达水平高低与44 例结直肠癌患者临床资料相匹配，通过分析总生存期（overall survival，OS）研究表达上调或下调的circRNA与患者预后之间的联系。

1.7 统计学方法

采用SPSS 23.0 软件进行数据分析。测量circRNA 表达量的RT-qPCR 实验重复3 次，最终结果用x±s 表示。结直肠癌组与癌旁正常组之间circRNA 表达水平的比较采用配对t检验。P<0.05认为差异有统计学意义。

2 结果

2.1 结直肠癌组织中circRNA的表达分析

经高通量测序后共检测出20 924个circRNA，其中有9 796个circRNA 在circBase数据库中，第6对样本预测表达circRNA数量最多（图1A、B）。相较于正常的结直肠癌旁组织，在结直肠癌组织中异常显著表达的circRNA共有373 个（︱log2FC︱>1，P<0.05），其中上调显著的有243个，下调显著的有130个。将有显著差异表达的circRNA数据制成火山图（图1C）和聚类分析热图（图1D），在RNA测序结果中筛选上调与下调最显著的5个circRNA，具体信息见表2。其中hsa_circ_0079534 表达上调最多，FC 为14.24；hsa_circ_0020093下调最明显，FC为−8.27。

2.2 结直肠癌组织中circRNA的功能分析

图1 结直肠癌组织及癌旁组织中差异表达circRNA的可视化分析Fig 1 Visual analysis of the differential expression of circRNAs in colorectal cancer and adjacent tissues

表2 circRNA测序结果中表达差异最显著的10种circRNATab 2 Ten circRNAs with the most significant expression in the sequencing result

通过GO 富集分析与KEGG 通路分析预测差异表达circRNA 的功能。部分差异表达的circRNA 对应1 个GO注释，分别根据生物过程、细胞定位、分子功能3个部分对应将最富集的前10条筛选出来。差异显著circRNA在生物过程预测显示较为富集的是细胞排泄、细胞内源性pH的调节、对cAMP合成过程的负向调控以及对草酸运输调控等；其细胞定位主要在生殖核、纺锤体中央区、细胞质膜以及刷状缘膜等处；差异表达的circRNA分子功能分析提示结直肠癌发病机制与部分转运蛋白活性密切相关，如硫酸盐跨膜转运蛋白与次级活性硫酸盐跨膜转运蛋白活性、草酸跨膜转运蛋白活性、G蛋白偶联受体活性等（图2A）。KEGG通路分析表明，在circRNA相关通路中最富集的通路是对肿瘤微环境中蛋白聚糖的调控，同时circRNA调控的血管内皮生长因子（vascular endothelial growth factor，VEGF）信号通路、趋化因子信号通路、红细胞白血病病毒癌基因同源物（v-erb-b2 avian erythroblastic leukemia viral oncogene homolog，ErbB）信号通路与结直肠癌发生息息相关，circRNA相关的白细胞经内皮迁移也可能参与了结直肠癌的致病过程（图2B）。

图2 差异表达circRNA 的GO分析和KEGG分析Fig 2 GO analysis and KEGG analysis of differentially expressed circRNAs

2.3 RT-qPCR 检 测 结 果 与circRNA-miRNA-mRNA 网 络分析

在所有差异表达的circRNA 中挑选出表达差异最显著的4 个circRNA，分别为hsa_circ_0023984、hsa_circ_0008192、hsa_circ_0020093 与hsa_circ_0069922。前两者在测序结果中预测表达上调，后两者预测表达下调。我们在44 对结直肠癌及癌旁组织样本中使用RT-qPCR 技术进一步验证circRNA 测序结果的准确性。用于RT-qPCR验证的44 例患者年龄为35～84 岁，男性26 例，女性18例；患者分期依据国际抗癌联盟第8 版结直肠癌TNM 分期标准：Ⅰ期3 例，Ⅱ期19 例，Ⅲ期20 例，Ⅳ期2 例。结果（图3A）显示：hsa_circ_0023984在结直肠癌肿瘤组织与癌旁组织中的相对表达量分别为0.179 6±0.446 6 与0.012 2±0.017 4 （t=2.480，P=0.017）；hsa_circ_0008192在结直肠癌肿瘤组织与癌旁组织中的相对表达量分别为0.089 4±0.186 7 与0.032 6±0.052 4（t=2.190，P=0.034）；hsa_circ_0020093 在结直肠癌肿瘤组织与癌旁组织中的相对表达量分别为0.057 8±0.146 1 与0.457 6±0.625 5（t=4.725，P<0.000 1）；hsa_circ_0069922 在结直肠癌肿瘤组织与癌旁组织中的相对表达量分别为0.030 9±0.030 9 与0.409 5±0.279 4（t=8.645，P<0.000 1）。通过RT-qPCR 证实了circRNA测序结果的准确性。

利用CircInteractome、MicroRNA、TargetScan 等网站对上调明显的hsa_circ_0023984 和hsa_circ_0008192 及下调显著的hsa_circ_0020093 和hsa_circ_0069922 进行下游miRNA 及最终调控的mRNA 预测。根据预测结果优先选择相关系数最高的4 个miRNA 和20 个mRNA 构建circRNA-miRNA-mRNA 互作网络，通过Cytoscape 软件将其可视化为核心网络图（图3B）。该图显示含锌指和BTB 结构域蛋白1（zinc finger and BTB domain containing 1， ZBTB1）、 SEC23 同 源 物A （SEC23 homolog A，SEC23A）、脆性X 智力低下蛋白翻译调节蛋白1（FMRP translational regulator 1，FMR1）、烟酰胺磷酸核糖转移酶（nicotinamide phosphoribosyltransferase，NAMPT）、具有植物同源结构域的锌指蛋白20 样蛋白1（PHD finger protein 20 like 1，PHF20L1）等为靶向关系最多的mRNA原始基因。

2.4 差异表达circRNA与患者预后的关系

结合患者的临床资料，通过Kaplan-Meier 生存分析显示，hsa_circ_0008192 在患者肿瘤组织中高表达时其OS 较差（P=0.026 2），同时对于hsa_circ_0023984 高表达的患者OS 也趋于变短（P=0.039 8）（图4A、B）。hsa_circ_0020093 在患者肿瘤组织中高表达时其OS 较好（P=0.000 1，图4C），对于同样表达下调的hsa_circ_0000111 也有类似的结果（P=0.002 2，图4D）。我们的研究结果提示，差异表达的circRNA 可能是判断结直肠癌患者预后的一个独立且有价值的指标。

3 讨论

尽管当前结直肠癌总的发病率与病死率持续下降，但进展期的结直肠癌依旧给公共卫生体系带来极大的负担；进一步研究结直肠癌的发病机制对早期识别结直肠癌与指导治疗有着重要的意义［13］。circRNA 作为单链共价闭合的RNA 分子，具有与线性RNA 完全不同的结构，其3' 末端和5' 末端连接在一起构成了特殊的环状结构［14］。circRNA 的稳定性很大程度由其结构决定。由于circRNA 缺乏类似线性RNA 的3'尾和5'帽子结构，其在很大程度避免了脱帽反应、脱腺苷化反应等相关的线性RNA 降解方式，因此细胞质中circRNA 的半衰期可达48 h以上而线性RNA半衰期均仅为10 h［15-16］。与此同时，circRNA 在不同组织及不同发育阶段存在特异性表达［9］。目前越来越多研究［12，17］发现，circRNA 与肿瘤的发生与促进肿瘤的进展密切相关。在肿瘤组织中差异表达的circRNA 如胃癌中高表达circMLLT10 的患者预后明显更差，提示circRNA 可作为判断肿瘤患者预后的生物标志物［18］。由于circRNA 所具有的稳定性、细胞类型或发育阶段特异性、结构保守性，使circRNA 可能成为新型生物诊断标志物［9］。因此，研究circRNA 在结直肠癌中的作用意义重大。

本研究通过高通量测序获得了结直肠癌circRNA 表达数据库，从2 万多个circRNA 中筛选了具有显著表达差异的373 个circRNA。在对表达差异的circRNA 进行功能分析时，我们通过GO分析和KEGG通路注释分析预测了这些circRNA 的潜在功能。在GO 分析中显示，在生物学过程中最富集的是细胞排泄与对细胞内源性pH 的调节，分别各有6 种circRNA 与此过程相关。在细胞定位分析中发现差异表达的circRNA 大多都分布在细胞质膜上，这与circRNA主要位于细胞质或储存于外泌体中相一致［19］。分子功能分析提示结直肠癌的发病机制与部分转运蛋白活性密切相关；目前已有研究已经证实了在结直肠中circRNA 与转运蛋白如G 蛋白偶联受体之间存在联系［20］。将这些差异表达的circRNA 进行KEGG 信号通路富集分析发现，共有173 条通路与结直肠癌发生机制密切相关，其中包括了VEGF 信号通路、趋化因子信号通路、ErbB信号通路等。对于结直肠癌中的VEGF 信号通路，Li等［21］发 现circCCT3 通 过 充 当miR-613 的 海 绵 来 调 节VEGF 信号，从而促进结直肠癌的远处转移。此外Zhong等［22］发现circMYLK 可充当miR-29a 的竞争内源性RNA（competing endogenous RNA，ceRNA），通过激活血管内皮生长因子A （vascular endothelial growth factor A，VEGFA）/血管内皮生长因子受体2（vascular endothelial growth factor receptor 2，VEGFR2）和下游Ras/细胞外信号调节蛋白激酶（extracellular signal-regulated protein kinase，ERK）信号通路来促进上皮间质转换，从而导致膀胱癌的进展。

图3 circRNA差异表达及circRNA-miRNA-mRNA网络图Fig 3 Differential expression of circRNA and circRNA-miRNA-mRNA network plot

图4 4种差异表达circRNA的临床预后分析Fig 4 Clinical prognosis analysis of four differentially expressed circRNAs

此外，我们在circRNA 表达数据库挑选了显著上调和显著下调共4个circRNA在44对新鲜的结直肠癌组织及其癌旁组织中进行了RT-qPCR分析，验证结果与在测序结果中验证基本一致。结合患者的临床资料，初步研究了高表达与低表达的circRNA与患者预后的关系。结果发现结直肠癌组织中高表达的hsa_circ_0008192 和hsa_circ_0023984 患者OS 较差，而高表达的hsa_circ_0020093 和hsa_circ_0069922 患者OS 较好。但考虑患者的例数较少，我们准备在后续实验中扩大样本量继续验证。

在目前已知的circRNA 中，仅有少部分同时含内部核糖体进入位点（internal ribosome entry site，IRES）和小开放阅读框（small open reading frame，sORF）结构才可翻译并编码多肽。如在肝癌中高表达的circβ-catenin 可以通过编码多肽β-catenin-370aa 来减弱β-catenin 磷酸化，从而激活Wnt/β-catenin 通路，促进肝癌的进展［23］。在结直肠癌中，差异表达的circLgr4 通过编码多肽Lgr4-aa 来激活Wnt/β-catenin 通路经上皮间质转换促进结直肠癌转移［24］。大多数circRNA通过与miRNA形成ceRNA，负向调控miRNA 的功能使其靶基因表达失调，参与并促进肿瘤发生和进展［21-22，25］。因此，我们对挑选的4 个circRNA进行了circRNA-miRNA-mRNA 的靶向预测，发现了ZBTB1、SEC23A、FMR1、NAMPT、PHF20L1 等为靶向关系最多的mRNA原始基因。

在这些靶向基因中，由烟酰胺磷酸核糖转移酶（nicotinamide phosphoribosyltransferase，NAMPT）基因编码的蛋白是烟酰胺腺嘌呤二核苷酸（NAD）生物合成过程中的关键酶。它可在组织或细胞水平影响机体代谢、细胞增殖分化与凋亡，尤其是衰老过程。NAMPT 在人类多种恶性肿瘤中上调，包括乳腺癌、前列腺癌、甲状腺癌、胃癌和多种造血系统恶性肿瘤［26-27］。与此同时，在一些恶性肿瘤如胃癌、甲状腺癌和前列腺癌中，NAMPT的表达与肿瘤侵袭和化学治疗耐药相关［28］。有文献［29］也报道了在结直肠癌中NAMPT 明显上调，且NAMPT 通过控制干细胞信号途径增加了癌症起始细胞的数量，从而导致了结直肠癌的发生。通过对circRNA 网络的分析，hsa-miR-576-5p 与hsa-miR-1228 与NAMPT 关 系 密 切，因此我们推测hsa_circ_0008192 或hsa_circ_0023984 高表达患者OS 总体变短可能与circ0023984-miR-576-5p-NAMPT或circ0008192-miR1228-NAMPT调控轴相关。

综上，本研究利用高通量测序筛选了差异表达的circRNA并对相应功能进行分析，用RT-qPCR验证了测序结果的准确性，结合患者资料分析了预后，此外构建了部分与结直肠癌密切相关的circRNA-miRNA-mRNA 调控轴。在后续研究中，我们将从体内外方向共同验证和完善本研究预测的circRNA-miRNA 靶向结合网络，旨在进一步研究结直肠癌发病机制，从而为结直肠癌的早期诊断与治疗提供新的研究思路。