跟着SCI学转录组分析|整合全外显子组和转录组测序揭示磨玻璃结节型肺腺癌从癌前病变到浸润癌的动态演变

肺癌早期检测可显著提升治疗成功率。随着肺癌筛查普及和高分辨率 CT 应用，肺部磨玻璃结节（GGNs）检出率大幅增加，其中包含大量原位腺癌（AIS）、微浸润腺癌（MIA）和早期浸润性腺癌（IAC）。AIS 作为 IAC 的癌前病变，是研究肺癌起始阶段的理想模型，但其分子特征及向 IAC 的演变机制仍不明确。尽管近年肺癌动态进展的分子机制成为研究热点，部分患者仍缺乏常见驱动突变。今天我们要解读学习的这篇文献，通过对35例早期LUAD患者的GGNs进行全外显子组测序（WES）和转录组测序（RNA-Seq），旨在明确突变特征、肿瘤异质性和信号通路改变，为早期LUAD的临床诊断和手术干预提供理论依据。

标题：Integrated whole-exome and bulk transcriptome sequencing delineates the dynamic evolution from preneoplasia to invasive lung adenocarcinoma featured with ground-glass nodules 链接：https://doi.org/10.1002/cam4.7383

今天我们重点学习该文中关于RNA-Seq的研究思路和分析方法。

摘要

本研究通过全外显子组测序（WES）和转录组测序（RNA-Seq）分析 35 例磨玻璃结节（GGNs）型肺腺癌（LUAD）患者样本，揭示了从原位腺癌（AIS）到微浸润腺癌（MIA）再到浸润性腺癌（IAC）的基因组动态演变。发现 EGFR（49%）和 RBM10（14%）突变与病理进展显著相关，鉴定出 466 个渐进上调的差异表达基因（DEGs），涉及细胞迁移、侵袭通路（如 GPR143、CCR9）。EGFR/RBM10 共突变组显示胆固醇代谢通路激活，免疫微环境分析显示 NK/CD8⁺ T 细胞减少、Treg/B 细胞增加。研究构建了 CNV-mRNA 调控网络（如 FAM83A、MAL2），并通过 CMap 分析筛选出潜在药物 W13，为早期 LUAD 精准治疗提供依据。

今天我们主要学习文中关于转录组测序（RNA-Seq）的研究思路和分析方法。

方法

本研究RNA提取使用TRIzol试剂，经Agilent 2100 Bioanalyzer评估质量后，用VAHTS Universal V6 RNA - Seq Library Prep Kit构建文库并在Illumina NovaSeq 6000平台测序。 RNA - Seq分析，用fastp处理原始数据，HISAT2比对，HTSeq count确定基因表达，用DESeq2确定差异表达基因（DEGs），GO、KEGG和Reactome进行富集分析（GSEA），STRING数据库和Cytoscape软件构建蛋白质 - 蛋白质相互作用（PPI）网络。用R包limma获得差异表达基因，进行CMap分析；用R包IOBR进行免疫浸润分析；用SPSS软件进行统计分析。

注：如果你不会写代码亦或者觉得安装软件配置环境太麻烦，也可以选择使用在线网站完成分析，比如Galaxy云平台（网址：usegalaxy.cn）。

结果

35例患者的临床信息

共收集35例肺GGNs患者的70个样本，包括12例AIS、13例MIA、10例IAC和35例癌旁组织。经三位病理学家确诊为LUAD。队列中有18例纯磨玻璃结节（pGGN）患者和17例混合磨玻璃结节（mGGN）患者，27名女性，8名男性，中位年龄55岁。病理亚型与GGO类型、结节大小显著相关，与性别、吸烟史和结节位置无关。

基因富集分析

对SMGs进行富集分析，GO富集分析显示AIS的生物学过程与MIA和IAC完全不同；KEGG富集分析表明AIS的SMGs主要富集在间隙连接、磷脂酶D信号通路等，而MIA和IAC主要富集在粘着斑和MAPK信号通路；Reactome富集分析显示AIS与MIA、IAC的SMGs富集信号通路不同（图3A - C）。

