外周血基因表达谱探测新型冠状病毒肺炎潜在诊断标志物

张思嘉， 张 顺， 蔡 挺

(中国科学院大学宁波生命与健康产业研究院， 宁波315100)

新型冠状病毒(2019-nCoV)流感与其他冠状病毒肺炎或季节性流感相比，虽然多数流行特征和临床学特征相类似，但其发病例数更多、传播更快、重症率和病死率更高，且目前尚无有效药物治疗[1-3]。新型冠状病毒肺炎(COVID-19)早期临床及影像表现与多种呼吸系统疾病，包括间质性肺炎、肺损伤、细菌性肺炎等有诸多相似之处[4-5]，主要表现为发热、干咳、气促、外周血白细胞降低及胸部CT改变。COVID-19患者起病更急，在1～2周出现呼吸困难，严重者肺部病变进展迅速，并累及心、肾器官最终引发多脏器功能衰竭[6-8]。

目前对COVID-19患者的筛查与确诊主要基于流行病学史、临床表现及核酸检测的综合分析。流感、普通肺炎和新型冠状病毒肺炎的流行病特征和临床特征相似，为疑似病例诊断带来较大挑战。基于患者的临床，目前RT-PCR核酸检测作为COVID-19检验金标准，具有时间短、操作简便等优势，但无法提供临床分型等信息，从而影响患者后期的分级诊疗。Shu等[9]基于机器学习的方法，通过对血浆蛋白质组学分析，共识别11种蛋白质作为预测COVID-19患者临床结局的生物标志物，用于密切监测和评估COVID-19患者的病情发展并为其临床治疗提供及时建议。Shen等[10]通过血清样本的蛋白质组学和代谢组学数据，确定22种蛋白质7种代谢产物用于COVID-19重症患者的早期识别。Levi等[11]发现COVID-19 重症病例常出现凝血异常和炎症反应并反映在多种血液生物标志物的含量变化上，如炎症反应可导致TNF-α、IL-1、IL-6 的水平显著升高，凝血异常表现为D-二聚体水平升高、血小板减少和凝血酶原时间增加。

为探究COVID-19与季节性流感、细菌性肺炎患者临床表现差异的原因，筛选潜在COVID-19诊断标志物，本研究基于早期COVID-19、季节性流感和细菌性肺炎患者的外周血基因表达谱，采用偏最小二乘分析、lasso回归、权重基因共表达网络分析等生物信息学手段，用于筛选鉴别早期COVID-19患者的血浆诊断生物标志物，同时针对早期COVID-19患者血浆中表达异常的基因进行进一步分析，旨在为后续治疗提供基础依据。

1 材料与方法

1.1 基因芯片数据集信息

通过GEO数据库(http:∥www.ncbi.nlm.nih.gov/geo/)下载高通量测序数据集GSE161731。该数据包括19例早期COVID-19患者外周血样本、19例健康人群样本、24例细菌性肺炎患者样本、60例其他冠状病毒肺炎患者样本、17例季节性流感人群样本。利用R软件对数据进行标准化预处理。

1.2 潜在差异基因筛选

采用R软件的limma软件包，针对全部疾病样本(n=120)和健康控制样本(n=19)进行差异分析，差异基因的筛选标准为以Benjamini-Hochberg方法对假设检验的P值进行校正，得到错误发现率(FDR)，满足FDR<0.01的基因作为潜在差异基因，并通过热图对前100个潜在差异基因的表达进行可视化。

1.3 偏最小二乘分析组间差异性

以疾病类型作为分类变量，健康人群的分组名称为Healthy，早期COVID-19患者的分组名称为COVID-19，季节性流感人群的分组名称为Influenza，其他冠状病毒肺炎患者分组名称Other COV，细菌性肺炎患者分组名称为Bacteria。利用1.2节中获得差异基因的表达矩阵建立偏最小二乘分析模型(partial least-square method, PLS)。通过样本得分图对比各组基因表达的差异情况，利用载荷图筛选差异较显著的基因进行后续分析，筛选标准为变量投影重要度(variable importance in projection, VIP)大于1。

1.4 Lasso回归对潜在COVID-19诊断标志物的筛选

