張思嘉， 張 順， 蔡 挺

(中國(guó)科學(xué)院大學(xué)寧波生命與健康產(chǎn)業(yè)研究院， 寧波315100)

新型冠狀病毒(2019-nCoV)流感與其他冠狀病毒肺炎或季節(jié)性流感相比，雖然多數(shù)流行特征和臨床學(xué)特征相類(lèi)似，但其發(fā)病例數(shù)更多、傳播更快、重癥率和病死率更高，且目前尚無(wú)有效藥物治療[1-3]。新型冠狀病毒肺炎(COVID-19)早期臨床及影像表現(xiàn)與多種呼吸系統(tǒng)疾病，包括間質(zhì)性肺炎、肺損傷、細(xì)菌性肺炎等有諸多相似之處[4-5]，主要表現(xiàn)為發(fā)熱、干咳、氣促、外周血白細(xì)胞降低及胸部CT改變。COVID-19患者起病更急，在1～2周出現(xiàn)呼吸困難，嚴(yán)重者肺部病變進(jìn)展迅速，并累及心、腎器官最終引發(fā)多臟器功能衰竭[6-8]。

目前對(duì)COVID-19患者的篩查與確診主要基于流行病學(xué)史、臨床表現(xiàn)及核酸檢測(cè)的綜合分析。流感、普通肺炎和新型冠狀病毒肺炎的流行病特征和臨床特征相似，為疑似病例診斷帶來(lái)較大挑戰(zhàn)。基于患者的臨床，目前RT-PCR核酸檢測(cè)作為COVID-19檢驗(yàn)金標(biāo)準(zhǔn)，具有時(shí)間短、操作簡(jiǎn)便等優(yōu)勢(shì)，但無(wú)法提供臨床分型等信息，從而影響患者后期的分級(jí)診療。Shu等[9]基于機(jī)器學(xué)習(xí)的方法，通過(guò)對(duì)血漿蛋白質(zhì)組學(xué)分析，共識(shí)別11種蛋白質(zhì)作為預(yù)測(cè)COVID-19患者臨床結(jié)局的生物標(biāo)志物，用于密切監(jiān)測(cè)和評(píng)估COVID-19患者的病情發(fā)展并為其臨床治療提供及時(shí)建議。Shen等[10]通過(guò)血清樣本的蛋白質(zhì)組學(xué)和代謝組學(xué)數(shù)據(jù)，確定22種蛋白質(zhì)7種代謝產(chǎn)物用于COVID-19重癥患者的早期識(shí)別。Levi等[11]發(fā)現(xiàn)COVID-19 重癥病例常出現(xiàn)凝血異常和炎癥反應(yīng)并反映在多種血液生物標(biāo)志物的含量變化上，如炎癥反應(yīng)可導(dǎo)致TNF-α、IL-1、IL-6 的水平顯著升高，凝血異常表現(xiàn)為D-二聚體水平升高、血小板減少和凝血酶原時(shí)間增加。

為探究COVID-19與季節(jié)性流感、細(xì)菌性肺炎患者臨床表現(xiàn)差異的原因，篩選潛在COVID-19診斷標(biāo)志物，本研究基于早期COVID-19、季節(jié)性流感和細(xì)菌性肺炎患者的外周血基因表達(dá)譜，采用偏最小二乘分析、lasso回歸、權(quán)重基因共表達(dá)網(wǎng)絡(luò)分析等生物信息學(xué)手段，用于篩選鑒別早期COVID-19患者的血漿診斷生物標(biāo)志物，同時(shí)針對(duì)早期COVID-19患者血漿中表達(dá)異常的基因進(jìn)行進(jìn)一步分析，旨在為后續(xù)治療提供基礎(chǔ)依據(jù)。

1 材料與方法

1.1 基因芯片數(shù)據(jù)集信息