图3与SMGs相关的富集分析

图3 与SMGs相关的富集分析注：(A) AIS中SMGs的GO、KEGG和Reactome富集结果。(B) MIA中SMGs的GO、KEGG和Reactome富集结果。(C) IAC中SMGs的GO、KEGG和Reactome富集结果

RNA测序和基因表达谱分析

肿瘤组织与正常组织的转录组差异：通过层次聚类分析发现肿瘤组织和正常组织共有1503个DEGs ，其中864个上调，639个下调（图4A）。不同病理亚型与正常组织的DEGs存在差异，AIS中显著上调的基因有ARHGAP40、HABP2等；MIA中有HABP2、ABCA4等；IAC中有GLB1L3、HABP2等（图4B - D）。功能富集分析显示，AIS主要富集在细胞粘附、免疫反应等；MIA主要富集在细胞粘附、细胞 - 细胞信号传导等；IAC主要富集在细胞粘附、血管生成等（图4B - D）。

图4 肿瘤组织和正常组织之间的转录组景观

注：(A) log2倍数变化>1且p < 0.05的差异表达基因（DEGs）热图。对肿瘤组织和正常组织中鉴定出的DEGs进行层次聚类分析。表达水平用颜色表示，蓝色表示低表达，粉橙色表示高表达。(B) log2倍数变化>1且p < 0.05的DEGs火山图，以及AIS中DEGs的富集结果。(C) log2倍数变化>1且p < 0.05的DEGs火山图，以及MIA中DEGs的富集结果。(D) log2倍数变化>1且p < 0.05的DEGs火山图，以及IAC中DEGs的富集结果

三种病理亚型之间的转录组差异：用R包edgeR和DEseq鉴定出AIS中有1124个DEGs ，MIA中有1310个，IAC中有1895个（|log2FoldChange| > 1且FDR < 0.01），其中597个为三种病理亚型共有（图5A）。对这些共有DEGs进行Reactome富集分析，发现其富集在“GPCR配体结合”“细胞连接组织”等致癌通路（图5B）。构建共有DEGs的PPI网络，鉴定出ECM - 受体相互作用、细胞因子 - 细胞因子相互作用等关键信号通路（图5C）。用MCODE插件分析得到两个高分聚类，其中都包含参与细胞因子和细胞因子受体相互作用的蛋白质（图5D、5E）。

图5 病理亚型DEGs一致性比较注：(A) 病理亚型之间DEGs的维恩图。(B) 病理亚型之间的Reactome富集分析。(C) 共有DEGs的蛋白质 - 蛋白质相互作用（PPI）网络。(D, E) MCODE插件聚类分析结果。节点表示蛋白质，边表示相互作用

为深入了解LUAD进展的动态机制，进行STEM分析，发现有539个基因在AIS、MIA和IAC中逐渐增加（图6A）。对这些基因绘制热图，显示其在肿瘤样本中的表达模式（图6B）。排除肿瘤组织中下调的基因后，剩余466个基因可能是LUAD动态演变的关键基因。Reactome和Wikipathway富集分析表明，这些关键基因显著富集在细胞迁移、侵袭等相关通路（图6C）。最终确定63个在AIS、MIA和IAC中逐渐增加的基因，强调其在从AIS到IAC动态演变中的重要性（图6D）。

图6 恶性转化中上调的DEGs及相关生物学过程注：(A) STEM分析的DEGs趋势图。(B) 不同组织学亚型肿瘤组织和正常组织之间上调DEGs的热图。(C) 上调DEGs的Wikipathway和Reactome富集分析。(D) 逐渐上升的上调DEGs热图

磨玻璃结节进展过程中CNV - mRNA相关性分析

