2.1   类风湿关节炎差异表达基因筛选结果   通过R函数包“limma”和“affy”分析类风湿关节炎组和正常组的基因转录组数据，并根据平均绝对折叠变化法获得的log FC Cut Off及-log10 (P.Value) < 0.05，筛选差异基因总计807个，其中上调基因数量为439个，下调基因数量为368个。将差异基因进行火山图可视化(图1A)，并对上调及下调基因中表达量最高的各30位基因进行热图可视化 (图1B)。"

2.2   “WGCNA”分析结果   “WGCNA”可作为一种初步筛选重要性基因的强大方法，通过对标准化后的矩阵数据进一步处理，发现数据集中无不良样本及基因被删除，然后对所有样本聚类，发现GSM389724与GSM389732两例样本离群，当“cutHeight” =45时，切除离群样本。通过计算不同软阈值的“scale free topology model fit”，继而构建无标度网络(图2A)。随后，集合不同表达水平基因的基因模块将被分离出来，通过聚类分析和动态剪切树法识别模块，并合并相似的模块(图2B)。计算模块特征基因与临床表型(正常对照和类风湿关节炎)之间的相关性，生成相关性矩阵和P值矩阵，继而绘制模块与临床表型的相关性热图(图2C) 。结果显示绿松石色模块(r=±0.66，P=6×10-5)、粉色模块(r=±0.54，P=0.002)意义最显著。提取上述模块基因，其中绿松石色模块含基因970个，粉色模块含基因238个，共有1 208个模块特征基因被认为与临床表型最相关。"

2.3   类风湿关节炎-线粒体自噬相关基因筛选及“GO”和“KEGG”分析结果   从“MitoCarta3.0”数据库中下载人类线粒体特征表达基因，数据集包含1 136个基因，以807个类风湿关节炎差异表达基因基因、1 208个“WGCNA”特征基因和人类线粒体1 136个特征表达基因取交集，获得类风湿关节炎-线粒体自噬相关基因总计53个(图3A)。对上述交集基因进行“GO”和“KEGG”功能富集分析，“GO”富集分析排行前五位的可视化结果显示生物过程主要与有氧呼吸、细胞呼吸、电子传递链、有氧电子传递链、ATP合成偶联电子传递有关，细胞组分主要与含蛋白线粒体复合物、线粒体内膜、线粒体基质、呼吸链复合物、线粒体呼吸体相关，分子功能主要与电子传递活性、NADH脱氢酶(泛醌)活性、NADH脱氢酶(醌)活性、NADH脱氢酶活性、NAD(P)H脱氢酶(醌)活性有关(图3B)。“KEGG”功能富集分析结果包含4条细胞过程通路，富集程度由高到低依次为线粒体自噬、过氧化物酶体代谢、细胞衰老、坏死性凋亡(图3C)。"

2.4   “随机森林”机器学习算法筛选核心诊断基因   采用R函数包“randomForest”的“RF”算法对53个类风湿关节炎-线粒体自噬相关基因进行进一步筛选并进行随机森林模型误差率曲线图(图4A)和基因重要性排名图(图4B)可视化，获得30个相关性较高的基因，当重要性评分>0.6时，获得“RF”算法下相关性最佳的11个基因，分别为SLIRP、NDUFB3、PUS1、DNAJA3、PDHA1、PDHB、PTRH2、C12orf65、AKR7A2、UQCRQ、ECHS1。"

2.5   “LASSO”机器学习算法筛选核心诊断基因   通过R函数包“tidyverse”、“broom”和“glmnet”对53个类风湿关节炎-线粒体自噬相关基因进行“LASSO”回归分析，绘制出“LASSO”回归模型交叉验证误差曲线图(图5A)和“LASSO”回归系数路径图(图5B)。根据交叉验证误差曲线图可知，“LASSO”筛选的特征基因数量在6-9之间，输出基因筛选结果，获得9个相关性最高的基因，包括DNAJA3、C12orf65、AKR7A2、HADH、PDHB、BLOC1S1、COX7C、SLC25A5、NDUFS5。"

2.6   “GMM”机器学习算法同级验证核心互作基因   将“RF”算法与“LASSO”算法筛选结果取交集，进行Venn图可视化(图6A)，取得4个交集基因分别为DNAJA3、C12orf65、AKR7A2、PDHB。通过R函数包“mclust”和“SimDesign”对“RF”和“LASSO”回归分析筛选出的特征基因同时分别进行“GMM”独立聚类验证，“RF”和“LASSO”的基因筛选结果分别被归聚为9个类别，验证“RF”结果下创建了2 047个训练模型(图6B)；验证“LASSO”结果下创建了511个训练模型(图6C)。输出两侧最佳模型结果，包括DNAJA3、PUS1、AKR7A2、PTRH2、ECHS1、PDHB、C12orf65、SLC25A5、COX7C共9个基因，而DNAJA3、C12orf65、AKR7A2、PDHB均位于“GMM”算法下最优训练模型基因结果内，因此被筛选作为4个RA-Mitophagy核心互作基因。"

2.7   预测模型构建及验证结果   在二分类逻辑回归分析模式下，对4个RA-Mitophagy核心互作基因进行诊断列线图模型构建，列出各个基因的分数区间，以用来评估基因表达水平在不同分数下的风险贡献，并通过相加每个基因的分数计算出总分，通过总分评估相对应的风险概率。总分越高，患者患得类风湿关节炎的风险越大(图7A)。同时基于Bootstrap方法，重复1 000次，进行模型校正曲线绘制(图7B)，以进行内部验证，结果发现校正曲线与理想曲线贴合度高(均方误差 = 0.005 15)，证明模型的构建精度良好。进一步构建DCA (图7C)，发现预测模型曲线与所有样本全部患病的曲线之间离散度较高，表明模型的临床效用较高。"

