SARS-CoV-2病毒感染潛在關(guān)鍵分子生物標(biāo)志物及免疫浸潤特征分析

于敏，王敏，魏延煥，劉毅毅

1.日照市人民醫(yī)院急診醫(yī)學(xué)科，山東 日照276800；

2.海軍軍醫(yī)大學(xué)第三附屬醫(yī)院檢驗科，上海200438

新型冠狀病毒肺炎（corona virus disease 2019，COVID-19）是一種由嚴重急性呼吸綜合征冠狀病毒2（severe acute respiratory syndrome-coronavirus 2，SARS-CoV-2）感染造成的嚴重傳染病［1］。目前，已在世界范圍內(nèi)持續(xù)流行時間超過2年，致使人類生命健康及經(jīng)濟發(fā)展面臨嚴重威脅［2］。隨著對SARS-CoV-2病毒相關(guān)研究的不斷深入及治療方法的不斷改進，感染患者的治愈率持續(xù)攀升，病死率不斷降低［3］。研究表明，COVID-19發(fā)病率持續(xù)增長是由SARS-CoV-2不斷突變產(chǎn)生新的高毒毒株引起的，且病毒刺突蛋白的突變與當(dāng)前疫苗接種的有效性降低有關(guān)［4］。SARS-CoV-2毒株發(fā)生突變后感染力和病毒載量顯著增加［5］。因此，篩選SARSCoV-2病毒與人體免疫細胞相關(guān)的關(guān)鍵分子生物標(biāo)志物對確定免疫治療的潛在靶標(biāo)至關(guān)重要。

基因表達綜合（gene expression omnibus，GEO）數(shù)據(jù)庫是目前最大的基因組測序數(shù)據(jù)庫之一，包含來自COVID-19患者免疫細胞的詳細基因組測序信息［6］。加權(quán)基因共表達網(wǎng)絡(luò)分析（weighted gene co-expression network analysis，WGCNA）是一種廣泛使用的生物信息學(xué)工具，可用于識別具有高協(xié)同變異的基因集［7］。GSE152418數(shù)據(jù)集［8］是基于SARS-CoV-2病毒感染與機體免疫系統(tǒng)之間相互作用研究較少的背景下開展的，通過對外周血單核細胞的轉(zhuǎn)錄組學(xué)分析探究COVID-19患者的免疫反應(yīng)情況。本研究從GEO數(shù)據(jù)庫下載GSE152418數(shù)據(jù)集，綜合運用生物信息學(xué)方法篩選SARS-CoV-2病毒感染潛在的關(guān)鍵分子生物標(biāo)志物并分析其免疫浸潤特征，旨在為開發(fā)新的免疫治療的潛在靶標(biāo)提供依據(jù)。

1 材料與方法

1.1 資料來源

從GEO數(shù)據(jù)庫下載GSE152418數(shù)據(jù)集（https：//www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE152418），測序平臺為GPL24676 Illumina NovaSeq 6000（Homo sapiens）。數(shù)據(jù)集共包含34例樣本，樣本來源于外周血單核細胞，其中COVID-19患者17例，健康對照17例。樣本數(shù)據(jù)由美國亞特蘭大市耶克斯國家靈長類動物研究中心提供。

1.2 差異表達分析

對數(shù)據(jù)矩陣中低表達基因進行過濾后，用“DESeq2”包［9］進行差異分析。高表達基因確定標(biāo)準(zhǔn)為log2FC≥1，P＜0.05；低表達基因確定標(biāo)準(zhǔn)為log2FC≤-1，P＜0.05。保存差 異基因（differential genes，DEGs）用于后續(xù)分析，并用熱圖和火山圖進行可視化。

1.3 加權(quán)基因共表達網(wǎng)絡(luò)分析

