基于全基因組測序技術(shù)篩選兒童金黃色葡萄球菌的致病相關(guān)分子標(biāo)志物

陳建余 王旭麟 李文宇 陳敏琪 周俊立 姚振江 傅錦堅 葉小華

（1.廣東藥科大學(xué)公共衛(wèi)生學(xué)院，廣東廣州 510310；2.柳州市婦幼保健院，廣西柳州 545001）

金黃色葡萄球菌（Staphylococcus aureus，以下簡稱“金葡菌”）是兒童感染的重要致病菌之一［1］。金葡菌的多重耐藥率較高，且社區(qū)和醫(yī)院中的金葡菌交叉?zhèn)鞑ワL(fēng)險較高，對人類健康造成極大威脅，給全球帶來了沉重的疾病經(jīng)濟(jì)負(fù)擔(dān)，已成為國際上日益關(guān)注的臨床問題和公共衛(wèi)生問題［1］。既往研究表明，兒童人群中的金葡菌鼻腔定植現(xiàn)象較為嚴(yán)重且明顯高于成人，在嬰幼兒時期的定植率高達(dá)45%［2］；同時，兒童身體各方面尚未發(fā)育完善，對病原體的抵抗力較低，兒童金葡菌感染率明顯高于成人，且一旦感染則更容易發(fā)展為菌血癥、骨髓炎、壞死性肺炎等侵襲性疾?。?］，提示兒童是金葡菌定植和感染的重點(diǎn)高危人群。既往研究主要局限于研究醫(yī)院感染人群金葡菌的耐藥性、分子特征，而關(guān)于健康兒童鼻腔定植金葡菌的分子流行病學(xué)研究仍缺乏系統(tǒng)報告，且金葡菌致病相關(guān)標(biāo)志物仍未闡明。金葡菌感染菌株與定植菌株的分子分型具有一定的相似性，僅通過克隆群、多位點(diǎn)序列分型等傳統(tǒng)分子分型方法難以區(qū)分菌株間的微小遺傳變異，有必要通過高通量全基因組測序技術(shù)闡明兩類金葡菌間全基因組水平上的潛在遺傳差異。金葡菌的致病力強(qiáng)弱與其自身含有的毒力因子高度相關(guān)，其中與金葡菌感染有關(guān)的毒力因子包括腸毒素、溶血素、中毒休克綜合征毒素、殺白細(xì)胞毒素等。2020 年研究［4］表明，中毒休克綜合征毒素基因tsst-1、表皮剝脫毒素基因eta和etb的表達(dá)可增強(qiáng)金葡菌的毒力及致病力，fnbA、cna、sdrE、sej、eta、hlg和ica基因在感染菌株中比定植菌株更加普遍，因此這些毒力因子為揭示金葡菌的致病相關(guān)標(biāo)志物提供重要思路。隨機(jī)森林是一種基于決策樹的非參數(shù)模型，可以有效處理基因組數(shù)據(jù)的高維、高相關(guān)性等問題，其特點(diǎn)是可以對變量的重要性進(jìn)行排序和識別致病相關(guān)因子，為揭示金葡菌致病相關(guān)標(biāo)志物提供新策略［5］。本研究對512 株金葡菌（272株感染菌株和240株定植菌株）進(jìn)行全基因組測序分析，全面檢測金葡菌的耐藥基因和毒素基因，運(yùn)用隨機(jī)森林篩選金葡菌致病相關(guān)分子標(biāo)志物，為追溯高致病性金葡菌和開展精準(zhǔn)的靶向干預(yù)提供遺傳學(xué)證據(jù)。

1 資料與方法

1.1 研究對象

