引用本文： 徐賀龍, 李燕, 任程暉, 宋曉靜, 張磊, 李汛. 利用加權基因共表網絡識別肝細胞癌干細胞相關基因的研究. 中國普外基礎與臨床雜志, 2020, 27(5): 560-568. doi: 10.7507/1007-9424.202001067 復制
肝細胞癌（hepatocellular carcinoma，HCC）是最常見的原發性肝惡性腫瘤，居全球腫瘤發病率的第 6 位，腫瘤相關性死亡率位居第 3 位[1]。雖然早期肝癌根治性切除后的 5 年生存率可達 70%，但累積復發率卻高達 100%[2]。對于晚期肝癌患者，化療是目前重要的治療手段，但耐藥性一直是影響其療效的重要因素，因而晚期肝癌的預后仍然很差[3]。近年來，腫瘤干細胞（cancer stem cells，CSCs）理論為腫瘤的診斷和治療提供了新的思路。研究[4]證實，CSCs 具有自我更新、分化和腫瘤發生的顯著特征，可導致各種腫瘤的轉移、化療耐藥和復發。雖然 CSCs 可能在腫瘤的進展和治療中發揮關鍵作用，但其內在的干細胞特性在分子水平上的調控機制仍不清楚。以往的研究主要集中在單個生物標志物在 CSCs 中的應用。隨著測序技術的進步，使高分辨率的全基因組測序成為可能，這些資源構成了許多數據庫，其中在腫瘤研究領域的癌癥基因組圖譜（TCGA，https://portal.gdc.cancer.gov）能使我們對 33 種腫瘤類型的近 12 000 個樣本進行廣泛的分析。基于 mRNA 表達的腫瘤干細胞干性指數（mRNA stemness indices，mRNAsi）是樣本干細胞特性的量化表示（范圍在 0～1 之間）。該指數模型是由 Malta 等[5]基于祖細胞數據集的干細胞樣本、利用機器學習算法建立的（工作流程見
1 資料與方法
1.1 數據下載及整理
基于 TCGA 中經病理學檢查確診的 HCC 患者的基因表達數據（具體選擇參數如下：Project Id is TCGA-LIHC，and Data Category is Transcriptome Profiling，and Data Type is Gene Expression Quantification），找出 424 例樣本。其中，HCC 組織樣本 374 例，正常肝組織樣本 50 例。使用 Perl 語言（
1.2 mRNAsi 與病理學分級及預后的相關性
mRNAsi 是描述腫瘤細胞和干細胞之間相似程度的指數，可以認為是對 CSCs 干性的量化。通過 Perl 語言合并腳本，分別將病理學分級矩陣、生存信息矩陣與 mRNAsi 模型中的基因 ID 進行配對，整合成新的矩陣。以腫瘤樣本 mRNAsi 值的中位數分成高 mRNAsi 組和低 mRNAsi 組，使用 R 語言中的“survival”和“survminer”包進行組間差異的生存分析。使用 R 語言中的“beeswarm”包進行病理學分級相關性的可視化分析，采用 Kruskal-Wallis 檢驗分析各組間 mRNAsi 值的差異。
1.3 HCC 差異表達基因分析
FPKM（expected number of Fragments Per Kilobase of transcript sequence per Millions base pairs sequenced）值用來評估樣本的基因表達水平，刪除在 2 組樣本中表達水平均小于 1 的基因，重復基因取表達平均值。使用 R 語言中的“Limma”包來篩選 HCC 組織和正常肝組織樣本之間的差異表達基因，截取標準為：|log2FC|>1，P<0.05，其中 FC 為樣本表達量的差異倍數。
1.4 WGCNA 網絡構建
WGCNA 由“WGCNA”包構建。首先對過濾后的差異基因表達水平進行分析，利用函數 hclust 將樣本中的表達數據進行樣本聚類分析，刪除離散樣本（圖 1a）。再使用函數 pick soft threshold 選擇軟閾值，軟閾值是基于近似無標度網絡的一種準則，使構建的網絡符合冪律分布，這種處理方式強化了強相關，弱化了弱相關，使得相關性數值更符合無標度網絡特征，更具有生物學意義。設置相關系數 R2為 0.8，使用 pick soft threshold 函數自動篩選 8 為最佳軟閾值，同時保證了較高的連通性（圖 1b 和圖 1c）。通過對模塊的進一步分析，計算模塊間的差異，并構建模塊樹狀圖。設定模塊中最小基因數為 50，設置 0.25 為切割高度的閾值，將相似度較高的模塊合并。
圖1
示 WGCNA 網絡構建及軟閾值的自動選取
a：樣本聚類樹分析圖，縱坐標為樣本的樹枝高度，高度越高代表該樣本的基因表達量越大，對 374 例樣本進行聚類樹分析，刪除高度大于 10 000 的離群樣本；b：通過無尺度網絡指數確定軟閾值，
1.5 確定關鍵模塊，鑒定 hub 基因并篩選關鍵基因
為確定關鍵模塊，需計算模塊的顯著性。首先通過計算基因顯著性（gene significance，GS）來評估基因與樣本間的相關性。將 GS 定義為基因表達與 mRNAsi 進行線性回歸時 P 值的 log10 轉換值（GS=lgP）。將模塊顯著性定義為模塊中所有基因的平均 GS，從而計算出模塊與樣本性狀間的相關性，選擇相關性最高的模塊為關鍵模塊。確定關鍵模塊后，計算模塊內每個基因的 GS 和模塊隸屬度（module membership，MM），即模塊自身基因與基因表達譜的相關性，并根據 GS 和 MM 設置 hub 基因的篩選閾值。本研究中 hub 基因的截取標準設置為 GS>0.38，MM>0.80。利用在線數據庫 DAVID[7]（
1.6 關鍵基因差異分析及相關性分析
為進一步分析關鍵基因表達情況，使用 R 語言中的“ggpubr”包行差異分析，并基于基因表達水平計算關鍵基因和上游基因間的共表達關系，以確定這些關系在轉錄水平的強弱。使用 R 語言中的“corrplot”包計算基因間的 Pearson 相關性。
1.7 關鍵基因的生存分析
利用在線數據庫 Kaplan-Meier plotter[9]（
1.8 統計學方法
采用 R 軟件進行生信分析和統計學分析。檢驗水準 α=0.05。
2 結果
2.1 mRNAsi 與臨床病理學特征及預后的相關性
正常肝組織與 HCC 組織間 mRNAsi 值的差異有統計學意義（P<0.001），見圖 2a。本研究結果還表明，mRNAsi 與高病理學等級之間存在正相關關系（P=0.001），見圖 2b。Kaplan-Meier 生存曲線結果顯示，高 mRNAsi 組的死亡率高于低 mRNAsi 組，差異有統計學意義（P=0.006），見圖 2c。
圖2
示 mRNAsi 在 HCC 及正常肝組織中的差異及其與病理學分級和預后的相關性、HCC 組織差異基因的火山圖、WGCNA 基因聚類樹、關鍵模塊的顯著性及顯著性較高模塊的鑒定結果
a：正常肝組織（50 個樣品）和 HCC 組織（374 個樣品）之間 mRNAsi 值的差異；b：mRNAsi 值與 HCC 病理學分級相關性分析的箱線圖；c：Kaplan-Meier 生存曲線圖，可見高 mRNAsi 組的死亡率高于低 mRNAsi 組；d：HCC 差異表達基因火山圖，fdr 為假陽性率，綠色圓點代表下調基因，紅色圓點代表上調基因，黑色圓點代表無顯著差異基因；e：HCC 的共表達基因模塊聚類樹圖，橫坐標 22 個不同顏色模塊分別代表一個基因集模塊，縱坐標為基因占比，樹狀圖的分支對應于 22 個不同的基因模塊，樹狀圖上的每一片葉子都對應一個基因，相似的基因被聚類到相同顏色的模塊中；f：基因模塊和臨床特征之間的相關性熱圖，該模塊內數字為模塊中基因的平均 mRNAsi 得分，顏色越紅即與 mRNAsi 和 EREG-mRNAsi 正相關越強，顏色越藍即與 mRNAsi 和 EREG-mRNAsi 負相關越強，并注釋相應的
2.2 HCC 組織與正常肝組織差異表達基因篩選結果
由于人類基因組數據庫龐大，且部分基因表達不具有顯著差異，而正常肝組織的 mRNAsi 值與腫瘤組織明顯不同，故本研究首先進行基因清洗以篩選出差異基因。過濾得到差異表達倍數關系>1 且有統計學意義（P<0.05）的基因。最終篩選了 7 668 個差異基因，其中 7 273 個上調，395 個下調（圖 2d）。
2.3 WGCNA 關鍵模塊及 hub 基因的確定
使用 R 語言中的“WGCNA”包將表達模式相似的基因分配到一個模塊中，形成模塊樹狀圖（圖 2e），得到 22 個基因模塊。然后，計算 22 個模塊內基因與 mRNAsi 和表觀遺傳調控的 mRNAsi（EREG-mRNAsi）的相關性（GS），以及基因與模塊的相關性（MM）。結果表明，綠黃色模塊與 mRNAsi 呈現顯著負相關且相關系數最大（圖 2f），對截取的基因進一步分析后發現并無明顯 mRNAsi 相關性及 LSCS 特性（圖 2g）。而藍色模塊與 HCC 的 mRNAsi 的正相關性最高（r=0.28，P=0.004），表明藍色模塊中的基因可能與 LSCS 密切相關，從藍色模塊篩選 hub 基因（圖 2h），共得到 24 個 Hub 基因：CDCA8、SGOL1、KIF15、GINS1、CDC20、TTK、CKAP2L、ZWINT、NCAPG、EXO1、RAD51AP1、NDC80、CCNB1、MCM2、KIF20A、MCM10、RFC4、CHAF1A、EME1、SPAG5、CDC45、TEDC2、CDCA5 和 OIP5。
2.4 hub 基因的 GO 和 KEGG 通路富集分析
篩選藍色模塊中 GS>0.38，MM>0.80 的基因作為本研究的 hub 基因，并將篩選得到的 24 個 hub 基因進行 GO 和 KEGG 通路分析，以鑒定相關的生物學功能和信號傳導途徑。GO 分析的結果表明，這些基因主要參與的生物過程有：有絲分裂姐妹染色單體分離、核染色體分離、有絲分裂核分裂、染色體分離、DNA 構象變化等；主要的細胞組成有：染色體區域、染色體著絲粒區域、紡錘體等（圖 3a）。此外，KEGG 通路分析發現，這些基因集主要富集于細胞周期、DNA 錯配修復系統、卵母細胞減數分裂和 DNA 復制信號通路（圖 3b）。
圖3
示藍色模塊的 24 個 hub 基因富集分析柱狀圖、hub 基因的 PPI 網絡圖及基因鄰邊節點數，以及關鍵基因的表達差異箱線圖及基因相關性熱圖
a 和 b：GO 分析結果（a）和 KEGG 通路富集分析結果（b），BP 為生物學過程，CC 為細胞組分，直方圖顏色越紅表示統計學差異越大；c：hub 基因 PPI 網絡，圓圈為基因名稱，線條表示基因間的聯系，隱藏了網絡中斷開的節點；d：直方圖示 PPI 網絡中 24 個中樞基因的節點數，其中 2 個斷開節點基因被隱藏，橫坐標為該基因在 PPI 網絡中的鄰邊節點數；e：HCC 中 10 個關鍵基因的差異表達箱線圖，HCC 組織和正常肝組織比較，***
2.5 hub 基因 PPI 網絡的可視化分析
通過 STRING 在線數據庫獲得了 24 個 hub 基因編碼的 PPI 網絡（圖 3c），設置顯示最高置信度（最低互作得分>0.9）。結果顯示，PPI 網絡共由 24 個節點和 44 條邊組成，平均局部聚類度：0.695（P<0.001）。除 CKAP2L2 和 TEDC2 外，其余 22 個基因編碼的蛋白產物均參與了 PPI 網絡的構建，且有多個蛋白位于網絡的中心，與周圍蛋白互作關系密切。使用 R 語言統計結果，并繪制柱狀圖（圖 3d），結果提示，hub 基因之間存在廣泛而牢固的關系，本研究選擇圖 3d 中鄰接節點數最多的前 10 個基因作為關鍵基因，分別為 CCNB1、CDC20、CDCA8、NDC80、KIF20A、CDC45、KIF15、MCM2、TTK 和 NCAPG。
2.6 關鍵基因差異分析及基因共表達分析
使用 R 語言完成關鍵基因的表達分析，并繪制箱線圖（圖 3e），可見 10 個關鍵基因在 HCC 中均顯著高表達（P<0.001）。繪制其共表達相關性熱圖（圖 3f），結果表明基因間相關系數均>0.7，提示關鍵基因相關性較高，上下游作用聯系緊密。
2.7 在線 Kaplan-Meier 生存曲線分析及驗證
使用在線 Kaplan-Meier Plotter 數據庫，以 HCC 組織中關鍵基因表達量的中位數作為分界值將基因表達量分為低表達組和高表達組。共分析了 364 例 HCC 患者的 OS，結果如圖 4a–圖 4j 所示，10 個關鍵基因高表達組均具有不良預后，差異具有統計學意義（P<0.001）。
圖4
示以在線數據庫 Kaplan-Meier plotter 和 HCCDB18 數據集分析 10 個關鍵基因對 HCC 患者預后的影響
a–j：分別示在線數據庫 Kaplan-Meier plotter 分析的 10 個關鍵基因 KIF15（56992）、CDC45（8318）、TTK（7272）、KIF20A（10112）、NDC80（10403）、CDCA8（55143）、CCNB1（891）、MCM2（4171）、CDC20（991）和 NCAPG（64151），括號內數字為基因的數字 ID； k–t：分別示 HCCDB18 數據集分析的 10 個關鍵基因 KIF15、CDC45、TTK、KIF20A、NDC80、CDCA8、CCNB1、MCM2、CDC20和 NCAPG
進一步使用 HCC 綜合分子數據庫中的數據集 HCCDB18 對 HCC 患者的生存分析結果進行驗證。共下載 203 例 HCC 患者的預后數據，為了減少與疾病無關的死亡原因的影響，本數據集中將長于 5 年的生存時間默認為 5 年，并將相應患者的狀態設置為存活。結果如圖 4k-圖 4t 所示，10 個關鍵基因高表達組均與不良預后有關（P<0.001）。
3 討論
基于 WGCNA 的方法對惡性腫瘤的相關基因進行挖掘在近年來逐漸成為腫瘤研究的重要課題。晚期肝癌患者因失去手術機會更加依靠精確的分子靶向治療，基于此背景，如何在精準醫療及基因水平治療的理念和模式下尋找潛在的靶向治療標志物，制定個體化的治療方案顯得尤為重要。
CSCs 是惡性腫瘤無限增殖和復發的源泉，近年來已成為腫瘤領域的研究熱點。CSCs 與腫瘤細胞增殖、侵襲轉移、耐藥等多種惡性表型密切相關[11-12]。本研究中，首先下載了由 Malta 等提供的 HCC 的 mRNAsi，并基于 TCGA 數據庫中的 374 例 HCC 樣本及 50 例正常肝組織樣本進行 mRNAsi 差異的驗證，發現 HCC 樣本中 mRNAsi 值顯著高于正常肝組織樣本，mRNAsi 值和腫瘤病理學分級成正相關；生存曲線分析結果顯示，高 mRNAsi 組的存活率較低 mRNAsi 組顯著下降，這些結果均符合 CSCs 特性，與預期結果一致。基于此，本研究根據肝癌差異基因與 mRNAsi 的關系構建了 WGCNA 網絡，發現藍色模塊與 mRNAsi 的相關性最高，截取的 24 個 hub 基因主要與干細胞自我更新和增殖特征有關。WGCNA 利用分子間的表達相關系數來衡量他們的共表達關系，同一模塊中的分子表達模式相似，表達模式相似的分子可能參與同一生物學過程或通路，然而基因間的相互作用還需要依托蛋白質組學、表觀遺傳組學等綜合判斷，因此限于 WGCNA 的局限性，關鍵基因的挑選需要參考蛋白組學和轉錄組學的結果。本研究中 hub 基因蛋白互作關系緊密，選取 PPI 網絡中鄰邊節點數最多的 10 個基因作為關鍵基因進一步分析。結果關鍵基因在腫瘤組織中的表達量明顯高于正常組織，在轉錄水平上存在很強的共表達關系，且 Kaplan-Meier 生存曲線顯示關鍵基因高表達組的生存率顯著下降并在 HCCDB18 數據集中得到了驗證。既往也發現這些基因的差異表達與腫瘤狀態密切相關。
KIF20A 和 KIF15 是驅動蛋白家族成員，該家族由一組共享高度保守的運動域蛋白質組成[13]。有研究[14-16]表明，KIF20A 在紡錘體組裝和染色體分離中起到重要作用，其在多種腫瘤中高表達且與預后差有關。Shi 等[17]研究發現，KIF20A 是 Hh 信號傳導的重要下游靶基因，此外 Gli 家族的鋅指蛋白是 Shh 信號傳導的介質，被廣泛認為是胚胎癌的致癌基因；Gli2-KIF20A 軸對于 HCC 細胞的增殖、侵襲和腫瘤復發有著至關重要的作用。KIF15 的過表達增加了 CSCs 相關的生物標志物和干性相關的轉錄因子水平，過表達 KIF15 的細胞在連續高比例（1∶20）傳代過程中仍維持類器官分化能力，此外 KIF15 能促進 LSCS 自我更新，并通過 HCC 細胞中活性氧的失衡促進肝癌組織中 LSCS 的干性維持[18]。
GINS1 是一種由 SLD5、PSF1、PSF2 和 PSF3 組成的四聚體復合物，該復合物在 DNA 復制過程中起著至關重要的作用。本研究中 GINS 作為 hub 基因與關鍵基因間有著緊密的聯系，其過表達與肝癌的預后不良相關，這與 Lian 等[19]的研究一致。Nagahama 等[20]證明了 PSF1 是在干細胞和祖細胞中廣泛表達的 DNA 復制因子，具有致瘤和轉移特性。GINS1 四聚體聯合 MCM2-7 與 CDC45 形成 DNA 解旋酶復合物（CDC45-MCM2-7-GINS 復合體）參與 DNA 復制，其參與 DNA 解鏈的機制研究被認為是腫瘤藥物研發的新方向[21-22]。
CCNB1 是細胞周期蛋白 B1（Cyclin B1），是細胞周期家族的成員之一，對于控制細胞周期的 G2/M（有絲分裂）過渡至關重要。近年來，越來越多的研究表明，CCNB1 在肝癌、非小細胞肺癌、結直腸癌、乳腺癌、胃癌等多種腫瘤中均呈高表達[23-24]。最近的研究[25]表明，CCNB1 基因參與了結直腸癌 CSCs 細胞周期的調控。Lee 等[26]通過實驗證實，頭頸部鱗狀細胞癌中 CCNB1 的過表達引起胚胎干細胞關鍵蛋白表達異常，并出現 CSCs 特性，導致腫瘤復發和預后不良。此外，細胞周期蛋白家族中多個成員均被證實參與 CSCs 細胞周期的調控。
NDC80 編碼蛋白質的功能是組織和穩定微管-線粒體相互作用，該蛋白是染色體分離所必需的結構。Shen 等[27]通過綜合分析多個 GEO 數據集發現，NDC80 在 HCC 的發生和復發過程中起到核心作用。近期有研究[28]表明，NDC80 參與前列腺癌干細胞的干性維持，高表達的 NDC80 具有更強的腫瘤細胞成球性，并且參與了多西他賽耐藥性。這些報道均直接或間接證實了本研究的結論。
TTK 編碼具有磷酸化酪氨酸、絲氨酸和蘇氨酸能力的雙重特異性蛋白激酶，這種蛋白質對于有絲分裂過程中染色體的著絲粒對齊至關重要，并且是中心體復制所必需的。已經發現它是有絲分裂期間作用于染色體精確分離的關鍵檢查點蛋白，當該蛋白不能降解并產生過量的中心體而導致異常的有絲分裂紡錘體堆積時，可能會導致肝癌的發生[29]。關于 TTK 基因與 CSCs 的生物學功能的研究已有很多發現，如 TTK 基因在從食管癌細胞系分離的 CSCs 樣細胞群及通過乙醛脫氫酶 1（ALDH1）分選的多發性骨髓瘤干細胞中均過表達[30-31]。有研究[32]表明，TTK 基因的敲低可降低 CD44 陽性的三陰性乳腺癌中 CSCs 的數量。最新研究[33]發現，TTK 幾乎在所有具有 CSCs 特性的癌組織中過表達，小分子化合物 TTK-27f 可有效抑制其作為致癌蛋白激酶的功能，具有強大的類藥物特性并優先抑制 CSCs，因此可能是 CSCs 靶向治療的新策略。
NCAPG 編碼凝縮蛋白復合物的一個亞基，負責有絲分裂和減數分裂過程中染色體的凝結和穩定。研究[34-35]發現，NCAPG 通過磷脂酰肌醇 3 激酶（PI3K）/蛋白激酶 B（Akt）信號通路信號促進 HCC 的增殖，NCAPG 的高表達與肝癌患者整體生存率和無病生存率顯著相關，NCAPG 可預測手術切除后 HCC 的早期復發并作為肝癌患者的預后生物標志物。此外，一項基于生物信息分析的研究[36]發現，NCAPG 可作為膀胱癌干細胞特征的生物標志物，但是該基因在 CSCs 中的生物學功能及作用機制尚未闡明，鑒定 NCAPG 的蛋白質底物及其信號傳導通路將更好地了解該基因在 CSCs 中的關鍵作用。
本研究通過基因相關性分析發現 10 個關鍵基因之間高度相關，且蛋白互作關系緊密。除 KIF15 外，其余關鍵基因雖沒有直接證據表明與 LSCS 的特性有關，但以往的研究結果間接表明，關鍵基因可以促進 LSCS 的干性維持且在 HCC 的復發和化療耐藥性中起重要作用。本研究中 hub 基因的 GO 分析和 KEGG 通路分析結果表明，這些基因主要集中在與細胞有絲分裂相關的功能上并參與細胞周期和 DNA 復制，這也間接反映了這些基因與干細胞的自我更新和增殖特性有關。
綜上所述，本研究結果表明，mRNAsi 在 HCC 中起著重要作用，基于 WGCNA 篩選了與 LSCS 相關的 10 個關鍵基因可能在維持 LSCS 特性方面發揮重要作用，這些基因可能是抑制 LSCS 特性的治療靶點。然而，這些結論僅僅是基于 mRNAsi 模型的分析結論，因此還需進一步研究來驗證這些發現。
重要聲明
利益沖突聲明：本文全體作者閱讀并理解了《中國普外基礎與臨床雜志》的政策聲明，我們沒有相互競爭的利益。
作者貢獻聲明：徐賀龍負責研究設計、數據處理和論文撰寫；李燕參與研究設計；任程暉參與數據處理；宋曉靜和張磊參與論文撰寫；李汛負責論文審改。
肝細胞癌（hepatocellular carcinoma，HCC）是最常見的原發性肝惡性腫瘤，居全球腫瘤發病率的第 6 位，腫瘤相關性死亡率位居第 3 位[1]。雖然早期肝癌根治性切除后的 5 年生存率可達 70%，但累積復發率卻高達 100%[2]。對于晚期肝癌患者，化療是目前重要的治療手段，但耐藥性一直是影響其療效的重要因素，因而晚期肝癌的預后仍然很差[3]。近年來，腫瘤干細胞（cancer stem cells，CSCs）理論為腫瘤的診斷和治療提供了新的思路。研究[4]證實，CSCs 具有自我更新、分化和腫瘤發生的顯著特征，可導致各種腫瘤的轉移、化療耐藥和復發。雖然 CSCs 可能在腫瘤的進展和治療中發揮關鍵作用，但其內在的干細胞特性在分子水平上的調控機制仍不清楚。以往的研究主要集中在單個生物標志物在 CSCs 中的應用。隨著測序技術的進步，使高分辨率的全基因組測序成為可能，這些資源構成了許多數據庫，其中在腫瘤研究領域的癌癥基因組圖譜（TCGA，https://portal.gdc.cancer.gov）能使我們對 33 種腫瘤類型的近 12 000 個樣本進行廣泛的分析。基于 mRNA 表達的腫瘤干細胞干性指數（mRNA stemness indices，mRNAsi）是樣本干細胞特性的量化表示（范圍在 0～1 之間）。該指數模型是由 Malta 等[5]基于祖細胞數據集的干細胞樣本、利用機器學習算法建立的（工作流程見
1 資料與方法
1.1 數據下載及整理
基于 TCGA 中經病理學檢查確診的 HCC 患者的基因表達數據（具體選擇參數如下：Project Id is TCGA-LIHC，and Data Category is Transcriptome Profiling，and Data Type is Gene Expression Quantification），找出 424 例樣本。其中，HCC 組織樣本 374 例，正常肝組織樣本 50 例。使用 Perl 語言（
1.2 mRNAsi 與病理學分級及預后的相關性
mRNAsi 是描述腫瘤細胞和干細胞之間相似程度的指數，可以認為是對 CSCs 干性的量化。通過 Perl 語言合并腳本，分別將病理學分級矩陣、生存信息矩陣與 mRNAsi 模型中的基因 ID 進行配對，整合成新的矩陣。以腫瘤樣本 mRNAsi 值的中位數分成高 mRNAsi 組和低 mRNAsi 組，使用 R 語言中的“survival”和“survminer”包進行組間差異的生存分析。使用 R 語言中的“beeswarm”包進行病理學分級相關性的可視化分析，采用 Kruskal-Wallis 檢驗分析各組間 mRNAsi 值的差異。
1.3 HCC 差異表達基因分析
FPKM（expected number of Fragments Per Kilobase of transcript sequence per Millions base pairs sequenced）值用來評估樣本的基因表達水平，刪除在 2 組樣本中表達水平均小于 1 的基因，重復基因取表達平均值。使用 R 語言中的“Limma”包來篩選 HCC 組織和正常肝組織樣本之間的差異表達基因，截取標準為：|log2FC|>1，P<0.05，其中 FC 為樣本表達量的差異倍數。
1.4 WGCNA 網絡構建
WGCNA 由“WGCNA”包構建。首先對過濾后的差異基因表達水平進行分析，利用函數 hclust 將樣本中的表達數據進行樣本聚類分析，刪除離散樣本（圖 1a）。再使用函數 pick soft threshold 選擇軟閾值，軟閾值是基于近似無標度網絡的一種準則，使構建的網絡符合冪律分布，這種處理方式強化了強相關，弱化了弱相關，使得相關性數值更符合無標度網絡特征，更具有生物學意義。設置相關系數 R2為 0.8，使用 pick soft threshold 函數自動篩選 8 為最佳軟閾值，同時保證了較高的連通性（圖 1b 和圖 1c）。通過對模塊的進一步分析，計算模塊間的差異，并構建模塊樹狀圖。設定模塊中最小基因數為 50，設置 0.25 為切割高度的閾值，將相似度較高的模塊合并。
圖1
示 WGCNA 網絡構建及軟閾值的自動選取
a：樣本聚類樹分析圖，縱坐標為樣本的樹枝高度，高度越高代表該樣本的基因表達量越大，對 374 例樣本進行聚類樹分析，刪除高度大于 10 000 的離群樣本；b：通過無尺度網絡指數確定軟閾值，
1.5 確定關鍵模塊，鑒定 hub 基因并篩選關鍵基因
為確定關鍵模塊，需計算模塊的顯著性。首先通過計算基因顯著性（gene significance，GS）來評估基因與樣本間的相關性。將 GS 定義為基因表達與 mRNAsi 進行線性回歸時 P 值的 log10 轉換值（GS=lgP）。將模塊顯著性定義為模塊中所有基因的平均 GS，從而計算出模塊與樣本性狀間的相關性，選擇相關性最高的模塊為關鍵模塊。確定關鍵模塊后，計算模塊內每個基因的 GS 和模塊隸屬度（module membership，MM），即模塊自身基因與基因表達譜的相關性，并根據 GS 和 MM 設置 hub 基因的篩選閾值。本研究中 hub 基因的截取標準設置為 GS>0.38，MM>0.80。利用在線數據庫 DAVID[7]（
1.6 關鍵基因差異分析及相關性分析
為進一步分析關鍵基因表達情況，使用 R 語言中的“ggpubr”包行差異分析，并基于基因表達水平計算關鍵基因和上游基因間的共表達關系，以確定這些關系在轉錄水平的強弱。使用 R 語言中的“corrplot”包計算基因間的 Pearson 相關性。
1.7 關鍵基因的生存分析
利用在線數據庫 Kaplan-Meier plotter[9]（
1.8 統計學方法
采用 R 軟件進行生信分析和統計學分析。檢驗水準 α=0.05。
2 結果
2.1 mRNAsi 與臨床病理學特征及預后的相關性
正常肝組織與 HCC 組織間 mRNAsi 值的差異有統計學意義（P<0.001），見圖 2a。本研究結果還表明，mRNAsi 與高病理學等級之間存在正相關關系（P=0.001），見圖 2b。Kaplan-Meier 生存曲線結果顯示，高 mRNAsi 組的死亡率高于低 mRNAsi 組，差異有統計學意義（P=0.006），見圖 2c。
圖2
示 mRNAsi 在 HCC 及正常肝組織中的差異及其與病理學分級和預后的相關性、HCC 組織差異基因的火山圖、WGCNA 基因聚類樹、關鍵模塊的顯著性及顯著性較高模塊的鑒定結果
a：正常肝組織（50 個樣品）和 HCC 組織（374 個樣品）之間 mRNAsi 值的差異；b：mRNAsi 值與 HCC 病理學分級相關性分析的箱線圖；c：Kaplan-Meier 生存曲線圖，可見高 mRNAsi 組的死亡率高于低 mRNAsi 組；d：HCC 差異表達基因火山圖，fdr 為假陽性率，綠色圓點代表下調基因，紅色圓點代表上調基因，黑色圓點代表無顯著差異基因；e：HCC 的共表達基因模塊聚類樹圖，橫坐標 22 個不同顏色模塊分別代表一個基因集模塊，縱坐標為基因占比，樹狀圖的分支對應于 22 個不同的基因模塊，樹狀圖上的每一片葉子都對應一個基因，相似的基因被聚類到相同顏色的模塊中；f：基因模塊和臨床特征之間的相關性熱圖，該模塊內數字為模塊中基因的平均 mRNAsi 得分，顏色越紅即與 mRNAsi 和 EREG-mRNAsi 正相關越強，顏色越藍即與 mRNAsi 和 EREG-mRNAsi 負相關越強，并注釋相應的
2.2 HCC 組織與正常肝組織差異表達基因篩選結果
由于人類基因組數據庫龐大，且部分基因表達不具有顯著差異，而正常肝組織的 mRNAsi 值與腫瘤組織明顯不同，故本研究首先進行基因清洗以篩選出差異基因。過濾得到差異表達倍數關系>1 且有統計學意義（P<0.05）的基因。最終篩選了 7 668 個差異基因，其中 7 273 個上調，395 個下調（圖 2d）。
2.3 WGCNA 關鍵模塊及 hub 基因的確定
使用 R 語言中的“WGCNA”包將表達模式相似的基因分配到一個模塊中，形成模塊樹狀圖（圖 2e），得到 22 個基因模塊。然后，計算 22 個模塊內基因與 mRNAsi 和表觀遺傳調控的 mRNAsi（EREG-mRNAsi）的相關性（GS），以及基因與模塊的相關性（MM）。結果表明，綠黃色模塊與 mRNAsi 呈現顯著負相關且相關系數最大（圖 2f），對截取的基因進一步分析后發現并無明顯 mRNAsi 相關性及 LSCS 特性（圖 2g）。而藍色模塊與 HCC 的 mRNAsi 的正相關性最高（r=0.28，P=0.004），表明藍色模塊中的基因可能與 LSCS 密切相關，從藍色模塊篩選 hub 基因（圖 2h），共得到 24 個 Hub 基因：CDCA8、SGOL1、KIF15、GINS1、CDC20、TTK、CKAP2L、ZWINT、NCAPG、EXO1、RAD51AP1、NDC80、CCNB1、MCM2、KIF20A、MCM10、RFC4、CHAF1A、EME1、SPAG5、CDC45、TEDC2、CDCA5 和 OIP5。
2.4 hub 基因的 GO 和 KEGG 通路富集分析
篩選藍色模塊中 GS>0.38，MM>0.80 的基因作為本研究的 hub 基因，并將篩選得到的 24 個 hub 基因進行 GO 和 KEGG 通路分析，以鑒定相關的生物學功能和信號傳導途徑。GO 分析的結果表明，這些基因主要參與的生物過程有：有絲分裂姐妹染色單體分離、核染色體分離、有絲分裂核分裂、染色體分離、DNA 構象變化等；主要的細胞組成有：染色體區域、染色體著絲粒區域、紡錘體等（圖 3a）。此外，KEGG 通路分析發現，這些基因集主要富集于細胞周期、DNA 錯配修復系統、卵母細胞減數分裂和 DNA 復制信號通路（圖 3b）。
圖3
示藍色模塊的 24 個 hub 基因富集分析柱狀圖、hub 基因的 PPI 網絡圖及基因鄰邊節點數，以及關鍵基因的表達差異箱線圖及基因相關性熱圖
a 和 b：GO 分析結果（a）和 KEGG 通路富集分析結果（b），BP 為生物學過程，CC 為細胞組分，直方圖顏色越紅表示統計學差異越大；c：hub 基因 PPI 網絡，圓圈為基因名稱，線條表示基因間的聯系，隱藏了網絡中斷開的節點；d：直方圖示 PPI 網絡中 24 個中樞基因的節點數，其中 2 個斷開節點基因被隱藏，橫坐標為該基因在 PPI 網絡中的鄰邊節點數；e：HCC 中 10 個關鍵基因的差異表達箱線圖，HCC 組織和正常肝組織比較，***
2.5 hub 基因 PPI 網絡的可視化分析
通過 STRING 在線數據庫獲得了 24 個 hub 基因編碼的 PPI 網絡（圖 3c），設置顯示最高置信度（最低互作得分>0.9）。結果顯示，PPI 網絡共由 24 個節點和 44 條邊組成，平均局部聚類度：0.695（P<0.001）。除 CKAP2L2 和 TEDC2 外，其余 22 個基因編碼的蛋白產物均參與了 PPI 網絡的構建，且有多個蛋白位于網絡的中心，與周圍蛋白互作關系密切。使用 R 語言統計結果，并繪制柱狀圖（圖 3d），結果提示，hub 基因之間存在廣泛而牢固的關系，本研究選擇圖 3d 中鄰接節點數最多的前 10 個基因作為關鍵基因，分別為 CCNB1、CDC20、CDCA8、NDC80、KIF20A、CDC45、KIF15、MCM2、TTK 和 NCAPG。
2.6 關鍵基因差異分析及基因共表達分析
使用 R 語言完成關鍵基因的表達分析，并繪制箱線圖（圖 3e），可見 10 個關鍵基因在 HCC 中均顯著高表達（P<0.001）。繪制其共表達相關性熱圖（圖 3f），結果表明基因間相關系數均>0.7，提示關鍵基因相關性較高，上下游作用聯系緊密。
2.7 在線 Kaplan-Meier 生存曲線分析及驗證
使用在線 Kaplan-Meier Plotter 數據庫，以 HCC 組織中關鍵基因表達量的中位數作為分界值將基因表達量分為低表達組和高表達組。共分析了 364 例 HCC 患者的 OS，結果如圖 4a–圖 4j 所示，10 個關鍵基因高表達組均具有不良預后，差異具有統計學意義（P<0.001）。
圖4
示以在線數據庫 Kaplan-Meier plotter 和 HCCDB18 數據集分析 10 個關鍵基因對 HCC 患者預后的影響
a–j：分別示在線數據庫 Kaplan-Meier plotter 分析的 10 個關鍵基因 KIF15（56992）、CDC45（8318）、TTK（7272）、KIF20A（10112）、NDC80（10403）、CDCA8（55143）、CCNB1（891）、MCM2（4171）、CDC20（991）和 NCAPG（64151），括號內數字為基因的數字 ID； k–t：分別示 HCCDB18 數據集分析的 10 個關鍵基因 KIF15、CDC45、TTK、KIF20A、NDC80、CDCA8、CCNB1、MCM2、CDC20和 NCAPG
進一步使用 HCC 綜合分子數據庫中的數據集 HCCDB18 對 HCC 患者的生存分析結果進行驗證。共下載 203 例 HCC 患者的預后數據，為了減少與疾病無關的死亡原因的影響，本數據集中將長于 5 年的生存時間默認為 5 年，并將相應患者的狀態設置為存活。結果如圖 4k-圖 4t 所示，10 個關鍵基因高表達組均與不良預后有關（P<0.001）。
3 討論
基于 WGCNA 的方法對惡性腫瘤的相關基因進行挖掘在近年來逐漸成為腫瘤研究的重要課題。晚期肝癌患者因失去手術機會更加依靠精確的分子靶向治療，基于此背景，如何在精準醫療及基因水平治療的理念和模式下尋找潛在的靶向治療標志物，制定個體化的治療方案顯得尤為重要。
CSCs 是惡性腫瘤無限增殖和復發的源泉，近年來已成為腫瘤領域的研究熱點。CSCs 與腫瘤細胞增殖、侵襲轉移、耐藥等多種惡性表型密切相關[11-12]。本研究中，首先下載了由 Malta 等提供的 HCC 的 mRNAsi，并基于 TCGA 數據庫中的 374 例 HCC 樣本及 50 例正常肝組織樣本進行 mRNAsi 差異的驗證，發現 HCC 樣本中 mRNAsi 值顯著高于正常肝組織樣本，mRNAsi 值和腫瘤病理學分級成正相關；生存曲線分析結果顯示，高 mRNAsi 組的存活率較低 mRNAsi 組顯著下降，這些結果均符合 CSCs 特性，與預期結果一致。基于此，本研究根據肝癌差異基因與 mRNAsi 的關系構建了 WGCNA 網絡，發現藍色模塊與 mRNAsi 的相關性最高，截取的 24 個 hub 基因主要與干細胞自我更新和增殖特征有關。WGCNA 利用分子間的表達相關系數來衡量他們的共表達關系，同一模塊中的分子表達模式相似，表達模式相似的分子可能參與同一生物學過程或通路，然而基因間的相互作用還需要依托蛋白質組學、表觀遺傳組學等綜合判斷，因此限于 WGCNA 的局限性，關鍵基因的挑選需要參考蛋白組學和轉錄組學的結果。本研究中 hub 基因蛋白互作關系緊密，選取 PPI 網絡中鄰邊節點數最多的 10 個基因作為關鍵基因進一步分析。結果關鍵基因在腫瘤組織中的表達量明顯高于正常組織，在轉錄水平上存在很強的共表達關系，且 Kaplan-Meier 生存曲線顯示關鍵基因高表達組的生存率顯著下降并在 HCCDB18 數據集中得到了驗證。既往也發現這些基因的差異表達與腫瘤狀態密切相關。
KIF20A 和 KIF15 是驅動蛋白家族成員，該家族由一組共享高度保守的運動域蛋白質組成[13]。有研究[14-16]表明，KIF20A 在紡錘體組裝和染色體分離中起到重要作用，其在多種腫瘤中高表達且與預后差有關。Shi 等[17]研究發現，KIF20A 是 Hh 信號傳導的重要下游靶基因，此外 Gli 家族的鋅指蛋白是 Shh 信號傳導的介質，被廣泛認為是胚胎癌的致癌基因；Gli2-KIF20A 軸對于 HCC 細胞的增殖、侵襲和腫瘤復發有著至關重要的作用。KIF15 的過表達增加了 CSCs 相關的生物標志物和干性相關的轉錄因子水平，過表達 KIF15 的細胞在連續高比例（1∶20）傳代過程中仍維持類器官分化能力，此外 KIF15 能促進 LSCS 自我更新，并通過 HCC 細胞中活性氧的失衡促進肝癌組織中 LSCS 的干性維持[18]。
GINS1 是一種由 SLD5、PSF1、PSF2 和 PSF3 組成的四聚體復合物，該復合物在 DNA 復制過程中起著至關重要的作用。本研究中 GINS 作為 hub 基因與關鍵基因間有著緊密的聯系，其過表達與肝癌的預后不良相關，這與 Lian 等[19]的研究一致。Nagahama 等[20]證明了 PSF1 是在干細胞和祖細胞中廣泛表達的 DNA 復制因子，具有致瘤和轉移特性。GINS1 四聚體聯合 MCM2-7 與 CDC45 形成 DNA 解旋酶復合物（CDC45-MCM2-7-GINS 復合體）參與 DNA 復制，其參與 DNA 解鏈的機制研究被認為是腫瘤藥物研發的新方向[21-22]。
CCNB1 是細胞周期蛋白 B1（Cyclin B1），是細胞周期家族的成員之一，對于控制細胞周期的 G2/M（有絲分裂）過渡至關重要。近年來，越來越多的研究表明，CCNB1 在肝癌、非小細胞肺癌、結直腸癌、乳腺癌、胃癌等多種腫瘤中均呈高表達[23-24]。最近的研究[25]表明，CCNB1 基因參與了結直腸癌 CSCs 細胞周期的調控。Lee 等[26]通過實驗證實，頭頸部鱗狀細胞癌中 CCNB1 的過表達引起胚胎干細胞關鍵蛋白表達異常，并出現 CSCs 特性，導致腫瘤復發和預后不良。此外，細胞周期蛋白家族中多個成員均被證實參與 CSCs 細胞周期的調控。
NDC80 編碼蛋白質的功能是組織和穩定微管-線粒體相互作用，該蛋白是染色體分離所必需的結構。Shen 等[27]通過綜合分析多個 GEO 數據集發現，NDC80 在 HCC 的發生和復發過程中起到核心作用。近期有研究[28]表明，NDC80 參與前列腺癌干細胞的干性維持，高表達的 NDC80 具有更強的腫瘤細胞成球性，并且參與了多西他賽耐藥性。這些報道均直接或間接證實了本研究的結論。
TTK 編碼具有磷酸化酪氨酸、絲氨酸和蘇氨酸能力的雙重特異性蛋白激酶，這種蛋白質對于有絲分裂過程中染色體的著絲粒對齊至關重要，并且是中心體復制所必需的。已經發現它是有絲分裂期間作用于染色體精確分離的關鍵檢查點蛋白，當該蛋白不能降解并產生過量的中心體而導致異常的有絲分裂紡錘體堆積時，可能會導致肝癌的發生[29]。關于 TTK 基因與 CSCs 的生物學功能的研究已有很多發現，如 TTK 基因在從食管癌細胞系分離的 CSCs 樣細胞群及通過乙醛脫氫酶 1（ALDH1）分選的多發性骨髓瘤干細胞中均過表達[30-31]。有研究[32]表明，TTK 基因的敲低可降低 CD44 陽性的三陰性乳腺癌中 CSCs 的數量。最新研究[33]發現，TTK 幾乎在所有具有 CSCs 特性的癌組織中過表達，小分子化合物 TTK-27f 可有效抑制其作為致癌蛋白激酶的功能，具有強大的類藥物特性并優先抑制 CSCs，因此可能是 CSCs 靶向治療的新策略。
NCAPG 編碼凝縮蛋白復合物的一個亞基，負責有絲分裂和減數分裂過程中染色體的凝結和穩定。研究[34-35]發現，NCAPG 通過磷脂酰肌醇 3 激酶（PI3K）/蛋白激酶 B（Akt）信號通路信號促進 HCC 的增殖，NCAPG 的高表達與肝癌患者整體生存率和無病生存率顯著相關，NCAPG 可預測手術切除后 HCC 的早期復發并作為肝癌患者的預后生物標志物。此外，一項基于生物信息分析的研究[36]發現，NCAPG 可作為膀胱癌干細胞特征的生物標志物，但是該基因在 CSCs 中的生物學功能及作用機制尚未闡明，鑒定 NCAPG 的蛋白質底物及其信號傳導通路將更好地了解該基因在 CSCs 中的關鍵作用。
本研究通過基因相關性分析發現 10 個關鍵基因之間高度相關，且蛋白互作關系緊密。除 KIF15 外，其余關鍵基因雖沒有直接證據表明與 LSCS 的特性有關，但以往的研究結果間接表明，關鍵基因可以促進 LSCS 的干性維持且在 HCC 的復發和化療耐藥性中起重要作用。本研究中 hub 基因的 GO 分析和 KEGG 通路分析結果表明，這些基因主要集中在與細胞有絲分裂相關的功能上并參與細胞周期和 DNA 復制，這也間接反映了這些基因與干細胞的自我更新和增殖特性有關。
綜上所述，本研究結果表明，mRNAsi 在 HCC 中起著重要作用，基于 WGCNA 篩選了與 LSCS 相關的 10 個關鍵基因可能在維持 LSCS 特性方面發揮重要作用，這些基因可能是抑制 LSCS 特性的治療靶點。然而，這些結論僅僅是基于 mRNAsi 模型的分析結論，因此還需進一步研究來驗證這些發現。
重要聲明
利益沖突聲明：本文全體作者閱讀并理解了《中國普外基礎與臨床雜志》的政策聲明，我們沒有相互競爭的利益。
作者貢獻聲明：徐賀龍負責研究設計、數據處理和論文撰寫；李燕參與研究設計；任程暉參與數據處理；宋曉靜和張磊參與論文撰寫；李汛負責論文審改。
