近两年，以“Bioinformatics analysis”为关键词进行检索，大部分文章都是关于预后模型的构建。看了这些文章后，发现思路大同小异。特殊一些的，结合了免疫，自噬，铁死亡等相关基因进行研究。如何在大众的文章中独树一帜呢？是个值得思考的问题。就在这时，一个老朋友，ceRNA，进入我们的视线，看似过气的它能否与预后模型擦出新的火花？让我们来看看吧。
本次分享的文献题目是“Identification of Hub lncRNAs Along With lncRNA-miRNA-mRNA Network for Effective Diagnosis and Prognosis of Papillary Thyroid Cancer”（Hub lncRNA和lncRNA-miRNA-mRNA网络的鉴定对甲状腺乳头状癌的有效诊断和预后的帮助），发表于
“Frontiers in pharmacology”，影响因子为5.810分。
分析流程图
从GEO数据库中获得4个GEO数据集，提取lncRNA数据；TCGA数据库中提取miRNA和mRNA数据。经过差异分析，筛选出差异lncRNA，miRNA和mRNA。预测lncRNA-miRNA以及mRNA-miRNA关系对，构建ceRNA网络。基于该网络中的RNAs，构建总生存期（OS）和无病生存期（DFS）两个预后模型，随后进行验证。
主要结果
Hub lncRNAs的识别和验证
四个数据集获得的差异lncRNA数量：GSE29265，8个差异lncRNA（3上调和5下调）；GSE3678，6个差异lncRNA（3上调和3下调）；GSE33630，9个差异lncRNA（5上调和4下调）；GSE3467:1下调。RRA方法整合上述的差异基因，共获得5个hub lncRNAs。
这5个hub lncRNAs是SLC26A4-AS1, RNF157-AS1, NR2F1-AS1, ST7-AS1,和MIR31HG。随后，这些基因的表达水平在四个GEO数据集以及GEPIA中进行了验证，结果表明NR2F1-AS1和MIR31HG在乳头甲状腺癌中上调，而其余的lncRNAs在癌症组织中下调。
为了进一步验证这5个hub lncRNAs用于标记乳头甲状腺癌的潜力，我们基于4个GEO数据集和TCGA进行ROC曲线分析，验证这5个lncRNAs的诊断性能。结果表明对于每个lncRNA，AUC值均大于0.75。因此，这5个hub lncRNAs都是可靠的，并且是标记乳头甲状腺癌（PTC）的潜在生物标志物。
lncRNA-mRNA共表达网络
进一步，筛选了与hub lncRNAs共表达的mRNAs。通过Pearson相关性分析，研究了TCGA中5个lncRNAs和2716个差异mRNAs之间的相互作用。构建的lncRNA-mRNA网络包括5个lncRNAs和647个mRNAs。出图中可以看出，SLC26A4-AS1，RNF157-AS1和ST7-AS1与较多的mRNAs具有相互作用，这可能表明三个lncRNA之间通过共享共同的mRNA而存在协同作用。但NR2F1-AS1单独具有最多的相互作用mRNA, MIR31HG具有最少的相互作用mRNA。
为了探究该网络中mRNA的生物功能，我们随后进行了功能富集分析。这些mRNA被分成常见mRNA和特异mRNA。如果一个mRNA与两个以上（包括两个）lncRNAs有相关性的话，定义为常见mRNA；而只与一个lncRNA相关的mRNA被认为是特异mRNA。
常见的mRNA参与甲状腺激素的生成（GO-BP）和动力蛋白中间链结合功能（GO-MF）。 它们通常与甲状腺激素合成（KEGG）和甲状腺疾病有关（DO）。
特异性mRNA主要参与轴突发生和跨膜受体蛋白酪氨酸激酶活性。它们可能与p53信号通路和乳头状癌有关。
lncRNA-miRNA-mRNA ceRNA网络
通过预测miRNA，得到了lncRNA-miRNA以及miRNA-mRNA关系对，整合后构建lncRNA-miRNA-mRNA ceRNA网络。该网络包括17miRNAs，5个lncRNAs和68个mRNAs。
随后，对该网络中的mRNA的功能进行预测。结果指出这些mRNA富集在皮肤形态发生和皮质类固醇（GO-BP）以及血小板衍生生长因子结合和细胞外基质（GO-MF）。另外，它们与甲状旁腺功能亢进（KEGG）和甲状旁腺疾病（DO）有关。
预后signature的构建和生存分析
为了识别具有预后特征的潜在RNAs，我们对网络中5个lncRNAs，17个miRNAs和68个mRNAs进行了单变量Cox比例风险回归分析。P<0.05作为阈值筛选与OS以及DFS有关的基因。共有9个mRNAs与DFS有关；4个miRNAs和9个mRNAs与患者的OS显著相关。为了进一步筛选构建预后模型的最优组合，随后进行了逐步多因素Cox回归分析。获得了OS相关的预后模型（PTC-mi1m4），包括了1个miRNA: hsa-miR-1305和4个mRNAs: PAPSS2, SCD, ID4, 和DCBLD2)。DFS相关的预后模型（PTC-m3）包括3个mRNA: TMEM184A, TOB1, 和FADS6。
为了评估这些模型的预测性能，分别对OS和DFS的预后模型进行了验证。
OS相关预后模型：
根据risk score的中位数，将所有的患者分为高风险和低风险组。KM分析曲线显示出低风险组的患者比高风险的患者具有更高的生存率。时间依赖的ROC曲线表明5年和10年的AUC值为0.781和0.823，C指数为0.755。这些结果表明该模型对预测乳头甲状腺癌的OS具有较强的预测能力。
为了进一步评估了该模型的独立预测能力，我们纳入了一些临床指标（年龄，性别和肿瘤分期）。单因素Cox回归分析显示，风险评分、年龄和分期与PTC患者OS相关，但多因素Cox回归分析显示，风险评分和年龄是PTC患者OS的独立预后指标。
DFS相关预后模型：
使用同样的方法对DFS模型的预测性能和独立预测能力进行了评估。结果表明该模型对预测PTC的DFS具有较强的预测能力；并且，风险评分和分期是PTC患者DFS的独立预后指标。
PTC患者的免疫景观
肿瘤浸润免疫细胞比例和水平的差异可能代表了不同个体的内在特征，因此我们研究了乳头状甲状腺癌患者的特异性免疫特征。运用CIBERSORT对样本中22个免疫细胞的比例进行评估。我们观察到静息记忆CD4 T细胞的比例最高，而中性粒细胞的比例很低，提示这两种免疫细胞可能在PTC肿瘤的发生发展中起重要作用。
随后，分别对OS和DFS预后模型中高风险和低风险组免疫细胞比例的差异进行了分析。
在OS模型中，与低危患者相比，高危患者中静息记忆CD4 T细胞，巨噬细胞M0和激活的树突状细胞的比例显著升高。另外，CD8 T细胞，滤泡辅助性T细胞和调节性T细胞（treg）比例在高危患者中比例较低。随后，对特征性mRNA和免疫细胞进行了Pearson相关分析。结果表明活化树突状细胞，巨噬细胞M0和treg不仅在高，低危组存在差异，还与基因表达水平密切相关。
对这3中免疫细胞进行了单因素和多因素Cox回归分析。结果表明活化树突状细胞与PTC患者的OS显著相关，另外，它的浸润程度在高危患者中较高。
在DFS模型中，与低危患者相比，高危患者静息树突状细胞，巨噬细胞M0，静息肥大细胞和treg的比例较高；而嗜酸性粒细胞，巨噬细胞M1，活化肥大细胞，单核细胞和静息记忆CD4 T细胞的比例较低。随后，对基因和免疫细胞浸润进行了Pearson相关分析。结果表明巨噬细胞M0，嗜酸性粒细胞，treg和静息记忆CD4 T细胞不仅与高低风险组有关，还与预后的特征基因有关。
随后，对差异的4种免疫细胞进行单因素和多因素Cox回归分析，发现巨噬细胞M0与PTC的DFS密切相关。
应用CMap分析确定治疗药物
发现新的有效药物可能改善PTC患者的预后。在我们的两个预后模型中，获得了与PTC预后相关的7个特征mRNA。靶向这些基因的药物有望在PTC的治疗中发挥巨大的潜力。除2个没有GPL96探针ID的外，其余5个mRNA: PAPSS2, TOB1, ID4, SCD, 和DCBLD2，上传至CMap web工具，以筛选逆转这5个hub基因表达的化合物。如果连接性得分为负，则表明该化合物抑制了查询基因的表达。因此，我们确定连接性得分接近-1的前三种生物活性化合物为潜在的PTC治疗药物。共获得了3种化合物，piolitazone，benserazide和SB-203580。通过文献检索，这3个生物活性化合物均可通过靶向5个特征mRNA调控PTC相关的生物过程或通路。
亮点总结
本文有以下几点创新之处。涉及了多个数据库整合分析（GEO和TCGA），而不是单一数据库。疾病的诊断和预后都进行了分析。另外，基于ceRNA网络构建了OS和DFS两个预后模型，相对比较全面。在免疫浸润分析中，还添加了免疫细胞和基因/预后之间的相关性。最后，还做了药物预测分析。
整体的分析有逻辑性，每个步骤都不难，可以复刻，这样能发到5分+也很不错！小伙伴们，如果你们也准备研究预后模型而不知道从哪入手的话，可以参考本次分享的文献哦。
另外，欢迎关注我们，我们会定期分享文献解读的视频以及写作技巧等干货。如果有感兴趣的内容，可以在评论区留言提醒小编哦。