肿瘤免疫 · 实操教程Day05.K-M Plotter教程

《肿瘤免疫 · 实操教程——可以操作的肿瘤免疫学》以推文为基础，以肿瘤免疫和纯生信操作为特色，并配有视频讲解、文献资料和期刊推荐。该教程完全免费，只需要每天完成相应内容。

该课程主要包括：①肿瘤免疫前6讲理论和相应公共数据平台的实操，即理论+实操；②从绘图角度分享生信论文复现，解读思路；③结合已经发表生信论文，推荐期刊，并分享投稿策略。第一周安排如下：

Day01、生信入门总动员

Day02、定义肿瘤，肿瘤相关表型

Day03、肿瘤微环境TME

Day04、GEPIA入门教程☆☆

Day05、K-M Plotter使用教程

Day06、生信论文解读

Day07、生信期刊投稿时的要求

生存分析是生信论文中经常出现的表型，也就是说基因在正常和肿瘤组织中表达的差异，与生存参数密切相关。比如A基因在肿瘤中表达明显上调，生存参数显著下降，这就是明确的相关关系。当然，这种相关关系缺乏严谨的证据，仅仅能证明存在相关关系。至于A基因是否参与调控患者的生存，这是不能保证的。如果A基因与患者的生存密切相关，且参与调控患者的生存，就认为A基因是driver gene；如果A基因与患者的生存切相关，但不参与调控患者的生存，就认为A基因是passenger gene。

生存分析是非常重要的表型。其中，最权威的是生存分析数据库是Kaplan-Meier (K-M) Plotter，能够评估 21 种肿瘤中54000个基因与生存参数的关系。此外，PrognoScan是生存分析信息比较全面的数据库之一。

PrognoScan（网址： http://dna00.bio.kyutech.ac.jp/PrognoScan/index.html ）的页面设置相对简单，操作也不复杂。PrognoScan的优势在于数据来源是GEO，因为TCGA的生存分析工具已经非常多。比如GEPIA生存分析的数据来源是TCGA+GTEx数据库，UALCAN、cBioportal、OncoLnc和TIMER等，都来源于TCGA。PrognoScan分析生存参数时，登录官网，在对话框中输入基因名称，点击submit即可。关键是数据整理。

然后出现下表。点击cox p value一栏，使得数据按照cox p从小到大的顺序排列，以获得具有显著差异的生存分析数据。

那么，PrognoScan数据库中的p value，我们该用哪一个，是corrected p还是cox p？分别代表什么意思？如果我们已经筛选到差异表达的基因（蛋白），并已经计算了表达量和p value。如果差异基因（或蛋白）太多了，如何筛选。其中最为关键的是需要对p value进行校正。（以下内容属于拓展，可略过）

基本概念：

零假设：在随机条件下的分布。

p值：在零假设下，观测到某一特定实验结果的概率称为p值。

假阳性：得到了阳性结果，但这个阳性结果是假的。

假阴性：得到了阴性结果，但这个阴性结果是假的。

对单个基因进行统计检验（单次检验）时，假设采用p < 0.05，通常认为这个基因在两个（组）样本中的表达是有显著差异的，但是仍旧有5%的概率，这个基因并不是差异基因。

但是，当两个（组）样本中有10000个基因采用同样的方式进行统计检验（多次检验）时，这个时候就有一个问题，单次犯错的概率为0.05， 进行10000次检验的话，那么就有0.05*10000=500 个基因的差异被错误估计了。为了解决多次检验带来的问题，我们需要对多次检验进行校正。那如何校正呢？多重检验矫正！（在此介绍两种方法）

Bonferroni校正法：假设进行N次检验，p值的筛选阈值应设定为p/N。比如，进行10000次检验的话，如果设定p值为0.05，校正的筛选p值则为0.000005。只有p < 0.000005的基因才是显著性差异基因。该方法简单粗暴，但过于严格，导致筛选到的差异基因很少，甚至找不到差异的基因。