与正常组织相比，三种肿瘤亚型均观察到大规模CNVs ，ANKRD46、PABPC1等基因显著扩增且从AIS到IAC逐渐增加（图7A）。整合DEGs和CNVs分析发现，466个在病理亚型中逐渐上调的基因中，238个基因的拷贝数与相应RNA水平显著相关，且多数CNV - mRNA对与GPCR配体结合、PI3K - Akt - mTOR信号通路等正相关（图7B）。FAM83A、MAL2等基因的mRNA表达水平与CNVs显著相关（图7C、7D）。Western Blot实验证实，随着病理进展，这些显著上调基因的表达增加（图7E）。

图7 mRNA表达与CNV状态的相关性注：(A) 前30个具有CNVs基因的概述。(B) mRNA表达与CNVs之间的相关性。(C) mRNA与CNVs之间相关性得分最高的前十个基因。(D) mRNA与CNVs之间相关性的散点图。(E) 12例不同病理阶段患者肿瘤组织中前4

3.6 EGFR/RBM10共突变与EGFR突变患者的比较

我们观察到RBM10突变常与EGFR突变同时出现。为进一步了解肺癌突变的机制，我们将17例EGFR突变患者分为两组：EGFR/RBM10共突变组和仅EGFR突变组。比较两组的mRNA表达谱，发现EGFR/RBM10共突变组中DPP4、BCAT1、FNDC9、VCX2、PGLYRP3和NPY的mRNA表达水平显著上调（图8A）。基因集富集分析（GSEA）表明，共突变组在血红素代谢（富集分数[NES]= -1.485，p = 0.0235）和胆固醇稳态（NES = -1.5362，p = 0.0267）信号通路中显著富集（图8B）。我们使用癌症基因组图谱（TCGA）数据库中的临床数据比较了两组的总生存率，结果显示共突变组的生存率略低于仅EGFR突变组，但差异无统计学意义（图8C）。此外，我们通过连通图（CMap）分析探索了可能靶向EGFR/RBM10共突变肺癌的潜在药物。分别分析了我们自己患者队列和TCGA-LUAD数据库中的前500个差异表达基因，并使用极端和（XSum）方法将EGFR/RBM10共突变相关的基因特征与CMap基因特征进行匹配。每种化合物的Xsum分数如图8D所示。W13是一种有效的MsbA抑制剂，在两个数据集的结果中均同时出现且分数明显较低，这表明W13可能对携带EGFR/RBM10共突变的肺癌细胞具有抑制作用。

图8 mRNA表达与EGFR/RBM10共突变的相关性注：(A) EGFR/RBM10共突变与仅EGFR突变之间差异表达基因（DEGs）的火山图。(B) DEGs的GSEA富集分析。(C) 通过TCGA数据库比较EGFR/RBM10共突变与仅EGFR突变的总生存率。(D) 使用极端和（XSum）方法的CMap分析结果。右图展示了逆转效力排名前十的化合物

3.7 不同病理亚型的免疫细胞微环境

为进一步研究免疫微环境与LUAD侵袭进展之间的关系，我们使用六种计算免疫微环境评分的方法，对聚类为正常、AIS、MIA和IAC组的免疫细胞表达谱进行了分析。由于解析算法的异质性，我们总结了六种方法中最一致的趋势。如图9所示，在CIBERSORT和MCPcounter的结果中，随着肺结节的进展，自然杀伤（NK）细胞的浸润显著减少（图9B）。同样，TIMER和CIBERSORT分析表明，恶性结节组织中的CD8⁺ T细胞减少（图9C）。此外，在CIBERSORT和xCell的分析中，随着肺结节的进展，调节性T细胞（Tregs）逐渐增加（图9D）。有趣的是，根据quanTIseq和xCel估计的分数，肿瘤组织中的B细胞数量增加（图9E）。

图9 恶性肺结节不同病理阶段的免疫细胞特征

注：(A) 热图显示基于恶性肺结节不同病理阶段RNA-Seq数据，用六种方法计算的免疫浸润分数。(B) 各病理亚型NK细胞免疫分数的散点图。(C) 各病理亚型CD8⁺ T细胞免疫分数的散点图。(D) 各病理亚型Tregs免疫分数的散点图。(E) 各病理亚型B细胞免疫分数的散点图