本研究采用橫斷面研究設(shè)計。（1）感染菌株來源：收集2014—2017 年在柳州市婦幼保健院就診的感染兒童臨床標(biāo)本。感染兒童的納入標(biāo)準(zhǔn)：年齡≤7 歲；家長簽署知情同意書；有臨床感染癥狀（如咳嗽、聽診異常、呼吸困難、發(fā)熱超過38℃等），且從感染部位（包括痰液、膿液、血液、支氣管肺泡灌洗液等）采集的臨床樣本中分離和鑒定出符合標(biāo)準(zhǔn)的陽性金葡菌，最終確定為金葡菌感染的患兒。排除標(biāo)準(zhǔn)：同一患兒的不同部位采集標(biāo)本，取主要發(fā)病部位的標(biāo)本；同一患兒對于同一疾病多次采樣時，只取首次標(biāo)本。感染兒童的信息采集包括：年齡、性別、病案號等基本信息；標(biāo)本來源、采樣時間、疾病診斷等臨床信息。（2）定植菌株來源：采用分層整群抽樣抽取2018 年4—6 月期間廣西柳州市6 所幼兒園，開展健康幼兒園兒童的鼻咽拭子采樣。健康兒童的納入標(biāo)準(zhǔn)：年齡≤7 歲的健康兒童，且家長簽署知情同意書。排除標(biāo)準(zhǔn)：調(diào)查前7 d 內(nèi)患有肺炎、支氣管炎等急性感染性疾病的兒童，以及患有血友病等易出血疾病導(dǎo)致不便采樣的兒童。共調(diào)查1 702名健康兒童，檢出277株金葡菌，其中240株成功進(jìn)行全基因組測序分析。本研究獲得柳州市婦幼保健院醫(yī)學(xué)倫理委員會（批件號：2018-085）和廣東藥科大學(xué)倫理委員會（批件號：2015-20）批準(zhǔn)。

1.2 菌株分離鑒定

菌株分離鑒定試驗(yàn)包括：甘露醇氯化鈉瓊脂培養(yǎng)基顯色反應(yīng)、革蘭染色、純化培養(yǎng)、溶血試驗(yàn)及血漿凝固酶試驗(yàn)，依次進(jìn)行上述試驗(yàn)，結(jié)果陽性則鑒定為金葡菌。

1.3 全基因組測序分析

經(jīng)菌液制備、試劑盒提取DNA、質(zhì)量控制后獲得金葡菌DNA，經(jīng)分光光度檢測儀檢測濃度和純度及格的基因組DNA 進(jìn)行全基因組測序，采用Illumina Hiseq 2000 測序儀對基因組DNA 進(jìn)行雙端測序，使用FastQC 進(jìn)行原始測序數(shù)據(jù)的質(zhì)控［6］，采用SPAdes 3.12.0 軟件進(jìn)行序列拼接。使用基因組流行病學(xué)中心（Center for Genomic Epidemiology,CGE） 網(wǎng)站（http://www.genomicepidemiology.org/）中的ResFinder 4.1 （https://cge.cbs.dtu.dk/services/ResFinder/）進(jìn)行序列比對，對15 類2 208 個耐藥基因進(jìn)行檢測［7］。使用CGE 中VirulenceFinder 2.0（https://cge.cbs.dtu.dk/services/VirulenceFinder/） 進(jìn)行序列比對獲得菌株的毒素基因，主要包括免疫逃逸基因、胞外酶基因、溶血毒素基因、腸毒素基因、表皮剝脫素基因、殺白細(xì)胞毒素基因、中毒休克綜合征毒素基因。以S. aureussubsp.aureusMRSA252（accession number: NC_002952）為參考基因組，使用Snippy 軟件（https://github.com/tseemann/snippy/）對菌株進(jìn)行核心單核苷酸多態(tài)性鑒定，使用Gubbins軟件去除基因重組或水平基因轉(zhuǎn)移相關(guān)區(qū)域［8］。使用Fasttree 軟件構(gòu)建基于核心基因組的系統(tǒng)進(jìn)化樹，采用廣義時間可逆替代（general time reversible, GTR）模型［9］，估計方法為極大似然估計法，Boostrap 值設(shè)置為100。使用Chiplot 網(wǎng)站（http://www.evolgenius.info/）對系統(tǒng)進(jìn)化樹進(jìn)行可視化和美化。

1.4 統(tǒng)計學(xué)分析

采用EpiData 3.0 軟件進(jìn)行數(shù)據(jù)的雙人雙錄入，并通過一致性檢驗(yàn)保證數(shù)據(jù)的準(zhǔn)確性。采用Stata 16.0 和R 4.1.2 軟件對數(shù)據(jù)進(jìn)行統(tǒng)計學(xué)分析。計數(shù)資料以頻數(shù)和百分率（%）表示，采用Pearsonχ2檢驗(yàn)或Fisher確切概率法比較金葡菌感染菌株與定植菌株的耐藥基因、毒素基因攜帶率，從而初步篩選致病相關(guān)因子。P＜0.05 為差異有統(tǒng)計學(xué)意義。采用單因素logistic 回歸或確切l(wèi)ogistic 回歸模型估計OR及95%CI。運(yùn)用機(jī)器學(xué)習(xí)方法中的隨機(jī)森林進(jìn)一步篩選重要的致病相關(guān)特征變量，模型以菌株類型（1=感染菌株，0=定植菌株）作為結(jié)局變量，以分子特征（毒素基因和耐藥基因）作為自變量。隨機(jī)森林分析中，通過randomForestSRC 程序包中隨機(jī)森林的變量捕獲法（variable selection using random forests, VSURF）對重要特征變量進(jìn)行篩選［10］；采用randomForest程序包對最終模型中的變量進(jìn)行重要性評分，特征變量的重要性以平均基尼指數(shù)減少量（mean decrease in the Gini, MDG）排序，從而判斷各個變量在模型中的重要性。采用組內(nèi)回代和十折交叉驗(yàn)證方法評價最終模型的預(yù)測效果，評價指標(biāo)包括正確率、靈敏度、特異度、陽性預(yù)測值、陰性預(yù)測值、受試者操作特征曲線（receiver operating characteristic curve, ROC曲線）。