采用R软件中的glmnet包，对1.3节筛选出的显著差异基因进行Lasso回归算法分析，回归模型类型为multinomial。以最终建立的回归模型中权重系数不为0的变量作为COVID-19患者的潜在诊断标志物，并采用R软件中的pROC包，针对Lasso回归模型筛选出的潜在诊断标志物绘制受试者工作特征曲线(receiver operating characteristic, ROC)，用于评估标志物对早期COVID-19患者的诊断性能，以ROC曲线下方的面积大小(area under curve, AUC)大于0.8为标准，确定用于区分COVID-19患者与其他组样本的最终诊断标志物。

1.5 基于WGCNA探索不同疾病组的关键基因及通路

利用R软件的 WGCNA 软件包对显著差异基因(P<0.01且VIP>1)进行共表达网络的构建[9]：首先构建基因共表达相似性矩阵，确定软阈值后将相似性矩阵转换为邻接矩阵。再将邻接矩阵转换成拓扑矩阵，并采用拓扑覆盖法对差异基因进行层次聚类，按照混合动态剪切树的方法确定基因模块，同时绘制树状图并对基因模块进行可视化，其中每个模块的最小基因数目为30。最后将基因模块与疾病类型进行关联分析，用于寻找不同疾病对应的显著相关基因模块。同时利用KEGG数据库(kyoto encyclopedia of genes and genomes, KEGG)，对不同疾病分组对应的显著相关基因模块进行功能和通路富集分析。

2 结果与分析

2.1 数据处理及差异表达分析

首先对GSE161731数据集进行数据预处理：若缺失值占全部样本20%以上，则该基因将被过滤；若同一基因对应多个探针，则使用该基因对应所有探针表达值的平均值；采用R语言中的scale函数对数据中心化和标准化。最终获得包含139例样本及12 880种化合物的表达矩阵。利用R语言edge R包筛选差异表达基因，基于FDR<0.01的标准，共筛选出1 397个潜在差异表达基因。通过火山图(图1)对潜在差异基因进行可视化，并用热图显示差异最明显的前100个基因(logFC绝对值前100)各样本中的表达情况(图2)。

图1 疾病组与健康控制组外周血基因表达谱差异火山图Figure 1 Volcano plot of peripheral blood gene expression profile between disease groups and healthy controls

每个小方格代表一个基因，颜色表示基因表达量大小，表达量越大方格色度越偏向红色，表达量越小方格色度越偏向绿色。图2 健康控制组与4种疾病组外周血基因表达谱差异分析热图Figure 2 Heatmap of peripheral blood gene expression profile between 4 disease groups and healthy controls

2.2 偏最小二乘法分析各组基因表达谱的差异性

利用偏最小二乘回归分析对外周血基因表达谱数据进行模式识别分析，得到样本得分图(图3)和变量载荷图(图 4)。其中样本得分图的X轴与Y轴分别代表第一主成分(PC1)、第二主成分(PC2)，第一主成分可解释数目最多的原变量，第二主成分其次。图3所示，细菌性肺炎与其他组样本点明显分离，聚集于第三象限，说明细菌性肺炎和其他疾病组的基因轮廓存在显著差异。其他冠状病毒肺炎(Other COV)样本与早期COVID-19样本区分度较差，均位于原点附近，说明二者的基因轮廓高度相似。变量载荷(图4)及变量投影重要度(图5)中红色标记的化合物为显著差异基因(VIP>1)，显著差异基因主要分布于距离原点较远的各象限内。通过基因的VIP值排序，筛选出样本分组贡献度前10的基因为IFI27、SPATS2L、SERPING1、CMPK2、BATF2、SEPTIN4、RNASE1、CD177、MYBL2和CCNB2。

图3 健康控制组与4种疾病组外周血基因轮廓PLS样本得分图Figure 3 The PLS scores plots of healthy control group and 4 disease groups gene profile

图4 健康控制组与4种疾病组外周血基因轮廓PLS变量载荷图Figure 4 The PLS loading plots of healthy control group and 4 disease groups gene profile

图5 各基因在PLS-DA模型中的变量投影重要度Figure 5 Variable importance of each gene in the PLS model

2.3 Lasso回归对潜在COVID-19诊断标志物的筛选

根据差异分析与偏最小二乘分析的结果，共有509种基因同时满足FDR<0.01及VIP>1，利用Lasso Logistic回归模型对509个显著差异基因进行变量选择。该模型通过交叉验证获得最优λ为0.006，最终确定的COVID-19诊断标志物共有17种，包括MYBL2、PKMYT1、HJURP、TCN2、TTC24、ESPL1、GZMK、RPA3、ATP5F1E、CD2、ZAP70、IL15RA、NFYB、CIAO2A、PLRG1、GPATCH11和MRPL33。通过受试者工作特征曲线(ROC)计算以上17种标志物对早期COVID-19患者的诊断能力，并评价每种标志物对区分早期COVID-19患者与其他相似疾病的诊断价值(表1)。

