Scoring model based on the signature of non-m6A-related neoantigen-coding lncRNAs assists in immune microenvironment analysis and TCR-neoantigen pair selection in gliomas
基于非M6A相关新抗原编码lncRNAs特征的评分模型有助于胶质瘤的免疫微环境分析和TCR-新抗原配对选择
发表期刊：J Transl Med
发表日期：2022 Oct 29
影响因子：8.440
DOI: 10.1186/s12967-022-03713-z
一、研究背景
胶质瘤是最普遍的脑肿瘤，5年生存率不超过3%。导致胶质瘤难治的一个因素是相对较低的免疫反应，即"免疫冷"或免疫抑制的反应。在大多数情况下，癌症产生于基因突变或DNA损伤的积累。新抗原来源于这些非同义的基因突变，包括单核苷酸变异（SNV）、染色体缺失和插入、基因融合和替代剪接。由于所呈现的新抗原可能会引起T细胞介导的抗肿瘤免疫，专门针对肿瘤细胞，因此它们被认为是有前途的免疫治疗目标。
lncRNAs的翻译可以被许多因素调节，包括小的或短的开放阅读框架（smORF或sORF）、eIF4E和N6-甲基腺苷（m6A）修饰。eIF4E是一种翻译起始因子，在磷酸化后与RNA的5′帽弱结合，从而诱导抑制mRNA的翻译，促进lncRNA和核糖体之间的相互作用。此外，m6A修饰已被证明可以影响mRNA的翻译。研究表明，m6A修饰的位点可作为环状RNA的翻译起始点。
二、材料与方法
1、数据来源
1）mRNA测序数据从TCGA数据库下载：169个多形性胶质瘤（GBM）样本和529个低级别胶质瘤（LGG）样本
2）验证数据从中国胶质瘤基因组图谱（CGGA）数据库下载："mRNAseq_325 "和 "mRNAseq_693"，在此分别标为CCGA325和CCGA693；在CGGA693数据集中包括249个GBM样本和444个LGG样本
3）单细胞测序（scSeq）数据来自GEO数据库：GSE84465，包含3589个细胞
4）还包括来自CGGA数据库的另一个数据集，包含6148个细胞
5）GSE129671被用于SCENIC分析
6）细胞的注释与一项同样使用GSE84465数据集的研究中描述的一样，经过轻微的修改
7）TCR测序数据从GEO数据库获得，两个数据集GSE79338和GSE188620被用于分析
8）LN229、U118MG、A172、U251MG和Jurkat细胞
9）患者和组织：WHO II级（n = 3）、III级（n = 3）和IV级（n = 3）的石蜡包埋胶质瘤组织；正常脑组织（n = 3）来自颅脑外伤患者；WHO II级（n = 4）和IV级（n = 4）的冷冻胶质瘤组织
2、分析流程
三、实验结果
01 - 非m6A修饰在高等级胶质瘤中被激活，NAS预测胶质瘤患者的预后
首先，收集TCGA数据集中非M6A相关调节因子的表达水平。非m6A修饰根据其水平分为以下三类：1）高到中等水平的修饰，有数百到数千个修饰位点（Ψ和m5C）；2）超低水平的修饰，有少数修饰位点（m1A）；3）未知水平的修饰，需要进一步确认（如N4-乙酰胞苷[ac4C]，2'-O-甲基化[Nm]，和7-甲基鸟苷[m7G]）。
接下来，前两类的修饰，包括Ψ、m5C和m1A，被选为要研究的主要非m6A修饰。通过相关性分析并在R中用“igraph”包构建相关网络选择与神经胶质瘤患者总生存时间显著相关的lncRNA，根据TransLnc，确定了潜在的肽编码lncRNAs，其中的肽可以被MHC I呈递，因为MHC I是内源性抗原的主要贡献者，包括肿瘤衍生的抗原。总共发现了13个lncRNAs，其中5个的HR大于1，8个的HR小于1。
建立新抗原活化评分（NAS）模型的过程见图1A。根据非M6A相关的新抗原编码lncRNAs的表达进行共识聚类。来自TCGA的胶质瘤样本被分为两个聚类，以获得每个聚类中最高的相关性，主成分分析（PCA）显示了聚类1和聚类2之间的分布差异。作者确定了这两个簇中选定的非m6A调节因子的表达，结果表明，大多数writers和readers在簇2中被上调（图1B）。对于erasers，FTO和TET2在群集1中明显升高，而ALKBH1和ALKBH3的水平在两个群集之间没有明显差异（图1B）。这表明与群集1相比，群集2中的非m6A修饰状态增强。在群集2中，HR大于1的lncRNAs表达水平明显高于群集1，而HR小于1的lncRNAs表达水平明显较低（图1C）。然后用SVM算法学习群集1和群集2的基因表达特征，并在验证数据集中重现聚类模型，在CGGA325数据集中观察到类似的结果。此外，集群模型有效地预测了所有胶质瘤和LGG的预后。在该模型中，群集2表现出比群集1明显更差的预后，而在GBM方面，群集之间没有明显差异。然而，在ROC分析中，集群模型的预后效果并不理想，非M6A相关的新抗原集群模型的AUC值为0.72（图1G），甚至低于年龄（AUC=0.82）和grade（AUC=0.82）。
作者假设聚类模型的准确性受到其对TCGA样本二元分类的限制。为了提高非M6A相关的新抗原聚类模型的准确性，在聚类1和聚类2之间的DEGs基础上建立了一个NAS模型。该公式包括基因表达水平及其由PCA产生的权重，所以它提供了一个量化的参数。根据对聚类模型的分析，聚类2比聚类1显示出更差的预后，与聚类1相比，所有HR>1的具有明显预后效果的DEGs在聚类2中都有所升高。HR>1也表明一个基因的高表达与预后恶化有关。另外，所有HR<1的具有明显预后效果的DEG在聚类2中都下调。因此，该公式基本上反映了RNA-seq样本的基因表达模式与聚类2之间的相似程度。与集群模型相比，NAS模型更精确地区分了不同层次的胶质瘤亚型，并与选定的调节器或lncRNAs的表达高度相关（图1B，C）。在CGGA325和CGGA693数据集中计算NAS时，也得到了类似的结果。此外，NAS模型在预测TCGA（图1D）、CGGA325（图1E）和CGGA694（图1F）数据集中GBM患者的预后时表现良好，因为高NAS组呈现较短的平均生存时间。在这三个数据集中的所有胶质瘤和LGG组都有类似的结果。非M6A相关NAS模型的AUC值为0.88，远远高于集群模型，年龄和等级的值也是如此（图1G）。
为了确定非m6A相关NAS模型和m6A相关NAS模型之间的差异，以类似的方式构建m6A相关NAS模型。m6A聚类模型的AUC值为0.59，而NAS模型为0.87（图1G）。然后将所有3个数据集（TCGA、CGGA325和CGGA693）的数据合并，再次计算所有模型的AUC值。结果表明，非m6A相关NAS模型的AUC值为0.76，比m6A相关NAS模型（AUC=0.66）大。这表明在所有三个数据集中，非m6A相关模型比m6A相关模型有更好的预后准确性。因此，非m6A相关的NAS模型被选为进一步研究的对象。
02 - 高NAS与胶质瘤的侵袭性亚型相关联
为了确定NAS与胶质瘤亚型之间的关系，采用NMF聚类法，将样本按照Verhaak分类分为三个聚类。结果显示，在三个亚型中，间质（MES）亚型的平均NAS最高，而俯卧神经（PN）亚型在所有三个数据集中的平均NAS最低（图2A-C）。MES亚型是最具侵略性的亚型，与不良的生存结果有关，而PN亚型显示出最低的侵略性水平。因此，NAS与胶质瘤的侵袭性呈正相关，正如Verhaak分类法使用大量测序数据所预测的那样。
GSE84465数据集中胶质瘤细胞的scSeq数据被聚类和注释（图S6A），每个聚类中的标记基因被确定（图S6C）。使用同样的程序在CGGA scSeq数据集上再现了结果（图S6B，D）。此外，还使用SVM将细胞分为两个聚类。基于聚类模型，使用R中的 "DEsingle "包探索DEGs，并计算NAS。聚类2与t分布随机邻居嵌入（t-SNE）还原图中NAS相对较高的细胞广泛重叠（图2D，E），Mann-Whitney检验证实了这些结果（图2F）。此外，还发现在GSE84465和CGGA数据集中，低肿瘤细胞的NAS明显低于高肿瘤和炎症相关的胶质瘤细胞（图2G）。然后将GSE84465数据集中的细胞分为四个先前建立的单细胞水平的分子亚型，然后通过分析剪接和未剪接RNA的丰度来计算GSE84465数据集中胶质瘤细胞的RNA速度，以便对胶质瘤细胞的进化过程进行分析。第2组中的大多数细胞处于其革命过程的末端（图2H）。同时，MES和少突胶质细胞祖细胞样（OPC）亚型位于革命途径的末端（图2I），表明大多数胶质瘤细胞随着时间的推移变得更具侵略性。NAS随着胶质瘤细胞的发展而增加（图2J），表明NAS反映了胶质瘤细胞革命的阶段。这些结果也通过伪时间分析得到证实（图S6I）。此外，细胞轨迹分析显示，在GSE84465（图S6E，G）和CGGA数据集（图S6F，H）中，OPCs和高新生胶质瘤细胞处于细胞轨迹的顶端，而低新生细胞处于细胞轨迹的上游。肿瘤性高的胶质瘤细胞比肿瘤性低的细胞有更高的NAS，这一事实表明，当这种发展伴随着更高的NAS时，胶质瘤细胞会发展成更具侵略性。这些数据共同表明，NAS和胶质瘤的侵略性之间存在正相关关系。
03 - 较高的NAS胶质瘤与较高的免疫浸润水平相关联
为了确定参与NAS水平差异的生物功能，作者对低或高NAS的TCGA样本的DEGs进行富集分析。GSVA富集分析显示，在GO和KEGG途径中，高NAS的样本在T细胞介导的免疫、NK细胞介导的细胞毒性以及抗原处理和表达方面显示出较高的富集分数（图3A，B）。GO富集分析在T细胞相关功能、抗原处理和呈现方面有非常相似的结果（图3C），以及在KEGG数据库中富集的重要途径，包括免疫相关途径，如NK细胞介导的细胞毒性、趋化因子信号途径和抗原处理和呈现（图3D）。scSeq数据的GSEA富集分析显示，高NAS组的抗肿瘤免疫因子（图3E）。然而，在这个数据集中，细胞-细胞粘附力被下调。对CGGA数据集中scSeq数据的分析也显示，高NAS组的T细胞相关功能得到促进（图3F）。总之，NAS和非M6A相关的新抗原编码lncRNAs与T细胞相关的免疫和抗原处理和表达有关。
为了确定NAS和免疫浸润之间的关系，采用了R中的 "ESTIMATE "包来评估TCGA数据集中样本的免疫浸润。NAS与免疫评分呈正相关，与肿瘤纯度呈负相关（图4A），表明较高的NAS意味着较高的免疫浸润水平。对CGGA325和CGGA693数据集的样本分析显示了类似的趋势。之后，用CIBERSORT详细分析了浸润的免疫细胞的变化，结果表明，虽然CD8 + T细胞和T辅助细胞的比例升高，但调节性T（Treg）细胞的比例增加，而激活的NK细胞的比例下降（图4B，C）。鉴于较高的NAS提示较差的生存结果和更具侵略性的肿瘤，分析了不抑制胶质瘤细胞的较高免疫浸润的可能机制。当胶质瘤等级上升时，PD-L1的表达升高，表明免疫抑制程度更高（图4D）。还发现四个lncRNAs（即AC060766.4、AC0738962、LEF-AS1和LINC00893）的表达与PD-L1呈正相关（图4E）。因此，很明显，高NAS组的免疫浸润高于低NAS组，PD-L1表达的增加可能抑制了T细胞介导的免疫的激活。
用IHC检测了NAS相关基因和PD-L1的表达。在NAS中，PC1+PC2是影响最终结果的主要指标。检查了NAS与TMSB10、VIM和PD-L1的表达之间的相关性（图5A）。结果显示，NAS与TMSB10、VIM和PD-L1的表达明显相关。然后用IHC检测TMSB10、VIM和PD-L1的表达。结果发现，更高级别的胶质瘤表现出更高的TMSB10、VIM和PD-L1水平（图5B，C）。这些结果表明，与NAS和PD-L1正相关的基因在较高等级的胶质瘤中表达较多，这也表明较高的NAS可能与较多的PD-L1表达有关。
04 - 高NAS组的T细胞阳性调节因子的异常表达可能导致T细胞功能失调
作者分析了T细胞的功能，它们是抗肿瘤免疫的直接执行者。与Ca2+通量（AHNAK和CALML3）、DNA修复（ZNF830）和自噬（HOMER1）相关的基因在TCGA数据集的高NAS组中被明显下调（图6A）。对CGGA325和CGGA693数据集的分析也取得了类似的结果，其他调节器在高NAS组中被上调。
此外，TCGA样本的SNV数据显示，这33个基因的TMB在高NAS组高于低NAS组，没有统计学意义，AHNAK对高和低NAS组贡献了大部分的突变（图6B，C）。同样，第2组的样本比第1组的样本显示出更多与这33个基因有关的突变负担。高等级胶质瘤的T细胞可能是对Ca2+通量、DNA修复、自噬和醛酮代谢的干扰。在这些基因的功能中，由AHNAK介导的Ca2+通量被强调，因为该基因在研究的33个基因中显示了最多的突变。然后，作者试图通过体外功能试验来验证Ca2+在T细胞功能中的作用。细胞内Ca2+螯合剂BAPTA-AM以0、10、20和40μM的浓度作用于Jurkat细胞48小时，然后检测细胞内钙，在BAPTA-AM组中钙明显下降（图6D）。Jurkat细胞的增殖也受到BAPTA-AM处理的明显抑制（图6E）。对于共培养试验，结果表明，IFN-γ的分泌在BAPTA-AM组也受到抑制（图6F）。BAPTA-AM组中剩余的LN229细胞的增加也表明激活的Jurkat细胞的功能受到抑制（图6G，H）。综上所述，在高NAS样本中，异常的Ca2+通量可能在T细胞介导的胶质瘤生长抑制失败中发挥了重要作用。
05 - 高NAS胶质瘤与参与干性的转录因子有关
为了确定TCGA、CGGA325和CGGA693数据集的上游调节转录因子网络，采用了X2K。将低NAS组和高NAS组的DEGs导入X2K，然后得到TCGA（图7A）、CGGA325（图7B）和CGGA693（图7C）数据集的前20个上游调节转录因子。注意到前20个转录因子中的一些在所有三个数据集中是共同的，如SUZ12、REST、EZH2、SMAD4和AR。SUZ12、REST和EZH2基因有助于癌细胞的干性，它们在高NAS组中表现出更高的表达水平（图7D）。这些结果表明，在RNA-seq数据中，较高的NAS与促进干性转录因子的表达有关。
在scSeq数据方面，应用pySCENIC构建转录因子的调控网络。计算了GSE84465数据集中的差异性激活转录因子（图7E）。前5个被激活的转录因子被证明，干性相关的基因，在较高的NAS细胞中被激活，而在较低的NAS细胞中，前5个被激活的转录因子中没有观察到此类基因。然而，在CGGA scSeq数据中，得到了相反的结果。作者在数据集GSE129671中进行了另一次验证。在五个最活跃的转录因子中，高NAS组的MYC结果为促进细胞干性。
为了研究胶质瘤标本中干性相关和NAS相关基因的表达，作者应用lasso来确定NAS的主要促成因素。通过合并TCGA、CGGA325和CGGA693数据集，确定了五个基因，简化的NAS可计算为：NAS=0.22864885*（COL5A2表达水平）+0.1083532*（PVP1表达水平）+0.07381116*（CHI3L2表达水平）+0.03597545*（SERPINE1表达水平）+0.02452755*（SOCS3表达水平）。得出的NAS与所有三个数据集的NAS高度相关（图8A）。然后在4个II级和4个IV级胶质瘤中用qRT-PCR检测这五个基因。结果显示，所有五个基因在IV级胶质瘤中都有所升高，COL5A2和PVR1显示出统计学意义（图8B），表明IV级胶质瘤的NAS更高。而上述分析所确定的五个干性相关基因也被检测到，包括EZH2、SUZ12、REST、SOX10和MYC。结果显示，只有EZH2在IV级胶质瘤中明显升高（图8C）。这五个基因在病人衍生的胶质母细胞瘤干细胞和分化细胞中也被检测到。它表明EZH2和SOX10在干样细胞中明显升高（图8D）。这些结果表明，大多数样本和数据集的高NAS组存在较高的干性。而EZH2是检测中最明显升高的干性相关基因。
作者应用R语言中的"celltalker "包来分析细胞间的交流，以确定高NAS组中T细胞功能被抑制的潜在机制。在低NAS组中，低肿瘤细胞在大多数重要的相互作用中是活跃的，通过ADAM12和ITGA9/SDC4途径与树突状细胞和T细胞相互作用，两者在总的相互作用图中没有点（图9A）。相反，在高NAS组，大多数相互作用发生在T细胞和OPC或炎症相关胶质瘤细胞之间（图9B，C）。ADAM12和ITGA9/SDC4途径是这些细胞之间的两个主要途径。
此外，还分析了T细胞与四种类型的胶质瘤细胞之间的重要相互作用，结果显示，高NAS组的T细胞与低肿瘤性、炎症相关胶质瘤细胞和OPC之间的相互作用比低NAS组多。此外，在CGGA数据集中，T细胞没有参与低NAS组的大多数重要的相互作用，而在高NAS组，它们与其他免疫细胞相互作用。此外，关于显著的相互作用，在高NAS组中，只有一些T细胞和炎症相关胶质瘤细胞之间的相互作用是显著的，表明炎症相关胶质瘤细胞参与了T细胞介导的免疫。
值得注意的是，T细胞与NAS相对高于上述的高肿瘤性胶质瘤细胞之间没有明显的相互作用。因此推断较少的T细胞结合和细胞-细胞粘附可能是基本机制。结果显示，与其他三种胶质瘤细胞类型相比，IFNGR1和JAK1的表达明显较高，但JAK2的表达没有明显差异（图9D-F）。此外，IFNGR1被发现与NAS有明显的负相关（图8G），而JAK1和JAK2则没有。还发现在CGGA数据集中，低肿瘤细胞显示出较高的IFNGR1表达，但JAK1在低肿瘤细胞中的表达略有下降，IFNGR1和NAS之间存在明显的负相关。这表明T细胞和胶质瘤细胞之间的相互作用增加，而高肿瘤细胞可能通过下调IFNGR1以减少T细胞的结合来逃避这种相互作用。
06 - 预测具有两种模式的TCR会与非m6A相关lncRNA的新抗原结合
为了确定一些基于非m6A相关新抗原模型的胶质瘤治疗的可能方法，筛选了已发表的TCR测序数据集，以探索可能与13个选定的非m6A相关lncRNAs编码的肽结合的TCR克隆型。来自GSE79338的LGG和GBM样本，被用来识别GBM和LGG的独特TCR克隆型。在寻找MHC限制性肽抗原的过程中，使用GLIPH2算法将TCR克隆型聚类为CDR3模式，在GBM和LGG样本中发现了许多在正常组织中找不到的独特的TCR CDR3模式。从这些TCR CDR3模式中，提取了52种LGG和GBM样本中的常见模式（图10A）。通过比较GBM和LGG中的模式的频率，确定了两组之间差异最大的10个模式（图10A）。其中五个氨基酸模式在GBM模式中被上调，从而暗示这些模式可能广泛存在于胶质瘤患者中，并可能与胶质瘤的进展有关。接下来，应用DLpTCR算法来评估13个选定的lncRNA编码的MHC I呈递肽与上述GBM或LGG样本中5个选定模式的TCR克隆型之间的识别和结合可能性。MNEQ和HDEQ是具有最大结合概率的模式（图10B）。
为了探索5种选定的模式，检查了GSE188620数据集中四个胶质瘤患者在肿瘤细胞裂解液接种前后的胶质瘤组织的scTCR数据。在接种疫苗后的样本中发现了独特的克隆型（。然而，对于5个选定的模式，虽然HDEQ和RNKQ表现出极低的表达量，但其他三个模式的表达量却明显较高。此外，接种疫苗后，%GSTDTQYF和MNEQ的总表达量明显升高，同时FGEQ也有轻微的、统计学上不明显的下降（图10C）。患者1、3和4中%GSTDTQYF的表达量升高，但只有患者4中发现有统计学意义（图10D）。患者1和2显示出MNEQ表达的升高（图10D）。还确定了疫苗接种后扩大的克隆型（图10E），这可能包括新抗原反应性TCR克隆型。确定了扩大的TCR克隆型和肽的结合可能性，有44个TCR-肽对的结合概率超过0.98；这表明潜在的新抗原反应性TCR克隆型和选定的肽之间有很高的结合可能性。总之，作者确定13个选定的lncRNAs可能编码的肽与潜在的新抗原反应性TCR克隆型有很高的结合概率，这些克隆型在胶质瘤组织中广泛存在。
四、结论
作者建立了一个基于非M6A相关的新抗原编码lncRNAs和NAS的预后模型，发现它们与T细胞免疫、抗原处理和表达以及免疫浸润呈正相关关系。还筛选了由所选lncRNAs翻译的普遍靶向的新抗原的可能TCR克隆型。这项研究详细说明了non-m6A修饰在lncRNA编码的多肽中的重要作用，它提供了TCR克隆型，可用于未来潜在的CAR-T研究和治疗。