FDR（False Discovery Rate）校正法：FDR错误控制法是Benjamini于1995年提出的一种方法，基本原理是通过控制FDR值来决定p值的值域。相对Bonferroni来说，FDR用比较温和的方法对p值进行校正，试图在假阳性和假阴性间达到平衡，将假/真阳性比例控制到一定范围之内。 那么，怎么用p值来估算FDR呢？人们设计出了几种不同的估算模型。其中使用最多的是Benjamini and Hochberg法，简称BH法，分两步完成，流程如下： 

1）假设共有m个候选基因，每个基因对应的p值从小到大排列分别是p(1),p(2),…,p(m) ；

2）若想控制FDR不超过q，则只需找到最大的正整数i，使得 p(i)≤ (i*q)/m；然后，挑选对应p(1),p(2),…,p(i)的基因做为差异表达基因，这样就能从统计学上保证FDR不超过q。

根据肿瘤类型，cox p value的结果，点击探针Probe ID（蓝色划线，可以点击进去）。出现类似下图结果，可以下载，文章需要的是kaplan-Meier plot数据。截图后在PPT中标注相应信息即可，非常简单。

在生信论文中，我们一直强调差异表达和生存分析。对于差异表达，只要有p < 0.05证明即可。但是，这种差异表达对生存或者其他临床表型如复发、转移等到底有多大效用，或直白说，某个基因高表达或低表达，对表型是好还是不好的预示呢？这就要用到风险比。

在生信论文31文献中，作者在Figure2中通过森林图展示PDCD1和CTLA4等免疫检查点与多种肿瘤预后的相关性，一目了然。那么，森林图代表什么意思呢？我们简要介绍森林图（风险比）的含义。

风险比（Hazard ratio，HR）是生存分析资料中用于评估某种因素导致死亡/缓解/复发等风险改变的倍数。风险函数指危险率函数、条件死亡率、瞬时死亡率，主要通过COX回归分析得出，可用于临床治疗性研究，也可用于队列研究。计算公式为：

HR=暴露组的风险函数h1(t)/非暴露组的风险函数h2(t)，t指在相同时间点。

在大型临床研究中，HR是非常重要的评估指标之一，多以森林图（Forest plot）形式展示。森林图是指在平面直角坐标系中，以一条垂直的无效线 (横坐标刻度为1或0) 为中心，用许多条平行于横轴的线段描述每个被纳入研究的效应值和95%可信区间 (confidence interval，CI)，用菱形或其它图形描述多个研究合并的效应量及可信区间。

来自Lancet oncology上发表的题为Atezolizumab in combination with carboplatin plusnab-paclitaxel chemotherapy compared with chemotherapy alone as first-linetreatment for metastatic non-squamous non-small-cell lung cancer (IMpower130): a multicentre, randomised, open-label, phase 3 trial 的研究，根据基线患者的特征，Figure 4展示意向治疗野生型人群中总体生存的危险比森林图。无效线左边是有利于阿特珠单抗（免疫治疗）联合化疗，右边是有利于化疗组。该森林图表明，除肝转移亚组外，其余亚组均可在免疫治疗联合化疗组获益。风险比越接近1，治疗措施或干预对预后影响的权重越小；越偏离1，影响越大。

在生信分析论文中，HR比较有用，也比较常用。既可以用于生存分析参数的评估，也可以用于差异表达，年龄，淋巴转移，分化阶段和病例阶段对生存的影响。具体说起来，就是待评估因素（比如差异表达，比如年龄）对生存的影响，p < 0.05和可信度区间用于确定待评估因素的效用。至于是有利还是不利，则需视实际情况做出分析。

本次系列学习中，我们重点分享K-M Plotter（网址： http://kmplot.com/analysis/ ）的使用教程。该数据库研究最多的肿瘤包括乳腺癌（6234例）、卵巢癌（2190例）、肺癌（3452例）和胃癌（1440例），包括基因芯片和RNA-seq来源的高通量数据，数据来源包括GEO、EGA和TCGA，主要用途是基于荟萃分析（meta-analysis）发现和评估生存标志物的价值。miRNA系统包括20种不同肿瘤类型的11000个样本。