通過(guò)GEO數(shù)據(jù)庫(kù)(http:∥www.ncbi.nlm.nih.gov/geo/)下載高通量測(cè)序數(shù)據(jù)集GSE161731。該數(shù)據(jù)包括19例早期COVID-19患者外周血樣本、19例健康人群樣本、24例細(xì)菌性肺炎患者樣本、60例其他冠狀病毒肺炎患者樣本、17例季節(jié)性流感人群樣本。利用R軟件對(duì)數(shù)據(jù)進(jìn)行標(biāo)準(zhǔn)化預(yù)處理。

1.2 潛在差異基因篩選

采用R軟件的limma軟件包，針對(duì)全部疾病樣本(n=120)和健康控制樣本(n=19)進(jìn)行差異分析，差異基因的篩選標(biāo)準(zhǔn)為以Benjamini-Hochberg方法對(duì)假設(shè)檢驗(yàn)的P值進(jìn)行校正，得到錯(cuò)誤發(fā)現(xiàn)率(FDR)，滿(mǎn)足FDR<0.01的基因作為潛在差異基因，并通過(guò)熱圖對(duì)前100個(gè)潛在差異基因的表達(dá)進(jìn)行可視化。

1.3 偏最小二乘分析組間差異性

以疾病類(lèi)型作為分類(lèi)變量，健康人群的分組名稱(chēng)為Healthy，早期COVID-19患者的分組名稱(chēng)為COVID-19，季節(jié)性流感人群的分組名稱(chēng)為Influenza，其他冠狀病毒肺炎患者分組名稱(chēng)Other COV，細(xì)菌性肺炎患者分組名稱(chēng)為Bacteria。利用1.2節(jié)中獲得差異基因的表達(dá)矩陣建立偏最小二乘分析模型(partial least-square method, PLS)。通過(guò)樣本得分圖對(duì)比各組基因表達(dá)的差異情況，利用載荷圖篩選差異較顯著的基因進(jìn)行后續(xù)分析，篩選標(biāo)準(zhǔn)為變量投影重要度(variable importance in projection, VIP)大于1。

1.4 Lasso回歸對(duì)潛在COVID-19診斷標(biāo)志物的篩選

采用R軟件中的glmnet包，對(duì)1.3節(jié)篩選出的顯著差異基因進(jìn)行Lasso回歸算法分析，回歸模型類(lèi)型為multinomial。以最終建立的回歸模型中權(quán)重系數(shù)不為0的變量作為COVID-19患者的潛在診斷標(biāo)志物，并采用R軟件中的pROC包，針對(duì)Lasso回歸模型篩選出的潛在診斷標(biāo)志物繪制受試者工作特征曲線(xiàn)(receiver operating characteristic, ROC)，用于評(píng)估標(biāo)志物對(duì)早期COVID-19患者的診斷性能，以ROC曲線(xiàn)下方的面積大小(area under curve, AUC)大于0.8為標(biāo)準(zhǔn)，確定用于區(qū)分COVID-19患者與其他組樣本的最終診斷標(biāo)志物。

1.5 基于WGCNA探索不同疾病組的關(guān)鍵基因及通路

利用R軟件的 WGCNA 軟件包對(duì)顯著差異基因(P<0.01且VIP>1)進(jìn)行共表達(dá)網(wǎng)絡(luò)的構(gòu)建[9]：首先構(gòu)建基因共表達(dá)相似性矩陣，確定軟閾值后將相似性矩陣轉(zhuǎn)換為鄰接矩陣。再將鄰接矩陣轉(zhuǎn)換成拓?fù)渚仃?，并采用拓?fù)涓采w法對(duì)差異基因進(jìn)行層次聚類(lèi)，按照混合動(dòng)態(tài)剪切樹(shù)的方法確定基因模塊，同時(shí)繪制樹(shù)狀圖并對(duì)基因模塊進(jìn)行可視化，其中每個(gè)模塊的最小基因數(shù)目為30。最后將基因模塊與疾病類(lèi)型進(jìn)行關(guān)聯(lián)分析，用于尋找不同疾病對(duì)應(yīng)的顯著相關(guān)基因模塊。同時(shí)利用KEGG數(shù)據(jù)庫(kù)(kyoto encyclopedia of genes and genomes, KEGG)，對(duì)不同疾病分組對(duì)應(yīng)的顯著相關(guān)基因模塊進(jìn)行功能和通路富集分析。