选取诊断性能排名前3的基因绘制小提琴图(图6)。结果显示该基因在各样本组的表达差异。其中：PKMYT1在区分健康人群与早期COVID-19患者(AUC=0.972)、季节性流感与早期COVID-19患者(AUC=0.861)、其他冠状病毒肺炎与早期COVID-19患者(AUC=0.873)均具有较强的诊断力；PLRG1在区分季节性流感与早期COVID-19患者(AUC=0.848)，其他冠状病毒肺炎与早期COVID-19患者(AUC=0.846)具有较强的诊断力；GZMK在区分细菌性肺炎与早期COVID-19患者(AUC=0.956)，季节性流感与早期COVID-19患者(AUC=0.845)具有较强的诊断能力；MYBL2在区分健康人群与早期COVID-19患者(AUC=0.936)，其他冠状病毒肺炎与早期COVID-19患者(AUC=0.876)具有较强的诊断能力。此外CD2(AUC=0.989)与ZAP70(AUC=0.956)亦属于判断细菌性肺炎与早期COVID-19患者的重要诊断标志物，GPATCH11可有效区别健康人群与早期COVID-19患者(AUC=0.953)。

图6 用于区分早期COVID-19患者的高诊断能力基因的表达Figure 6 Expression of gene with significant diagnostic ability for patients with early-phase COVID-19

2.4 基于WGCNA探索不同疾病组的关键基因及通路

基于筛选得到的509个显著差异基因在4个不同疾病组的表达数据，利用加权基因共表达网络分析算法(WGCNA)构建共表达模块。参考Shen等[10]的方法，首先按照无尺度网络的标准，确定软阈值(power值)为17(图7)；再根据最佳power值进行基因模块划分。最终得到6个有效模块(图8)，每个模块包含基因数目大于30。

图7 共表达模块软阈值的确定Figure 7 Analysis of network topology for various soft-thresholding powers

图8 联合加权共表达网络:基因层次聚类树及基因模块Figure 8 WGCNA: clustering dendrogram of genes with assigned modules

基因模块与临床性状的关联分析，主要用于确定与临床性状显著相关的基因模块。本研究以不同疾病为临床性状，通过计算基因模块与临床性状的相关性系数和P值，确定每种疾病的特异性基因模块。如图9所示，根据疾病组与基因模块相关性的显著性检验结果(P值)和相关性系数的大小，确定每个疾病组对应的特异性基因模块：Bacteria(细菌性肺炎)特异性基因模块为Blue；Cov Other(其他冠状病毒肺炎)特异性基因模块为Yellow；早期COVID-19特异性基因模块为Turquoise；Influenza(季节性流感)特异性基因模块为Brown。

方格颜色表示相关系数大小，红色为正相关，蓝色为负相关；方格首行数字表示相关系数，次行数字表示P值。图9 基因模块与不同疾病分组的相关分析Figure 9 Associations between gene modules and different diseases

2.5 基因模块功能分析

对分析得到的6个基因模块分别进行KEGG富集分析，结果显示不同基因模块具有特异性的代谢通路、生物学过程。例如，Yellow模块显著与T细胞受体信号传导途径相关，而Turquoise、Grey模块则分别富集到氧化磷化与细胞周期。疾病富集分析表明：Yellow、Turquoise模块与COVID-19密切相关，其中Yellow模块包含的COVID-19特异基因共13种(表2)，属于Turquoise模块的COVID-19特异基因共27种，上述基因多用于编码核糖体蛋白；Blue模块则富集到多种细菌感染性疾病；Green模块则与神经系统性疾病相关。

3 讨论与结论

对新型冠状病毒肺炎，做到早发现、早诊断、早隔离、早治疗极为重要[12]。目前RT-PCR核酸检测是病原学诊断的金标准，具有时间短、操作简便等优势。但受采样质量的影响，核酸采样仍存在局限性，有报道指出鼻咽拭子和口咽拭子阳性率不高易出现假阴性[13-14]。本研究针对早期COVID-19患者，以及与COVID-19具有相似临床症状的3类患者的外周血进行分析，以期在分子层面做到迅速从季节性流感、其他冠状病毒肺炎、细菌性肺炎区分新型冠状肺炎患者，从而对接诊患者做到有效的隔离和对症治疗，减少不必要的资源浪费，同时针对COVID-19起病急、发展快、临床症状严重的原因做出一定探究。