一般认为，K-M Plotter适用于生存分析、临床特征的相关分析和免疫分析。但是，随着数据库的迭代，除了以上功能以外，K-M Plotter还提供ROC Plotter统计分析、基因突变和表达分析、多重测试矫正和差异分析功能。

K-M Plotter是如何工作的呢？后台数据库是手动管理的。从GEO和TCGA下载基因表达数据、无复发和总体生存信息。数据库由PostgreSQL服务器处理，该服务器同时集成基因表达和临床数据。为了分析特定基因的预后价值，根据所提出的生物标志物的不同分位数表达，将患者样本分为两组。采用Kaplan-Meier生存图对两组患者进行比较，计算95%置信区间的危险率和log rank P值。数据库和临床数据定期监督和扩展。

使用哪个基因ID？KM-plot识别70,632个基因符号(包括HUGO基因命名委员会批准的官方基因符号、以前的符号和别名——所有这些都列在结果页面中)。由于不同的名称可能重叠，我们建议交叉检查所选基因的身份。

如何进行多基因分析？点击“使用多个基因”按钮，输入多个基因。同时对所有这些生物标志物进行分析(默认设置)，或者使用基因的平均表达。为此，勾选“使用所选探测的平均值表达式”单选按钮。最多允许65个基因。

同时，根据网址推荐，还可以对化疗效果、基因突变在肿瘤中的作用做进一步分析。K-M plotter数据库是进行生存分析最权威的网站，数据最全，尤其在乳腺癌、卵巢癌、肺癌和胃癌，包括肝癌的研究中，具有极大的优势。尽管，象GEPIA数据库和UALCAN数据库也可以对生存率做分析，但是资料的全面性和权威性上，还稍差一些。在已发表的论文中，生存分析的数据多来自K-M Plotter数据库。

1.生存分析

因为生存分析是该数据库最重要的特色，我们首先分享如何通过KM plotter数据库，直接获得生存曲线。我们以FTH1基因为例，分析FTH1对肾癌存活率的影响。

首先选择Pan-cancer RNA-seq；在gene symbol位置输入基因名称，FTH1；survival有OS和RFS，随访时间最长是240个月；数据库中只有两种肾癌的亚型，KIRC和KIRP；亚型分析包括，病理阶段stage，性别gender，种族race，分级grade，突变负担Mutation burden；细胞类型包括各种免疫细胞，可以先选择all。

然后，点击画图，draw Kaplan-Meier Plotter即可。

得出基因表达与KIRC的关系，FTH1基因高表达和低表达显著影响KIRC患者的全生存率（OS），同样步骤再做RFS。

基因表达与KIRP的关系，FTH1基因高表达和低表达显著影响KIRC患者的全生存率（OS），同样步骤再做RFS。

根据结果，截图后在PPT中进行编辑，即获得用于论文发表的数据。

2.相关分析

基因表达相关性分析和与肿瘤分期之间的相关性，是GEPIA数据库的特色，而KM plotter数据库提供全面的相关性分析，主要是肿瘤类型、临床特征与生存率的相关性。肿瘤类型包括肿瘤分期、病理类型、分化等，临床特征包括性别、胃肠穿孔（perforation）有无（疾病专有特征）、治疗方式和HER2状态等。

依次选择，即可获得相应p值，病例数等信息，用三线表的形式在PPT中编辑，即可获得用于发表的图片。如性别gender处，female和male，可以查看基因差异表达在男性和女性肿瘤患者中，生存率是否有明显不同。其余分析，可以类推。

3.免疫分析

分享这么多，不能没有实例，我们以公众号分享的文章为例，阐述K-M plotter数据库生存分析和相关分析的使用，非常值得我们借鉴在自己的论文中使用，而且跟免疫关系、治疗等关系密切，对于免疫相关的基因，尤其适用。