2 結(jié)果與分析

2.1 數(shù)據(jù)處理及差異表達(dá)分析

首先對(duì)GSE161731數(shù)據(jù)集進(jìn)行數(shù)據(jù)預(yù)處理：若缺失值占全部樣本20%以上，則該基因?qū)⒈贿^(guò)濾；若同一基因?qū)?yīng)多個(gè)探針，則使用該基因?qū)?yīng)所有探針表達(dá)值的平均值；采用R語(yǔ)言中的scale函數(shù)對(duì)數(shù)據(jù)中心化和標(biāo)準(zhǔn)化。最終獲得包含139例樣本及12 880種化合物的表達(dá)矩陣。利用R語(yǔ)言edge R包篩選差異表達(dá)基因，基于FDR<0.01的標(biāo)準(zhǔn)，共篩選出1 397個(gè)潛在差異表達(dá)基因。通過(guò)火山圖(圖1)對(duì)潛在差異基因進(jìn)行可視化，并用熱圖顯示差異最明顯的前100個(gè)基因(logFC絕對(duì)值前100)各樣本中的表達(dá)情況(圖2)。

圖1 疾病組與健康控制組外周血基因表達(dá)譜差異火山圖Figure 1 Volcano plot of peripheral blood gene expression profile between disease groups and healthy controls

每個(gè)小方格代表一個(gè)基因，顏色表示基因表達(dá)量大小，表達(dá)量越大方格色度越偏向紅色，表達(dá)量越小方格色度越偏向綠色。圖2 健康控制組與4種疾病組外周血基因表達(dá)譜差異分析熱圖Figure 2 Heatmap of peripheral blood gene expression profile between 4 disease groups and healthy controls

2.2 偏最小二乘法分析各組基因表達(dá)譜的差異性

利用偏最小二乘回歸分析對(duì)外周血基因表達(dá)譜數(shù)據(jù)進(jìn)行模式識(shí)別分析，得到樣本得分圖(圖3)和變量載荷圖(圖 4)。其中樣本得分圖的X軸與Y軸分別代表第一主成分(PC1)、第二主成分(PC2)，第一主成分可解釋數(shù)目最多的原變量，第二主成分其次。圖3所示，細(xì)菌性肺炎與其他組樣本點(diǎn)明顯分離，聚集于第三象限，說(shuō)明細(xì)菌性肺炎和其他疾病組的基因輪廓存在顯著差異。其他冠狀病毒肺炎(Other COV)樣本與早期COVID-19樣本區(qū)分度較差，均位于原點(diǎn)附近，說(shuō)明二者的基因輪廓高度相似。變量載荷(圖4)及變量投影重要度(圖5)中紅色標(biāo)記的化合物為顯著差異基因(VIP>1)，顯著差異基因主要分布于距離原點(diǎn)較遠(yuǎn)的各象限內(nèi)。通過(guò)基因的VIP值排序，篩選出樣本分組貢獻(xiàn)度前10的基因?yàn)镮FI27、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回歸對(duì)潛在COVID-19診斷標(biāo)志物的篩選

根據(jù)差異分析與偏最小二乘分析的結(jié)果，共有509種基因同時(shí)滿(mǎn)足FDR<0.01及VIP>1，利用Lasso Logistic回歸模型對(duì)509個(gè)顯著差異基因進(jìn)行變量選擇。該模型通過(guò)交叉驗(yàn)證獲得最優(yōu)λ為0.006，最終確定的COVID-19診斷標(biāo)志物共有17種，包括MYBL2、PKMYT1、HJURP、TCN2、TTC24、ESPL1、GZMK、RPA3、ATP5F1E、CD2、ZAP70、IL15RA、NFYB、CIAO2A、PLRG1、GPATCH11和MRPL33。通過(guò)受試者工作特征曲線(xiàn)(ROC)計(jì)算以上17種標(biāo)志物對(duì)早期COVID-19患者的診斷能力，并評(píng)價(jià)每種標(biāo)志物對(duì)區(qū)分早期COVID-19患者與其他相似疾病的診斷價(jià)值(表1)。

