宫颈癌差异表达基因中内含子–mRNA匹配模式的计算分析

《Biology》：Intron–mRNA Matching Patterns in Differentially Expressed Genes in Cervical Cancer

　　

宫颈癌（CC）以广泛的转录组失调为特征，然而在差异表达基因（DEGs）中，mRNA与其对应内含子之间的序列匹配模式仍缺乏充分表征。本研究旨在描述CC差异表达基因参考序列中经计算预测的内含子–mRNA匹配特征。研究人员利用基因表达综合数据库（GEO）和癌症基因组图谱（TCGA）数据集识别DEGs，并从Ensembl人类参考基因组中检索相应的mRNA与内含子序列。采用反向互补结合Smith–Waterman局部比对，分析预测匹配片段在全长mRNA以及翻译起始区和终止区内的位置分布，并进一步按匹配片段的鸟嘌呤和胞嘧啶（GC）含量进行描述性分层分析。在标准化全长mRNA坐标上，高GC片段在靠近5′端的某些区段表现出较高的匹配频率，而低GC片段的匹配频率向3′端升高。在翻译起始区内，预测匹配频率在AUG密码子上游总体高于下游，低GC组在起始密码子附近呈现局部峰值。在翻译终止区内，低GC组在终止密码子下游呈现明显峰值，而高GC组未出现类似明显的下游峰值，且匹配频率的位置波动较大。总体而言，预测匹配分布在所分析的各区域和GC组之间存在差异。本研究在该CC差异表达基因参考序列集合中提供了内含子–mRNA匹配模式的序列水平描述；其生物学意义尚需进一步验证。

研究背景方面，内含子曾被视为前体mRNA剪接产生的非功能性副产物，如今已被认识为真核基因表达的重要调控因子，参与转录效率、核输出和翻译等多个RNA加工与转录后调控环节，并与基因组演化、蛋白质组多样化和基因组稳定性维持相关。内含子与mRNA内互补序列区域的潜在相互作用是内含子功能的重要方面，此类相互作用可能参与转录本加工、稳定性和翻译的调控，其破坏与异常剪接、非编码RNA活性失调及疾病相关转录组改变有关。既往研究采用反向互补序列匹配结合Smith–Waterman局部比对方法，将mRNA与对应内含子序列间相互作用概率最高的区域定义为"最优匹配片段"，并系统刻画了这些片段的序列特征及匹配强度沿mRNA序列的分布模式，在包括人类、小鼠、大鼠、蜜蜂和河豚等27个物种以及染色体1基因组序列的多物种验证中揭示了内含子–mRNA匹配特征的跨物种保守性。既往研究还提示低GC预测匹配区更常见于3′非翻译区（3′UTR），而高GC预测区优先富集于5′非翻译区（5′UTR），且编码区内的互补相互作用弱于非翻译区。然而，宫颈癌差异表达基因（DEGs）参考mRNA与其对应内含子之间的匹配特征仍缺乏充分表征，包括按片段GC含量分层后观察到的变异。在此背景下，研究人员利用GEO和TCGA转录组数据集识别宫颈癌差异表达基因，从Ensembl人类参考基因组而非组织特异性转录本组装中检索用于序列分析的mRNA和内含子序列，对富含AU和相对富含GC的片段组进行比较，系统描述预测的内含子–mRNA匹配在转录本不同位置和不同GC含量组间的分布差异，为后续采用适当序列对照和功能实验进行评估提供描述性序列水平信息。

在关键技术方法方面，研究人员从GEO数据库获取GSE138080、GSE39001、GSE46857、GSE63514和GSE63678五个数据集（共166例样本，含111例肿瘤样本和55例正常样本），并从TCGA获取304例CC样本和3例正常宫颈组织样本的表达数据。采用R语言的limma包进行背景校正和标准化，使用sva包的ComBat函数校正数据集间批次效应，以调整p值<0.05且|logFC|>1为标准识别DEGs，取GEO与TCGA交集得到84个重叠DEGs。从Ensembl数据库检索序列，经质量控制后保留84个基因的486条内含子序列和84条mRNA序列。采用Smith–Waterman局部比对算法（EDNAFULL评分矩阵，空位开放罚分50.0，空位延伸罚分5.0）对反向互补的内含子序列与对应mRNA进行比对，识别最优匹配区域，并按GC含量≤30%和≥50%分为低GC组和高GC组。以人类染色体1基因组基因序列构建背景数据集进行对照分析。