2 結(jié)果

2.1 基本情況

有272株金葡菌（即感染菌株組）來源于醫(yī)院感染兒童，年齡為0～7 歲，中位數(shù)為2 個月，172株（63.2%）來源于男性、100株（36.8%）來源于女性。有240 株鼻腔定植金葡菌（即定植菌株組）來源于健康兒童，年齡為3～7 歲，中位數(shù)為5 歲，124 株（51.7%）來源于男性、116 株（48.3%）來源于女性。

2.2 定植菌株組與感染菌株組的毒素基因攜帶情況

金葡菌感染菌株組的seb、sep、splA、splB、splE、edinC、lukD、lukE、lukF-PV、lukS-PV、eta和etb基因攜帶率高于定植菌株組（均P＜0.05）；但感染菌株組的sec、sec3、seg、seh、sei、sel、sem、sen、seo和seu基因攜帶率低于定植菌株組（均P＜0.05），其他毒素基因攜帶率比較差異無統(tǒng)計學(xué)意義（P＞0.05）。見表1。

表1 不同類型金葡菌的毒素基因攜帶情況比較 ［株（%）］

2.3 定植菌株組與感染菌株組的耐藥基因攜帶情況

金葡菌感染菌株組的lnuA、aadD、tetK和dfrG基因攜帶率明顯高于定植菌株組（P＜0.05）。但感染菌株組的blaTEm-1A基因攜帶率低于定植菌株（P＜0.05）。其余耐藥基因攜帶率比較差異均無統(tǒng)計學(xué)意義（P＞0.05）。見表2。

表2 不同類型金葡菌的耐藥基因攜帶情況比較 ［株（%）］

2.4 隨機(jī)森林篩選金葡菌致病相關(guān)分子特征