论文发表在AGING-US（预警期刊），论文题目——ACE2 correlated with immune infiltration serves as a prognostic biomarker in endometrial carcinoma and renal papillary cell carcinoma: implication for COVID-19。该文是分析新冠病毒受体ACE2差异表达与肿瘤免疫浸润的分析，提出ACE2对子宫内膜癌（endometrial carcinoma）和肾乳头状细胞癌（renal papillary cell carcinoma）的预后有生物标志物价值。虽然AGING-US被质疑灌水，甚至进入预警期刊，但是有的文章还是可以的，我们关注其中的思路。

作者对差异表达与存活率的相关性做分析。其实，对于为什么分析这几种肿瘤，作者没有很好的展示推理过程。我们说过，差异表达≠有临床意义，而对于差异表达不显著但低表达组和高表达组恰好有临床意义的肿瘤，也没必要作进一步的分析。有且只有如下两个条件时，生信分析才有价值：

①基因差异表达（p<0.05）；

②有临床意义（与存活率、病理分期或其他基因相关，或影响免疫浸润）。

然后，作者对差异表达与子宫内膜癌和肾乳头状细胞癌中免疫浸润的相关性进行分析。免疫浸润分析是对机制探究的一种。当然，机制探究还是要建立在表达有差异，有临床意义的基础上。结合TIMER分析，作者在分别设定免疫细胞decreased和enriched的条件下，再用K-M plotter深度分析各种免疫细胞富集或减少条件下的生存参数，巧妙地把TIMER和K-M Plotter数据库结合起来，无缝衔接，完美！推荐这种方法！

其他可以做生存分析的数据库：

GEPIA（国人之光，不断升级、迭代）

http://gepia.cancer-pku.cn/detail.php?gene=&clicktag=survival

UALCAN（甲基化是特色）

http://ualcan.path.uab.edu/

Oncolnc数据库（连mRNA, miRNA, or lncRNA也可以做生存分析）

http://www.oncolnc.org/

cBioPortal（组学分析神器也能做生存分析）

https://www.cbioportal.org/

当然，其他数据库，比如GENT2、TIMER、THPA和LinkedOmics等也可以做生存分析，特色或主要功能不是做生存分析，我们就不展开讨论了。

4.ROC绘图。ROC绘图能够使用3,104 名乳腺癌患者和2,369名卵巢癌患者的转录组水平数据将基因表达和对治疗的反应联系起来。生物标志物可以是预测性的或预后性的。预测性标记可预测特定治疗的益处，有助于选择一种特定的治疗方法而不是另一种治疗方法。预后标志物可预测疾病的自然病程（生存），与治疗无关。它可以表明需要进一步治疗，但无助于确定哪种治疗方法。ROC Plotter 是第一个用于预测性生物标志物的在线转录组级验证工具。

目前，ROC绘图可识别70,632 基因符号，包括HUGO基因命名委员会批准的官方基因符号、以前的符号和别名。所有这些都列在结果页面中。由于不同的名称可能重叠，建议交叉检查所选基因的身份。目前，该网站已经推出乳腺癌、卵巢癌、胶质瘤和结直肠癌的ROC绘图功能，其他功能尚在开发中。

比如我们以CD274基因在卵巢癌中的ROC绘图为例，简要展示其使用方法。首先，选择绿色的"ROC Plotter for ovarian cancer"，点击进入如下界面。输入基因名称CD274，根据目的，选择疗效，治疗方法、过滤条件和设置（选择No outlier和Signature即可。）

然后，点击绿色的“Calculate”选项，即可得到ROC的分析结果。可以下载数据库分析的数据PDF版本，也可以将原始数据到处用于绘图。

关于K-M Plotter（初级版）的使用就全部结束了。当然还有其他功能，感兴趣的果友可以通过官网具体了解。