2.8   外部患者样本的工作特征曲线验证   将GSE97779与GSE55235数据集去批次合并获得外部患者样本的独立验证集，通过对比去批次前后的PCA(图8A-B)，发现去批次效果良好。在独立外部验证集基因表达数据中进行核心互作基因的ROC分析(表2)，并进行可视化(图8C)。结果发现所有核心互作基因的AUC均大于0.7 (P < 0.01)，表明诊断效能良好，证明4个RA-Mitophagy核心互作基因在实际应用中对类风湿关节炎诊断的预测将具有良好的稳定性与准确性。"

"

2.9   “CIBERSORT”免疫浸润分析结果   通过“CIBERSORT”算法，利用“MsigDB”数据库的22种免疫细胞“gmt”数据和类风湿关节炎差异基因表达矩阵数据进行免疫浸润分析，进行免疫细胞堆积比例图可视化(图9A)及免疫细胞相关性热图可视化(图9B)。可以直观地看出，呈高正相关性的免疫细胞对有活化的记忆性CD4 T细胞与γδT细胞、活化的记忆性CD4 T细胞与静息态肥大细胞、M1型巨噬细胞与嗜酸性粒细胞、γδT细胞与静息态肥大细胞。其中，γδT细胞与静息态肥大细胞、活化的记忆性CD4 T细胞、静息态树突细胞、M0型巨噬细胞、中性粒细胞均具有较高正相关性。呈高负相关性的免疫细胞对有γδT细胞与静息态 NK 细胞、活化的记忆性CD4 T细胞与静息态 NK 细胞、活化的肥大细胞与静息态树突细胞。其中静息态 NK细胞与静息态肥大细胞、活化的记忆性CD4 T细胞、γδT细胞、静息态树突细胞、M0型巨噬细胞、中性粒细胞多种免疫细胞均有高负相关性。结果表明，在类风湿关节炎线粒体自噬过程中，γδT细胞与活化的记忆性CD4 T细胞的免疫状态较为活跃。"

2.10   “ssGSEA”单样本—免疫细胞亚群间关联性分析结果   通过R包“GSVA”相关函数初步分析各样本和免疫细胞亚群间的相关性，鉴定不同样本上相关免疫细胞的浸润丰度，继而进行免疫浸润热图可视化(图10A)及表达差异分组箱式图可视化(图10B)。免疫浸润热图显示不同样本的免疫细胞表达中，类风湿关节炎组γδT细胞、巨噬细胞的表达显著高于正常组，正常组调节性T细胞的表达显著高于类风湿关节炎组。表达差异分组箱式图显示类风湿关节炎组与正常组间活化的CD8 T细胞、CD56明亮自然杀伤细胞、嗜酸性粒细胞、γδT细胞、Ⅰ型辅助性T细胞、Ⅱ型辅助性T细胞的表达均具有显著性差异(P < 0.001)。"

2.11   类风湿关节炎-线粒体自噬核心互作基因—免疫细胞亚群间关联性分析结果  首先，通过R包“ggcorrplot”分析核心互作基因与免疫细胞亚群间的相关性，根据4个核心互作基因的表达矩阵数据与“CIBERSORT”免疫浸润结果数据，利用Shapiro-Wilk检验与封装好的Spearman相关性分析函数进行相关计算，获得基因与免疫细胞亚群间的关联性数据，继而进行核心互作基因—免疫细胞关联性“棒棒糖”图可视化(图11A-D)。图中所示，DNAJA3、C12orf65、PDHB与记忆B细胞间均具有强正相关性，DNAJA3、AKR7A2、PDHB与M0型巨噬细胞均呈强负相关性，此外，C12orf65还与活化的记忆性CD4 T细胞呈强负相关性，AKR7A2与静息态记忆性CD4 T细胞呈强正相关性。 其次，通过R包“linkET”进行汇总性可视化，绘制核心互作基因—免疫细胞亚群间关联性网络矩阵图(图11E)。可以直观地发现，多核心互作基因共同作用时，其中以记忆B细胞、活化的记忆性CD4 T细胞、调节性T细胞、M0型巨噬细胞、中性粒细胞的参与最为重要。同时，可视化网络矩阵图也显示，在类风湿关节炎线粒体自噬过程中，活化的记忆性CD4 T细胞与γδT细胞、静息态肥大细胞之间的正相关表达作用最强，单核细胞和中性粒细胞之间的负相关表达作用最强。 "

2.12   基于“ssGSEA”的核心互作基因相关生物学功能分析结果   “ssGSEA”算法下共获得821种与RA-Mitophagy核心互作基因相关的“GO”生物学功能(|cor| > 0.8，P < 0.001)，其中生物过程数量为600，细胞组分数量为72，分子功能数量为149，并筛选出与每个核心互作基因正相关性和负相关性最强的前三位生物学功能进行展示(表3)。同时获得48条与RA-Mitophagy核心互作基因相关的“KEGG”信号通路(|cor| > 0.8，P < 0.001)，同样筛选出与每个核心互作基因正相关性和负相关性最强的前三位信号通路进行展示(表4)。最后，逐条筛选上述全部获得的821种“GO”生物学功能和48条“KEGG”信号通路，发现有6种生物学功能与线粒体自噬密切相关 (表5)。"

"

"