主讲老师第一课：文献套路解析
在正式的文章复现之前，我们先来看一下整个文章的分析思路。
第一板块，带大家整体的浏览一下全文，与三十六策中基础科研的五恒量三变量相互呼应。某老师总结出生信文中的“挑”、“圈”、“联”、“靠”四字逻辑套路，这四个呢，分别是指差异表达、功能聚类、交互网络，以及临床意义。掌握其中奥妙的临床小伙伴们，从此就可以对干、湿套路的研究驾轻就熟了~~~文献就是在群里发布的——范例文献 1 和范例文献 2第一篇范例文献首先分享的第一篇内容是基于基因突变与免疫治疗分析的文献，于 2020 年 8 月发表于 Frontiers in Genetics 杂志上。
题 目 为 ATM Mutations Benefifit Bladder Cancer Patients Treated With ImmuneCheckpoint Inhibitors by Acting on the Tumor Immune Microenvironment. 杂志的影响因子为 3.258 分先来看下文章的背景内容：DNA 修复基因的缺陷在肿瘤的发生，发展，以及治疗反应和预后中起着重要作用。在膀胱癌中，常常发现 DNA 修复相关基因的突变，并且其突变与对顺铂类新辅助化疗的良好应答反应密切相关。同时，有研究报道，在铂类药物治疗的晚期膀胱癌中，DNA 损伤反应基因的突变与患者相对较好的预后相关。另外，DNA 修复机制的缺陷与对 ICI（Immune Checkpoint Inhibitors）的反应增加有关。因此，基于前期的这些结果表明，DNA 损伤与修复途径可能在 ICI治疗中发挥重要作用再来看看今天讲解的主变量分子，ATM 基因。ATM，Ataxia telangiectasia mutated，翻 译过来是“共济失调毛细血管扩张突变基因”，ATM 最早发现于毛细血管扩张性共济失调症患者。ATM 基因是 DNA 修复系统的核心成员，可以在识别 DNA 双链断裂后激活 ATM 诱导的同源重组修复途径。近年来，抑制 DDR 途径已成为许多癌症中的常见治疗策略，并且已知 ATM 信号传导途径在乳腺癌，生殖细胞癌和其他癌症的发展中起重要作用。同时，有研究证实，ATM、RB1 和 FANCC 基因的突变可以预测新顺铂化疗在膀胱癌中的疗效。肿瘤免疫治疗，其中以肿瘤免疫检查点抑制剂（Checkpoint inhibitors）的临床研究最为成熟和充分，应用最为广泛。免疫检查点药物主要分为两大类，一类是以 PD-1 为代表的抑制剂，另一类是激活剂，激活剂目前还处于临床研究阶段下面，我们一起来看一下整篇文章的思路
1.数据来源
根据实验方法描述，我们可以看到，作者一共使用了两个不同来源的数据集，第一个是一个免疫治疗的数据集，其在 2019 年发表于 Nat. Genet.杂志上，题目是Tumor mutational load predicts survival after immunotherapy across multiplecancer types. 在该数据集中，包含了 210 例膀胱癌患者的临床信息和基因突变信息，而且，值得一提的是，所有患者均经过了抗 PD-1/PD-L1 的免疫治疗。对于这个数据集，大家可以标记一下，后期在进行自己的研究项目时可以参考使用。根据ATM 基因的突变情况，作者将所有患者分成了 ATM 突变和非突变组，分别命名为ATM-MT 和 ATM-WT。而另一个数据集则是大家所熟悉的 TCGA 数据库中膀胱癌患者，在此，作者通过 TCGAbiolinks 包分别下载了 TCGA-BLCA 患者的突变数据（somatic mutation）和生存相关数据。这里注意一点，在下载 TCGA 突变数据时，需要下载体细胞突变数据来进行后续分析，同样的，对 TCGA-BLCA 患者根据 ATM基因的突变状态进行分组。当然，如果看过文献的同学可以发现，其中还存在另外一个数据集，即来自于 GDSC 数据库的相关数据，这个数据集是由一系列肿瘤细胞表达和抗癌药物敏感性数据所组成，了解了整体的数据来源和简单分组后，接下来进入正餐部分
二、ATM 突变与预后相关性
生信分析，挑圈联靠，四大步骤。由于作者直接指出了“主变量分子 ATM 基因”，因此，直接省略了“挑”的过程。
首先，作者从第四步“靠”来引入分析结果显示，在免疫治疗的膀胱癌 ICI-Treated BLCA 数据集中，经过免疫治疗后，ATM 突变组患者的预后明显优于未突变组，P = 0.007。在图 B-D 中，作者分析了ATM 基因突变对 TCGA-BLCA 患者的预后影响，结果提示，ATM 基因突变对 BLCA患者的总生存率，无进展生存率，以及无病生存率均无显著影响。因此，作者提出，ATM 基因突变可能会影响 BLCA 患者对抗 PD-1/PD-L1 免疫治疗的反应性，进而影响患者的预后。提出预后相关性现象后，作者接着分析了 ATM 基因突变对 BLCA患者整体突变水平的影响。
三、ATM 基因突变全景图
在这其中，一共包含两方面的内容，单核苷酸突变（SNP）和拷贝数变异（CNV）。
与普通的突变瀑布图不一样，作者分别展示了 ATM-MT 和 ATM-WT 两组患者中最常见突变基因的突变水平，同时，图中还展示了 ATM 突变与肿瘤突变负荷 TMB，微卫星不稳定性评分 MSI，以及不同的临床病理特征之间的相关性。基因突变可以直接的影响氨基酸的翻译过程，而重要位置氨基酸位点的改变往往可以影响蛋白的表达和功能。随后，根据 SNP 的改变位点，结合 ATM 蛋白的结构域，去除突变但氨基酸未改变的位点，将其余突变后氨基酸的改变情况，使用棒棒糖图（Lollipopplot）展示在整个蛋白结构示意图中在基因组水平，除了 SNP 水平的改变，拷贝数 CNV 的变化同样是一个十分重要的改变因素。在图 D 中，通过 GISTIC 2.0 分析，作者整体分析了 TCGA-BLCA 数据集中患者拷贝数变异水平，随后分别分析了 ATM-MT 和 ATM-WT 两组患者中 CNV 的变化水平，结果显示，无论 ATM 的状态如何，ATM-MT 组和 ATM-WT 组都具有较大的拷贝数差异，与 ATM-WT 组相比，我们发现 ATM-MT 中有更多 CNV 变异。
ATM-MT 患者的显着扩增区域包括 1q21.3、1q23.3、8q22.3 和 11q13.3，其中11q13.3 跨越 CCND1（cyclin D1）基因，而 ATM-MT 患者的显著缺失包括 9p21.3，包括抑癌基因 CDKN2A 和 CDKN2B
四、突变与不同生物学特征
接着，作者从不同的角度来阐述 ATM 基因突变对 BLCA 患者不同生物学特征的影响。
1.免疫学特征
考虑到 ATM-MT 患者能从免疫治疗过程中获益，因此作者首先分析了 ATM 基因突变对不同免疫学特征的影响，在这分析过程中，又包括免疫基因表达和免疫细胞浸润水平。
在免疫基因水平，作者分别根据功能来源和细胞来源，挑选了一系列免疫相关基因，并使用热图展示了这些基因在 ATM-MT 组和 ATM-WT 组之间的差异水平。结果显示，多种免疫相关基因的表达在两组间存在显著差异，同时，肿瘤突变负荷 TMB值在免疫治疗数据集，TCGA-BLCA 数据集，以及 GDSC 数据库来源的 BLCA 细胞系中均存在显著差异，而且肿瘤新生抗原水平在 TCGA-BLCA 数据集中也存在显著差异。肿瘤突变负荷和肿瘤新生抗原水平的差异，可能可以解释为什么 BLCA 患者能从免疫治疗过程中获益的原因。
随后，使用 CIBERSORT 算法评估了 BLCA 患者中 22 种免疫细胞浸润的相对丰度，经比较显示，活化的树突状细胞的浸润水平在 ATM-MT 组中显著上调。接着，作者对 TCGA-BLCA 数据集和 GDSC 数据集进行了 GSEA 富集分析。结果显示，抗原刺激的免疫炎症反应，IFN 途径和巨噬细胞活化调节途径在 ATM-MT 患者中显著富集。胰岛素样生长因子受体信号通路与血管生成途径在突变患者中呈下降趋势。
2.药物敏感性
除了免疫学特征外，突变带来的药物敏感性变化也是一个重要特征。然后，基于GDSC 数据库中不同细胞系对不同化疗药物和小分子抗癌药物的 IC50，作者分析了ATM 基因突变对这些抗癌药物的敏感性改变。结果发现，ATM 突变可能导致膀胱癌对 29 种不同的化疗药物和小分子抗癌药物的敏感性增加，其中包括顺铂（P<0.05），BMS-536924（一种 IGF-1R- 5 -抑制剂（P <0.05）），莫替沙尼（一种VEGFR 抑制剂）（P <0.05）和 WHI-P97（一种 JAK 抑制剂）（P <0.05）等。此外，结合前面 GSEA 分析中的分析显示的 ATM-MT 组中胰岛素样生长因子受体表达显著下调的现象，作者发现，BMS-536924 在抑制 IGF-1R 方面可能与 ATM 突变具有协同作用
3.DNA 损伤修复通路
在肿瘤发生发展过程中，DNA 的损伤与修复过程往往伴随其中，且发挥着重要的作用，而且，DNA 的损伤与修复，可以引起氨基酸的突变，并产生肿瘤新生抗原，出现新的治疗靶点。因此，作者从 MSigDB 数据库选取了 8 条重要的 DNA 损伤与修复通路，从整体水平分析了这 8 条通路的富集水平在三个不同数据集中 ATM-MT和 ATM-WT 组之间的差异。
结果显示，免疫治疗队列中具有 ATM 突变患者的同源重组（HR）（P <0.0001），MMR（P <0.001），DSB（P <0.0001）和 FA（P <0.0001）存在明显的改变。然而，SSB（单链 DNA 结合）通路在三个数据集中均无差异。
五、总结
总结一下整篇文章的内容，综合前面分析得到的结果，作者给出了 ATM 基因突变在膀胱癌中增强免疫检查点抑制剂效果的可能机制模型图。
作者推断，ATM 突变可能通过下调胰岛素样生长因子受体途径的活性来改善膀胱癌患者的预后。ATM-MT 膀胱癌患者可能会受益于 ICI 和 IGF-1R 抑制剂的联合应用。此外，ATM-MT 组对 VEGFR 抑制剂和 JAK 抑制剂的敏感性增强，这与 GSEA的结果一致，同样表明 ATM 突变与 VEGFR 抑制剂和 JAK 抑制剂具有协同作用。
当然，对于文章而言，内容分析到此为止了，实际上，我们还是可以补充很多其他相关的内容进来，挑圈联靠，他也远远没有把整个流程走完，我们再来看看另外一篇免疫相关的文献内容。
第二篇范例文献第二篇分享的这篇文章是 2020 年 4 月发表在 Frontiers in Oncology 上的文章，文章的题目是 BTK Has Potential to Be a Prognostic Factor for Lung Adenocarcinomaand an Indicator for Tumor Microenvironment Remodeling: A Study Based on TCGAData Mining，影响因子为 4.8 分。
BTK 基因可能成为肺腺癌的预后因素和肿瘤微环境重塑的 biomarker。根据题目内容，我们可以看到，作者基于 TCGA 数据库中肺腺癌患者的表达情况，用生信分析的方法探讨了 BTK 基因对肺腺癌患者预后和肿瘤免疫微环境两个维度的影响。
初步浏览全文内容，一共用了 1 个数据集（TCGA-LUAD 项目，包含 497 例肿瘤样本和 54 例正常样本），共有 8 个 Figures。其中，Figure 1 总体陈述了文章的总体
分析流程图：
下面我们来快速浏览一下全文的结果，看看作者是如何花式凑数据的
1.靠（临床意义）
首先，作者总体评估了 TCGA-LUAD 患者的免疫浸润评分，包括免疫评分ImmuneScore，基质评分 StromalScore，以及总体评分 ESTIMAScore。基于各自的中位值，将所有患者平均分成高和低评分两组，分析了两组之间患者的总体生存率。
结果显示，具有高评分的患者预- 8 -后优于低评分的患者，尤其是对于ImmuneScore 和 ESTIMAScore，两组患者之间具有显著性差异。因此，作者提出免疫微环境的改变能影响肺腺癌患者的生存预后。在这里，两篇文献还是存在一定的共同之处的。
随后，作者进一步对这三个免疫相关评分进行了临床病理特征的相关性分析：包括Stage 分期，以及 T 分期，M 分期，N 分期。
2.挑（差异表达）&圈（功能聚类）
临床现象已经提出来了以后，接下来就需要寻找特征基因，去解释这一现象。
首先，作者基于免疫评分和基质评分，分别将所有患者平均分成了两组，根据|log(Fold change)|>1.0 和 q<0.05，各自进行了差异分析，分别得到了显著差异表达基因，也就是 DEGs（differentially expressed genes）。随后使用 Venn 图，结果显示，在两组中，共同显著上调基因共 317 个，共同显著下调基因为 62 个，最终得到 379 个 DEGs。
接着，作者对这 379 个 DEGs 分别进行了 GO 和 KEGG 功能富集分析，得到了这些差异基因可能调节的相关生物学功能和通路，并展示了其中最显著富集的四个 GO项目和四条 KEGG 通路。当然，在我们自己后续的分析过程中，可以增加两个Table，用来呈现富集结果，以增加数据量。
3.联（交互网络）
随后，作者将 379 个 DEGs 纳入 STRING 数据库中，基于高可信度（interactionconfidence value >0.95），构建了蛋白-蛋白相互作用网络，这就是我们在生信分析过程中常常听到的蛋白-蛋白互作 PPI 网络。接着，作者呈现其中 PPI 网络中的top30 基因节点。
与此同时，作者将 379 个 DEGs 基因纳入了单因素 Cox 回归分析，根据 P<0.05，以森林图的形式展示了对预后影响有意义的基因，其中 HR 值小于 1 的为保护性因素，而大于 1 的为危险性因素。最终，作者使用 Venn 图分析最终得到 PPI 网络中top30 基因和预后差异基因之间的交集基因，共 2 个，分别为 BTK 和 CCR2。
文章的结果分析到这里，第一板块的故事告一段落。简单总结一下，作者从临床现象出发，提出问题，通过一整套完整的“挑圈联靠”操作，完成了主变量分子 BTK的筛选，接下来的故事舞台就交给了基因 BTK。
4.单基因分析套路
随后，作者进行了一系列的单基因分析套路，在表达差异中，无论是 BTK 基因在癌和癌旁的表达，还是在配对组织中的表达分析，均提示 BTK 在肿瘤组织中表达降低。而且，BTK 基因高表达的患者表现出较好的生存预后。同时，作者也分析了 BTK 基因的表达与不同的临床病理特征之间的相关性接着，在富集分析中，为避免与前面呈现结果产生重复，在此作者使用了 GSEA 富集分析。基于两种数据集，分析了高和低 BTK 基因表达的肺腺癌患者之间通路富集水平的差异。
5.免疫相关分析
不忘初心，从免疫特征出发，最终还是回归免疫细胞浸润特征。但是与之前不同的是，在前期的分析过程中作者是对患者整体免疫环境的评分，在此则是基于CIBERSORT 算法，细分成 22 种不同的免疫细胞浸润水平。
首先，作者使用 barplot 图总体展示了 TCGA-LUAD 样本中 22 种免疫细胞浸润水平，并计算了这些细胞之间的相关系数，整体呈现了肺癌患者组织免疫细胞浸润的情况。
展现完整体的结果之后，作者就 BTK 基因的表达水平进行分组。
通过小提琴图呈现了两组之间不同免疫细胞浸润水平及其差异，其中 10 种免疫细胞的浸润水平在两组之间呈显著性差异（P <0.05）。接着，作者进行了相关性分析，分别计算例 Pearson 相关系数和 P 值。根据 P<0.05，结果显示其中 12 种免疫细胞水平与 BTK 表达之间呈显著相关性。最终，Venn 图显示，其中 8 种免疫细胞无论在差异分析中，还是相关性分析中均表现为与 BTK 基因显著相关性。到此呢，整个文章的结果介绍就基本结束了。
回顾一下，整篇文章仍然在总则“挑圈联靠”的范围。通过免疫浸润现象，结合临床特征的相关性，差异基因及其富集分析，PPI 网络的构建，单因素 Cox 分析，以及单基因的分析套路模式，讲述了 BTK 基因可能是影响肺腺癌患者预后和免疫微环境的靶点。整体而言，虽然结果部分使用的都是常见的图形展示方式，但对于我们初学者而言是一个很好的模仿示范材料——End