選取診斷性能排名前3的基因繪制小提琴圖(圖6)。結(jié)果顯示該基因在各樣本組的表達(dá)差異。其中：PKMYT1在區(qū)分健康人群與早期COVID-19患者(AUC=0.972)、季節(jié)性流感與早期COVID-19患者(AUC=0.861)、其他冠狀病毒肺炎與早期COVID-19患者(AUC=0.873)均具有較強(qiáng)的診斷力；PLRG1在區(qū)分季節(jié)性流感與早期COVID-19患者(AUC=0.848)，其他冠狀病毒肺炎與早期COVID-19患者(AUC=0.846)具有較強(qiáng)的診斷力；GZMK在區(qū)分細(xì)菌性肺炎與早期COVID-19患者(AUC=0.956)，季節(jié)性流感與早期COVID-19患者(AUC=0.845)具有較強(qiáng)的診斷能力；MYBL2在區(qū)分健康人群與早期COVID-19患者(AUC=0.936)，其他冠狀病毒肺炎與早期COVID-19患者(AUC=0.876)具有較強(qiáng)的診斷能力。此外CD2(AUC=0.989)與ZAP70(AUC=0.956)亦屬于判斷細(xì)菌性肺炎與早期COVID-19患者的重要診斷標(biāo)志物，GPATCH11可有效區(qū)別健康人群與早期COVID-19患者(AUC=0.953)。

圖6 用于區(qū)分早期COVID-19患者的高診斷能力基因的表達(dá)Figure 6 Expression of gene with significant diagnostic ability for patients with early-phase COVID-19

2.4 基于WGCNA探索不同疾病組的關(guān)鍵基因及通路

基于篩選得到的509個(gè)顯著差異基因在4個(gè)不同疾病組的表達(dá)數(shù)據(jù)，利用加權(quán)基因共表達(dá)網(wǎng)絡(luò)分析算法(WGCNA)構(gòu)建共表達(dá)模塊。參考Shen等[10]的方法，首先按照無(wú)尺度網(wǎng)絡(luò)的標(biāo)準(zhǔn)，確定軟閾值(power值)為17(圖7)；再根據(jù)最佳power值進(jìn)行基因模塊劃分。最終得到6個(gè)有效模塊(圖8)，每個(gè)模塊包含基因數(shù)目大于30。

圖7 共表達(dá)模塊軟閾值的確定Figure 7 Analysis of network topology for various soft-thresholding powers

圖8 聯(lián)合加權(quán)共表達(dá)網(wǎng)絡(luò):基因?qū)哟尉垲?lèi)樹(shù)及基因模塊Figure 8 WGCNA: clustering dendrogram of genes with assigned modules

基因模塊與臨床性狀的關(guān)聯(lián)分析，主要用于確定與臨床性狀顯著相關(guān)的基因模塊。本研究以不同疾病為臨床性狀，通過(guò)計(jì)算基因模塊與臨床性狀的相關(guān)性系數(shù)和P值，確定每種疾病的特異性基因模塊。如圖9所示，根據(jù)疾病組與基因模塊相關(guān)性的顯著性檢驗(yàn)結(jié)果(P值)和相關(guān)性系數(shù)的大小，確定每個(gè)疾病組對(duì)應(yīng)的特異性基因模塊：Bacteria(細(xì)菌性肺炎)特異性基因模塊為Blue；Cov Other(其他冠狀病毒肺炎)特異性基因模塊為Yellow；早期COVID-19特異性基因模塊為T(mén)urquoise；Influenza(季節(jié)性流感)特異性基因模塊為Brown。

方格顏色表示相關(guān)系數(shù)大小，紅色為正相關(guān)，藍(lán)色為負(fù)相關(guān)；方格首行數(shù)字表示相關(guān)系數(shù)，次行數(shù)字表示P值。圖9 基因模塊與不同疾病分組的相關(guān)分析Figure 9 Associations between gene modules and different diseases

2.5 基因模塊功能分析