用“WGCNA”包［10］進行WGCNA分析。①數(shù)據(jù)進行歸一化及對數(shù)化處理，選擇方差＞25%的基因進行分析；②棄除離群樣本；③確定最佳鄰接函數(shù)參數(shù)（軟閾值），構(gòu)建鄰接矩陣，后轉(zhuǎn)為拓撲重疊矩陣，計算基因間相異度；④使用動態(tài)剪切樹劃分基因模塊，合并相關(guān)性系數(shù)大于0.8（即相異性系數(shù)小于0.2）的模塊；⑤計算模塊與臨床特征的相關(guān)性。選擇與臨床特征最相關(guān)的模塊計算基因顯著性（gene significance，GS）和模塊隸屬度（module membership，MM）。最后，用韋恩圖提取模塊基因與DEGs的共同基因（common genes，CGs）用于后續(xù)分析。

1.4 基因功能及信號通路富集分析

用“ClusterProfiler”包［11］對1.3步 驟 得 到 的CGs進行GO功能及KEGG信號通路富集分析。

1.5 蛋白互作網(wǎng)絡(luò)分析

利用String數(shù)據(jù)庫（https：//string-db.org/）構(gòu)建蛋白 互 作（protein-protein interaction，PPI）網(wǎng)絡(luò)，最小交互分數(shù)設(shè)置為0.9。然后將結(jié)果在Cytoscape3.8.2中進行可視化，利用插件CytoHub-ba篩選關(guān)鍵基因（Hub）。

1.6 miRNA-轉(zhuǎn)錄因子-mRNA網(wǎng)絡(luò)分析

分 別 利 用miRTarBase（https：//mirtarbase.cuhk.edu.cn/～miRTarBase）、Starbase（https：//starbase.sysu.edu.cn/starbase2）和Targetscan（https：//www.targetscan.org/vert_72）預(yù)測調(diào)控Hub基因的miRNA。利用韋恩圖取3個數(shù)據(jù)庫預(yù)測的共同miRNA，以提高預(yù)測的準(zhǔn)確性。利用Enrichr（https：//maayanlab.cloud/Enrichr）數(shù)據(jù)庫預(yù)測調(diào)控Hub基因的轉(zhuǎn)錄因子（TF）。再用Cytoscape3.8.2構(gòu)建可視化miRNA-TF-mRNA調(diào)控網(wǎng)絡(luò)。

1.7 免疫浸潤特征分析

用CIBERSORT算法［12］分析樣本中浸潤的免疫細胞構(gòu)成比例。用“Vioplot”包［13］分析COVID-19與健康對照樣本間免疫細胞的水平及差異。用“Corrplot”包［12］繪制相關(guān)性熱圖。

2 結(jié)果與分析

2.1 差異表達基因分析結(jié)果

DEGs的可視化見圖1，共得到2 049個DEGs，其中上調(diào)及下調(diào)基因分別有1 873個及176個。

圖1 差異基因的可視化圖Fig.1 Visualization of differential genes

2.2 WGCNA結(jié)果

