大规模-高通量-多组学-居然只能把食管癌分成两个亚型

最近在朋友圈刷到了2021年8月份发表在Nature Communications杂志的文章，标题是：《Large-scale and high-resolution mass spectrometry-based proteomics profiling defines molecular subtypes of esophageal cancer for therapeutic targeting》，我就下载看了看，研究者们对124对食管癌患者的肿瘤和配对癌旁组织进行了蛋白质组学和磷酸蛋白质组学分析，揭示了EC中失调的蛋白质和磷酸化位点，并根据蛋白质组分析将EC分为S1和S2两个分子亚型。

• 其中：S2 subtype char- acterized by the upregulation of spliceosomal and ribosomal proteins, and being more aggressive

• 而且 根据 ELOA and SCAF4, 两个基因，可以构建诊断和预后模型

• 并且针对 S2 subtype, and three candidate drugs

首先看看实验设计

可以看到，有两个技术：25 groups( 125 samples )  for TMT proteomics, and 31 paired samples were subjected to lable-free phosphoproteomics.

简单的层次聚类，并且辅助热图以及PCA的可视化，就可以看到正常样品和肿瘤样品是泾渭分明：

蛋白质组学数据本质上也是表达量矩阵

所以前面的质量控制也是三张图，我很早之前就在生信技能树的教程：《你确定你的差异基因找对了吗？》提到过，必须要对你的转录水平的全局表达矩阵做好质量控制，最好是看到标准3张图：

• 左边的热图，说明我们实验的两个分组，normal和npc的很多基因表达量是有明显差异的
• 中间的PCA图，说明我们的normal和npc两个分组非常明显的差异
• 右边的层次聚类也是如此，说明我们的normal和npc两个分组非常明显的差异

如果分组在3张图里面体现不出来，实际上后续差异分析是有风险的。这个时候需要根据你自己不合格的3张图，仔细探索哪些样本是离群点，自行查询中间过程可能的问题所在，或者检查是否有其它混杂因素，都是会影响我们的差异分析结果的生物学解释。

然后因为有proteins and phosphoproteins 两个技术，所以这个文章有两个表达量矩阵，各种都可以在正常样品和肿瘤样品，就有两次火山图，两个差异分析的对比：

这些我们都反复演示过。目前简单的差异分析流程，基本上转录组测序技术和芯片技术拿到的表达量矩阵后续分析大同小异，哪怕是本文的基于蛋白质组学技术拿到的表达量矩阵，也是如此。公众号推文在：

• 解读GEO数据存放规律及下载，一文就够
• 解读SRA数据库规律一文就够
• 从GEO数据库下载得到表达矩阵 一文就够
• GSEA分析一文就够（单机版+R语言版）
• 根据分组信息做差异分析- 这个一文不够的
• 差异分析得到的结果注释一文就够

使用R包ConsensusClusterPlus进行一致性聚类

前面的蛋白质组学也是矩阵，所以很容易使用R包ConsensusClusterPlus进行一致性聚类，这里作者对全部的肿瘤样品定下来了2个分组，然后依据两个分组继续对proteins and phosphoproteins 两个技术的矩阵进行差异分析：

可以看到，两个分组，有生存差异，一般来说，这就是普通临床研究的终点了。

首先呢，让我们仔细看看上面的火山图（都是S2相对于S1的差异哦，因为S2的预后差的组 ）：

• 在 6468 proteins that were quantifiable in the 124 paired samples  矩阵里面的差异分析，是  984 and 1012 proteins were significantly upregulated and downregulated ，如果提高阈值，比如  1.5-fold,  是 137 and 93 个上下调基因。

• 在  61,471 phosphorylation sites quantifiable in at least half of the samples ，分组是 15 patients belonged to S1 and 16 belonged to the S2 subtype。差异分析结果是 1446 and 1415 sites were significantly increased and decreased, respectively, in the S2 samples

构建诊断模型和预后模型

这里的诊断模型就是  SVM model implemented in the R package mlr ，如果有数据，代码实现起来并不难，就是需要一些摸索。

最后入选的两个基因是；

• SCAF4, also known as splicing fac- tor, arginine/serine-rich 15, belongs to the splicing factor SR family.
• ELOA, also known as elongin A, is a component of the SIII complex, which activates RNA polymerase II elongation by suppressing transient pausing of the polymerase

如下所示：