研究结果方面，第一，不同GC含量下匹配频率沿mRNA序列的分布：将靶mRNA标准化为100个核苷酸长度并计算各标准化位置的匹配频率，总体分布在约20 nt和80 nt处出现两个明显局部极小值，分别大致对应5′UTR和3′UTR；低GC组匹配频率谱在3′UTR区域呈现显著峰值，高GC组则在5′UTR区域表现出更强富集，表明预测的内含子–mRNA匹配位置分布在不同GC含量条件下存在差异。第二，不同GC含量下翻译起始区匹配频率分布：以AUG密码子为位置0构建包含上游30个核苷酸和下游60个核苷酸的分析窗口，总体分布中AUG上游匹配频率普遍高于下游；分层后低GC组在起始密码子附近呈现明显局部峰值，高GC组AUG上游匹配频率明显高于下游且下游区域波动更明显。第三，不同GC含量下翻译终止位点匹配频率分布：以终止密码子（UAA、UAG、UGA）为位置0构建包含上游30个核苷酸、完整终止密码子及下游90个核苷酸的分析窗口，总体分布显示匹配频率从上游向下游呈升高趋势；低GC组在下游区域快速升高并在终止密码子后紧邻区间达到显著峰值后逐渐下降，高GC组下游匹配频率略低于上游水平且围绕终止密码子呈现多个局部波动。第四，匹配频率差异分析：由于同一基因可产生多个内含子且mRNA相邻位置可被同一比对片段覆盖，内含子水平和位置水平观察不作为独立重复，匹配谱以描述性方式总结；全长mRNA序列匹配频率低GC组低于高GC组，翻译起始区匹配频率低GC组低于高GC组，而翻译终止区匹配频率低GC组高于高GC组，表明GC含量组间差异随所分析转录本区域而变化；翻译起始区的组间差异最大且效应量很大，全长mRNA和翻译终止区的比较效应量较小且组间差异方向相反。

讨论部分总结：既往研究表明内含子与对应mRNA间经计算预测的互补序列匹配呈非随机分布，并可能随转录本位置和序列组成而变化；本研究将这一描述性分析扩展至宫颈癌数据集中差异表达基因的参考序列，并与人类染色体1背景数据集进行比较。尽管两个基因集共享若干总体趋势，宫颈癌相关转录本集表现出独特的区域模式。在翻译起始位点附近，两者均显示AUG上游匹配频率高于下游且低GC片段在起始密码子附近局部富集，但染色体1背景中AUG周围的转变相对平缓，而宫颈癌DEGs在AUG边界处呈现更明显的特征差异，尤其在高GC组。在翻译终止位点附近，两组低GC预测匹配频率在终止密码子上游均较低并向下游升高，但宫颈癌DEGs高GC片段的下游匹配频率更低且位置变异更大。全长mRNA分析进一步显示高GC预测匹配片段优先分布于5′UTR而低GC片段更常见于3′UTR，与既往报道基本一致。研究存在若干局限：基于公共转录组数据集和计算序列分析，预测匹配特征未经实验验证；分析仅限于GEO和TCGA重叠DEGs，染色体1仅作为背景参考而非完整正常基因组基线；观察到的差异应主要解释为该特定DEG集的固有位置和统计序列特征，而非物理相互作用或功能调控的直接证据。未来需采用全基因组蛋白编码背景、构建基于基因长度、GC含量、内含子数目和UTR长度严格匹配的随机基因集并进行多次置换检验，结合更广泛基因组参考、独立临床队列、RNA结构分析和功能验证，以确定这些计算特征是否与RNA调控过程相关及其作为宫颈癌相关基因转录本分子特征的潜在价值。

研究结论翻译：本研究从序列位置和GC含量角度系统刻画了宫颈癌DEGs的mRNA与其对应内含子序列间的匹配频率特征。结果表明，内含子–mRNA匹配在不同功能转录本区域间存在差异，并在该DEG集内表现位置异质性。高GC预测匹配片段优先分布于5′UTR，而低GC预测匹配片段更常见于3′UTR。在翻译起始和终止位点周围也观察到不同的位置模式。与染色体1背景参考相比，该宫颈癌相关转录本集在AUG边界周围呈现更明显的特征差异，并在翻译起始区和终止区表现高、低GC片段的不对称分布。这些发现与既往研究基本一致，即内含子与mRNA间的互补区域沿转录本不均匀分布并受序列组成影响。基于这些探索性计算发现，未来研究可采用组成保留序列对照和独立基因集进一步检验观察到的内含子–mRNA匹配模式的可重复性，整合公共转录组数据集和独立临床队列数据以探索这些序列水平特征与内含子保留等RNA加工事件间的潜在关联，并可通过结合靶向体内RNA结构探测（如DMS-MaPseq）与RNA–RNA相互作用图谱（如LIGR-seq）的进一步计算分析和实验研究，刻画局部RNA二级结构并探究涉及选定候选区域的潜在RNA–RNA相互作用，从而为后续评估预测匹配模式的生物学相关性提供依据。