瘢痕疙瘩為傷口皮膚結締組織過度增生引起的良性皮膚腫瘤。精準預測創傷者瘢痕疙瘩患病風險并及時做出早期診斷，對深度管理瘢痕疙瘩和控制其發展具有重大意義。本研究對高通量基因表達綜合（GEO）數據庫中的4個瘢痕疙瘩數據集進行分析，篩選出瘢痕疙瘩的診斷標志物，并建立列線圖預測模型。首先，通過加權基因共表達網絡分析（WGCNA）、差異表達分析和蛋白質互作網絡中心性算法，篩選出37個核心蛋白質編碼基因。隨后，利用最小絕對值收斂和選擇算子（LASSO）以及支持向量機?遞歸特征消除（SVM-RFE）兩種機器學習算法，從中篩選出4個最具預測能力的瘢痕疙瘩診斷標志物，分別為肝細胞生長因子（HGF）、多配體蛋白聚糖4（SDC4）、外核苷酸焦磷酸酶/磷酸二酯酶2（ENPP2）和Rho家族三磷酸鳥苷酶3（RND3），并通過單基因的基因集富集分析（GSEA）探索可能涉及的生物途徑。最后，對診斷標志物進行單因素與多因素邏輯回歸分析，并構建列線圖預測模型。經內外部驗證發現，該模型校準曲線貼近理想曲線，決策曲線優于其他策略，接受者操作特征曲線下面積高于對照模型（最佳截斷值為0.588），表明該模型具有較高的校準度、臨床收益率以及預測能力，有望為臨床診斷提供有效先期手段。
引用本文： 李政宇, 田保華, 梁海霞. 基于加權基因共表達網絡分析和機器學習的瘢痕疙瘩列線圖預測模型. 生物醫學工程學雜志, 2023, 40(4): 725-735. doi: 10.7507/1001-5515.202212048 復制
0 引言
瘢痕疙瘩是一種良性皮膚纖維增生性腫瘤[1-2]，其主要特征是損傷愈合后，局部皮膚過度生長，形成瘢痕疙瘩[3]。瘢痕疙瘩的發生率在4.5%～16.0%之間，且可持續增長多年，很少會自然消退[4]。盡管研究人員對瘢痕疙瘩的發病機制進行了一些探索，但目前仍缺乏明確的病理診斷標準和根治方案[5]。目前，較為有效的治療方式是盡早判斷瘢痕疙瘩的發病可能，在創傷后盡快干預，并進行深度的瘢痕疙瘩管理[6]，可有效延緩或控制瘢痕疙瘩的無限發展。在臨床實踐中，瘢痕疙瘩的診斷主要依賴于醫生的主觀判斷[7]，若患者錯誤進行手術切除治療，可能會導致更大范圍的瘢痕疙瘩形成[8]。因此，為了能盡早診斷瘢痕疙瘩，降低誤診率，避免瘢痕疙瘩的過度發展，急需建立一種有效的預測模型。
列線圖預測模型的設計初衷是為了對數學函數或公式進行近似的圖形化計算，通過圖形方式表達多因素邏輯回歸結果[9]。列線圖預測模型已廣泛應用于各種疾病的臨床風險預測，例如前列腺癌和糖尿病患者術中皮膚壓力性損傷的風險預測[10-11]。目前，關于瘢痕疙瘩發病預測的研究在國內外都相對有限，相關的預測模型也尚未建立。近年來，許多研究開始探索基因作為臨床預測因素的可能性[12-13]。同時，生物信息學技術也在飛速發展，眾多研究者嘗試運用各種方法，如差異表達分析、最小絕對值收斂和選擇算子（least absolute shrinkage and selection operator，LASSO）、支持向量機?遞歸特征消除（support vector machine-recursive feature elimination，SVM-RFE）[14-16] 、加權基因共表達網絡分析（weighted gene co-expression network analysis，WGCNA）以及蛋白質互作網絡中心性算法等[17]，盡可能地挖掘疾病的潛在生物標志物。在此過程中，Bi等[18]通過差異表達分析與蛋白質互作網絡尋找到關鍵的蛋白質集群。隨后，Li等[19]依據WGCNA、差異表達分析與蛋白質互作網絡中心性算法篩選出重要蛋白質的編碼基因。同時，Yin等[20]利用差異表達分析結合機器學習算法（SVM-RFE和LASSO）挑選出瘢痕疙瘩的特征基因。盡管上述方法已展開相關應用，但遺憾的是，瘢痕疙瘩的預測特征均未與瘢痕疙瘩疾病進行關聯，也未建立基于上述方法的列線圖預測模型。
本文借助上述方法，為瘢痕疙瘩的預測因素提出了一種創新的深層基因信息挖掘模式。首先，結合WGCNA、差異表達分析和蛋白質互作網絡中心性算法，分析基因的瘢痕疙瘩關聯性、組間表達水平以及它們在蛋白質互作網絡中的重要性，從而確定核心蛋白質編碼基因。隨后，利用機器學習算法（SVM-RFE和LASSO），識別出最具預測能力的診斷標志物，并對其進行單基因的基因集富集分析（gene set enrichment analysis，GSEA）。然后，通過多因素邏輯回歸分析，建立了列線圖預測模型。從內部和外部兩方面對模型進行驗證，評估它的準確性、校準度以及臨床效用，并將其與缺失單項算法的對照模型進行比對。如圖1研究流程圖所示，其成果或將為臨床上瘢痕疙瘩的診斷和治療提供新的思路。
圖1
研究流程圖
Figure1.
Flowchart of the research
1 方法
1.1 數據獲取和差異表達分析
美國國立生物技術信息中心的高通量基因表達綜合（gene expression omnibus，GEO）數據庫（網址為：http: //www.ncbi.nlm.nih.gov/geo/）是一個全球性公開資源庫。該庫收錄并整理了世界各地研究人員上傳的高通量基因組數據，并向公眾提供免費下載和使用的服務。在GEO數據庫中檢索并下載了瘢痕疙瘩相關的GEO系列（GEO series，GSE）數據，包括GSE121618[21]（5例瘢痕疙瘩和6例正常皮膚）和GSE145725[22]（9例瘢痕疙瘩和10例正常皮膚）。首先，利用數據處理可視化軟件R 4.2.0（Lucent Inc，美國）的替代變量分析（sva）R程序包消除數據集的批次效應。然后，使用微陣列數據線性模型（limma）R程序包[23]，以|log2倍數變化| >1和檢驗水準為0.05（即P < 0.05表示差異具有統計學意義）為標準，對數據集基因進行篩選。最后，利用 繪圖語法（ggplot2）R程序包，將這些基因展現為火山圖。
1.2 加權基因共表達網絡分析篩選關鍵基因
根據數據準確性和樣本數量的需求，選擇GSE145725數據集進行WGCNA，并設置標準差閾值為0.1，剔除波動過小基因。使用WGCNA R程序包對樣本進行聚類，并檢測數據集是否存在離群樣本[24]。在建立無尺度網絡的過程中，將無尺度網絡模型指數R2>0.85作為標準，計算出軟閾值，并構建拓撲重疊矩陣。用模塊的特征向量來計算模塊關聯度，并用層次聚類方法將相似模塊合并，設置剪切線高度為0.65。然后，繪制了模塊?特征關聯圖，以分析模塊特征向量與疾病的相關性以及它們之間差異是否具有統計學意義。通過比較基因顯著性（gene significance，GS）和模塊成員關系（module membership，MM）的中位數，篩選出了高GS與MM模塊基因。最后，將這些基因與差異表達基因相交，得到瘢痕疙瘩的關鍵基因。
1.3 構建蛋白互作網絡并篩選核心蛋白質編碼基因
使用預測蛋白質?蛋白質相互作用的相鄰基因重復實例搜索工具（search tool for recurring instances of neighboring genes，STRING）數據庫（網址為：https://string-db.org/）[25]構建關鍵基因的蛋白質互作網絡。接著，利用可視化分子相互作用網絡軟件Cytoscape 3.8.2（國際開源開發者聯盟，美國）[26]的網絡中心性分析插件（cytoNCA）對蛋白質互作網絡的介數中心性和度中心性進行計算，選擇兩種算法排名前50的交集作為核心蛋白質編碼基因。
1.4 診斷標志物的鑒定與使用單基因基因集富集分析探索其生物途經
利用廣義線性模型（glmnet）R程序包[27]對數據集進行LASSO回歸分析，包括模型擬合、數據預測、10折交叉驗證和提取模型系數。通過選取使誤差最小的非負調節參數λ，提取出模型系數不為0的預測特征。接下來，使用分類與回歸訓練（caret）R程序包 [28]進行SVM?RFE分析，通過反向特征選擇輔助函數、10折交叉驗證以及帶有徑向基核函數的支持向量機識別最優特征集。兩種機器算法的結果相交，即得到瘢痕疙瘩的診斷標志物。接著，進行單基因GSEA。根據單基因表達量的中位數，將疾病組分為低表達和高表達兩個亞組，并對兩個亞組進行差異表達分析。然后，基于倍數變化對數據集基因進行排序，利用生物信息學富集分析（clusterProfiler）R程序包基于京都基因和基因組百科全書（kyoto encyclopedia of genes and genomes，KEGG）數據庫進行GSEA，并繪制山脊圖來展示診斷標志物可能涉及的生物途經。
1.5 列線圖預測模型的構建與驗證
為了增加樣本數量，整合GSE121618、GSE145725、GSE92566[29]和GSE44270[30]（共35例瘢痕疙瘩和33例正常皮膚）四個數據集。對這些數據集消除批次效應和進行標準化處理。利用caret R程序包進行分層比例抽樣，將其分為訓練組和驗證組。在訓練集數據上，進行診斷標志物的單因素和多因素邏輯回歸分析。接著，使用回歸模型策略（rms）R程序包建立列線圖預測模型，并在訓練組數據上進行了內部驗證。繪制校準圖，以評價模型的校準度。接下來，使用風險模型決策分析（rmda）R程序包，對模型進行決策曲線分析（decision curve analysis，DCA），從而評估模型的臨床收益。該過程需繪制接受者操作特征（receiver operating characteristic curve，ROC）曲線，故使用部分ROC（partial ROC，pROC）R程序包，比較該模型與缺失單項算法對照模型的優劣。然后，利用報告ROC（reportROC）R程序包，進行約登指數計算，確定模型的最佳截斷值。最后，利用DCA和ROC曲線，對驗證組數據進行分析，以便評估模型在外部驗證中的表現。
1.6 統計分析
采用R軟件 4.2.0和統計處理軟件SPSS 13.0（IBM Corp.，美國）進行統計分析。首先，進行了正態分布檢驗，結果顯示所得數據符合正態分布。對于瘢痕疙瘩和正常皮膚的組間差異，采用獨立樣本t檢驗進行比較。使用皮爾遜相關系數構建了加權基因共表達網絡。對診斷標志物進行了單因素與多因素回歸分析，并使用χ2檢驗計算了優勢比和95%置信區間。為了驗證模型的性能，繪制了校準曲線以觀察預測結果與實際情況的擬合程度，并使用霍斯默·萊梅肖（Hosmer-Lemeshow）擬合優度檢驗評價模型的校準度。最后，使用約登指數確定模型的最佳截斷值，并通過繪制ROC曲線和DCA評價模型的準確性和臨床效用。P<0.05，表示差異具有統計學意義。
2 實驗結果與討論
2.1 差異表達分析和加權基因共表達網絡分析篩選的瘢痕疙瘩關鍵基因
首先，將GSE121618和GSE145725兩個數據集整合，并進行差異表達分析。如圖2差異表達基因的火山圖所示，在正常皮膚和瘢痕疙瘩之間存在247個差異表達基因，其中包括109個上調（紅點）和138個下調基因（藍點）。
圖2
差異表達分析和WGCNA的結果以及瘢痕疙瘩的關鍵基因
Figure2.
Results from differential expression analysis and WGCNA, and key genes in keloids
隨后，針對GSE145725數據集進行了WGCNA，以識別出與瘢痕疙瘩疾病相關的基因。依據最佳軟閾值β=12構建了無尺度網絡，并進行了基因的層次聚類。如圖2基因樹圖和動態樹剪切圖所示，共識別出了60個模塊。在合并相似模塊并保留剪切樹高度低于0.65的模塊后，最終得到了5個模塊。如圖2模塊?特征關聯圖所示，深綠松石和綠色模塊的模塊特征基因與瘢痕疙瘩疾病相關（相關性分別為0.88、?0.59，P<0.01）。深綠松石模塊包含3 191個基因，綠色模塊則包含1 926個基因。此外，深綠松石和綠色模塊的GS與MM之間也存在相關性（相關性系數r分別為0.96、?0.75，P<0.01），如圖2深綠松石模塊和綠色模塊的GS與MM圖所示。在疾病相關模塊中，共1 840個基因的GS和MM值大于模塊的中位數，包括深綠松石模塊的1 323個基因和綠色模塊的517個基因。最后，通過將差異表達基因和高GS與MM模塊基因取交集，如圖2瘢痕疙瘩關鍵基因的韋恩圖所示，得出了136個與瘢痕疙瘩疾病密切相關的關鍵基因。
2.2 診斷標志物的鑒定及單基因基因集富集分析結果
為了篩選出蛋白質互作網絡中的核心蛋白質編碼基因，本文基于STRING數據庫構建了關鍵基因的蛋白質互作網絡（節點數：136，邊數：779，平均節點度：11.5，平均局部聚類系數：0.305），如圖3關鍵基因的蛋白質互作網絡圖所示。使用度中心性和介數中心性作為篩選指標，取兩者的交集，以鑒定在蛋白質互作網絡中具有信息傳遞能力且處于中心地位的核心蛋白質。通過篩選，得到了介數中心性和度中心性排名前50的基因的交集，共有37個核心蛋白質編碼基因，如圖3 核心蛋白質編碼基因的韋恩圖所示。接著，采用了兩種機器學習算法，即SVM-RFE和LASSO，來計算對瘢痕疙瘩最具預測能力的診斷標志物。如圖3 SVM-RFE交叉驗證圖和LASSO回歸變量篩選圖所示，以最小化均方根誤差為目標進行交叉驗證，在SVM-RFE算法中篩選出了10個基因，分別是肝細胞生長因子（hepatocyte growth factor，HGF）、多配體蛋白聚糖4（syndecan-4，SDC4）、外核苷酸焦磷酸酶/磷酸二酯酶2（ectonucleotide pyrophosphatase/phosphodiesterase 2，ENPP2）、Rho家族三磷酸鳥苷酶3（Rho family guanosine triphosphatase 3，RND3）、第1類D型肌動蛋白、C3肉毒桿菌毒素底物2、成對框基因9、肌肉結節蛋白同源序列1、前列腺素D2合成酶、腎上腺素受體β2。選擇最小化10折交叉驗證誤差的最佳非負調節參數λ，LASSO算法篩選出了6個基因，包括ENPP2、HGF、SDC4、RND3、同源框B7、易洛魁族同源框蛋白5。最后，綜合LASSO和SVM-RFE算法的結果，成功地篩選出了4個最具預測能力的基因，分別為HGF、SDC4、ENPP2和RND3，如圖3瘢痕疙瘩診斷標志物的韋恩圖所示。
圖3
瘢痕疙瘩診斷標志物篩選
Figure3.
Screening for diagnostic markers of keloids
為了更深入地理解診斷標志物在瘢痕疙瘩中的作用，依據診斷標志物的表達量中位數將瘢痕疙瘩樣本分為兩組進行單基因GSEA。如圖4所示，在ENPP2高表達亞組中，發現核苷酸結合寡聚化結構域（nucleotide-binding oligomerization domain，NOD）樣信號通路、過氧化物酶體增殖物激活受體（peroxisome proliferator-activated receptor，PPAR）信號通路、兩面神激酶?信號傳導及轉錄激活蛋白（Janus kinase-signal transducer and activator of transcription，JAK?STAT）信號通路和核因子活化B細胞κ輕鏈增強子（nuclear factor kappa-light-chain-enhancer of activated B cells，NF-kappa B）信號通路等相關通路的富集。而在ENPP2低表達亞組中，N?聚糖的生物合成過程更為活躍。ENPP2可以產生磷脂酸，它通過同源G蛋白偶聯受體，影響多種生物學功能，包括應激上皮細胞產生促炎信號、激活轉化生長因子信號和積累成纖維細胞[31]。通過核糖核酸干擾（RNA interference，RNAi）的小鼠實驗，研究人員發現ENPP2受到無翼/整合（wingless/integrated，Wnt）信號蛋白/β?連環蛋白信號通路的調控，這是調節肌細胞分化的核心機制之一[32]。此外，ENPP2還可通過N?聚糖的生物合成通路，改變細胞外基質的組成和結構[33]。
圖4
單基因GSEA結果的山脊圖
Figure4.
Ridge plot of single-gene GSEA results
對于HGF高表達亞組，觀察到JAK?STAT信號通路、奇異（Toll）樣受體信號通路、叉頭盒O蛋白（forkhead box O，FoxO）信號通路、NOD樣信號通路以及NF-kappa B信號通路等相關通路的富集。HGF的通路主要集中在TGF-β、腫瘤壞死因子（tumor necrosis factor，TNF）和趨化因子等信號通路，失調可能導致成纖維細胞的活化和細胞外基質蛋白的過量產生[34]。HGF在瘢痕疙瘩的發展過程中主要通過抑制膠原合成、調控細胞外基質的周轉以及增強基質金屬蛋白酶的表達來發揮作用[35-36]。
在RND3高表達亞組中，發現白細胞介素17信號通路、NF-kappa B信號通路、轉化生長因子?β信號通路等相關通路活性增強。而在RND3低表達亞組中，單磷酸腺苷酸活化蛋白激酶（adenosine monophosphate-activated protein kinase，AMPK）信號通路、環磷酸鳥苷（cyclic guanosine monophosphate，cGMP）?蛋白激酶G（protein kinase G，PKG）信號通路和血管平滑肌收縮顯得更為活躍。RND3能夠調節張力纖維的收縮和肌動蛋白絲的伸長，這些過程是控制細胞骨架的重組以及細胞遷移的關鍵環節[37-38]。此外，RND3被發現能夠降低特發性肺纖維化中纖連蛋白、I型膠原蛋白和α-平滑肌肌動蛋白的表達[39]。
在SDC4高表達亞組中，黏液素型O?聚糖生物合成、細胞周期、視黃酸誘導基因蛋白I（retinoic acid-inducible gene-I，RIG?I）樣受體信號通路、細胞黏附分子、細胞外基質?受體相互作用、黏著斑等生物學功能更為活躍。而在SDC4的低表達亞組中發現，TNF信號通路和白細胞介素17信號通路的富集。SDC4是一種細胞表面糖蛋白，它通過與多個配基結合，以共受體的形式來調節細胞黏附、遷移、信號轉導以及細胞增殖等生物過程[40-41]。SDC4也可能間接影響細胞周期進程，從而調控細胞增殖[42]。這些診斷標志物在免疫反應、細胞增殖、分化、遷移以及細胞外基質的組成和結構調控等多個生物過程和信號通路中發揮作用。
2.3 基于診斷標志物的列線圖預測模型
瘢痕疙瘩是皮膚損傷后傷口異常愈合的極端表現[43]，無論接受何種形式的治療，都面臨極高的復發率[1]。目前尚未發現針對瘢痕疙瘩易感人群的生物標志物，這導致手術后發病風險增大[44]。研究發現，瘢痕疙瘩患者比普通人更容易患上癌癥，尤其是皮膚癌[45]。此外，具有瘢痕疙瘩家族史的人群，其患病風險明顯增加且病情更嚴重[46]。不同瘢痕疙瘩患者的易感位點各不相同，這表明瘢痕疙瘩可能符合寡基因遺傳模式[47]。盡管對瘢痕疙瘩的臨床診斷主要依賴病理檢查，但手術操作可能加速瘢痕疙瘩的進展[48]。因此，構建基于診斷標志物的列線圖預測模型可以為醫生的診斷提供輔助依據，減少侵入性診斷技術帶來的危害，具有較高的臨床應用價值。同時，不同人種或不同地區的瘢痕疙瘩發生率存在差異，甚至同一種族在不同地區也存在發病率的差異[4, 6, 49]。因此，制定適應各個區域人群的早期瘢痕疙瘩診斷和治療策略將對管理瘢痕疙瘩起到積極的作用。
在本研究中，采用訓練組數據，以診斷標志物為自變量，病理檢查結果是否為瘢痕疙瘩為因變量，進行了單因素和多因素邏輯回歸分析。通過分析，得到了四個變量的回歸分析參數，包括回歸系數、優勢比、置信區間和P值，詳細參數如表1所示。單因素回歸分析結果表明，ENPP2是瘢痕疙瘩發生的危險因素，而HGF、RND3和SDC4則表現為保護因素，這一發現與多因素回歸分析的結果保持一致。依據多因素回歸分析結果，構建了一個線性邏輯模型，如式（1）所示：
表1
訓練組各變量的單因素與多因素回歸分析
Table1.
Single factor and multi-factor logistic regression analysis of each variable in the training group
 |
