今天给大家带来的是3分 学习笔记。文末阅读原文可获取笔记原文。
这篇通过相似性网络融合(SNF)和一致性聚类(CC)划分整合多组学数据对软组织肉瘤(STS)进行聚类。并筛选预后不良和预后良好亚群之间差异表达的lncRNAs、miRNAs和mrna，以此构建ceRNA网络进行分析。
题目：综合性聚类揭示一种预后不良的软组织肉瘤新亚型
一.研究背景
软组织肉瘤(STS)是一组恶性肿瘤。目前大多数的STS临床研究都是基于一种或几种特定的组织学类型，而由于STS罕见且组织学亚型太多（目前约有50种），所以患者数量有限，临床试验进展缓慢。采用基于基因组学的基因分型方法，将具有相似基因组特征的肿瘤归入同一类别，从而建立新的分类系统。这有利于分子靶点的筛选和未来临床试验的开展，有助于开发更合理、特异和有效的治疗方法。
二.分析流程
三.结果解读
1.样本详情
从TCGA数据库中下载STS的RNASeq、miRNAseq、体细胞突变数据和患者的临床信息（表1），共247例。
表1 训练组患者的临床信息
2.一致性聚类
使用Cancer Subtypes包中的 "ExecuteSNF.CC "函数对样本进行聚类。该函数联合使用相似性网络融合(SNF)和一致性聚类(CC)的方法整合多组学数据进行聚类。
- 首先，SNF分别利用RNAseq数据和miRNA数据计算患者之间的相似度。
- 然后，整合RNAseq和miRNA这两类数据类型计算出的患者之间的相似度，根据整合的相似度构建融合的患者相似度矩阵，作为CC的样本距离，进行聚类。
通过累计分布函数曲线下面积相对变化量来确定最佳聚类数。当聚类数(k)为3时，面积变化最大(图1A)，因此STS被确定聚为3个类。图1C用热图展现了样本距离。结果显示，聚类内的样本彼此接近，而不同聚类之间的样本彼此相距较远。
还使用轮廓宽度（Silhouette Width）来判断聚类效果。轮廓宽度结合了内聚度和分离度两种因素进行判断，其取值范围为[-1,1]。将所有点的轮廓宽度求平均，就是该聚类结果总的轮廓宽度。取值越接近1则说明聚类效果越好；相反，取值越接近-1则说明聚类效果越差。图1B结果表明，每个聚类的轮廓宽度和总轮廓宽度都接近1，聚类效果良好。
图1 对STS进行聚类
3.三个亚群的特点
为了比较聚类到的亚群与组织学亚型的关系，统计了每个亚群中组织学亚型的分布情况（图2A）：
- 对于聚类到的亚群：
- C1包含各种的组织学亚型，其中dedifferentiated liposarcoma最多（36.29%）。
- C2主要包含dedifferentiated liposarcoma和leiomyosarcoma(LMS)这两种组织学亚型（共占74%）。
- C3主要包含leiomyosarcoma(LMS)（99%）。
- 对于组织学亚型：
- dedifferentiated liposarcoma在C1和C2中占比很高。
- leiomyosarcoma(LMS)主要存在于C3中，在C1和C2中也有一些。
- Myxofibrosarcoma, pleomorphic MFH/UPS, UPS主要存在于C1中。
图2A 聚类到的三个亚群中肿瘤组织学亚型的分布
以聚类结果分组进行生存分析，结果显示，C2患者总体生存率显著低于C1和C3患者(P = 0.02)(图2B)。但考虑到不同聚类中患者组织学亚型的构成不同（C3主要为LMS亚型），为了验证总体生存期差异是由不同聚类而不是组织学亚型引起的，进一步针对C2和C3中的LMS亚型患者进行了生存分析。结果表明，C2中LMS患者的总体生存率显著低于C3中的LMS患者(P = 0.0015)(图2C)。
图2A 聚类到的三个亚群生存分析结果
接下来将每个聚类分组分别对其它两个聚类分组进行差异表达分析（ |log2 (FC)| > 2.0 , P < 0.01）。每个聚类的两次差异表达分析中均高表达的基因被定义为亚群特异性基因。并将两次差异表达分析中FC之和最大的基因定义为亚群的标记基因（图3），分别是：
- C1：LncRNA LINC01133, mRNA DKK1, miRNA has-miR-511-5p。
- C2：LncRNA MEG3, mRNA BMP7, miRNA has-miR-483-3p。
- C3：LncRNA HAND2-AS1, mRNA MYH11, miRNA has-miR-133a-3p。
以往研究表明，TP53和RB1是STS中最常见的突变基因，利用卡方检验分析比较这两个基因的突变率在不同亚群之间的差异。结果表明，TP53和RB1在C1（TP53：17.0%，RB1：4.0%）和C3（TP53：15.0%，RB1：6.9%）的突变率显著高于在C2（TP53：2.0%，RB1：0.0%）中的突变率（P＜0.05）。
图3 每个亚群中表达差异最大的lncRNAs、mRNAs、miRNAs和基因突变的热图
4.构建ceRNA网络并进行功能富集分析
上一步分析中可知C2患者预后差，但这并不是由抑癌基因TP53和RB1的突变所驱动的。尝试基于C2和C1&C3之间差异表达基因（lncRNAs、miRNAs和mRNAs），利用GDCRNATools包构建ceRNA网络。为了深入了解ceRNA网络的功能，还进行了GO富集分析。这样就可能有助于了解C2的生物学特征。
- C2和C1&C3的差异分析中，发现了100个差异表达的LncRNAs，152个差异表达的miRNAs和1663个差异表达的mRNAs。
- 建立了一个具有91个节点(4个lncRNA、8个miRNA和79个mRNA)和167条边的ceRNA网络（图4）。竞争性的基因对通过miRNAs连接，形成ceRNA调控网络。
- GO富集分析表明，ceRNA调控网络主要在发育和cellular processes中发挥作用（图5）。
图4 构建的ceRNA调控网络
图5 ceRNA的GO富集分析结果
5.筛选hub网络并深入研究
接下来筛选hub ceRNA调控网络。对上一步构建的ceRNA网络的每个节点进行Cox回归分析，并以每个节点的中位表达水平将样本分为两组生存分析。如果一个竞争基因对(lncRNA-miRNA-mRNA)的所有成分都都有显著的生存差异(Cox P值<0.05，logRank P值<0.05)，则选择该基因对作为hub ceRNA网络的一部分。
- 在hub ceRNA调控网络中，LncRNA KCNQ1OT1、JARID2、CDK6、DNMT3A、TET1竞争性地结合了hsa-miR-29c-3p（图6）。
图6 hub ceRNA调控网络
随后进一步研究hub ceRNA调控网络中的这6个基因：
- KCNQ1OT1和mRNAs（JARID2、CDK6、DNMT3A、TET1）在C2中的表达水平明显高于C1、C3（图7A）。
- KCNQ1OT1与has-miR-29c-3p的表达呈负相关，与mRNAs的表达呈正相关；mRNAs均与hsa-miR-29c-3p的表达呈负相关（图7B）。
- 生存分析方面（图7C）：
- KCNQ1OT1和mRNAs的高表达与患者预后较差有关；hsa-miR-39c-3p的高表达则表示预后较好。
- 对mRNAs进行主成分分析（PCA），根据第一主成分（PC1）的中位表达水平将样本分为两组，进行生存分析。结果显示PC1表达水平越高预后越差。
- 将这6个基因表达值经Cox回归系数加权求和后构建风险评分评分公式，并计算每个样本的风险评分。根据风险评分中位值分为高低组，进行生存分析。结果显示，高风险评分意味着较差的总体生存率。
图7 hub ceRNA调控网络特点
小结
从TCGA数据库中获取STSs的RNAseq和miRNAseq数据，通过相似性网络融合(SNF)和一致性聚类(CC)划分整合多组学数据对样本进行聚类，划分成了三个亚群。分析各个亚群临床特点，进行生存分析，找到预后较差的亚群。进行突变分析，发现预后较差的亚群抑癌基因突变反而低，说明预后差并非由抑癌基因突变导致。随后筛选预后不良和预后良好亚群之间差异表达的lncRNAs、miRNAs和mrna，并构建ceRNA网络。进行了功能富集分析，并从构建的ceRNA网络中提取了一个hub网络。hub网络中的基因可以为未来STS的诊治提供参考靶点。