基于核心單核苷酸多態(tài)性的系統(tǒng)進(jìn)化樹（圖1）表明：感染菌株與定植菌株在系統(tǒng)進(jìn)化樹上的分布沒有明顯差異，提示兩者可能具有相似的遺傳背景，需進(jìn)一步揭示致病相關(guān)分子特征。以菌株類型（1=感染菌株，0=定植菌株）為因變量，建立隨機(jī)森林模型篩選兒童金葡菌致病相關(guān)的重要特征變量（69 個），采用VSURF 法最終篩選出16 個特征變量。對變量篩選前后的模型進(jìn)行預(yù)測效果評價，結(jié)果表明：變量篩選前后預(yù)測模型的交叉驗(yàn)證正確率分別為69%、68%（表3）；變量篩選前模型的曲線下面積（area under the curve,AUC）稍高于篩選后（0.75 vs 0.70）（圖2A～B）。最終篩選出的16 個特征變量為腸毒素基因（sem、sep、ser、sea）、表皮剝脫毒素基因etb、胞外酶編碼基因splE、β-內(nèi)酰胺類耐藥基因（mecA、blaZ、blaTEm-1A）、大環(huán)內(nèi)酯類耐藥基因（ermA、ermB、lnuA）、氯霉素類耐藥基因cat(pC233)、氨基糖苷類耐藥基因［aph(3')-Ⅲ、ant(9)-Ⅰa、ant(6)-Ⅰa］；變量重要性排序結(jié)果顯示，前5個最重要的特征變量分別為sem、etb、splE、sep、ser（圖2C）。使用隨機(jī)森林模型進(jìn)一步對16 個重要致病相關(guān)標(biāo)志物進(jìn)行風(fēng)險預(yù)測，結(jié)果顯示：有10 個致病相關(guān)標(biāo)志物位于對角線上方，提示攜帶這些致病相關(guān)標(biāo)志物會增加金葡菌的致病風(fēng)險；而有6個致病相關(guān)標(biāo)志物位于對角線下方，提示攜帶這些致病標(biāo)志物會降低金葡菌的致病風(fēng)險（圖2D）。

圖1 512 株金葡菌的系統(tǒng)進(jìn)化樹 從內(nèi)環(huán)到外環(huán)的變量分別為菌株來源（紅色為感染菌株，藍(lán)色為定植菌株）、毒力基因數(shù)、耐藥基因數(shù)、16個致病相關(guān)特征變量（深紅表示有，淺紅表示無）。

圖2 隨機(jī)森林模型的擬合效果和風(fēng)險得分圖 A：變量篩選前模型的ROC曲線。B：變量篩選后模型的ROC曲線。C：最終模型中16個特征變量的重要性排序。D：最終模型中16個特征變量的風(fēng)險得分圖。

表3 基于分子特征的隨機(jī)森林模型預(yù)測效果評價

3 討論

關(guān)于金葡菌的耐藥機(jī)制研究較多，常見于以下幾類耐藥機(jī)制［11-14］：（1）酶降解抗生素，如β-內(nèi)酰胺酶；（2）產(chǎn)生拮抗物質(zhì)，減少抗生素的有效濃度；（3）改變細(xì)胞靶標(biāo)，使抗生素不能正常發(fā)揮作用；（4）主動外排泵作用。blaZ基因受blaR1-BlaI 系統(tǒng)的調(diào)節(jié)控制，blaR1 蛋白受到β-內(nèi)酰胺類抗生素的刺激后，導(dǎo)致抑制蛋白BlaI水解脫離結(jié)合位點(diǎn)，編碼產(chǎn)生β-內(nèi)酰胺酶，繼而破壞β-內(nèi)酰胺從而使金葡菌對青霉素類抗生素耐藥［11］。mecA編碼PBP2a 蛋白，產(chǎn)生過量的PBP 蛋白消耗β-內(nèi)酰胺類藥物，從而導(dǎo)致金葡菌對β-內(nèi)酰胺類藥物耐藥。既往研究［12］顯示，雖然blaTEm基因編碼的β-內(nèi)酰胺酶主要通過絲氨酸殘基滅活β-內(nèi)酰胺環(huán)產(chǎn)生耐藥性，但是blaTEm在多重耐藥菌株和敏感菌株中均被檢出且差異無統(tǒng)計學(xué)意義，提示blaTEm可能不是導(dǎo)致其廣泛耐藥的主要原因；值得注意的是，本研究的隨機(jī)森林預(yù)測模型揭示，blaTEm-1A基因是金葡菌致病相關(guān)標(biāo)志物，提示其可能影響菌株的致病能力。erm編碼核糖體甲基化酶催化金葡菌23S rRNA 發(fā)生甲基化反應(yīng)使核糖體靶位點(diǎn)發(fā)生變化，進(jìn)而阻止大環(huán)內(nèi)酯類藥物與核糖體靶位結(jié)合，從而導(dǎo)致金葡菌對大環(huán)內(nèi)酯類藥物耐藥［13］。lnuA基因編碼核苷酸轉(zhuǎn)移酶，從而介導(dǎo)對林可酰胺類抗生素的耐藥性［14］。

細(xì)菌的毒力因子在感染性疾病中起著關(guān)鍵作用，它是由多種毒素基因決定的。金葡菌具有數(shù)十種毒素因子，包括腸毒素、表皮剝脫毒素、胞外酶等。本研究隨機(jī)森林結(jié)果顯示，金葡菌致病相關(guān)的毒素基因有腸毒素基因（sem、sep、ser和sea）、表皮剝脫毒素基因（etb）、胞外酶編碼基因（splE）。葡萄球菌腸毒素在食源性中毒中最常見，它不僅可引起食物中毒，還可引起中毒性休克綜合征等嚴(yán)重的侵襲性疾病［15］。本研究中，感染菌株組sep基因的攜帶率明顯高于定植菌株組，但是sem基因的攜帶率明顯低于定植菌株組，提示感染菌株與定植菌株存在毒力差異，與中國臺灣的一項研究結(jié)果一致［16］。表皮剝脫毒素A 和表皮剝脫毒素B是導(dǎo)致大多數(shù)人類葡萄球菌燙傷樣皮膚綜合征的重要原因［17］；本研究中，感染菌株組較定植菌株組的eta（7.4% vs 2.9%）、etb（9.2% vs 2.5%）基因攜帶率高，以上結(jié)果提示表皮剝脫毒素與菌株的致病能力密切相關(guān)。splE編碼絲氨酸蛋白酶樣蛋白，可以在健康人，尤其是哮喘患者中誘導(dǎo)Ⅱ型超敏反應(yīng)［18］；西班牙一項研究提示，splE基因可作為金葡菌菌血癥的致病相關(guān)標(biāo)志物［19］。以上結(jié)果提示不同類型疾病可能存在致病相關(guān)標(biāo)志物的差異，因此后續(xù)研究可針對不同疾病類型進(jìn)一步探討潛在的特異性標(biāo)志物。

基因組數(shù)據(jù)具有“高維度、小樣本”的數(shù)據(jù)特點(diǎn)，且存在“非線性、高相關(guān)”的復(fù)雜關(guān)聯(lián)，使得傳統(tǒng)分析方法不再適用，這對基因組關(guān)聯(lián)分析提出了新挑戰(zhàn)。隨機(jī)森林是利用多棵決策樹對樣本進(jìn)行訓(xùn)練和預(yù)測的機(jī)器學(xué)習(xí)算法，能有效處理線性、非線性、具有交互作用的復(fù)雜數(shù)據(jù)，具有抗噪聲、防止過擬合、不受共線性影響的重要特征。隨機(jī)森林不僅具有良好的預(yù)測能力，且具有精度高、穩(wěn)定性好、易操作等優(yōu)點(diǎn)，因此已成為處理醫(yī)學(xué)領(lǐng)域基因組數(shù)據(jù)的精確機(jī)器學(xué)習(xí)方法之一。為了獲得更準(zhǔn)確的風(fēng)險預(yù)測模型，本研究以金葡菌的常見致病相關(guān)基因（32 個耐藥基因和37 個毒素基因）為預(yù)測變量、菌株類型（感染和定植菌株）為因變量建立隨機(jī)森林模型，最終模型的16個預(yù)測因子分別為腸毒素基因（sem、sep、ser和sea）、表皮剝脫毒素基因etb、胞外酶編碼基因splE、β-內(nèi)酰胺類耐藥基因（mecA、blaZ和blaTEm-1A）、大環(huán)內(nèi)酯類耐藥基因（ermA、ermB和lnuA）、氯霉素類耐藥基因cat(pC233)、氨基糖苷類耐藥基因［aph(3')-III、ant(9)-Ⅰa 和ant(6)-Ⅰa］，最終模型的交叉驗(yàn)證正確率為68%，AUC 為0.70，提示模型的擬合效果較好。隨機(jī)森林模型揭示：前5 個最重要的預(yù)測變量為sem、etb、splE、sep、ser；其中sep（OR=3.97）、etb（OR=3.95）、ser（OR=1.68）和splE（OR=1.68）明顯提高金葡菌的致病風(fēng)險，提示腸毒素、表皮剝脫毒素和胞外酶是金葡菌的重要致病相關(guān)標(biāo)志物。以上結(jié)果揭示這些致病相關(guān)分子標(biāo)志物有較大的潛力預(yù)測金葡菌致病株。

綜上，本研究使用高通量全基因組測序技術(shù)全面檢測金葡菌基因組的耐藥基因和毒素基因，采用隨機(jī)森林篩選出金葡菌的16 個致病相關(guān)標(biāo)志物（6 個毒素基因和10 個耐藥基因），且模型預(yù)測效果較優(yōu)，為追溯高致病性金葡菌和開展精準(zhǔn)的靶向干預(yù)提供遺傳學(xué)證據(jù)。本研究也存在一定局限性：（1）本研究提示致病相關(guān)特征變量與結(jié)局之間存在統(tǒng)計學(xué)關(guān)聯(lián)，但篩選的致病相關(guān)標(biāo)志物與疾病狀態(tài)之間不一定存在因果關(guān)系，且仍未能深入闡明潛在的生物學(xué)功能與作用機(jī)制，因此今后可以結(jié)合生物學(xué)功能富集分析、通路富集分析對其潛在作用機(jī)制進(jìn)行深入研究。（2）本研究于2014—2018 年間對金葡菌感染患兒和健康兒童進(jìn)行金葡菌采樣，菌株的采集時間較早；但是，本研究主要從基因組水平上探究金葡菌的致病相關(guān)標(biāo)志物，且既往研究表明菌株分子特征與菌株的采樣時間無關(guān)［20］。今后可進(jìn)一步開展多中心調(diào)查來增大樣本量和補(bǔ)充新菌株。

利益沖突聲明：所有作者均聲明不存在利益沖突。