其中，式（1）是邏輯回歸模型的對數形式，即對數幾率（logit）形式，主要用于計算事件發生和不發生的比值的自然對數。該式包含了截距和各個預測變量的回歸系數（截距：1.86，ENPP2：2.203，HGF：?1.301，RND3：?0.339，SDC4：?1.667）。
將式（1）的左式表示為變量y，采用logit變換將線性邏輯模型y轉化為概率值θ，如式（2）所示：
 |
其中，e 是自然對數的底數，約為2.71828。概率值θ表示瘢痕疙瘩發生的概率，線性邏輯模型y是預測變量（在本例中，這些變量為ENPP2、HGF、RND3和 SDC4）的加權總和。當預測變量的值已知，可以通過模型計算出疾病發生的概率θ。
如圖5訓練組和驗證組的ROC圖所示，比較了四個模型（模型1：本文模型；模型2：WGCNA + 蛋白質互作網絡中心性算法+機器學習；模型3：差異表達分析 + WGCNA + 機器學習；模型4：差異表達分析 + 蛋白質互作網絡中心性算法 + 機器學習）的曲線下面積。對照實驗結果表明，任何單個環節的缺失都會降低模型的準確性。相比之下，本文模型展現出較高的診斷能力。如圖5列線圖預測模型所示，它不僅提供了預測事件發生的概率，還闡明了預測變量影響事件發生概率的方式和程度。例如，某創傷初診患者，其樣本數據中診斷標志物RND3、SDC4、ENPP2和HGF的表達量log2值分別為5.16、5.19、4.85和1.12。通過將這些值投射到頂部標尺，得到每個變量的對應分數，將它們相加得到總分數167。然后再將總分數向下投射，即可預測該患者病理診斷為瘢痕疙瘩的概率為66.7%。通過計算約登指數，發現模型的最佳截斷值為0.588；在此閾值下，模型的特異度為58.8%，敏感度為100.0%。這意味著當預測疾病概率大于58.8%時，建議進行病理檢查；反之，應進行主動監測。此后，對列線圖預測模型的校準曲線進行了Hosmer-Lemeshow擬合優度檢驗，如圖5列線圖預測模型的校準圖所示，結果表明校準曲線與理想曲線的差異不具有統計學意義（χ2 = 5.368，P = 0.718），該模型的預測概率與實際概率接近，具有較高的校準度。此外，如圖5訓練組和驗證組的DCA圖所示，不論是在訓練組還是驗證組中，列線圖預測模型的臨床凈收益都優于“單診斷標志物預測”、“全部處理”以及“不作處理”策略，顯示出良好的臨床收益率。經過內部和外部兩方面的驗證，該模型展現了較強的適應性，可用于瘢痕疙瘩的風險預測。
圖5
列線圖預測模型建立與驗證
Figure5.
Nomogram model establishment and validation
3 結論
本文提出了一種新穎的深層基因信息挖掘模式，旨在尋找瘢痕疙瘩的預測因素。該模式運用了WGCNA、差異表達分析、蛋白質互作網絡中心性算法以及機器學習，并基于診斷標志物建立了瘢痕疙瘩的列線圖預測模型。與傳統依賴于臨床因素的列線圖預測模型相比，本文的方法結合了新興的生物信息學技術，順應精準醫療的發展趨勢，并為瘢痕疙瘩預測研究提供了新的思路。然而，該方法仍存在改進空間，未來可考慮結合本地醫院的臨床數據以提高瘢痕疙瘩風險預測的準確性，并通過擴大樣本規模來進一步驗證其效用。此外，還需要進一步探究本文鑒定的診斷標志物在瘢痕疙瘩的分子機制和生物功能方面的作用。
重要聲明
利益沖突聲明：本文全體作者均聲明不存在利益沖突。
作者貢獻聲明：李政宇主要負責數據處理、數據分析、論文撰寫與修改；田保華主要負責數據分析指導、協助論文撰寫、修改文章關鍵內容；梁海霞主要負責協助論文撰寫、論文審閱修訂并提供基金支持。
0 引言
瘢痕疙瘩是一種良性皮膚纖維增生性腫瘤[1-2]，其主要特征是損傷愈合后，局部皮膚過度生長，形成瘢痕疙瘩[3]。瘢痕疙瘩的發生率在4.5%～16.0%之間，且可持續增長多年，很少會自然消退[4]。盡管研究人員對瘢痕疙瘩的發病機制進行了一些探索，但目前仍缺乏明確的病理診斷標準和根治方案[5]。目前，較為有效的治療方式是盡早判斷瘢痕疙瘩的發病可能，在創傷后盡快干預，并進行深度的瘢痕疙瘩管理[6]，可有效延緩或控制瘢痕疙瘩的無限發展。在臨床實踐中，瘢痕疙瘩的診斷主要依賴于醫生的主觀判斷[7]，若患者錯誤進行手術切除治療，可能會導致更大范圍的瘢痕疙瘩形成[8]。因此，為了能盡早診斷瘢痕疙瘩，降低誤診率，避免瘢痕疙瘩的過度發展，急需建立一種有效的預測模型。
列線圖預測模型的設計初衷是為了對數學函數或公式進行近似的圖形化計算，通過圖形方式表達多因素邏輯回歸結果[9]。列線圖預測模型已廣泛應用于各種疾病的臨床風險預測，例如前列腺癌和糖尿病患者術中皮膚壓力性損傷的風險預測[10-11]。目前，關于瘢痕疙瘩發病預測的研究在國內外都相對有限，相關的預測模型也尚未建立。近年來，許多研究開始探索基因作為臨床預測因素的可能性[12-13]。同時，生物信息學技術也在飛速發展，眾多研究者嘗試運用各種方法，如差異表達分析、最小絕對值收斂和選擇算子（least absolute shrinkage and selection operator，LASSO）、支持向量機?遞歸特征消除（support vector machine-recursive feature elimination，SVM-RFE）[14-16] 、加權基因共表達網絡分析（weighted gene co-expression network analysis，WGCNA）以及蛋白質互作網絡中心性算法等[17]，盡可能地挖掘疾病的潛在生物標志物。在此過程中，Bi等[18]通過差異表達分析與蛋白質互作網絡尋找到關鍵的蛋白質集群。隨后，Li等[19]依據WGCNA、差異表達分析與蛋白質互作網絡中心性算法篩選出重要蛋白質的編碼基因。同時，Yin等[20]利用差異表達分析結合機器學習算法（SVM-RFE和LASSO）挑選出瘢痕疙瘩的特征基因。盡管上述方法已展開相關應用，但遺憾的是，瘢痕疙瘩的預測特征均未與瘢痕疙瘩疾病進行關聯，也未建立基于上述方法的列線圖預測模型。
本文借助上述方法，為瘢痕疙瘩的預測因素提出了一種創新的深層基因信息挖掘模式。首先，結合WGCNA、差異表達分析和蛋白質互作網絡中心性算法，分析基因的瘢痕疙瘩關聯性、組間表達水平以及它們在蛋白質互作網絡中的重要性，從而確定核心蛋白質編碼基因。隨后，利用機器學習算法（SVM-RFE和LASSO），識別出最具預測能力的診斷標志物，并對其進行單基因的基因集富集分析（gene set enrichment analysis，GSEA）。然后，通過多因素邏輯回歸分析，建立了列線圖預測模型。從內部和外部兩方面對模型進行驗證，評估它的準確性、校準度以及臨床效用，并將其與缺失單項算法的對照模型進行比對。如圖1研究流程圖所示，其成果或將為臨床上瘢痕疙瘩的診斷和治療提供新的思路。
圖1
研究流程圖
Figure1.
Flowchart of the research
1 方法
1.1 數據獲取和差異表達分析
美國國立生物技術信息中心的高通量基因表達綜合（gene expression omnibus，GEO）數據庫（網址為：http: //www.ncbi.nlm.nih.gov/geo/）是一個全球性公開資源庫。該庫收錄并整理了世界各地研究人員上傳的高通量基因組數據，并向公眾提供免費下載和使用的服務。在GEO數據庫中檢索并下載了瘢痕疙瘩相關的GEO系列（GEO series，GSE）數據，包括GSE121618[21]（5例瘢痕疙瘩和6例正常皮膚）和GSE145725[22]（9例瘢痕疙瘩和10例正常皮膚）。首先，利用數據處理可視化軟件R 4.2.0（Lucent Inc，美國）的替代變量分析（sva）R程序包消除數據集的批次效應。然后，使用微陣列數據線性模型（limma）R程序包[23]，以|log2倍數變化| >1和檢驗水準為0.05（即P < 0.05表示差異具有統計學意義）為標準，對數據集基因進行篩選。最后，利用 繪圖語法（ggplot2）R程序包，將這些基因展現為火山圖。
1.2 加權基因共表達網絡分析篩選關鍵基因
根據數據準確性和樣本數量的需求，選擇GSE145725數據集進行WGCNA，并設置標準差閾值為0.1，剔除波動過小基因。使用WGCNA R程序包對樣本進行聚類，并檢測數據集是否存在離群樣本[24]。在建立無尺度網絡的過程中，將無尺度網絡模型指數R2>0.85作為標準，計算出軟閾值，并構建拓撲重疊矩陣。用模塊的特征向量來計算模塊關聯度，并用層次聚類方法將相似模塊合并，設置剪切線高度為0.65。然后，繪制了模塊?特征關聯圖，以分析模塊特征向量與疾病的相關性以及它們之間差異是否具有統計學意義。通過比較基因顯著性（gene significance，GS）和模塊成員關系（module membership，MM）的中位數，篩選出了高GS與MM模塊基因。最后，將這些基因與差異表達基因相交，得到瘢痕疙瘩的關鍵基因。
1.3 構建蛋白互作網絡并篩選核心蛋白質編碼基因
使用預測蛋白質?蛋白質相互作用的相鄰基因重復實例搜索工具（search tool for recurring instances of neighboring genes，STRING）數據庫（網址為：https://string-db.org/）[25]構建關鍵基因的蛋白質互作網絡。接著，利用可視化分子相互作用網絡軟件Cytoscape 3.8.2（國際開源開發者聯盟，美國）[26]的網絡中心性分析插件（cytoNCA）對蛋白質互作網絡的介數中心性和度中心性進行計算，選擇兩種算法排名前50的交集作為核心蛋白質編碼基因。
1.4 診斷標志物的鑒定與使用單基因基因集富集分析探索其生物途經
利用廣義線性模型（glmnet）R程序包[27]對數據集進行LASSO回歸分析，包括模型擬合、數據預測、10折交叉驗證和提取模型系數。通過選取使誤差最小的非負調節參數λ，提取出模型系數不為0的預測特征。接下來，使用分類與回歸訓練（caret）R程序包 [28]進行SVM?RFE分析，通過反向特征選擇輔助函數、10折交叉驗證以及帶有徑向基核函數的支持向量機識別最優特征集。兩種機器算法的結果相交，即得到瘢痕疙瘩的診斷標志物。接著，進行單基因GSEA。根據單基因表達量的中位數，將疾病組分為低表達和高表達兩個亞組，并對兩個亞組進行差異表達分析。然后，基于倍數變化對數據集基因進行排序，利用生物信息學富集分析（clusterProfiler）R程序包基于京都基因和基因組百科全書（kyoto encyclopedia of genes and genomes，KEGG）數據庫進行GSEA，并繪制山脊圖來展示診斷標志物可能涉及的生物途經。
1.5 列線圖預測模型的構建與驗證
為了增加樣本數量，整合GSE121618、GSE145725、GSE92566[29]和GSE44270[30]（共35例瘢痕疙瘩和33例正常皮膚）四個數據集。對這些數據集消除批次效應和進行標準化處理。利用caret R程序包進行分層比例抽樣，將其分為訓練組和驗證組。在訓練集數據上，進行診斷標志物的單因素和多因素邏輯回歸分析。接著，使用回歸模型策略（rms）R程序包建立列線圖預測模型，并在訓練組數據上進行了內部驗證。繪制校準圖，以評價模型的校準度。接下來，使用風險模型決策分析（rmda）R程序包，對模型進行決策曲線分析（decision curve analysis，DCA），從而評估模型的臨床收益。該過程需繪制接受者操作特征（receiver operating characteristic curve，ROC）曲線，故使用部分ROC（partial ROC，pROC）R程序包，比較該模型與缺失單項算法對照模型的優劣。然后，利用報告ROC（reportROC）R程序包，進行約登指數計算，確定模型的最佳截斷值。最后，利用DCA和ROC曲線，對驗證組數據進行分析，以便評估模型在外部驗證中的表現。
1.6 統計分析
采用R軟件 4.2.0和統計處理軟件SPSS 13.0（IBM Corp.，美國）進行統計分析。首先，進行了正態分布檢驗，結果顯示所得數據符合正態分布。對于瘢痕疙瘩和正常皮膚的組間差異，采用獨立樣本t檢驗進行比較。使用皮爾遜相關系數構建了加權基因共表達網絡。對診斷標志物進行了單因素與多因素回歸分析，并使用χ2檢驗計算了優勢比和95%置信區間。為了驗證模型的性能，繪制了校準曲線以觀察預測結果與實際情況的擬合程度，并使用霍斯默·萊梅肖（Hosmer-Lemeshow）擬合優度檢驗評價模型的校準度。最后，使用約登指數確定模型的最佳截斷值，并通過繪制ROC曲線和DCA評價模型的準確性和臨床效用。P<0.05，表示差異具有統計學意義。
2 實驗結果與討論
2.1 差異表達分析和加權基因共表達網絡分析篩選的瘢痕疙瘩關鍵基因
首先，將GSE121618和GSE145725兩個數據集整合，并進行差異表達分析。如圖2差異表達基因的火山圖所示，在正常皮膚和瘢痕疙瘩之間存在247個差異表達基因，其中包括109個上調（紅點）和138個下調基因（藍點）。
圖2
差異表達分析和WGCNA的結果以及瘢痕疙瘩的關鍵基因
Figure2.
Results from differential expression analysis and WGCNA, and key genes in keloids
隨后，針對GSE145725數據集進行了WGCNA，以識別出與瘢痕疙瘩疾病相關的基因。依據最佳軟閾值β=12構建了無尺度網絡，并進行了基因的層次聚類。如圖2基因樹圖和動態樹剪切圖所示，共識別出了60個模塊。在合并相似模塊并保留剪切樹高度低于0.65的模塊后，最終得到了5個模塊。如圖2模塊?特征關聯圖所示，深綠松石和綠色模塊的模塊特征基因與瘢痕疙瘩疾病相關（相關性分別為0.88、?0.59，P<0.01）。深綠松石模塊包含3 191個基因，綠色模塊則包含1 926個基因。此外，深綠松石和綠色模塊的GS與MM之間也存在相關性（相關性系數r分別為0.96、?0.75，P<0.01），如圖2深綠松石模塊和綠色模塊的GS與MM圖所示。在疾病相關模塊中，共1 840個基因的GS和MM值大于模塊的中位數，包括深綠松石模塊的1 323個基因和綠色模塊的517個基因。最后，通過將差異表達基因和高GS與MM模塊基因取交集，如圖2瘢痕疙瘩關鍵基因的韋恩圖所示，得出了136個與瘢痕疙瘩疾病密切相關的關鍵基因。
2.2 診斷標志物的鑒定及單基因基因集富集分析結果
為了篩選出蛋白質互作網絡中的核心蛋白質編碼基因，本文基于STRING數據庫構建了關鍵基因的蛋白質互作網絡（節點數：136，邊數：779，平均節點度：11.5，平均局部聚類系數：0.305），如圖3關鍵基因的蛋白質互作網絡圖所示。使用度中心性和介數中心性作為篩選指標，取兩者的交集，以鑒定在蛋白質互作網絡中具有信息傳遞能力且處于中心地位的核心蛋白質。通過篩選，得到了介數中心性和度中心性排名前50的基因的交集，共有37個核心蛋白質編碼基因，如圖3 核心蛋白質編碼基因的韋恩圖所示。接著，采用了兩種機器學習算法，即SVM-RFE和LASSO，來計算對瘢痕疙瘩最具預測能力的診斷標志物。如圖3 SVM-RFE交叉驗證圖和LASSO回歸變量篩選圖所示，以最小化均方根誤差為目標進行交叉驗證，在SVM-RFE算法中篩選出了10個基因，分別是肝細胞生長因子（hepatocyte growth factor，HGF）、多配體蛋白聚糖4（syndecan-4，SDC4）、外核苷酸焦磷酸酶/磷酸二酯酶2（ectonucleotide pyrophosphatase/phosphodiesterase 2，ENPP2）、Rho家族三磷酸鳥苷酶3（Rho family guanosine triphosphatase 3，RND3）、第1類D型肌動蛋白、C3肉毒桿菌毒素底物2、成對框基因9、肌肉結節蛋白同源序列1、前列腺素D2合成酶、腎上腺素受體β2。選擇最小化10折交叉驗證誤差的最佳非負調節參數λ，LASSO算法篩選出了6個基因，包括ENPP2、HGF、SDC4、RND3、同源框B7、易洛魁族同源框蛋白5。最后，綜合LASSO和SVM-RFE算法的結果，成功地篩選出了4個最具預測能力的基因，分別為HGF、SDC4、ENPP2和RND3，如圖3瘢痕疙瘩診斷標志物的韋恩圖所示。
圖3
瘢痕疙瘩診斷標志物篩選
Figure3.
Screening for diagnostic markers of keloids
為了更深入地理解診斷標志物在瘢痕疙瘩中的作用，依據診斷標志物的表達量中位數將瘢痕疙瘩樣本分為兩組進行單基因GSEA。如圖4所示，在ENPP2高表達亞組中，發現核苷酸結合寡聚化結構域（nucleotide-binding oligomerization domain，NOD）樣信號通路、過氧化物酶體增殖物激活受體（peroxisome proliferator-activated receptor，PPAR）信號通路、兩面神激酶?信號傳導及轉錄激活蛋白（Janus kinase-signal transducer and activator of transcription，JAK?STAT）信號通路和核因子活化B細胞κ輕鏈增強子（nuclear factor kappa-light-chain-enhancer of activated B cells，NF-kappa B）信號通路等相關通路的富集。而在ENPP2低表達亞組中，N?聚糖的生物合成過程更為活躍。ENPP2可以產生磷脂酸，它通過同源G蛋白偶聯受體，影響多種生物學功能，包括應激上皮細胞產生促炎信號、激活轉化生長因子信號和積累成纖維細胞[31]。通過核糖核酸干擾（RNA interference，RNAi）的小鼠實驗，研究人員發現ENPP2受到無翼/整合（wingless/integrated，Wnt）信號蛋白/β?連環蛋白信號通路的調控，這是調節肌細胞分化的核心機制之一[32]。此外，ENPP2還可通過N?聚糖的生物合成通路，改變細胞外基質的組成和結構[33]。
圖4
單基因GSEA結果的山脊圖
Figure4.
Ridge plot of single-gene GSEA results
對于HGF高表達亞組，觀察到JAK?STAT信號通路、奇異（Toll）樣受體信號通路、叉頭盒O蛋白（forkhead box O，FoxO）信號通路、NOD樣信號通路以及NF-kappa B信號通路等相關通路的富集。HGF的通路主要集中在TGF-β、腫瘤壞死因子（tumor necrosis factor，TNF）和趨化因子等信號通路，失調可能導致成纖維細胞的活化和細胞外基質蛋白的過量產生[34]。HGF在瘢痕疙瘩的發展過程中主要通過抑制膠原合成、調控細胞外基質的周轉以及增強基質金屬蛋白酶的表達來發揮作用[35-36]。
在RND3高表達亞組中，發現白細胞介素17信號通路、NF-kappa B信號通路、轉化生長因子?β信號通路等相關通路活性增強。而在RND3低表達亞組中，單磷酸腺苷酸活化蛋白激酶（adenosine monophosphate-activated protein kinase，AMPK）信號通路、環磷酸鳥苷（cyclic guanosine monophosphate，cGMP）?蛋白激酶G（protein kinase G，PKG）信號通路和血管平滑肌收縮顯得更為活躍。RND3能夠調節張力纖維的收縮和肌動蛋白絲的伸長，這些過程是控制細胞骨架的重組以及細胞遷移的關鍵環節[37-38]。此外，RND3被發現能夠降低特發性肺纖維化中纖連蛋白、I型膠原蛋白和α-平滑肌肌動蛋白的表達[39]。
在SDC4高表達亞組中，黏液素型O?聚糖生物合成、細胞周期、視黃酸誘導基因蛋白I（retinoic acid-inducible gene-I，RIG?I）樣受體信號通路、細胞黏附分子、細胞外基質?受體相互作用、黏著斑等生物學功能更為活躍。而在SDC4的低表達亞組中發現，TNF信號通路和白細胞介素17信號通路的富集。SDC4是一種細胞表面糖蛋白，它通過與多個配基結合，以共受體的形式來調節細胞黏附、遷移、信號轉導以及細胞增殖等生物過程[40-41]。SDC4也可能間接影響細胞周期進程，從而調控細胞增殖[42]。這些診斷標志物在免疫反應、細胞增殖、分化、遷移以及細胞外基質的組成和結構調控等多個生物過程和信號通路中發揮作用。
2.3 基于診斷標志物的列線圖預測模型
瘢痕疙瘩是皮膚損傷后傷口異常愈合的極端表現[43]，無論接受何種形式的治療，都面臨極高的復發率[1]。目前尚未發現針對瘢痕疙瘩易感人群的生物標志物，這導致手術后發病風險增大[44]。研究發現，瘢痕疙瘩患者比普通人更容易患上癌癥，尤其是皮膚癌[45]。此外，具有瘢痕疙瘩家族史的人群，其患病風險明顯增加且病情更嚴重[46]。不同瘢痕疙瘩患者的易感位點各不相同，這表明瘢痕疙瘩可能符合寡基因遺傳模式[47]。盡管對瘢痕疙瘩的臨床診斷主要依賴病理檢查，但手術操作可能加速瘢痕疙瘩的進展[48]。因此，構建基于診斷標志物的列線圖預測模型可以為醫生的診斷提供輔助依據，減少侵入性診斷技術帶來的危害，具有較高的臨床應用價值。同時，不同人種或不同地區的瘢痕疙瘩發生率存在差異，甚至同一種族在不同地區也存在發病率的差異[4, 6, 49]。因此，制定適應各個區域人群的早期瘢痕疙瘩診斷和治療策略將對管理瘢痕疙瘩起到積極的作用。
在本研究中，采用訓練組數據，以診斷標志物為自變量，病理檢查結果是否為瘢痕疙瘩為因變量，進行了單因素和多因素邏輯回歸分析。通過分析，得到了四個變量的回歸分析參數，包括回歸系數、優勢比、置信區間和P值，詳細參數如表1所示。單因素回歸分析結果表明，ENPP2是瘢痕疙瘩發生的危險因素，而HGF、RND3和SDC4則表現為保護因素，這一發現與多因素回歸分析的結果保持一致。依據多因素回歸分析結果，構建了一個線性邏輯模型，如式（1）所示：
表1
訓練組各變量的單因素與多因素回歸分析
Table1.
Single factor and multi-factor logistic regression analysis of each variable in the training group
|
其中，式（1）是邏輯回歸模型的對數形式，即對數幾率（logit）形式，主要用于計算事件發生和不發生的比值的自然對數。該式包含了截距和各個預測變量的回歸系數（截距：1.86，ENPP2：2.203，HGF：?1.301，RND3：?0.339，SDC4：?1.667）。
將式（1）的左式表示為變量y，采用logit變換將線性邏輯模型y轉化為概率值θ，如式（2）所示：
|
其中，e 是自然對數的底數，約為2.71828。概率值θ表示瘢痕疙瘩發生的概率，線性邏輯模型y是預測變量（在本例中，這些變量為ENPP2、HGF、RND3和 SDC4）的加權總和。當預測變量的值已知，可以通過模型計算出疾病發生的概率θ。
如圖5訓練組和驗證組的ROC圖所示，比較了四個模型（模型1：本文模型；模型2：WGCNA + 蛋白質互作網絡中心性算法+機器學習；模型3：差異表達分析 + WGCNA + 機器學習；模型4：差異表達分析 + 蛋白質互作網絡中心性算法 + 機器學習）的曲線下面積。對照實驗結果表明，任何單個環節的缺失都會降低模型的準確性。相比之下，本文模型展現出較高的診斷能力。如圖5列線圖預測模型所示，它不僅提供了預測事件發生的概率，還闡明了預測變量影響事件發生概率的方式和程度。例如，某創傷初診患者，其樣本數據中診斷標志物RND3、SDC4、ENPP2和HGF的表達量log2值分別為5.16、5.19、4.85和1.12。通過將這些值投射到頂部標尺，得到每個變量的對應分數，將它們相加得到總分數167。然后再將總分數向下投射，即可預測該患者病理診斷為瘢痕疙瘩的概率為66.7%。通過計算約登指數，發現模型的最佳截斷值為0.588；在此閾值下，模型的特異度為58.8%，敏感度為100.0%。這意味著當預測疾病概率大于58.8%時，建議進行病理檢查；反之，應進行主動監測。此后，對列線圖預測模型的校準曲線進行了Hosmer-Lemeshow擬合優度檢驗，如圖5列線圖預測模型的校準圖所示，結果表明校準曲線與理想曲線的差異不具有統計學意義（χ2 = 5.368，P = 0.718），該模型的預測概率與實際概率接近，具有較高的校準度。此外，如圖5訓練組和驗證組的DCA圖所示，不論是在訓練組還是驗證組中，列線圖預測模型的臨床凈收益都優于“單診斷標志物預測”、“全部處理”以及“不作處理”策略，顯示出良好的臨床收益率。經過內部和外部兩方面的驗證，該模型展現了較強的適應性，可用于瘢痕疙瘩的風險預測。
圖5
列線圖預測模型建立與驗證
Figure5.
Nomogram model establishment and validation
3 結論
本文提出了一種新穎的深層基因信息挖掘模式，旨在尋找瘢痕疙瘩的預測因素。該模式運用了WGCNA、差異表達分析、蛋白質互作網絡中心性算法以及機器學習，并基于診斷標志物建立了瘢痕疙瘩的列線圖預測模型。與傳統依賴于臨床因素的列線圖預測模型相比，本文的方法結合了新興的生物信息學技術，順應精準醫療的發展趨勢，并為瘢痕疙瘩預測研究提供了新的思路。然而，該方法仍存在改進空間，未來可考慮結合本地醫院的臨床數據以提高瘢痕疙瘩風險預測的準確性，并通過擴大樣本規模來進一步驗證其效用。此外，還需要進一步探究本文鑒定的診斷標志物在瘢痕疙瘩的分子機制和生物功能方面的作用。
重要聲明
利益沖突聲明：本文全體作者均聲明不存在利益沖突。
作者貢獻聲明：李政宇主要負責數據處理、數據分析、論文撰寫與修改；田保華主要負責數據分析指導、協助論文撰寫、修改文章關鍵內容；梁海霞主要負責協助論文撰寫、論文審閱修訂并提供基金支持。
