乳腺癌是近年來女性發病率和死亡率最高的惡性腫瘤，是一種影響人類健康的復雜疾病。研究表明，動態網絡生物標志物（DNB）可以有效識別乳腺癌等復雜疾病由正常狀態轉變為疾病狀態的臨界狀態。然而，傳統的 DNB 方法需要用到同一疾病狀態下多個樣本的數據，這通常在臨床診斷中是無法實現的。本文定量分析了人類乳腺癌細胞系（MCF-7）的時間序列數據，基于三維景觀動態網絡生物標志物（L-DNB）方法，找到時間序列中單個樣本的 DNB 模塊，然后構建綜合指標檢測乳腺癌細胞分化過程中的早期預警信號，從而確定其臨界狀態。本文研究結果對于乳腺癌的防治和早期診斷或具有重要意義，期待本文研究能夠為乳腺癌的相關研究提供參考。
引用本文： 趙宏倩, 高潔. 基于三維景觀動態網絡生物標志物識別乳腺癌細胞分化的臨界狀態. 生物醫學工程學雜志, 2020, 37(2): 304-310. doi: 10.7507/1001-5515.201908013 復制
引言
乳腺癌是發生在乳腺上皮組織的惡性腫瘤，目前已在全球女性惡性腫瘤的相關死亡原因中位居首位[1]。乳腺癌在中國的發病率年增幅速度是世界平均水平的 2 倍，到 2021 年將高達 250 萬，因此乳腺癌已成為嚴重威脅我國女性身心健康的惡性腫瘤之一[2]。乳腺癌具有發病率高且逐漸趨于年輕化的特點，現已成為全球公共衛生事業面臨的重大問題，因此采取相應措施來預防和阻斷乳腺癌的發生發展變得十分迫切。由于乳腺癌在初期階段往往沒有特異性癥狀，所以如何針對乳腺癌實現早期診斷，就成為一項具有挑戰性的任務。
復雜疾病的發展是一個動態的過程，通常由正常狀態逐漸積累小的病理變化，最終跨越臨界狀態轉變為疾病狀態。已有研究表明，在許多復雜疾病的正常狀態到疾病狀態發生急劇轉變之前，存在臨界狀態[3-4]。未到達臨界狀態時，疾病往往是可逆的，也就是可以向正常生理狀態轉變，而一旦跨越臨界狀態，許多復雜疾病就無法逆轉了，因此對于乳腺癌這一類復雜疾病，其臨界狀態的檢測和確定顯得尤為重要。2012 年，Chen 等[5]提出用動態網絡生物標志物（dynamic network biomarkers，DNB）的方法檢測疾病惡化的早期預警信號。DNB 方法將疾病的發展視為隨時間變化的非線性動力系統，而疾病的臨界狀態對應著動力系統中的分岔理論，這種無模型的數學方法已經得到了生物信息學界的廣泛認可。盡管 DNB 方法可以檢測復雜疾病的臨界狀態，但該方法在同一疾病狀態下需要多個樣本數據，從而限制了其臨床應用。為解決 DNB 方法需要大量樣本和計算成本高等問題，Liu 等[6]提出了一種利用單個樣本準確地識別出復雜疾病中的臨界狀態的方法。基于上述研究，2019 年，Liu 等[7]提出了一種全新的三維景觀動態網絡生物標志物（landscape dynamic network biomarkers，L-DNB）方法，能夠系統地找到調控疾病惡化的關鍵基因，而無需使用以往研究中的聚類算法或者其他啟發式程序。
乳腺癌是一種發展性疾病，其惡化過程主要是腫瘤細胞的分化或增殖，檢測腫瘤細胞分化的預警信號有利于乳腺癌及其并發癥的防治，極大減少了患者的治療周期[8]。編號為密歇根癌癥基金會-7（Michigan Cancer Foundation–7，MCF-7）的細胞系是人類乳腺癌細胞系，被歸類為浸潤性乳腺導管癌，盡管 MCF-7 細胞分化、增殖的潛在分子機制尚不清楚，但已有研究發現人表皮生長因子受體調節蛋白（heregulin，HRG）、表皮生長因子（epidermal growth factor，EGF）與誘導 MCF-7 細胞分化或增殖的關鍵轉變相關[9]。本文以 MCF-7 細胞為材料，定量分析 HRG 誘導 MCF-7 細胞分化過程的基因表達數據，基于 L-DNB 方法確定參與相關過程的關鍵基因，即 DNB 分子模塊，然后構建單個樣本的綜合指標（composite index，CI）以檢測 MCF-7 細胞分化的臨界狀態。通過本文研究，提供了一種新的途徑來探究細胞分化的潛在機制，可為生物醫學領域中乳腺癌的早期診斷和防治提供參考。
1 材料與方法
本文選取了來自美國國立生物技術信息中心創建的基因表達芯片數據庫（gene expression omnibus，GEO）（網址：https://www.ncbi.nlm.nih.gov/geo）中的一組公開數據集（編號：GSE13009）。該實驗中，使用 HRG 分子持續誘導兩組 MCF-7 細胞分化，在 10 min～72 h 之間共選取 17 個時間點進行采樣，另外有一個基線對照樣本（0 min 樣本）。本文選取 0 min～4 h 共 10 個時間點的樣本進行分析計算，其中 0 min 樣本和 10、15、20 min 的樣本視為參考樣本。
由于該數據集在每個時間點只有一個樣本，無法利用傳統手段計算出基因的相似性網絡，因此可使用樣本特異性網絡來反應每個樣本的特征[10-11]。首先，需要給出足量的參考樣本來表征正常時期基因之間的相關性。給定 n 個參考樣本，參考樣本數據中任意一對基因 x 和 y 之間的相關性可以使用皮爾森相關系數（Pearson correlation coefficient，PCC）表示，記為 
。接下來，將一個新樣本添加到參考樣本組后，計算任意一對基因 x 和 y 在加入樣本后新的 PCC 值，記為 
。因此，針對 n 個參考樣本中基因 x 和 y 的單樣本 PCC（single-sample Pearson correlation coefficient，sPCC）［以符號 
表示］定義如式（1）所示：
 |