共得到方差大于25%的基因4 743個。設(shè)置閾值h=120，剔除2個異常樣本，保留聚類1中的32個樣本。選擇最佳軟閾值β=4（R2=0.9）使基因表達關(guān)系符合無尺度網(wǎng)絡(luò)（圖2A）。通過動態(tài)剪切樹劃分基因模塊，設(shè)置每個模塊的最低基因數(shù)為50，合并相近的模塊（即相異性系數(shù)小于0.2）（圖2B）。共得到7個模塊，其中“土耳其藍色”模塊基因與SARS-CoV-2感染這一臨床特征相關(guān)性最高（r=0.91，P＜0.001）（圖2C）?！巴炼渌{色”模塊中MM和GS之間呈顯著正相關(guān)（r=0.96，P＜0.001）（圖2D）。韋恩圖結(jié)果共得到766個CGs（圖2E）。

2.3 CGs功能及信號通路富集分析結(jié)果

功能及信號通路富集結(jié)果（表1～2）顯示，CGs主要參與有絲分裂、微管結(jié)合、陽離子通道活性及卵母細胞減數(shù)分裂、細胞衰老、心肌病等。

表1 GO功能注釋及富集分析Table 1 GO functional annotation and enrichment analysis

2.4 PPI網(wǎng)絡(luò)分析結(jié)果

如圖3所示，利用String數(shù)據(jù)庫及Cytoscape3.8.2軟件構(gòu)建CGs的PPI網(wǎng)絡(luò)圖，然后利用CytoHubba插件篩選到Top10的Hub基因分別為CDK1、BUB1、CCNA2、CDC20、KIF11、BUB1B、CDCA8、TOP2A、CCNB2、KIF20A。

圖3 CGs的PPI網(wǎng)絡(luò)圖Fig.3 PPI network diagram of CGs

2.5 miRNA-TF-mRNA調(diào)控網(wǎng)絡(luò)分析

通過miRTarBase、Starbase、Targetscan及Enrichr數(shù)據(jù)庫分別預(yù)測調(diào)控Hub基因的miRNA和TF后，在Cytoscape3.8.2軟件中構(gòu)建miRNA-TFmRNA調(diào)控網(wǎng)絡(luò)，調(diào)控網(wǎng)絡(luò)中包含51個miRNA、5個TF和10個mRNA，其相互作用關(guān)系見圖4。

圖4 Hub基因的miRNA-TF-mRNA調(diào)控網(wǎng)絡(luò)圖Fig.4 miRNA-TF-mRNA regulatory network of Hub genes

2.6 免疫細胞浸潤分析

根據(jù)CIBERSORT算法預(yù)測COVID-19患者和健康對照組間的免疫細胞浸潤水平。樣本中22種免疫細胞的相對百分比見圖5。與健康對照組比較，COVID-19患者幼稚B細胞、嗜酸性粒細胞浸潤水平顯著降低（P＜0.05），漿細胞、活化肥大細胞浸潤水平顯著升高（P＜0.05）（圖6）。相關(guān)性分析結(jié)果顯示，活化肥大細胞與嗜酸性粒細胞（r=0.83，P＜0.05），靜息肥大細胞與單核細胞（r=0.77，P＜0.05）呈顯著正相關(guān)。靜息自然殺傷細胞與靜息肥大細胞（r=-0.69，P＜0.05）及單核細胞（r=-0.68，P＜0.05）呈顯著負相關(guān)（圖7）。

圖5 免疫細胞相對百分比圖Fig.5 Relative percentage of immune cells

圖6 兩組樣本間免疫細胞浸潤水平比較圖Fig.6 Comparison of immune cell infiltration levels between the two groups of samples

圖7 免疫浸潤細胞亞型間相關(guān)性分析圖Fig.7 Correlation analysis of immune infiltrating cell subtypes

表2 KEGG信號通路富集分析Table 2 KEGG signaling pathway enrichment analysis

3 討論

COVID-19已在全球多個國家和地區(qū)造成大流行。據(jù)世界衛(wèi)生組織（World Health Organization，WHO）官網(wǎng)統(tǒng)計數(shù)據(jù)顯示，截至2022年9月全球確診COVID-19病例已超過6.03億例，死亡人數(shù)超過648萬，但其感染的分子機制仍不完全明確［14］。本研究對包含17例COVID-19患者和17例健康對照樣本的高通量測序數(shù)據(jù)集進行分析，共篩選出2 049個DEGs，其中1 873個基因表達上調(diào)，176個基因表達下調(diào)。通過WGCNA分析得到7個模塊，其中“土耳其藍色”模塊與COVID-19相關(guān)性最高。利用韋恩圖得到766個CGs。WGCNA分析更側(cè)重于識別整個模塊中功能相似的基因，而不是單個基因的差異表達［15］。通過富集分析發(fā)現(xiàn)，CGs主要參與有絲分裂、微管結(jié)合、陽離子通道活性及卵母細胞減數(shù)分裂、細胞衰老、心肌病等，富集結(jié)果的差異性也可能解釋了COVID-19臨床表現(xiàn)和患者預(yù)后的異質(zhì)性［16］。

本研究進一步利用String數(shù)據(jù)庫對CGs進行PPI網(wǎng)絡(luò)分析，篩選到的前10位Hub基因分別為CDK1、BUB1、CCNA2、CDC20、KIF11、BUB1B、CDCA8、TOP2A、CCNB2、KIF20A。Hahn等［17］研究發(fā)現(xiàn)，細胞周期蛋白依賴性激酶1（CDK1）與抑制SARS-CoV-2病毒復(fù)制有關(guān)。Agrawal等［18］也發(fā)現(xiàn)BUB1是抗SARS-CoV-2治療的潛在靶點。Kim等［19］研究表明，Bcl-2抑制劑ABT-737介導(dǎo)的Bcl-2抑制與細胞周期蛋白A2（CCNA2）和細胞周期蛋白B1（CCNB1）的低表達有關(guān)，而Bcl-2抑制劑對SARS-CoV-2病毒的活性有抑制作用。此外，構(gòu)建的miRNA-TF-mRNA調(diào)控網(wǎng)絡(luò)中調(diào)控Hub基因的miRNA 51個、TF 5個。與多數(shù)病毒一樣，SARSCoV-2病毒RNA與宿主發(fā)生相互作用最關(guān)鍵的步驟是通過miRNA靶向調(diào)控RNA干擾宿主基因表達［20］。此外，5個TF蛋白是Hub基因的關(guān)鍵轉(zhuǎn)錄調(diào)節(jié)因子，與免疫缺陷及多種癌癥有關(guān)［21］。

免疫應(yīng)答在COVID-19發(fā)病機制中的關(guān)鍵作用越來越受到關(guān)注，包括固有免疫和適應(yīng)性免疫細胞浸潤特征［22］。本研究利用CIBERSORT算法評估了22種免疫細胞在樣本中的浸潤水平［12］。結(jié)果表明，COVID-19患者較健康對照樣本幼稚B細胞、嗜酸性粒細胞浸潤水平顯著降低，漿細胞、活化肥大細胞浸潤水平顯著升高。研究表明，B細胞的功能狀態(tài)及分化情況與COVID-19患者的預(yù)后有關(guān)［23］。Kuri-Cervantes等［24］研究發(fā)現(xiàn)，與健康對照及輕中癥COVID-19患者比較，重癥患者的外周血B細胞水平降低，但漿細胞水平升高，表明SARS-CoV-2病毒感染導(dǎo)致機體B細胞數(shù)量減少且分化異常，這可能也是COVID-19病情程度不一的原因之一。Zhang等［25］研究發(fā)現(xiàn)，COVID-19住院患者中大多數(shù)嗜酸性粒細胞低于正常水平，這有望作為診斷時的參考指標(biāo)之一。Nagashima等［26］研究發(fā)現(xiàn)，COVID-19患者中存在更多的活化肥大細胞，這些細胞與影響患者病情嚴重程度的血管通透性過高、水腫和彌漫性肺泡損傷事件直接相關(guān)。Sun等［27］利用GSE152418數(shù)據(jù)集及多中心腎移植受者隊列數(shù)據(jù)發(fā)現(xiàn)，受COVID-19影響，腎移植受者血液轉(zhuǎn)錄組結(jié)果顯示T細胞和適應(yīng)性免疫激活顯著減少。

本研究也具有一定的局限性：第一，本研究樣本數(shù)據(jù)是從GEO數(shù)據(jù)庫中獲取的，數(shù)據(jù)受到樣本量、臨床信息缺乏等因素影響，使得研究結(jié)論的可靠性受限；第二，對于免疫浸潤特征的分析尚不能證實免疫細胞與miRNA、TF、mRNA等標(biāo)志物之間的相互作用關(guān)系。因此，今后應(yīng)開展體內(nèi)外試驗研究，進一步驗證本研究結(jié)論的準(zhǔn)確性。

綜上所述，本研究通過WGCNA及PPI網(wǎng)絡(luò)分析篩選出10個Hub基因，并預(yù)測到調(diào)控Hub基因的5個TF及51個miRNA，且COVID-19患者與健康對照的免疫浸潤特征存在顯著差異，這些免疫細胞相關(guān)的分子標(biāo)志物可能作為COVID-19免疫治療的潛在靶標(biāo)。