模型效果确实好，因为是诊断模型。An SVM classification model with ELOA and SCAF4 as features was con- structed for subtype prediction， and an AUC of 0.976 was obtained for the 124 patients we performed proteomics analysis on (Cohort 1)

这两个基因恰好也是预后相关的基因，可以比较好的区分生存，其实主要是因为这个S1和S2本来就是预后有差异，那么它们的差异基因，本来就是区分S1和S2就同时区分了生存。生存分析是目前肿瘤等疾病研究领域的点睛之笔！我在生信技能树多次分享过生存分析的细节；

• 人人都可以学会生存分析（学徒数据挖掘）
• 学徒数据挖掘之谁说生存分析一定要按照表达量中位值或者平均值分组呢？
• 基因表达量高低分组的cox和连续变量cox回归计算的HR值差异太大?
• 学徒作业-两个基因突变联合看生存效应
• TCGA数据库里面你的基因生存分析不显著那就TMA吧
• 对“不同数据来源的生存分析比较”的补充说明
• 批量cox生存分析结果也可以火山图可视化
• 既然可以看感兴趣基因的生存情况，当然就可以批量做完全部基因的生存分析
• 多测试几个数据集生存效应应该是可以找到统计学显著的！
• 我不相信kmplot这个网页工具的结果(生存分析免费做)
• 为什么不用TCGA数据库来看感兴趣基因的生存情况
• 200块的代码我的学徒免费送给你，GSVA和生存分析
• 集思广益-生存分析可以随心所欲根据表达量分组吗
• 生存分析时间点问题
• 寻找生存分析的最佳基因表达分组阈值
• apply家族函数和for循环还是有区别的（批量生存分析出图bug）
• TCGA数据库生存分析的网页工具哪家强
• KM生存曲线经logRNA检验后也可以计算HR值

最后也不可免俗的使用公共数据集挖掘可能的药物作用

首先是tumor and non-tumor samples 差异分析，得到  189 upregulated and 271 downregulated proteins in tumor samples ，然后去Connectivity Map (CMAP)搜索上下调基因相关的药物。

然后是搜索这个S2亚型特异性的药物 ，使用 86 upregulated and 24 downregulated proteins in the S2 subtype ，得到  six candidate drugs，可以使用 six EC cell lines—KYSE30, KYSE150, KYSE450, TE1, TE3, and TE5来验证这些药物的作用。

两个蛋白质组学技术

前面反复提到了 proteins and phosphoproteins 两个技术， 中文一般来说叫做 TMT蛋白质组检测和  label-free磷酸化蛋白质组检测，大家可以自行搜索学习。

我们也有一下简单的基础介绍，目录是：

• 蛋白质组学第1期-认识基础概念
• 蛋白质组学第2期-认识蛋白质组学原始数据
• 蛋白质组学第3期-蛋白质组学的三大元素
• 蛋白质组学第4期 文章搜库过程复现
• 蛋白质组学第5期搜库软件之 MaxQuant 再介绍
• 蛋白质组学第6期 搜库软件之 MaxQuant 结果数据介绍
• 蛋白质组学第7期 复现文章数据- 预处理之Perseus 的使用
• 蛋白质组学第8期 文章复现之数据处理
• 蛋白质组学第9期 文章数据分析之差异蛋白筛选和功能分析
• 蛋白质组学第10期 定量方法介绍

这个蛋白质组学技术也是越来越正规了，可以看到原始数据，比如这个文章的就是：

• The raw files of proteome and phosphoproteome datasets can be obtained from PRIDE database (accession number PXD021701)  or iProX database (accession number IPX0002501000) .

但是从原始数据文件，处理到可以进行下游统计可视化的表达量矩阵，仍然是任重道远。

写在文末

我在《生信技能树》，《生信菜鸟团》，《单细胞天地》的大量推文教程里面共享的代码都是复制粘贴即可使用的， 有任何疑问欢迎留言讨论，也可以发邮件给我，详细描述你遇到的困难的前因后果给我，我的邮箱地址是 jmzeng1314@163.com

如果你确实觉得我的教程对你的科研课题有帮助，让你茅塞顿开，或者说你的课题大量使用我的技能，烦请日后在发表自己的成果的时候，加上一个简短的致谢，如下所示：

We thank Dr.Jianming Zeng(University of Macau), and all the members of his bioinformatics team, biotrainee, for generously sharing their experience and codes.

十年后我环游世界各地的高校以及科研院所（当然包括中国大陆）的时候，如果有这样的情谊，我会优先见你