對(duì)分析得到的6個(gè)基因模塊分別進(jìn)行KEGG富集分析，結(jié)果顯示不同基因模塊具有特異性的代謝通路、生物學(xué)過(guò)程。例如，Yellow模塊顯著與T細(xì)胞受體信號(hào)傳導(dǎo)途徑相關(guān)，而Turquoise、Grey模塊則分別富集到氧化磷化與細(xì)胞周期。疾病富集分析表明：Yellow、Turquoise模塊與COVID-19密切相關(guān)，其中Yellow模塊包含的COVID-19特異基因共13種(表2)，屬于Turquoise模塊的COVID-19特異基因共27種，上述基因多用于編碼核糖體蛋白；Blue模塊則富集到多種細(xì)菌感染性疾病；Green模塊則與神經(jīng)系統(tǒng)性疾病相關(guān)。

3 討論與結(jié)論

對(duì)新型冠狀病毒肺炎，做到早發(fā)現(xiàn)、早診斷、早隔離、早治療極為重要[12]。目前RT-PCR核酸檢測(cè)是病原學(xué)診斷的金標(biāo)準(zhǔn)，具有時(shí)間短、操作簡(jiǎn)便等優(yōu)勢(shì)。但受采樣質(zhì)量的影響，核酸采樣仍存在局限性，有報(bào)道指出鼻咽拭子和口咽拭子陽(yáng)性率不高易出現(xiàn)假陰性[13-14]。本研究針對(duì)早期COVID-19患者，以及與COVID-19具有相似臨床癥狀的3類(lèi)患者的外周血進(jìn)行分析，以期在分子層面做到迅速?gòu)募竟?jié)性流感、其他冠狀病毒肺炎、細(xì)菌性肺炎區(qū)分新型冠狀肺炎患者，從而對(duì)接診患者做到有效的隔離和對(duì)癥治療，減少不必要的資源浪費(fèi)，同時(shí)針對(duì)COVID-19起病急、發(fā)展快、臨床癥狀嚴(yán)重的原因做出一定探究。

本研究通過(guò)對(duì)比健康控制人群以及其他疾病患者的外周血mRNA表達(dá)水平，共鑒定出1 397個(gè)差異基因?；? 397個(gè)差異基因表達(dá)水平進(jìn)行偏最小二乘判別分析(PLS)，以VIP > 1作為篩選標(biāo)準(zhǔn)，共納入509個(gè)顯著差異基因。對(duì)顯著差異基因進(jìn)行Lasso回歸分析，最終篩選出COVID-19潛在診斷標(biāo)志物17種，包括MYBL2、PKMYT1、HJURP、TCN2、TTC24、ESPL1、GZMK、RPA3、ATP5F1E、CD2、ZAP70、IL15RA、NFYB、CIAO2A、PLRG1、GPATCH11和MRPL33。

在篩選出的診斷標(biāo)志物中，許多標(biāo)志物與內(nèi)質(zhì)網(wǎng)及高爾基體蛋白的合成和加工有關(guān)。這可能是由于冠狀病毒作為單鏈正股RNA病毒，缺乏增殖所需要的酶系統(tǒng)需要通過(guò)吸附于宿主細(xì)胞的核糖體，利用宿主細(xì)胞的能源系統(tǒng)翻譯出病毒復(fù)制酶和結(jié)構(gòu)蛋白，再利用復(fù)制酶以及RNA模板合成子代RNA，最后與蛋白組裝成為子代病毒顆粒完成病毒的復(fù)制[15]。如PLRG1作為CDC5L的核心組件，是PRP19-CDC5L復(fù)合物的組成部分，主要參與Pre-mRNA轉(zhuǎn)化為成熟mRNA 用于蛋白質(zhì)的生產(chǎn)[16-17]。此外SARS-Cov-2引起的COVID-19往往伴隨過(guò)激的免疫應(yīng)答，甚至可加重患者肺組織炎性損傷及其他器官損害。本研究篩選出的部分診斷標(biāo)志物在COVID-19早期患者中的過(guò)表達(dá)也印證了該觀點(diǎn)，如：GZMK編碼的Granzyme K(GZMK)被認(rèn)為在細(xì)胞免疫中起到保護(hù)宿主細(xì)胞的作用，通過(guò)裂解宿主細(xì)胞表面附著的抗原，防止宿主細(xì)胞成為T(mén)細(xì)胞或NK細(xì)胞的靶細(xì)胞，該基因表達(dá)量的顯著上升可能是機(jī)體防止過(guò)度免疫的一種保護(hù)機(jī)制[18]；CD2基因編碼的蛋白質(zhì)是T細(xì)胞上重要的功能標(biāo)志分子，既參與T細(xì)胞活化，又參與T細(xì)胞與其他細(xì)胞的粘連過(guò)程，起到介導(dǎo)信號(hào)傳遞和細(xì)胞激活的重要作用。研究指出CD2基因的過(guò)表達(dá)與多種自身免疫性疾病密切相關(guān)[19-21]；ZAP70編碼的蛋白酪氨酸激酶(PTK)主要用于激活T細(xì)胞早期活化信號(hào)的轉(zhuǎn)導(dǎo)，是一類(lèi)重要的信號(hào)轉(zhuǎn)導(dǎo)分子，同時(shí)在維持T細(xì)胞的免疫反應(yīng)和免疫耐受的平衡中起到重要作用，可防止免疫失調(diào)引起的對(duì)正常組織細(xì)胞的損害[22-24]。

