今天介绍的是一篇发表在我们的老朋友期刊——Cancer Biomarkers上的文章。一眼看去，映入眼帘的就是Machine learning algorithm and deep neural networks，脑海便浮现出一个想法，不简单！！！首先文章选题方向的是我们熟悉的肝癌基因，然后利用一些统计学和机器学习方法进行临床与基因表达的关联分析，筛选出了一些与肝癌具有显著相关性的基因，最后讲这些基因作为构建深度神经网络模型的训练集，搭建完成后再使用富集分析揭示关键基因所在的代谢途径。怎么样？是不是思路比较简单？并没有想象中的那么复杂。接下来让小果带你解读一下这篇文章。

题目：机器学习算法和深度神经网络确定了一个新的肝细胞癌亚型

杂志：Cancer Biomarkers

影响因子：IF=3.828

发表时间：2022年9月

在世界上最常见的癌症中，肝细胞癌(HCC)排名第五，也是导致癌症相关死亡的三大原因之一。在北美和欧洲，这些历史上较低的HCC发病率近年来有所增加。且由于疾病早期缺乏特异性特征，患者通常在疾病进展的晚期被诊断出来。因此，迫切需要一种准确预测HCC患者预后生存风险的模型，以指导临床治疗。

使用“Limma”包对GSE112790数据进行差异表达分析，结果如图1(A)所示，获得624个下调基因和794个上调基因。同时对TCGA-LIHC数据进行差异表达分析，如图1(B)所示，获得464个下调基因和2394个上调基因。且图1(C)展示了所有的差异基因热图。且从两次差异分析中挑选了重叠基因，如图1(D)所示。

图1 GSE112790和TCGA-LIHC数据差异表达分析

使用单变量COX回归模型来检验每个基因的预后影响。筛查阈值p < 0.05的基因，提示其预后效果有统计学意义，基于单因素COX结果，进一步考虑临床因素，包括性别、年龄、T分期、分期和分级分期。然后通过多元COX回归计算各基因的回归系数、风险比、相应的置信区间和p值。最终得到了546个符合条件的差异表达基因。

得到生存相关的关键基因后，利用无监督K-means距离中的PAM函数计算不同K值下的簇间轮廓宽度，其聚类效果如图2(A)所示。可以看到，当K=2时，K-means模型的聚类效果最佳。其空间分布如图2(B)所示。365例患者被分为两种新的亚型包含243例低危患者(C1)和122例高危患者(C2)。随后采用Kaplan-Meier方法分析两种亚型的生存率，K-M曲线如图2(C)所示(p < 0.0001)。该模型ROC曲线的具体参数图2(D)。最后比较了两种亚型中不同AJCC分级患者的比例，结果如图2(E)和2(F)所示，可以发现晚期癌症患者在高危亚型中的比例远高于晚期癌症患者在低危亚型中的比例。

图2 无监督K-means模型聚类分析

随机森林模型的决策树数量参数（ntree）由错误树图优化（图 3A），“mtry”使用默认参数。这模型的特征向量是生存相关基因，并且类别向量使用高风险和低风险亚型。经过反复训练，模型的出袋误差率为4.66%。将接近度矩阵转化为距离矩阵，计算每个MDS轴的重要性。从而绘制MDS图(图3B)。高危和低危患者的样本分布在MDS1轴的两端，相应的ROC曲线(图3C)也验证了随机森林模型的可靠性。根据基尼系数从高到低对基因进行排序，输出最重要的前50个基因(图3D)。

使用R包“glmnet”实现LASSO逻辑回归，族参数设置为“二项式”。不断增大惩罚(Lambda)，使模型中各基因的系数调整为趋于0，系数可视化如图4A和4B所示。交叉验证是优化模型以确保模型的均方误差满足要求的一个很好的选择(图4C)。每个模型对应的ROC曲线下面积存储在AUC中并输出(图4D)，其中AUC值均大于0.9。

图3 随机森林中候选基因的识别

图4 LASSO中候选基因的识别

利用Lasso模型中系数部位0的基因与随机森林筛选的TOP 50基因进行比较，最后如图5(A)所示，得到17个基因，将这17个基因作为DNN模型训练的特征向量，类别向量为无监督聚类得到的两组亚型。将具有完整临床数据的TCGA-LIHC数据集随机划分三分之一作为模型的预测集，其余三分之二样本作为训练集。图5(C-F)展示了训练过程中数据集的交叉熵损失。将DNN预测结果和K-Means聚类结果放入混淆矩阵结果如图5(B)。最终将预测集放入训练完成的KNN模型，利用预测结果绘制生存曲线，如图5(G)所示，可以看出该模型可以准确预测患者的亚型(p=0.00027)。

使用R包“enrichment go”和“enrichment kegg”对差异基因进行超几何分布检验分析。如图6所示为KEGG富集分析结果，如图7所示为GO富集分析结果。使用MSigDB数据库作为背景基因集，来展开GSEA分析。结果如图8(A)和8(B)展示了上调、下调基因所富集的8条通路。

图5  DNN预测模型的建立

图6  C1和C2之间的KEGG通路富集分析

图7  基因本体(GO)功能富集分析

图8  GSEA富集分析结果