大家好，今天和大家分享的是19年1月份在J. Cell. Physiol(IF：4.522) 上的一篇文章，"Identification of the molecular subgroups in coronary artery disease by gene expression profiles"，文章用生信方法研究冠状动脉疾病，通过一致性聚类根据基因表达谱对此病患者进行分型，并通过差异表达分析，WCGNA模型，通路分析等方法研究不同亚型间在临床特征和基因表达模式上的差异。
Identification of the molecular subgroups in coronary artery disease by gene expression profiles
通过基因表达谱识别冠状动脉疾病的分子亚型
一. 文章背景
冠状动脉疾病(Coronary artery disease，CAD)，是全球最常见的造成死亡的心血管疾病之一。随着高通量测序技术的发展，许多研究借此对CAD的分子机制展开研究。但是以往大量的研究都聚焦于CAD患者与正常患者间的差异而忽视了CAD患者间的差异。事实是CAD患者存在临床上的异质性。作者希望像研究肿瘤一样，根据CAD患者的基因表达谱进行分型，研究引起CAD的分子机制。
二.研究思路
三.结果解读
开始之前先看一下数据来源，作者一共分析了GEO数据库中两个数据集的样本的外周血芯片数据(nCAD=352，nControl=263)。两个数据集分别为GSE12288(n=222)，GSE20686(n=393)。其中GSE20686包含两个批次GSE20680，GSE20681)。两个数据集的两个芯片平台共同检测到的基因为11314个
1. 通过跨平台标准化去除批次效应
作者使用SVA包中的ComBat函数消除三组不同批次样本间的批次效应，并用消除批次前后的样本主成分分析图来验证。
A：消除批次效应前的主成分分析图，三个批次的样本分被区分开(不同颜色的点代表不同批次的样本)
B：消除批次效应后的主成分分析图，三个批次的样本混合在一起，说明批次效应被消除
图1. 对各样本的主成分分析
2. 对CAD样本进行一致性聚类
在消除了各样本基因表达谱间的批次效应后，作者想要根据基因表达谱对CAD样本进行聚类，划分CAD亚型。这里作者采用的是一致性聚类的方法(ConsensusClusterPlus包)。B图是不同聚类数目下各亚组的聚类一致性得分。当聚类数目为3时，各个亚组的聚类一致性得分均大于0.8，故作者决定将样本聚成三类。A图时聚类数目为3时的一致性矩阵热图。
作者通过以上方法将样本分为了3个亚组(n=68,183,97)
图2. 对CAD样本采取一致性聚类
3. 分析不同CAD亚型在临床特征上的差异
根据一致性聚类将CAD患者分为三个亚组后，作者对三组间在临床特征上的差异展开研究
A：三组在男性比例上的差异，只在Ⅱ和Ⅲ组间存在显著差异
B：三组在CAD指数(一种评价CAD患者等级的指标，这里只有GSE12288数据集中样本的数据)，Ⅰ组与其它两组相比，CAD指数显著偏高
C：三组在年龄上的差异，Ⅰ组的年龄与其它两组相比显著偏大
图3. 比较不同亚组之间在临床特征上的差异
此外，作者用双因素方差分析比较自己得到的CAD分类，年龄(两个主效应)以及两者的交互效应对CAD患者CAD指数的影响(表1)，发现自己得到的CAD分类是一个独立于年龄的显著效应(p=0.0369)，说明CAD患者根据基因表达谱进行分类不仅影响CAD患者的CAD指数，也暗示了CAD患者内在的生物学差异
表1. 方差分析的结果
4. 识别各亚组中的基因共表达模式
作者在三个CAD亚组间，每个CAD亚组与对照组之间进行基因差异表达分析(表2)。与以往比较某个基因在两组间的logFC不同，这里认定一个基因是差异表达基因的条件是该基因在在两个比较组间均值(mean)的差异要大于0.2(可能由于在CAD中基因差异表达倍数并不大所以选用此标准)且adjp<0.05
在三个CAD亚组间的差异表达分析中，Ⅰ组有2410个表达上调的DEGs，Ⅱ组中有40个表达上调的DEGs，Ⅲ组中有2217个表达上调的DEGs(表2第4列)
表2. 不同亚组间基因差异分析结果
上述通过组间比较得到的各组中表达上调的基因，需要排除在对照组中也高表达的可能
A-C：作者根据上文差异分析得到的三个基因集在各CAD亚组与对照组间进行GSEA分析，结果表明这些基因显著富集在各自的亚组中
D：WGCNA分析得到的6个共表达模块中的基因(纵轴)在各亚组以及对照组样本(横轴)中的表达量热图
图4. CAD亚组中特定的表达上调基因的表达模式
5. WCGNA
作者根据上文的得到的在三个CAD亚组间差异表达的共4667个基因对所有CAD样本进行WCGNA分析(基因加权共表达网络分析)。共得到6个共表达模块(图4.D)，表一中的第5列标出了各亚组间表达上调的基因对应的共表达模块
图4.D 共表达模块中基因在各样本中表达量热图
6. 对各共表达模块中的基因进行功能富集分析
在得到了6个共表达模块后，作者从MsigDB数据库中获取了全部的KEGG通路基因集，对共表达模块中的基因近行KEGG通路富集分析，选出了各模块富集结果最显著的通路进行展示(图5，颜色越深表示富集越显著)
7.分析模块特征基因和临床特征的相关性
最后，作者分析了6个共表达模块的模块特征基因与CAD患者年龄以及CAD指数的相关性
模块1，3，6与CAD患者的年龄以及CAD指数呈显著正相关
模块2，4与CAD患者的年龄以及CAD指数呈显著负相关
图. 6模块特征基因与临床特征的相关性热图
小结
本文作为一篇非肿瘤生信分析例子有很多地方值得我们学习。作者先对跨平台的数据集消除批次效应，以根据基因表达谱对CAD患者分型为目的，用一致性聚类的方法将CAD患者分为三个亚组，并比较了三组在临床特征上的差异。接着作者分析了三个亚组间各自的表达上调基因，据此对所有样本进行WCGNA分析，得到了六个共表达模块。最后作者对各模块基因集进行KEGG通路分析以及研究各模块特征基因与临床特征的相关性。
还是和往常一样，后台回复「4f」，即可获取今天小编为大家解读的文献。本期的分享就到这里啦，一起期待下一期的精彩分享吧~
▼▼是否遗漏了往期精彩生信解析没看呢？▼▼
关注科研菌回复"1"即可获取"28G全网最全科研绘图素材"、"谷歌上网助手"(梯子)和"科研常用软件合集"、“肿瘤/消化内科指南”、"考研复试大礼包"；回复2020，可获得2020年生信套路合集；回复“单细胞文献”，可获得单细胞文献合集(均无需转发/集赞)