本研究通過(guò)WGCNA技術(shù)建立權(quán)重基因共表達(dá)網(wǎng)絡(luò)模塊，并對(duì)基因模塊與疾病分組進(jìn)行相關(guān)分析，發(fā)現(xiàn)每種疾病的外周血轉(zhuǎn)錄組具有一定的特異性和協(xié)同性：特異性表現(xiàn)為不同疾病分組與不同基因模塊相關(guān)，而不同基因模塊具有特異性功能；系統(tǒng)性表現(xiàn)為同一模塊基因不僅在轉(zhuǎn)錄水平上具有相似模式，還在基因功能及相關(guān)生理活動(dòng)中存在密切聯(lián)系，并在信號(hào)通路之間相互調(diào)節(jié)、相互配合。KEGG通路富集結(jié)果顯示，COVID-19早期患者組的特異模塊Turquoise最主要的兩個(gè)信號(hào)通路為核糖體蛋白合成及氧化磷酸化，而病毒顆粒的快速增殖需要耗費(fèi)大量能量，即隨著蛋白合成、核酸復(fù)制速率的加快伴隨著細(xì)胞內(nèi)磷酸化水平的逐漸升高。這些基因在COVID-19早期患者的高表達(dá)間接說(shuō)明了與普通流感、病毒性肺炎相比早期COVID-19患者體內(nèi)SARS-Cov-2病毒增殖速率更快，病情發(fā)展更迅速。

通過(guò)對(duì)比細(xì)菌性肺炎、季節(jié)性流感等4類(lèi)人群與早期COVID-19患者的外周血mRNA的表達(dá)水平，共鑒定出17種潛在診斷標(biāo)志物，這些標(biāo)志物在調(diào)節(jié)細(xì)胞周期、氧化磷酸化、細(xì)胞免疫等方面發(fā)揮重要作用。此外，利用WGCNA技術(shù)建立權(quán)重基因共表達(dá)網(wǎng)絡(luò)模塊并篩選的COVID-19特異性基因模塊，發(fā)現(xiàn)早期COVID-19患者外周血中多種核糖體蛋白基因與氧化磷酸化相關(guān)mRNA的表達(dá)量顯著上升，可能與SARS-Cov-2快速增殖有關(guān)。本研究尚存在局限性和不足：案例數(shù)偏少，增大了抽樣誤差的風(fēng)險(xiǎn)；細(xì)菌性肺炎樣本與其他組樣本基因表達(dá)譜差異過(guò)大，可能會(huì)導(dǎo)致差異基因篩選的偏倚：篩選到的診斷標(biāo)志物的具體適用情況尚不明確，后期需引入其他樣本，采用qPCR等手段對(duì)潛在診斷標(biāo)志物的表達(dá)量進(jìn)行驗(yàn)證。