式（1）為差分方程，它表示當加入一個新樣本之后，原參考樣本組基因之間 PCC 值發生的擾動。假設把每一個基因視為網絡中的一個節點，
的值視為節點 x 和 y 之間的邊，可以得到每個樣本的基因表達網絡結構。由于樣本中 PCC 遵循正態分布，因此 sPCC 遵循差分正態分布。對每一個 sPCC 值進行雙側 Z-檢驗，只保留置信度水平高于 0.95 的值，反之則視為基因 x 和 y 之間不存在邊。至此，已經構建出一個完整的單樣本差分網絡，此網絡為一個加權無向圖。
為了獲得每個樣本中的 DNB 模塊，需要賦予每個基因 x 一個局部得分。DNB 其本質是一組相互作用的基因，當疾病發展到臨界狀態附近時，它們出現 3 個特征：① 組內分子的平均變異系數增大；② 組內分子的相關性增大；③ 組內分子與組外分子的相關性減小[5]。因此需要構建的基因局部得分指標實際上就是量化上述 DNB 的 3 個特征。
具體來說，對于 DNB 的第一個特征，通過改進 L-DNB 方法[7]，使用基因 x 以及其一階鄰域內所有的基因（sPCC 網絡中與基因 x 有邊相連的基因構成基因 x 的一階鄰域）表達值的平均變異系數（coefficient of variation，CV）（以符號 CV 表示）來表示，表達式如式（2）所示：
 |
其中，Nx 表示基因 x 的一階鄰域，nx 表示一階鄰域中基因的個數，σx 表示基因 x 表達值的標準差，μx 是基因 x 表達值的平均值，σy 表示基因 y 表達值的標準差，μy 是基因 y 表達值的平均值。
第二個特征用基因 x 和它的一階鄰域內的基因的平均 sPCC 值來量化（以符號 sPCCin 表示），表達式如式（3）所示，而對于第三個特征，則考慮基因 x 的一階鄰域和二階鄰域之間的 sPCC 平均值（以符號 sPCCout 表示）如式（4）所示：
 |
 |
其中，Mx 表示基因 x 的二階鄰域，mx 表示二階鄰域中基因的個數。因此，基于單樣本的基因 x 局部模塊得分指標（local score index，LSI）［以符號 LSI（x）表示］如式（5）所示：
 |
將每個樣本按照細胞分化進度進行分組，得到一個縱軸為 LSI 的三維景觀圖。此外，對于每個基因 x，LSI（x）可以被認為是該基因在不同時間點的局部得分，因此可以在該樣本中將此得分作為基因關鍵性排序的依據。
在從單樣本中的相應模塊獲得每個基因的 LSI 之后，所有基因按得分降序排列。排序列表中得分最高的 k 個基因就可以被視為該樣本的 DNB，并且可以通過對這 k 個基因進行求和來簡單地估計相應的全局得分。對于任意一個樣本，全局得分即 CI（以符號 CI 表示）如式（6）所示。CI 作為檢測即將到來的臨界狀態的特定標識符，其快速上升預示著臨界狀態出現，而在其他時間點幾乎沒有顯著波動。
 |