本研究通过对比健康控制人群以及其他疾病患者的外周血mRNA表达水平，共鉴定出1 397个差异基因。基于1 397个差异基因表达水平进行偏最小二乘判别分析(PLS)，以VIP > 1作为筛选标准，共纳入509个显著差异基因。对显著差异基因进行Lasso回归分析，最终筛选出COVID-19潜在诊断标志物17种，包括MYBL2、PKMYT1、HJURP、TCN2、TTC24、ESPL1、GZMK、RPA3、ATP5F1E、CD2、ZAP70、IL15RA、NFYB、CIAO2A、PLRG1、GPATCH11和MRPL33。

在筛选出的诊断标志物中，许多标志物与内质网及高尔基体蛋白的合成和加工有关。这可能是由于冠状病毒作为单链正股RNA病毒，缺乏增殖所需要的酶系统需要通过吸附于宿主细胞的核糖体，利用宿主细胞的能源系统翻译出病毒复制酶和结构蛋白，再利用复制酶以及RNA模板合成子代RNA，最后与蛋白组装成为子代病毒颗粒完成病毒的复制[15]。如PLRG1作为CDC5L的核心组件，是PRP19-CDC5L复合物的组成部分，主要参与Pre-mRNA转化为成熟mRNA 用于蛋白质的生产[16-17]。此外SARS-Cov-2引起的COVID-19往往伴随过激的免疫应答，甚至可加重患者肺组织炎性损伤及其他器官损害。本研究筛选出的部分诊断标志物在COVID-19早期患者中的过表达也印证了该观点，如：GZMK编码的Granzyme K(GZMK)被认为在细胞免疫中起到保护宿主细胞的作用，通过裂解宿主细胞表面附着的抗原，防止宿主细胞成为T细胞或NK细胞的靶细胞，该基因表达量的显著上升可能是机体防止过度免疫的一种保护机制[18]；CD2基因编码的蛋白质是T细胞上重要的功能标志分子，既参与T细胞活化，又参与T细胞与其他细胞的粘连过程，起到介导信号传递和细胞激活的重要作用。研究指出CD2基因的过表达与多种自身免疫性疾病密切相关[19-21]；ZAP70编码的蛋白酪氨酸激酶(PTK)主要用于激活T细胞早期活化信号的转导，是一类重要的信号转导分子，同时在维持T细胞的免疫反应和免疫耐受的平衡中起到重要作用，可防止免疫失调引起的对正常组织细胞的损害[22-24]。

本研究通过WGCNA技术建立权重基因共表达网络模块，并对基因模块与疾病分组进行相关分析，发现每种疾病的外周血转录组具有一定的特异性和协同性：特异性表现为不同疾病分组与不同基因模块相关，而不同基因模块具有特异性功能；系统性表现为同一模块基因不仅在转录水平上具有相似模式，还在基因功能及相关生理活动中存在密切联系，并在信号通路之间相互调节、相互配合。KEGG通路富集结果显示，COVID-19早期患者组的特异模块Turquoise最主要的两个信号通路为核糖体蛋白合成及氧化磷酸化，而病毒颗粒的快速增殖需要耗费大量能量，即随着蛋白合成、核酸复制速率的加快伴随着细胞内磷酸化水平的逐渐升高。这些基因在COVID-19早期患者的高表达间接说明了与普通流感、病毒性肺炎相比早期COVID-19患者体内SARS-Cov-2病毒增殖速率更快，病情发展更迅速。

通过对比细菌性肺炎、季节性流感等4类人群与早期COVID-19患者的外周血mRNA的表达水平，共鉴定出17种潜在诊断标志物，这些标志物在调节细胞周期、氧化磷酸化、细胞免疫等方面发挥重要作用。此外，利用WGCNA技术建立权重基因共表达网络模块并筛选的COVID-19特异性基因模块，发现早期COVID-19患者外周血中多种核糖体蛋白基因与氧化磷酸化相关mRNA的表达量显著上升，可能与SARS-Cov-2快速增殖有关。本研究尚存在局限性和不足：案例数偏少，增大了抽样误差的风险；细菌性肺炎样本与其他组样本基因表达谱差异过大，可能会导致差异基因筛选的偏倚：筛选到的诊断标志物的具体适用情况尚不明确，后期需引入其他样本，采用qPCR等手段对潜在诊断标志物的表达量进行验证。