全部的算法流程示意圖如圖 1 所示。
圖1
算法的流程示意圖
Figure1.
Schematic diagram of the algorithm
2 結果
2.1 局部得分的三維景觀圖
數據集中包含兩位患者的數據，因此將數據分為兩組，分別計算前 10 個時間點的局部得分并利用計算機程序設計語言 Python（3.7.1, Anaconda Inc.，美國）繪制三維景觀圖，如圖 2 所示。圖 2 左圖為第一例樣本局部得分的三維景觀圖，右圖則為第二例樣本。雖然基因數高達 2 200，三維景觀圖中的數據較為密集，但可以觀察到兩位患者的 MCF-7 細胞分化情況十分相似，其過程中都存在著一個明顯的“高峰”，實際上本文算法的目的就是試圖去量化這一“高峰”，“高峰”表明該時刻基因的局部模塊得分達到了峰值，即 HRG 誘導的 MCF-7 細胞分化過程中出現臨界狀態。
圖2
兩例樣本的三維景觀圖
Figure2.
Landscape figures of two samples
2.2 利用全局得分量化疾病臨界狀態
在 2.1 小節的基礎上，首先對每一個樣本中所有基因的局部得分進行排序。將排序后的結果投影至二維平面并繪制熱圖，如圖 3 所示。圖 3 左圖為第一例樣本 LSI 映射到二維平面的熱圖，顏色越接近紅色表示數值越大，右圖則為第二例樣本。本文旨在找出乳腺癌細胞分化過程中起關鍵調控作用的基因，因此只需在每個樣本中選取局部得分最高的基因作為 DNB 模塊。特別地，本文選取每個樣本排在前 20 的基因。
圖3
三維景觀圖映射到二維平面的熱圖
Figure3.
Landscape figures mapped to heat maps of a two-dimensional plane
由式（6）分別計算出兩例樣本 MCF-7 細胞分化的時間過程曲線，如圖 4 所示。圖 4 中藍色曲線為第一例樣本的 MCF-7 細胞分化過程圖，橙色曲線為第二例樣本。兩例樣本均在 0.75～1 h 出現 CI 值的突增，而其他時間點的 CI 值明顯較低，這意味著 0.75～1 h 為 HRG 誘導的 MCF-7 細胞分化的關鍵階段。實際上，在 HRG 刺激后 0.75 h，DNB 會發出指示性的預警信號，CI 達到峰值，對應為 MCF-7 細胞分化過程的臨界狀態。這與實驗結果一致，用 EGF 和 HRG 刺激的 MCF-7 細胞可產生長達 0.75 h 非常相似的早期轉錄譜，隨后的細胞表型在 3 h 后有所不同，這表明分化時間約為 3 h[9]。
圖4
兩例樣本 CI 值的時間過程曲線
Figure4.
Time course curves of composite index of two cases
現已知，激活子蛋白-1（activator?protein?1，AP-1）是細胞內的一個轉錄激活因子，Fos 原癌基因（Fos proto-oncogene，FOS）作為一類核蛋白轉錄因子，在調控細胞生長、分裂、增殖、分化乃至程序性死亡等方面具有重要的作用，Jun 原癌基因（Jun proto-oncogene，JUN）可以編碼一種與病毒蛋白高度相似的蛋白，在控制細胞是否惡變方面有重要的作用。AP-1 是由 FOS 和 JUN 基因家族編碼的蛋白質的二聚體組成，已廣泛參與細胞分化、增殖和轉化，而 FOS 蛋白家族（c-FOS，FOSB，FRA-1 和 FRA-2）與 JUN 蛋白家族（c-JUN，JUNB 和 JUND）形成異二聚體，調節各種啟動子中 12—氧—十四烷佛波醇—13—乙酸酯（12-O-tetradecanoyl-phorbol-13-acetate，TPA）反應元件（TPA responsive element，TRE）基因表達。實驗分析表明，盡管 c-JUN 的締合是短暫的，但經過 HRG 處理的 MCF-7 細胞中的 AP-1 復合物包含 c-JUN，c-FOS 和 FRA-1[9]。因此，本文的指標應用結果與實驗觀察相吻合，并成功地檢測到了即將發生的臨界轉變的預警信號。兩例樣本同樣取自于 MCF-7 細胞系，而它們分化過程中的 CI 曲線非常一致，這表現出 L-DNB 算法也具有良好的協同性。事實上，L-DNB 算法是一種適用于臨床的個性化算法，只要給出足量的參考樣本，即可對任意單樣本進行計算分析，且不需要任何的先驗知識。
2.3 通路分析與文獻挖掘
在選取的兩個病例共 20 個樣本中，一共有 241 個基因被鑒定為參與調控乳腺癌細胞分化過程的 DNB 模塊。為了驗證這些基因的生物學功能與乳腺癌的關聯性，本文從兩個方面進行分析。首先利用京都基因和基因型百科全書同源注釋在線數據庫（Kyoto Encyclopedia of Genes and Genomes Orthology Based Annotation System，KOBAS）（網址：http://kobas.cbi.pku.edu.cn）對這 241 個基因進行京都基因與基因組百科全書（Kyoto Encyclopedia of Genes and Genomes，KEGG）通路分析[12]，篩選出校正 P < 0.05 且與 MCF-7 細胞增殖分化密切相關的 11 條通路，如表 1 所示。
表1
KEGG 通路分析結果
Table1.
KEGG pathway analysis results
從表 1 中可以發現，病毒致癌作用和癌癥中的通路這兩條通路的校正 P 值都在 10-5數量級以下，并且路徑中的基因分別有 11 和 14 個。這驗證了所選出的 DNB 確實有很多處于癌癥相關的通路中，是參與誘導和調控乳腺癌細胞增殖與分化的關鍵因素。表 1 中，其他具有較小校正 P 值（較高置信度）的通路也與乳腺癌或者惡性腫瘤密切相關，例如磷脂酰肌醇 3 激酶-蛋白質絲氨酸蘇氨酸激酶（phosphatidylinositol 3 kinases-protein serine threonine kinase signaling pathway，PI3K-Akt）信號通路和絲裂原活化蛋白激酶（mitogen-activated protein kinase signaling pathway，MAPK）信號通路已經被證實與腫瘤細胞凋亡相關，在人表皮生長因子受體-2（human epidermal growth factor receptor-2，HER2）陽性乳腺腫瘤中，HER2 激活 PI3K-Akt 信號通路和 MAPK 信號通路，并刺激細胞生長、存活和分化；腫瘤抑制基因 p53 信號通路中的 p53 基因則是許多癌癥的重要生物標志物[13-14]。
另外，在本文確定的 241 個 DNB 中，有相當一部分已經被文獻證實與乳腺癌相關。例如有多項研究表明乳腺癌易感基因（breast cancer 2，BRCA2）伴侶及定位蛋白 2（partner and localizer of BRCA2，PALB2）基因的突變與早發性乳腺癌具有密切關聯[15-16]。谷氨酰-tRNA 酰胺轉移酶結合蛋白 3（glutamyl-tRNA amidotransferase binding protein 3，GATA3）是由乳房中的腔上皮細胞高度表達的轉錄激活因子，它被鑒定為乳腺癌的預后標志物[17]，在所有乳腺癌中，GATA3 基因上體細胞突變率超過 10%。而在一項同樣關于 MCF-7 細胞的研究中，生長因子受體結合蛋白 2（growth factor receptor bound protein 2，GRB2）基因有明顯的 mRNA 過表達現象[18]。在乳腺癌細胞分化過程中臨界狀態新出現的 DNB 有肌動蛋白 α4（actinin alpha 4，ACTN4）、CRK 樣原癌基因（CRK like proto-oncogene，CRKL）、普列克底物蛋白同源樣結構域家族 A 成員 1 蛋白（Pleckstrin homology like domain family A member 1，PHLDA1）、衣殼蛋白復合體亞基 ε（coat protein complex subunit epsilon，COPE）、重組人 DNA 結合抑制劑 1（inhibitor of DNA binding 1，ID1）等。其中，ACTN4 基因被鑒定為與癌細胞生物學功能相關并且在包括乳腺癌在內的許多人類上皮癌中高表達[19]；CRKL 基因在基質細胞衍生因子-1（stromal cell-derived factor-1，SDF-1）誘導的細胞外調節蛋白激酶（extracellular regulated protein kinases，Erk1/2）和 PI3K-Akt 信號通路中起調節作用，并進一步控制了乳腺癌細胞的侵襲和遷移[20]；PHLDA1 基因的下調是乳腺癌患者預后不良的重要預兆[21]。轉錄輔助抑制因子 C-末端結合蛋白 2（C-terminal binding protein 2，CTBP2）、M2-型丙酮酸激酶（pyruvate kinase isozymes M2，PKM2）、PALB2、CRKL、三磷酸鳥苷（guanosinetriphosphate，GTP）結合絲裂原誘導蛋白（GTP binding protein overexpressed in skeletal muscle，GEM）等為乳腺癌細胞分化過程中多個時間點的 DNB。其中，CTBP2 基因能夠調節細胞轉化，并且在乳腺癌細胞的增殖和存活中發揮著重要作用[22]；腫瘤細胞轉移是女性乳腺癌相關死亡原因之一，腫瘤干細胞在轉移中起重要作用，PKM2 基因在乳腺癌干細胞中有特定作用，并且在乳腺癌細胞自我更新中發揮作用[23]；PALB2 基因被認為是乳腺癌易感基因，PALB2 種系功能喪失突變會導致罹患乳腺癌的風險增加[24]。
根據已有研究成果，充分證明了用于檢測 MCF-7 細胞分化過程中臨界狀態的 DNB 模塊具有可靠性。同時，基于 DNB 理論，癌癥的演變往往不是某幾個基因的單獨作用，而是部分基因的相互作用，本文尋找到的 241 個基因能夠為新的乳腺癌相關基因的發現提供精確的線索和指引。
3 討論
隨著高通量測序技術的日益發展，大量基因表達微陣列數據涌現，許多數學方法使得乳腺癌的早期預警識別成為了可能。本文基于 L-DNB 方法，選取了 HRG 誘導 MCF-7 細胞分化過程的微陣列數據，確定了兩組病例惡性腫瘤演變過程中 MCF-7 細胞分化的臨界狀態。實際上，在 HRG 刺激后的 0.75～1 h 間，DNB 會發出指示性預警信號。
本文提供了一種新的途徑來探究細胞分化的潛在機制，從而有助于及時對乳腺癌進行防治。與傳統的生物標志物方法相比，本文的算法能夠解決臨床診斷中單個樣本的問題。算法最終構造了單個樣本的綜合指標，它可以作為一個普適性的特征應用于任何乳腺癌樣本或其他的復雜疾病。根據通路分析以及查閱文獻的結果，參與調控疾病過程的 DNB 模塊有很多已經被證實與乳腺癌相關，這也驗證了該算法的生物學可靠性。此外，算法中參考樣本的選取對結果有直接影響，理論上需要選取足量的正常樣本才能構建一個穩定的參考網絡，但是生物學實驗中的對照樣本數量往往不會很多，這可能帶來一些較大的噪聲，本課題組后期將關注如何消除參考樣本帶來的誤差。
利益沖突聲明：本文全體作者均聲明不存在利益沖突。
引言
乳腺癌是發生在乳腺上皮組織的惡性腫瘤，目前已在全球女性惡性腫瘤的相關死亡原因中位居首位[1]。乳腺癌在中國的發病率年增幅速度是世界平均水平的 2 倍，到 2021 年將高達 250 萬，因此乳腺癌已成為嚴重威脅我國女性身心健康的惡性腫瘤之一[2]。乳腺癌具有發病率高且逐漸趨于年輕化的特點，現已成為全球公共衛生事業面臨的重大問題，因此采取相應措施來預防和阻斷乳腺癌的發生發展變得十分迫切。由于乳腺癌在初期階段往往沒有特異性癥狀，所以如何針對乳腺癌實現早期診斷，就成為一項具有挑戰性的任務。
復雜疾病的發展是一個動態的過程，通常由正常狀態逐漸積累小的病理變化，最終跨越臨界狀態轉變為疾病狀態。已有研究表明，在許多復雜疾病的正常狀態到疾病狀態發生急劇轉變之前，存在臨界狀態[3-4]。未到達臨界狀態時，疾病往往是可逆的，也就是可以向正常生理狀態轉變，而一旦跨越臨界狀態，許多復雜疾病就無法逆轉了，因此對于乳腺癌這一類復雜疾病，其臨界狀態的檢測和確定顯得尤為重要。2012 年，Chen 等[5]提出用動態網絡生物標志物（dynamic network biomarkers，DNB）的方法檢測疾病惡化的早期預警信號。DNB 方法將疾病的發展視為隨時間變化的非線性動力系統，而疾病的臨界狀態對應著動力系統中的分岔理論，這種無模型的數學方法已經得到了生物信息學界的廣泛認可。盡管 DNB 方法可以檢測復雜疾病的臨界狀態，但該方法在同一疾病狀態下需要多個樣本數據，從而限制了其臨床應用。為解決 DNB 方法需要大量樣本和計算成本高等問題，Liu 等[6]提出了一種利用單個樣本準確地識別出復雜疾病中的臨界狀態的方法。基于上述研究，2019 年，Liu 等[7]提出了一種全新的三維景觀動態網絡生物標志物（landscape dynamic network biomarkers，L-DNB）方法，能夠系統地找到調控疾病惡化的關鍵基因，而無需使用以往研究中的聚類算法或者其他啟發式程序。
乳腺癌是一種發展性疾病，其惡化過程主要是腫瘤細胞的分化或增殖，檢測腫瘤細胞分化的預警信號有利于乳腺癌及其并發癥的防治，極大減少了患者的治療周期[8]。編號為密歇根癌癥基金會-7（Michigan Cancer Foundation–7，MCF-7）的細胞系是人類乳腺癌細胞系，被歸類為浸潤性乳腺導管癌，盡管 MCF-7 細胞分化、增殖的潛在分子機制尚不清楚，但已有研究發現人表皮生長因子受體調節蛋白（heregulin，HRG）、表皮生長因子（epidermal growth factor，EGF）與誘導 MCF-7 細胞分化或增殖的關鍵轉變相關[9]。本文以 MCF-7 細胞為材料，定量分析 HRG 誘導 MCF-7 細胞分化過程的基因表達數據，基于 L-DNB 方法確定參與相關過程的關鍵基因，即 DNB 分子模塊，然后構建單個樣本的綜合指標（composite index，CI）以檢測 MCF-7 細胞分化的臨界狀態。通過本文研究，提供了一種新的途徑來探究細胞分化的潛在機制，可為生物醫學領域中乳腺癌的早期診斷和防治提供參考。
1 材料與方法
本文選取了來自美國國立生物技術信息中心創建的基因表達芯片數據庫（gene expression omnibus，GEO）（網址：https://www.ncbi.nlm.nih.gov/geo）中的一組公開數據集（編號：GSE13009）。該實驗中，使用 HRG 分子持續誘導兩組 MCF-7 細胞分化，在 10 min～72 h 之間共選取 17 個時間點進行采樣，另外有一個基線對照樣本（0 min 樣本）。本文選取 0 min～4 h 共 10 個時間點的樣本進行分析計算，其中 0 min 樣本和 10、15、20 min 的樣本視為參考樣本。
由于該數據集在每個時間點只有一個樣本，無法利用傳統手段計算出基因的相似性網絡，因此可使用樣本特異性網絡來反應每個樣本的特征[10-11]。首先，需要給出足量的參考樣本來表征正常時期基因之間的相關性。給定 n 個參考樣本，參考樣本數據中任意一對基因 x 和 y 之間的相關性可以使用皮爾森相關系數（Pearson correlation coefficient，PCC）表示，記為 。接下來，將一個新樣本添加到參考樣本組后，計算任意一對基因 x 和 y 在加入樣本后新的 PCC 值，記為 。因此，針對 n 個參考樣本中基因 x 和 y 的單樣本 PCC（single-sample Pearson correlation coefficient，sPCC）［以符號 表示］定義如式（1）所示：
|
式（1）為差分方程，它表示當加入一個新樣本之后，原參考樣本組基因之間 PCC 值發生的擾動。假設把每一個基因視為網絡中的一個節點，的值視為節點 x 和 y 之間的邊，可以得到每個樣本的基因表達網絡結構。由于樣本中 PCC 遵循正態分布，因此 sPCC 遵循差分正態分布。對每一個 sPCC 值進行雙側 Z-檢驗，只保留置信度水平高于 0.95 的值，反之則視為基因 x 和 y 之間不存在邊。至此，已經構建出一個完整的單樣本差分網絡，此網絡為一個加權無向圖。
為了獲得每個樣本中的 DNB 模塊，需要賦予每個基因 x 一個局部得分。DNB 其本質是一組相互作用的基因，當疾病發展到臨界狀態附近時，它們出現 3 個特征：① 組內分子的平均變異系數增大；② 組內分子的相關性增大；③ 組內分子與組外分子的相關性減小[5]。因此需要構建的基因局部得分指標實際上就是量化上述 DNB 的 3 個特征。
具體來說，對于 DNB 的第一個特征，通過改進 L-DNB 方法[7]，使用基因 x 以及其一階鄰域內所有的基因（sPCC 網絡中與基因 x 有邊相連的基因構成基因 x 的一階鄰域）表達值的平均變異系數（coefficient of variation，CV）（以符號 CV 表示）來表示，表達式如式（2）所示：
|
其中，Nx 表示基因 x 的一階鄰域，nx 表示一階鄰域中基因的個數，σx 表示基因 x 表達值的標準差，μx 是基因 x 表達值的平均值，σy 表示基因 y 表達值的標準差，μy 是基因 y 表達值的平均值。
第二個特征用基因 x 和它的一階鄰域內的基因的平均 sPCC 值來量化（以符號 sPCCin 表示），表達式如式（3）所示，而對于第三個特征，則考慮基因 x 的一階鄰域和二階鄰域之間的 sPCC 平均值（以符號 sPCCout 表示）如式（4）所示：
|
|
其中，Mx 表示基因 x 的二階鄰域，mx 表示二階鄰域中基因的個數。因此，基于單樣本的基因 x 局部模塊得分指標（local score index，LSI）［以符號 LSI（x）表示］如式（5）所示：
|
將每個樣本按照細胞分化進度進行分組，得到一個縱軸為 LSI 的三維景觀圖。此外，對于每個基因 x，LSI（x）可以被認為是該基因在不同時間點的局部得分，因此可以在該樣本中將此得分作為基因關鍵性排序的依據。
在從單樣本中的相應模塊獲得每個基因的 LSI 之后，所有基因按得分降序排列。排序列表中得分最高的 k 個基因就可以被視為該樣本的 DNB，并且可以通過對這 k 個基因進行求和來簡單地估計相應的全局得分。對于任意一個樣本，全局得分即 CI（以符號 CI 表示）如式（6）所示。CI 作為檢測即將到來的臨界狀態的特定標識符，其快速上升預示著臨界狀態出現，而在其他時間點幾乎沒有顯著波動。
|
全部的算法流程示意圖如圖 1 所示。
圖1
算法的流程示意圖
Figure1.
Schematic diagram of the algorithm
2 結果
2.1 局部得分的三維景觀圖
數據集中包含兩位患者的數據，因此將數據分為兩組，分別計算前 10 個時間點的局部得分并利用計算機程序設計語言 Python（3.7.1, Anaconda Inc.，美國）繪制三維景觀圖，如圖 2 所示。圖 2 左圖為第一例樣本局部得分的三維景觀圖，右圖則為第二例樣本。雖然基因數高達 2 200，三維景觀圖中的數據較為密集，但可以觀察到兩位患者的 MCF-7 細胞分化情況十分相似，其過程中都存在著一個明顯的“高峰”，實際上本文算法的目的就是試圖去量化這一“高峰”，“高峰”表明該時刻基因的局部模塊得分達到了峰值，即 HRG 誘導的 MCF-7 細胞分化過程中出現臨界狀態。
圖2
兩例樣本的三維景觀圖
Figure2.
Landscape figures of two samples
2.2 利用全局得分量化疾病臨界狀態
在 2.1 小節的基礎上，首先對每一個樣本中所有基因的局部得分進行排序。將排序后的結果投影至二維平面并繪制熱圖，如圖 3 所示。圖 3 左圖為第一例樣本 LSI 映射到二維平面的熱圖，顏色越接近紅色表示數值越大，右圖則為第二例樣本。本文旨在找出乳腺癌細胞分化過程中起關鍵調控作用的基因，因此只需在每個樣本中選取局部得分最高的基因作為 DNB 模塊。特別地，本文選取每個樣本排在前 20 的基因。
圖3
三維景觀圖映射到二維平面的熱圖
Figure3.
Landscape figures mapped to heat maps of a two-dimensional plane
由式（6）分別計算出兩例樣本 MCF-7 細胞分化的時間過程曲線，如圖 4 所示。圖 4 中藍色曲線為第一例樣本的 MCF-7 細胞分化過程圖，橙色曲線為第二例樣本。兩例樣本均在 0.75～1 h 出現 CI 值的突增，而其他時間點的 CI 值明顯較低，這意味著 0.75～1 h 為 HRG 誘導的 MCF-7 細胞分化的關鍵階段。實際上，在 HRG 刺激后 0.75 h，DNB 會發出指示性的預警信號，CI 達到峰值，對應為 MCF-7 細胞分化過程的臨界狀態。這與實驗結果一致，用 EGF 和 HRG 刺激的 MCF-7 細胞可產生長達 0.75 h 非常相似的早期轉錄譜，隨后的細胞表型在 3 h 后有所不同，這表明分化時間約為 3 h[9]。
圖4
兩例樣本 CI 值的時間過程曲線
Figure4.
Time course curves of composite index of two cases
現已知，激活子蛋白-1（activator?protein?1，AP-1）是細胞內的一個轉錄激活因子，Fos 原癌基因（Fos proto-oncogene，FOS）作為一類核蛋白轉錄因子，在調控細胞生長、分裂、增殖、分化乃至程序性死亡等方面具有重要的作用，Jun 原癌基因（Jun proto-oncogene，JUN）可以編碼一種與病毒蛋白高度相似的蛋白，在控制細胞是否惡變方面有重要的作用。AP-1 是由 FOS 和 JUN 基因家族編碼的蛋白質的二聚體組成，已廣泛參與細胞分化、增殖和轉化，而 FOS 蛋白家族（c-FOS，FOSB，FRA-1 和 FRA-2）與 JUN 蛋白家族（c-JUN，JUNB 和 JUND）形成異二聚體，調節各種啟動子中 12—氧—十四烷佛波醇—13—乙酸酯（12-O-tetradecanoyl-phorbol-13-acetate，TPA）反應元件（TPA responsive element，TRE）基因表達。實驗分析表明，盡管 c-JUN 的締合是短暫的，但經過 HRG 處理的 MCF-7 細胞中的 AP-1 復合物包含 c-JUN，c-FOS 和 FRA-1[9]。因此，本文的指標應用結果與實驗觀察相吻合，并成功地檢測到了即將發生的臨界轉變的預警信號。兩例樣本同樣取自于 MCF-7 細胞系，而它們分化過程中的 CI 曲線非常一致，這表現出 L-DNB 算法也具有良好的協同性。事實上，L-DNB 算法是一種適用于臨床的個性化算法，只要給出足量的參考樣本，即可對任意單樣本進行計算分析，且不需要任何的先驗知識。
2.3 通路分析與文獻挖掘
在選取的兩個病例共 20 個樣本中，一共有 241 個基因被鑒定為參與調控乳腺癌細胞分化過程的 DNB 模塊。為了驗證這些基因的生物學功能與乳腺癌的關聯性，本文從兩個方面進行分析。首先利用京都基因和基因型百科全書同源注釋在線數據庫（Kyoto Encyclopedia of Genes and Genomes Orthology Based Annotation System，KOBAS）（網址：http://kobas.cbi.pku.edu.cn）對這 241 個基因進行京都基因與基因組百科全書（Kyoto Encyclopedia of Genes and Genomes，KEGG）通路分析[12]，篩選出校正 P < 0.05 且與 MCF-7 細胞增殖分化密切相關的 11 條通路，如表 1 所示。
表1
KEGG 通路分析結果
Table1.
KEGG pathway analysis results
從表 1 中可以發現，病毒致癌作用和癌癥中的通路這兩條通路的校正 P 值都在 10-5數量級以下，并且路徑中的基因分別有 11 和 14 個。這驗證了所選出的 DNB 確實有很多處于癌癥相關的通路中，是參與誘導和調控乳腺癌細胞增殖與分化的關鍵因素。表 1 中，其他具有較小校正 P 值（較高置信度）的通路也與乳腺癌或者惡性腫瘤密切相關，例如磷脂酰肌醇 3 激酶-蛋白質絲氨酸蘇氨酸激酶（phosphatidylinositol 3 kinases-protein serine threonine kinase signaling pathway，PI3K-Akt）信號通路和絲裂原活化蛋白激酶（mitogen-activated protein kinase signaling pathway，MAPK）信號通路已經被證實與腫瘤細胞凋亡相關，在人表皮生長因子受體-2（human epidermal growth factor receptor-2，HER2）陽性乳腺腫瘤中，HER2 激活 PI3K-Akt 信號通路和 MAPK 信號通路，并刺激細胞生長、存活和分化；腫瘤抑制基因 p53 信號通路中的 p53 基因則是許多癌癥的重要生物標志物[13-14]。
另外，在本文確定的 241 個 DNB 中，有相當一部分已經被文獻證實與乳腺癌相關。例如有多項研究表明乳腺癌易感基因（breast cancer 2，BRCA2）伴侶及定位蛋白 2（partner and localizer of BRCA2，PALB2）基因的突變與早發性乳腺癌具有密切關聯[15-16]。谷氨酰-tRNA 酰胺轉移酶結合蛋白 3（glutamyl-tRNA amidotransferase binding protein 3，GATA3）是由乳房中的腔上皮細胞高度表達的轉錄激活因子，它被鑒定為乳腺癌的預后標志物[17]，在所有乳腺癌中，GATA3 基因上體細胞突變率超過 10%。而在一項同樣關于 MCF-7 細胞的研究中，生長因子受體結合蛋白 2（growth factor receptor bound protein 2，GRB2）基因有明顯的 mRNA 過表達現象[18]。在乳腺癌細胞分化過程中臨界狀態新出現的 DNB 有肌動蛋白 α4（actinin alpha 4，ACTN4）、CRK 樣原癌基因（CRK like proto-oncogene，CRKL）、普列克底物蛋白同源樣結構域家族 A 成員 1 蛋白（Pleckstrin homology like domain family A member 1，PHLDA1）、衣殼蛋白復合體亞基 ε（coat protein complex subunit epsilon，COPE）、重組人 DNA 結合抑制劑 1（inhibitor of DNA binding 1，ID1）等。其中，ACTN4 基因被鑒定為與癌細胞生物學功能相關并且在包括乳腺癌在內的許多人類上皮癌中高表達[19]；CRKL 基因在基質細胞衍生因子-1（stromal cell-derived factor-1，SDF-1）誘導的細胞外調節蛋白激酶（extracellular regulated protein kinases，Erk1/2）和 PI3K-Akt 信號通路中起調節作用，并進一步控制了乳腺癌細胞的侵襲和遷移[20]；PHLDA1 基因的下調是乳腺癌患者預后不良的重要預兆[21]。轉錄輔助抑制因子 C-末端結合蛋白 2（C-terminal binding protein 2，CTBP2）、M2-型丙酮酸激酶（pyruvate kinase isozymes M2，PKM2）、PALB2、CRKL、三磷酸鳥苷（guanosinetriphosphate，GTP）結合絲裂原誘導蛋白（GTP binding protein overexpressed in skeletal muscle，GEM）等為乳腺癌細胞分化過程中多個時間點的 DNB。其中，CTBP2 基因能夠調節細胞轉化，并且在乳腺癌細胞的增殖和存活中發揮著重要作用[22]；腫瘤細胞轉移是女性乳腺癌相關死亡原因之一，腫瘤干細胞在轉移中起重要作用，PKM2 基因在乳腺癌干細胞中有特定作用，并且在乳腺癌細胞自我更新中發揮作用[23]；PALB2 基因被認為是乳腺癌易感基因，PALB2 種系功能喪失突變會導致罹患乳腺癌的風險增加[24]。
根據已有研究成果，充分證明了用于檢測 MCF-7 細胞分化過程中臨界狀態的 DNB 模塊具有可靠性。同時，基于 DNB 理論，癌癥的演變往往不是某幾個基因的單獨作用，而是部分基因的相互作用，本文尋找到的 241 個基因能夠為新的乳腺癌相關基因的發現提供精確的線索和指引。
3 討論
隨著高通量測序技術的日益發展，大量基因表達微陣列數據涌現，許多數學方法使得乳腺癌的早期預警識別成為了可能。本文基于 L-DNB 方法，選取了 HRG 誘導 MCF-7 細胞分化過程的微陣列數據，確定了兩組病例惡性腫瘤演變過程中 MCF-7 細胞分化的臨界狀態。實際上，在 HRG 刺激后的 0.75～1 h 間，DNB 會發出指示性預警信號。
本文提供了一種新的途徑來探究細胞分化的潛在機制，從而有助于及時對乳腺癌進行防治。與傳統的生物標志物方法相比，本文的算法能夠解決臨床診斷中單個樣本的問題。算法最終構造了單個樣本的綜合指標，它可以作為一個普適性的特征應用于任何乳腺癌樣本或其他的復雜疾病。根據通路分析以及查閱文獻的結果，參與調控疾病過程的 DNB 模塊有很多已經被證實與乳腺癌相關，這也驗證了該算法的生物學可靠性。此外，算法中參考樣本的選取對結果有直接影響，理論上需要選取足量的正常樣本才能構建一個穩定的參考網絡，但是生物學實驗中的對照樣本數量往往不會很多，這可能帶來一些較大的噪聲，本課題組后期將關注如何消除參考樣本帶來的誤差。
利益沖突聲明：本文全體作者均聲明不存在利益沖突。
