在基因—藥物多重網絡中利用模塊鑒別方法推斷出新的基因—藥物關聯，可以為已知的藥物確定新的治療靶向基因。本文基于腫瘤藥物敏感性基因組學（GDSC）數據庫中肺癌的基因表達數據和藥物響應數據，提出一種多重網絡算法，首先構造出不同細胞系下肺癌基因與藥物的異構網絡，然后利用基于圖熵的網絡模塊鑒別方法，在該異構網絡中進行模塊的鑒別，通過迭代收斂共鑒別出5個肺癌基因—藥物關聯模塊。與其它方法對比，該算法在運行時間、準確性和魯棒性方面均有更好的效果，且鑒別出的模塊具有明顯的生物學意義。本文的研究結果對于肺癌的用藥治療具有指導意義，并且可以為其它具有相同靶向基因的疾病治療提供參考。
引用本文： 周剛, 高潔. 基于多重網絡算法的肺癌基因—藥物模塊鑒別. 生物醫學工程學雜志, 2021, 38(6): 1111-1117. doi: 10.7507/1001-5515.202102006 復制
引言
新藥物的研發是一個具有風險、周期長且成本高的過程[1-2]。盡管最近幾十年醫藥公司在藥物研發上的投入已經增加了很多，但是被美國食品藥品監督管理局（food and drug administration，FDA）認可的新藥物數量卻不是很多。為了解決這個難題，對已經存在的藥物進行重新開發利用是一項非常有前景的工作。
惡性腫瘤嚴重影響人類生命健康，其中肺癌是發病率和死亡率增長最快，對人體健康和生命威脅最大的惡性腫瘤之一。據報道，近50年來許多國家肺癌的發病率和死亡率均明顯增高，其中男性的肺癌發病率和死亡率均占所有惡性腫瘤的第一位，而女性的肺癌發病率和死亡率均占所有惡性腫瘤第二位。惡性腫瘤的發生與基因突變密切相關，最近的研究表明，每一個腫瘤患者都擁有一個特殊的基因表達圖譜，但是患同一種惡性腫瘤的不同患者，他們的基因表達圖譜卻又具有一定的相似性，即相似的基因突變可能產生同一種疾病，反之亦然。研究表明，用于治愈同一種疾病的不同藥物具有相似的化學結構，且具有相似基因表達圖譜的細胞系對同一種藥物具有相似的化學反應[3-5]。因此，可以通過研究同一種惡性腫瘤的基因表達圖譜和治療該惡性腫瘤藥物之間的模塊化結構，鑒別新的基因—藥物關聯，發現新的藥物靶向基因，進一步將已經存在的藥物用于其它惡性腫瘤治療中。這將在惡性腫瘤的精準治療方面發揮巨大的作用，增加惡性腫瘤治愈的可能性，減少新藥物研發的成本。
高通量測序技術的發展產生了豐富的多組學數據，為惡性腫瘤致病機制的研究提供了極大的便利。不少學者已經提出了各種方法來對惡性腫瘤的多組學數據進行分析和研究，其中從網絡的角度整合多組學數據，對惡性腫瘤進行研究，是現在比較流行的一種分析方法。多重網絡可以充分整合某一惡性腫瘤的多組學數據，與單一網絡相比，多重網絡可以更全面地描述與惡性腫瘤有關的信息。2019年，Yuan等[6]通過整合乳腺癌的DNA甲基化數據、拷貝數變異數據和二者在基因調控網絡中的關聯，使用雙權相關系數衡量基因調控網絡中的生物分子關聯，并且利用非凸懲罰函數使網絡模型達到最優，最終在基因調控網絡中推斷出與乳腺癌發生相關的DNA甲基化位點和拷貝數變異位點。Valdeolivas等[7]利用重啟隨機游走的方法在由蛋白質、基因和疾病構成的多重異構網絡上進行研究，結果證明了整合的多組學數據可以有效提高實驗結果的準確性。Shang等[8]基于受約束的網頁排序（constrained PageRank，CPR）算法，將肝細胞癌的多組學數據整合在一個統一的網絡框架中，通過實驗發現，與單一數據集相比，使用多組學數據進行網絡構造時取得的實驗效果明顯優于單一數據集。Ma等[9]利用全基因組的DNA甲基化數據和基因表達數據，分別構造出基因差異共表達網絡和基因差異共甲基化網絡，然后在這兩個網絡上利用基于圖熵的方法發現表觀遺傳模塊，該模塊可以作為生物標志物用于疾病的診斷和治療。由此可見，在網絡中利用多組學數據進行生物分子關聯的鑒別是可行且具有巨大發展潛力的。而在基因—藥物關聯研究方面，Chen等[10]采用矩陣分解的稀疏網絡正則化偏最小二乘法（spare network-regularized partial least square，SNPLS）整合基因表達數據和藥物響應數據，通過計算得到了基因—藥物模塊。
基于此，本文考慮將多重網絡的方法用于基因表達數據和藥物響應數據的整合分析上，但由于這兩種數據之間的異構性，難以整合在一個網絡內，因此，考慮將這兩種數據放在二重異構網絡中，通過多重網絡的方法對異構數據進行處理。首先利用肺癌基因在不同細胞系中的基因表達數據和治療肺癌的藥物在不同細胞系中的藥物響應數據構造出基因表達矩陣和藥物響應矩陣，其次利用皮爾遜相關系數（Pearson correlation coefficient，PCC）構造出基因表達相似網絡和藥物響應相似網絡，再借助已知的基因—藥物關聯，將這兩個網絡連接起來，得到一個異構網絡，然后利用提出的一種基于圖熵的多重網絡模塊鑒別算法（multiplex networks algorithm based on graph entropy，MNA_GE），在該異構網絡中尋找高度關聯的基因—藥物模塊，并對這些模塊進行生物學意義分析和文獻驗證。模塊鑒別實驗處理流程如圖1所示。
圖1
模塊鑒別實驗處理流程
Figure1.
Module identification experimental process
1 材料和方法
1.1 數據集選取
公開使用的腫瘤藥物敏感性基因組學數據庫（genomics of drug sensitivity in cancer，GDSC）（https://www.cancerrxgene.org/）包含了大約20萬個藥物反應實驗產生的數據，這些實驗是在近1 000個人類惡性腫瘤細胞系上測試了266種藥物后得到的。GDSC數據庫主要包括基因組圖譜和藥物響應數據兩大類：① 基因組圖譜：表征了細胞系的基因表達水平，包括基因表達、拷貝數變異和DNA甲基化數據集等，本文使用的是基因表達數據集。② 藥物響應：主要指藥物作用在生物體上的特征反應，包括半數最大抑制濃度、劑量—反應曲線斜率和劑量—反應曲線下的面積（area under curve，AUC）。本文藥物響應數據用的是AUC值，AUC值越低，藥物越有效。
本文從GDSC數據庫下載的原始數據包括肺癌中8 046個基因在734個細胞系中的表達量，以及201種藥物在734個細胞系中用藥后的AUC值。當細胞系存在AUC缺失值時，則剔除該細胞系，最終獲得了肺癌基因在217個細胞系的基因表達矩陣和201種藥物在217個細胞系的藥物響應矩陣。
基因—藥物關聯網絡從公開使用的基因—藥物相互作用數據庫（the drug gene interaction database，DGIdb） (http://www.dgidb.org/)中下載，并過濾掉實驗數據集中未出現的基因。
1.2 相似網絡構建
如圖1所示，對于基因表達和藥物響應數據，本文使用PCC定義基因與基因、藥物與藥物之間的相似度，進一步構造出肺癌的基因表達相似網絡和藥物響應相似網絡。用Pa 和Pb分別表示基因a和基因b在細胞系中的表達量，則基因a與基因b之間的相似度被定義為如式（1）所示：
 |
其中，N代表該基因所有細胞系的數量。兩個基因的PCC值越高，說明兩個基因之間的相似度越高。
在計算出每兩個基因的PCC值后，對所有的PCC值取絕對值，得到基因表達相似矩陣。用
表示n個基因，ARR表示基因相似矩陣。同理可構造藥物相似矩陣，用
表示m種藥物，ADD表示藥物與藥物之間的相似矩陣。通過整合基因表達相似網絡、藥物響應相似網絡和已知的基因—藥物關聯，構造出一個二重異構網絡。在已知的基因—藥物關聯基礎上，利用模塊化鑒別方法進行基因—藥物模塊鑒別。基因—藥物二重異構網絡用一個二向圖G(R, D, E)表示，這里 
包含了所有在基因R和藥物D中存在的連接，則該二重異構網絡的鄰接矩陣A表示為如式（2）所示：
 |
其中，ARD是基因—藥物關聯網絡的鄰接矩陣。在基因—藥物連接中，當基因與藥物存在關聯時，用1表示；反之，當基因與藥物之間不存在關聯時，用0表示，最終獲得ARD矩陣。
1.3 網絡模塊鑒別
為了推斷出基因—藥物模塊，基于Ma等[9]在多層網絡中以圖熵為基礎的模塊鑒別算法，本文提出了MNA_GE算法。該算法包括三個主要步驟：種子節點最優排序、通過種子節點拓展進行模塊搜尋和候選模塊的優化，MNA_GE算法流程圖如圖2所示。
圖2
MNA_GE算法流程圖
Figure2.
MNA_GE algorithm flowchart
因為兩個網絡的規模存在較大差異，所以本文先在基因表達相似網絡中進行種子節點的搜尋，之后再拓展到多重網絡上進行模塊查找和優化。首先在基因表達相似網絡中使用基因的拓撲特征對基因進行最優排序。定義基因表達相似網絡的權重鄰接矩陣 
，這里函數
，fk(i)表示第i個基因在該網絡中的重要性。函數fk定義如式（3）所示：
 '/> |
其中，
表示節點的拓撲重要性，Y是代表先驗信息的向量，參數 
控制節點拓撲重要性和先驗信息之間的關聯權重，
在0~1范圍內進行選取。通過調節 
，可以獲得基因相似網絡中每一個節點的重要性分數。
是一個標準化鄰接矩陣，即 
，這里 
。每一個節點的重要性取決于與它相鄰節點的數量、節點之間的連接強度和它相鄰節點的重要性。式（3）的準確解為 
，但是計算矩陣的逆很耗時，因此使用Zhou等[11]提出的快速迭代算法獲得方程的最優解 
，如式（4）所示：
 '/> |
其中，t代表迭代次數。設 
；如果沒有先驗信息，則設 
。該算法的迭代原理與著名的網頁排序（PageRank）算法類似[12]。
在多次調整參數 
的情況下，發現當 
取0.4時，在網絡中獲取的種子節點數量最多，且最合理。最后，每一個節點周圍相鄰的節點都將被添加到該節點的重要性中。當隨著迭代的進行，
與 
之間不再發生改變時，迭代終止。通常，20次迭代就可以達到預期目標。對每一個fk(i)計算z-分數，計算公式如式（5）所示：
 |
其中，
表示所有 
的均值，
表示所有 
的標準差。利用計算出的 
對所有基因進行排序，
值越高，基因越重要。
從每一個種子節點開始，向種子節點添加基因或者藥物節點信息，它的添加導致目標函數圖熵的最大減少，直到目標函數的圖熵不再發生改變時添加終止。給定一個模塊C，網絡 
的熵函數定義為如式（6）所示：
 |
其中，
。模塊C在兩個網絡中的熵函數定義為如式（7）所示：
 |
模塊搜尋的步驟是使用式（6）拓展模塊。每次添加的節點使得H(C)達到最大程度的下降，直到H(C)不再發生改變。
在優化步驟中，包含節點數小于5的候選模塊被移除。為了合并重疊的功能模塊，如果兩個模塊的杰卡德系數大于0.5，則將這兩個模塊合并。
2 實驗結果
在執行的過程中，首先計算肺癌217個細胞系中8 046個基因的相似矩陣，得到基因表達相似矩陣P8 046 × 8 046，對得到的相似矩陣每個元素取絕對值。再計算217個細胞系下的201種藥物響應的相似矩陣，得到了藥物響應相似矩陣P201 × 201。基因表達相似矩陣的維數是8 046，而藥物響應相似矩陣的維數是201，因此為了將這兩個矩陣匹配在一起，需要把藥物響應相似矩陣的維數也拓展到8 046維。在矩陣拓展的過程中，要滿足兩個條件：兩個矩陣維數相同；添加的新矩陣對原矩陣的影響要盡可能的小。通過多次實驗，最終選取增加一個7 845維的單位矩陣。這種添加方法滿足實驗條件，且使得實驗效果最優。最終將藥物響應相似矩陣拓展到8 046 × 8?046維數大小，方便后續在兩個網絡上進行模塊搜尋。
尋找種子節點時，因為藥物響應相似網絡中具有較多的零值，所以把基因表達相似網絡作為基礎網絡。首先在基因表達相似網絡上尋找種子節點，確定種子節點后在多重網絡上進行種子拓展，尋找到合適的基因—藥物網絡模塊。最后進行模塊優化，除了前面的優化條件外，當模塊中只包含一個基因或者只包含一種藥物時，同樣去除該模塊。
2.1 MNA_GE方法的對比實驗
為了驗證方法的有效性，本文將MNA_GE方法與另外兩種方法進行對比：一種是Cho等[13]提出的基于圖匹配的重新加權隨機游走（reweighted random walks for graph matching，RRWM）方法，另一種是由Chen等[10]提出的SNPLS方法。
首先將三種方法都運用在相同模擬數據集上[14]，分別考慮在添加無關樣本的條件和不同噪聲水平下三種方法在模塊鑒別準確性方面的表現。實驗結果如圖3所示，通過向原始數據添加無關樣本來檢測三種方法在有無關樣本時的模塊鑒別能力。通過調整無關樣本的比例，本文提出的MNA_GE方法在多數情況下優于SNPLS方法，且一直優于RRWM方法。可能的原因是MNA_GE方法在確定種子節點時考慮了先驗信息，這將使得種子節點的排序更加具有準確性，進而在后續種子節點的拓展過程中更具有優勢。雖然增加了無關樣本，但是這對本文方法的影響要小于另外兩種方法。
圖3
準確性和魯棒性測試
Figure3.
Accuracy and robustness test
在比較方法的魯棒性時，通過調整模擬數據集中的噪聲比例，對方法的準確性進行檢驗。可以發現MNA_GE方法對比SNPLS方法有較小的優勢，且明顯優于RRWM方法。這可能是因為RRWM方法在遍歷圖的節點時，容易受到圖中某一個節點的影響，使得方法的準確性改變較大。
從整體來看， MNA_GE方法在樣本中存在無關樣本或噪聲的情況下，均可以取得較好的實驗效果，這證明了方法的有效性。在運行時間方面， MNA_GE方法也得到了大幅提升，明顯優于另外兩種方法。
2.2 肺癌基因—藥物模塊鑒別
在確定了方法的有效性后，將MNA_GE方法應用于肺癌數據集中，通過實驗，共得到了5個基因—藥物模塊。在這5個模塊中，每個模塊包含5種藥物，大約包含145個基因，藥物的具體名稱如表1第3列所示。
表1
模塊的生物功能分析
Table1.
Biological function analysis of the modules
藥物靶點包括受體、酶、離子通道、轉運體、免疫系統、基因等。通過GDSC數據庫查詢模塊中每種藥物的靶點，將查詢到的藥物靶點與鑒別出的5個基因—藥物模塊中的基因進行比對，保留在二者中都出現的靶向基因，如表1第4列所示。蛋白抑制劑莫利瑞司（molibresib）的靶點為含溴結構域3（bromodomain containing 3，BRD3）編碼基因。激酶抑制劑魯索替尼（ruxolitinib）的靶點為詹納斯激酶1（janus kinase1，JAK1）和詹納斯激酶2（janus kinase2，JAK2）。蛋白脫乙酰酶抑制劑伏立諾他（vorinostst）的靶點為組蛋白去乙酰化酶1（histone deacetylases 1，HDAC1）。選擇性激酶抑制劑馬賽替尼（masitinib）的靶點為貓科肉瘤病毒癌基因同源物（v-kit hardy-zuckerman 4 feline sarcoma viral oncogene homolog，KIT）和血小板衍生生長因子受體 
（platelet derived growth factor receptor alpha，PDGFRA）。可逆的表皮生長因子受體（epidermal growth factor receptor，EGFR）抑制劑培利替尼（pelitinib）和針對EGFR的單克隆抗體西妥昔單抗（erbitux）的藥物靶點均為EGFR。這些藥物靶點出現在基因—藥物模塊中，說明MNA_GE方法可以鑒別出已經驗證的藥物靶向基因。
2.2.1 GO富集分析
為了確定基因—藥物模塊的功能相關性，本文測試了每個模塊的基因是否含有特定的生物學功能或信號通路。使用用于注釋、可視化和集成發現的數據庫（the database for annotation, visualization and integrated discovery，DAVID）（https://david.ncifcrf.gov/）工具對篩選出的基因進行基因本體論（gene ontology，GO）富集分析和京都基因與基因組百科全書（Kyoto encyclopedia of genes and genomes，KEGG）富集分析，最后發現所有的基因—藥物模塊都具有至少4個重要的GO生物學功能，如表1所示。
蛋白質結合（相互作用）存在于機體每個細胞的生命活動過程中，生物學中的許多現象如DNA復制、轉錄、翻譯，細胞周期調控和信號傳導等均受蛋白質相互作用的調控，這些信號分子調節的失控將導致惡性腫瘤等疾病的發生。
細胞黏附決定細胞的極性，并且參與正常組織結構的維持[15]。黏附分子E—鈣黏附蛋白（E-cadherin）是具有黏附作用的膜表面糖蛋白，它使細胞與細胞、細胞與基質間形成黏附力，阻止腫瘤細胞從惡性腫瘤組織中脫離出去，抑制腫瘤細胞浸潤。在彌漫性浸潤性腫瘤中，E-cadherin和 
-連環蛋白（
-catenin）、
-連環蛋白（
-catenin）的基因發生了突變。在細胞癌變初期，
-catenin的酪氨酸磷酸化使E-cadherin細胞黏附系統失活，進而導致細胞黏附性降低。細胞間黏附性的降低使腫瘤細胞違背細胞內的正常秩序，導致組織結構的破壞，并且對腫瘤的轉移也會產生影響。由此可見，細胞黏附在細胞癌變的過程中起著重要作用。
2.2.2 KEGG富集分析
在對5個基因—藥物模塊進行KEGG富集分析時，發現每個模塊至少有32%以上的基因數在KEGG上有富集通路，其中出現次數最多的富集通路包括磷脂酰肌醇3激酶—蛋白質絲氨酸蘇氨酸激酶（phosphatidylinositol 3 kinases-protein serine threonine kinase signaling pathway，PI3K-Akt）信號通路和p53信號通路。PI3K-Akt信號通路在促進惡性腫瘤的侵襲和轉移方面具有重要作用[16]。PI3K-Akt信號通路依靠降低腫瘤細胞之間的黏附作用影響腫瘤細胞的轉移。在卵巢癌和前列腺癌中，磷脂酰肌醇3-激酶（phosphpoinositide 3-kinase，PI3K）介導會使E-cadherin表達降低，使得細胞間的黏附性降低，導致腫瘤細胞更容易發生轉移。此外，PI3K-Akt信號通路也可以調控生長因子受體的表達，使腫瘤細胞發生轉移。
p53腫瘤抑制基因可以調節細胞周期和避免細胞發生癌變，通常50%以上的惡性腫瘤患者中p53基因發生了突變和失活。受p53基因調控的p53信號通路在遺傳突變中同樣具有重要地位。通過KEGG富集分析可以發現這5個基因—藥物模塊中的基因在影響人類惡性腫瘤產生的信號通路上明顯富集，驗證了基因—藥物模塊具有明顯的生物學意義。
2.3 文獻驗證
如表1所示，模塊4和5的藥物多西他賽對非小細胞肺癌有較好的療效，這已經在臨床得到驗證。多西他賽是紫杉類藥物，它與微管蛋白結合的能力比紫杉醇高2倍。多西他賽與紫杉醇有不完全交叉耐藥性，紫杉醇治療對患者無效后，多西他賽治療對患者依然有效。這證明了通過MNA_GE方法鑒別出的基因—藥物模塊所包含的藥物可以治療肺癌疾病。
曲美替尼在5個模塊中出現了3次，西妥昔單抗在5個模塊中出現了2次，這兩個藥物均有較高的出現頻率。西妥昔單抗已經作為結直腸癌治療的重要分子靶向藥物在臨床上廣泛應用，但是半數以上的結直腸癌患者因鼠類肉瘤病毒原癌基因（kirsten rat sarcoma viral proto-oncogene，KRAS）、絲氨酸蘇氨酸激酶原癌基因（B-Raf proto-oncogene serine threonine kinase，BRAF）、PI3K基因突變而對西妥昔單抗耐藥。此外，絲裂原活化蛋白激酶激酶（mitogen-activated proteinkinase kinase，MEK）靶點抑制劑曲美替尼在大多數惡性腫瘤患者的治療初期可以產生良好的效果，但是很快就會由于各種耐藥機制而導致藥物療效減退甚至無效[17]。聯合用藥是降低耐藥率，提高疾病治療效果的一種較為理想的治療策略。已有的研究結果顯示，EGFR抑制劑西妥昔單抗聯合MEK靶點抑制劑曲美替尼可以很好地改善絲裂原活化蛋白激酶（mitogen-activated protein kinase，MAPK）通路抑制，提高KRAS、BRAF基因突變型人結直腸癌細胞系對西妥昔單抗的藥物敏感性，增加BRAF基因第600位密碼子纈氨酸（V）被谷氨酸（E）取代（BRAF V600E）突變型結直腸癌的反應率[18-19]。聯合用藥后細胞增殖較西妥昔單抗組受到的抑制更明顯，同時會進一步誘導KRAS、BRAF突變型結直腸癌細胞的凋亡。這個實驗證明，包含基因KRAS、BRAF與藥物曲美替尼、西妥昔單抗的模塊具有組合的優勢，這兩種藥物聯合用藥對兩個基因的作用效果要優于單一用藥。
3 討論
本文首先從GDSC數據庫下載了肺癌的基因表達數據和相應的藥物響應數據，利用PCC構建了基因表達相似網絡和藥物響應相似網絡，借助已知的基因—藥物關聯將這兩個相似網絡關聯起來，然后利用MNA_GE方法在該異構網絡中進行模塊的鑒別，得到5個肺癌基因—藥物關聯模塊。通過在增加無關樣本和噪聲的情況下，與另外兩種方法對比，驗證了MNA_GE方法具有更好的準確性和魯棒性。通過對得到的基因—藥物模塊進行生物功能富集分析和文獻挖掘，驗證了模塊的生物學意義。
現在對復雜疾病的研究中，對于基因—疾病關聯的研究是比較多的，但是由于藥物與惡性腫瘤基因之間的作用關系比較復雜，導致對基因—藥物關聯的研究相對較少。本文利用提出的多重網絡算法處理了難以整合的異構數據，研究基因和藥物之間的關聯，取得了一定的效果，后續可以在這方面做進一步的研究。由于肺癌的患病基數大，數據來源廣泛，因此本文選擇了肺癌進行研究。當其他疾病與肺癌有類似的數據集時，該研究方法也同樣適用于其他疾病。本文選擇的數據集中基因的數量遠大于藥物的數量，導致基因表達相似矩陣和藥物響應相似矩陣維數差異較大，對實驗結果的準確性會產生一定的影響，后續會加以改進。
利益沖突聲明：本文全體作者均聲明不存在利益沖突。
引言
新藥物的研發是一個具有風險、周期長且成本高的過程[1-2]。盡管最近幾十年醫藥公司在藥物研發上的投入已經增加了很多，但是被美國食品藥品監督管理局（food and drug administration，FDA）認可的新藥物數量卻不是很多。為了解決這個難題，對已經存在的藥物進行重新開發利用是一項非常有前景的工作。
惡性腫瘤嚴重影響人類生命健康，其中肺癌是發病率和死亡率增長最快，對人體健康和生命威脅最大的惡性腫瘤之一。據報道，近50年來許多國家肺癌的發病率和死亡率均明顯增高，其中男性的肺癌發病率和死亡率均占所有惡性腫瘤的第一位，而女性的肺癌發病率和死亡率均占所有惡性腫瘤第二位。惡性腫瘤的發生與基因突變密切相關，最近的研究表明，每一個腫瘤患者都擁有一個特殊的基因表達圖譜，但是患同一種惡性腫瘤的不同患者，他們的基因表達圖譜卻又具有一定的相似性，即相似的基因突變可能產生同一種疾病，反之亦然。研究表明，用于治愈同一種疾病的不同藥物具有相似的化學結構，且具有相似基因表達圖譜的細胞系對同一種藥物具有相似的化學反應[3-5]。因此，可以通過研究同一種惡性腫瘤的基因表達圖譜和治療該惡性腫瘤藥物之間的模塊化結構，鑒別新的基因—藥物關聯，發現新的藥物靶向基因，進一步將已經存在的藥物用于其它惡性腫瘤治療中。這將在惡性腫瘤的精準治療方面發揮巨大的作用，增加惡性腫瘤治愈的可能性，減少新藥物研發的成本。
高通量測序技術的發展產生了豐富的多組學數據，為惡性腫瘤致病機制的研究提供了極大的便利。不少學者已經提出了各種方法來對惡性腫瘤的多組學數據進行分析和研究，其中從網絡的角度整合多組學數據，對惡性腫瘤進行研究，是現在比較流行的一種分析方法。多重網絡可以充分整合某一惡性腫瘤的多組學數據，與單一網絡相比，多重網絡可以更全面地描述與惡性腫瘤有關的信息。2019年，Yuan等[6]通過整合乳腺癌的DNA甲基化數據、拷貝數變異數據和二者在基因調控網絡中的關聯，使用雙權相關系數衡量基因調控網絡中的生物分子關聯，并且利用非凸懲罰函數使網絡模型達到最優，最終在基因調控網絡中推斷出與乳腺癌發生相關的DNA甲基化位點和拷貝數變異位點。Valdeolivas等[7]利用重啟隨機游走的方法在由蛋白質、基因和疾病構成的多重異構網絡上進行研究，結果證明了整合的多組學數據可以有效提高實驗結果的準確性。Shang等[8]基于受約束的網頁排序（constrained PageRank，CPR）算法，將肝細胞癌的多組學數據整合在一個統一的網絡框架中，通過實驗發現，與單一數據集相比，使用多組學數據進行網絡構造時取得的實驗效果明顯優于單一數據集。Ma等[9]利用全基因組的DNA甲基化數據和基因表達數據，分別構造出基因差異共表達網絡和基因差異共甲基化網絡，然后在這兩個網絡上利用基于圖熵的方法發現表觀遺傳模塊，該模塊可以作為生物標志物用于疾病的診斷和治療。由此可見，在網絡中利用多組學數據進行生物分子關聯的鑒別是可行且具有巨大發展潛力的。而在基因—藥物關聯研究方面，Chen等[10]采用矩陣分解的稀疏網絡正則化偏最小二乘法（spare network-regularized partial least square，SNPLS）整合基因表達數據和藥物響應數據，通過計算得到了基因—藥物模塊。
基于此，本文考慮將多重網絡的方法用于基因表達數據和藥物響應數據的整合分析上，但由于這兩種數據之間的異構性，難以整合在一個網絡內，因此，考慮將這兩種數據放在二重異構網絡中，通過多重網絡的方法對異構數據進行處理。首先利用肺癌基因在不同細胞系中的基因表達數據和治療肺癌的藥物在不同細胞系中的藥物響應數據構造出基因表達矩陣和藥物響應矩陣，其次利用皮爾遜相關系數（Pearson correlation coefficient，PCC）構造出基因表達相似網絡和藥物響應相似網絡，再借助已知的基因—藥物關聯，將這兩個網絡連接起來，得到一個異構網絡，然后利用提出的一種基于圖熵的多重網絡模塊鑒別算法（multiplex networks algorithm based on graph entropy，MNA_GE），在該異構網絡中尋找高度關聯的基因—藥物模塊，并對這些模塊進行生物學意義分析和文獻驗證。模塊鑒別實驗處理流程如圖1所示。
圖1
模塊鑒別實驗處理流程
Figure1.
Module identification experimental process
1 材料和方法
1.1 數據集選取
公開使用的腫瘤藥物敏感性基因組學數據庫（genomics of drug sensitivity in cancer，GDSC）（https://www.cancerrxgene.org/）包含了大約20萬個藥物反應實驗產生的數據，這些實驗是在近1 000個人類惡性腫瘤細胞系上測試了266種藥物后得到的。GDSC數據庫主要包括基因組圖譜和藥物響應數據兩大類：① 基因組圖譜：表征了細胞系的基因表達水平，包括基因表達、拷貝數變異和DNA甲基化數據集等，本文使用的是基因表達數據集。② 藥物響應：主要指藥物作用在生物體上的特征反應，包括半數最大抑制濃度、劑量—反應曲線斜率和劑量—反應曲線下的面積（area under curve，AUC）。本文藥物響應數據用的是AUC值，AUC值越低，藥物越有效。
本文從GDSC數據庫下載的原始數據包括肺癌中8 046個基因在734個細胞系中的表達量，以及201種藥物在734個細胞系中用藥后的AUC值。當細胞系存在AUC缺失值時，則剔除該細胞系，最終獲得了肺癌基因在217個細胞系的基因表達矩陣和201種藥物在217個細胞系的藥物響應矩陣。
基因—藥物關聯網絡從公開使用的基因—藥物相互作用數據庫（the drug gene interaction database，DGIdb） (http://www.dgidb.org/)中下載，并過濾掉實驗數據集中未出現的基因。
1.2 相似網絡構建
如圖1所示，對于基因表達和藥物響應數據，本文使用PCC定義基因與基因、藥物與藥物之間的相似度，進一步構造出肺癌的基因表達相似網絡和藥物響應相似網絡。用Pa 和Pb分別表示基因a和基因b在細胞系中的表達量，則基因a與基因b之間的相似度被定義為如式（1）所示：
|
其中，N代表該基因所有細胞系的數量。兩個基因的PCC值越高，說明兩個基因之間的相似度越高。
在計算出每兩個基因的PCC值后，對所有的PCC值取絕對值，得到基因表達相似矩陣。用表示n個基因，ARR表示基因相似矩陣。同理可構造藥物相似矩陣，用表示m種藥物，ADD表示藥物與藥物之間的相似矩陣。通過整合基因表達相似網絡、藥物響應相似網絡和已知的基因—藥物關聯，構造出一個二重異構網絡。在已知的基因—藥物關聯基礎上，利用模塊化鑒別方法進行基因—藥物模塊鑒別。基因—藥物二重異構網絡用一個二向圖G(R, D, E)表示，這里 包含了所有在基因R和藥物D中存在的連接，則該二重異構網絡的鄰接矩陣A表示為如式（2）所示：
|
其中，ARD是基因—藥物關聯網絡的鄰接矩陣。在基因—藥物連接中，當基因與藥物存在關聯時，用1表示；反之，當基因與藥物之間不存在關聯時，用0表示，最終獲得ARD矩陣。
1.3 網絡模塊鑒別
為了推斷出基因—藥物模塊，基于Ma等[9]在多層網絡中以圖熵為基礎的模塊鑒別算法，本文提出了MNA_GE算法。該算法包括三個主要步驟：種子節點最優排序、通過種子節點拓展進行模塊搜尋和候選模塊的優化，MNA_GE算法流程圖如圖2所示。
圖2
MNA_GE算法流程圖
Figure2.
MNA_GE algorithm flowchart
因為兩個網絡的規模存在較大差異，所以本文先在基因表達相似網絡中進行種子節點的搜尋，之后再拓展到多重網絡上進行模塊查找和優化。首先在基因表達相似網絡中使用基因的拓撲特征對基因進行最優排序。定義基因表達相似網絡的權重鄰接矩陣 ，這里函數，fk(i)表示第i個基因在該網絡中的重要性。函數fk定義如式（3）所示：
| '/> |
其中， 表示節點的拓撲重要性，Y是代表先驗信息的向量，參數 控制節點拓撲重要性和先驗信息之間的關聯權重， 在0~1范圍內進行選取。通過調節 ，可以獲得基因相似網絡中每一個節點的重要性分數。 是一個標準化鄰接矩陣，即 ，這里 。每一個節點的重要性取決于與它相鄰節點的數量、節點之間的連接強度和它相鄰節點的重要性。式（3）的準確解為 ，但是計算矩陣的逆很耗時，因此使用Zhou等[11]提出的快速迭代算法獲得方程的最優解 ，如式（4）所示：
| '/> |
其中，t代表迭代次數。設 ；如果沒有先驗信息，則設 。該算法的迭代原理與著名的網頁排序（PageRank）算法類似[12]。
在多次調整參數 的情況下，發現當 取0.4時，在網絡中獲取的種子節點數量最多，且最合理。最后，每一個節點周圍相鄰的節點都將被添加到該節點的重要性中。當隨著迭代的進行， 與 之間不再發生改變時，迭代終止。通常，20次迭代就可以達到預期目標。對每一個fk(i)計算z-分數，計算公式如式（5）所示：
|
其中， 表示所有 的均值， 表示所有 的標準差。利用計算出的 對所有基因進行排序， 值越高，基因越重要。
從每一個種子節點開始，向種子節點添加基因或者藥物節點信息，它的添加導致目標函數圖熵的最大減少，直到目標函數的圖熵不再發生改變時添加終止。給定一個模塊C，網絡 的熵函數定義為如式（6）所示：
|
其中，。模塊C在兩個網絡中的熵函數定義為如式（7）所示：
|
模塊搜尋的步驟是使用式（6）拓展模塊。每次添加的節點使得H(C)達到最大程度的下降，直到H(C)不再發生改變。
在優化步驟中，包含節點數小于5的候選模塊被移除。為了合并重疊的功能模塊，如果兩個模塊的杰卡德系數大于0.5，則將這兩個模塊合并。
2 實驗結果
在執行的過程中，首先計算肺癌217個細胞系中8 046個基因的相似矩陣，得到基因表達相似矩陣P8 046 × 8 046，對得到的相似矩陣每個元素取絕對值。再計算217個細胞系下的201種藥物響應的相似矩陣，得到了藥物響應相似矩陣P201 × 201。基因表達相似矩陣的維數是8 046，而藥物響應相似矩陣的維數是201，因此為了將這兩個矩陣匹配在一起，需要把藥物響應相似矩陣的維數也拓展到8 046維。在矩陣拓展的過程中，要滿足兩個條件：兩個矩陣維數相同；添加的新矩陣對原矩陣的影響要盡可能的小。通過多次實驗，最終選取增加一個7 845維的單位矩陣。這種添加方法滿足實驗條件，且使得實驗效果最優。最終將藥物響應相似矩陣拓展到8 046 × 8?046維數大小，方便后續在兩個網絡上進行模塊搜尋。
尋找種子節點時，因為藥物響應相似網絡中具有較多的零值，所以把基因表達相似網絡作為基礎網絡。首先在基因表達相似網絡上尋找種子節點，確定種子節點后在多重網絡上進行種子拓展，尋找到合適的基因—藥物網絡模塊。最后進行模塊優化，除了前面的優化條件外，當模塊中只包含一個基因或者只包含一種藥物時，同樣去除該模塊。
2.1 MNA_GE方法的對比實驗
為了驗證方法的有效性，本文將MNA_GE方法與另外兩種方法進行對比：一種是Cho等[13]提出的基于圖匹配的重新加權隨機游走（reweighted random walks for graph matching，RRWM）方法，另一種是由Chen等[10]提出的SNPLS方法。
首先將三種方法都運用在相同模擬數據集上[14]，分別考慮在添加無關樣本的條件和不同噪聲水平下三種方法在模塊鑒別準確性方面的表現。實驗結果如圖3所示，通過向原始數據添加無關樣本來檢測三種方法在有無關樣本時的模塊鑒別能力。通過調整無關樣本的比例，本文提出的MNA_GE方法在多數情況下優于SNPLS方法，且一直優于RRWM方法。可能的原因是MNA_GE方法在確定種子節點時考慮了先驗信息，這將使得種子節點的排序更加具有準確性，進而在后續種子節點的拓展過程中更具有優勢。雖然增加了無關樣本，但是這對本文方法的影響要小于另外兩種方法。
圖3
準確性和魯棒性測試
Figure3.
Accuracy and robustness test
在比較方法的魯棒性時，通過調整模擬數據集中的噪聲比例，對方法的準確性進行檢驗。可以發現MNA_GE方法對比SNPLS方法有較小的優勢，且明顯優于RRWM方法。這可能是因為RRWM方法在遍歷圖的節點時，容易受到圖中某一個節點的影響，使得方法的準確性改變較大。
從整體來看， MNA_GE方法在樣本中存在無關樣本或噪聲的情況下，均可以取得較好的實驗效果，這證明了方法的有效性。在運行時間方面， MNA_GE方法也得到了大幅提升，明顯優于另外兩種方法。
2.2 肺癌基因—藥物模塊鑒別
在確定了方法的有效性后，將MNA_GE方法應用于肺癌數據集中，通過實驗，共得到了5個基因—藥物模塊。在這5個模塊中，每個模塊包含5種藥物，大約包含145個基因，藥物的具體名稱如表1第3列所示。
表1
模塊的生物功能分析
Table1.
Biological function analysis of the modules
藥物靶點包括受體、酶、離子通道、轉運體、免疫系統、基因等。通過GDSC數據庫查詢模塊中每種藥物的靶點，將查詢到的藥物靶點與鑒別出的5個基因—藥物模塊中的基因進行比對，保留在二者中都出現的靶向基因，如表1第4列所示。蛋白抑制劑莫利瑞司（molibresib）的靶點為含溴結構域3（bromodomain containing 3，BRD3）編碼基因。激酶抑制劑魯索替尼（ruxolitinib）的靶點為詹納斯激酶1（janus kinase1，JAK1）和詹納斯激酶2（janus kinase2，JAK2）。蛋白脫乙酰酶抑制劑伏立諾他（vorinostst）的靶點為組蛋白去乙酰化酶1（histone deacetylases 1，HDAC1）。選擇性激酶抑制劑馬賽替尼（masitinib）的靶點為貓科肉瘤病毒癌基因同源物（v-kit hardy-zuckerman 4 feline sarcoma viral oncogene homolog，KIT）和血小板衍生生長因子受體 （platelet derived growth factor receptor alpha，PDGFRA）。可逆的表皮生長因子受體（epidermal growth factor receptor，EGFR）抑制劑培利替尼（pelitinib）和針對EGFR的單克隆抗體西妥昔單抗（erbitux）的藥物靶點均為EGFR。這些藥物靶點出現在基因—藥物模塊中，說明MNA_GE方法可以鑒別出已經驗證的藥物靶向基因。
2.2.1 GO富集分析
為了確定基因—藥物模塊的功能相關性，本文測試了每個模塊的基因是否含有特定的生物學功能或信號通路。使用用于注釋、可視化和集成發現的數據庫（the database for annotation, visualization and integrated discovery，DAVID）（https://david.ncifcrf.gov/）工具對篩選出的基因進行基因本體論（gene ontology，GO）富集分析和京都基因與基因組百科全書（Kyoto encyclopedia of genes and genomes，KEGG）富集分析，最后發現所有的基因—藥物模塊都具有至少4個重要的GO生物學功能，如表1所示。
蛋白質結合（相互作用）存在于機體每個細胞的生命活動過程中，生物學中的許多現象如DNA復制、轉錄、翻譯，細胞周期調控和信號傳導等均受蛋白質相互作用的調控，這些信號分子調節的失控將導致惡性腫瘤等疾病的發生。
細胞黏附決定細胞的極性，并且參與正常組織結構的維持[15]。黏附分子E—鈣黏附蛋白（E-cadherin）是具有黏附作用的膜表面糖蛋白，它使細胞與細胞、細胞與基質間形成黏附力，阻止腫瘤細胞從惡性腫瘤組織中脫離出去，抑制腫瘤細胞浸潤。在彌漫性浸潤性腫瘤中，E-cadherin和 -連環蛋白（-catenin）、-連環蛋白（-catenin）的基因發生了突變。在細胞癌變初期，-catenin的酪氨酸磷酸化使E-cadherin細胞黏附系統失活，進而導致細胞黏附性降低。細胞間黏附性的降低使腫瘤細胞違背細胞內的正常秩序，導致組織結構的破壞，并且對腫瘤的轉移也會產生影響。由此可見，細胞黏附在細胞癌變的過程中起著重要作用。
2.2.2 KEGG富集分析
在對5個基因—藥物模塊進行KEGG富集分析時，發現每個模塊至少有32%以上的基因數在KEGG上有富集通路，其中出現次數最多的富集通路包括磷脂酰肌醇3激酶—蛋白質絲氨酸蘇氨酸激酶（phosphatidylinositol 3 kinases-protein serine threonine kinase signaling pathway，PI3K-Akt）信號通路和p53信號通路。PI3K-Akt信號通路在促進惡性腫瘤的侵襲和轉移方面具有重要作用[16]。PI3K-Akt信號通路依靠降低腫瘤細胞之間的黏附作用影響腫瘤細胞的轉移。在卵巢癌和前列腺癌中，磷脂酰肌醇3-激酶（phosphpoinositide 3-kinase，PI3K）介導會使E-cadherin表達降低，使得細胞間的黏附性降低，導致腫瘤細胞更容易發生轉移。此外，PI3K-Akt信號通路也可以調控生長因子受體的表達，使腫瘤細胞發生轉移。
p53腫瘤抑制基因可以調節細胞周期和避免細胞發生癌變，通常50%以上的惡性腫瘤患者中p53基因發生了突變和失活。受p53基因調控的p53信號通路在遺傳突變中同樣具有重要地位。通過KEGG富集分析可以發現這5個基因—藥物模塊中的基因在影響人類惡性腫瘤產生的信號通路上明顯富集，驗證了基因—藥物模塊具有明顯的生物學意義。
2.3 文獻驗證
如表1所示，模塊4和5的藥物多西他賽對非小細胞肺癌有較好的療效，這已經在臨床得到驗證。多西他賽是紫杉類藥物，它與微管蛋白結合的能力比紫杉醇高2倍。多西他賽與紫杉醇有不完全交叉耐藥性，紫杉醇治療對患者無效后，多西他賽治療對患者依然有效。這證明了通過MNA_GE方法鑒別出的基因—藥物模塊所包含的藥物可以治療肺癌疾病。
曲美替尼在5個模塊中出現了3次，西妥昔單抗在5個模塊中出現了2次，這兩個藥物均有較高的出現頻率。西妥昔單抗已經作為結直腸癌治療的重要分子靶向藥物在臨床上廣泛應用，但是半數以上的結直腸癌患者因鼠類肉瘤病毒原癌基因（kirsten rat sarcoma viral proto-oncogene，KRAS）、絲氨酸蘇氨酸激酶原癌基因（B-Raf proto-oncogene serine threonine kinase，BRAF）、PI3K基因突變而對西妥昔單抗耐藥。此外，絲裂原活化蛋白激酶激酶（mitogen-activated proteinkinase kinase，MEK）靶點抑制劑曲美替尼在大多數惡性腫瘤患者的治療初期可以產生良好的效果，但是很快就會由于各種耐藥機制而導致藥物療效減退甚至無效[17]。聯合用藥是降低耐藥率，提高疾病治療效果的一種較為理想的治療策略。已有的研究結果顯示，EGFR抑制劑西妥昔單抗聯合MEK靶點抑制劑曲美替尼可以很好地改善絲裂原活化蛋白激酶（mitogen-activated protein kinase，MAPK）通路抑制，提高KRAS、BRAF基因突變型人結直腸癌細胞系對西妥昔單抗的藥物敏感性，增加BRAF基因第600位密碼子纈氨酸（V）被谷氨酸（E）取代（BRAF V600E）突變型結直腸癌的反應率[18-19]。聯合用藥后細胞增殖較西妥昔單抗組受到的抑制更明顯，同時會進一步誘導KRAS、BRAF突變型結直腸癌細胞的凋亡。這個實驗證明，包含基因KRAS、BRAF與藥物曲美替尼、西妥昔單抗的模塊具有組合的優勢，這兩種藥物聯合用藥對兩個基因的作用效果要優于單一用藥。
3 討論
本文首先從GDSC數據庫下載了肺癌的基因表達數據和相應的藥物響應數據，利用PCC構建了基因表達相似網絡和藥物響應相似網絡，借助已知的基因—藥物關聯將這兩個相似網絡關聯起來，然后利用MNA_GE方法在該異構網絡中進行模塊的鑒別，得到5個肺癌基因—藥物關聯模塊。通過在增加無關樣本和噪聲的情況下，與另外兩種方法對比，驗證了MNA_GE方法具有更好的準確性和魯棒性。通過對得到的基因—藥物模塊進行生物功能富集分析和文獻挖掘，驗證了模塊的生物學意義。
現在對復雜疾病的研究中，對于基因—疾病關聯的研究是比較多的，但是由于藥物與惡性腫瘤基因之間的作用關系比較復雜，導致對基因—藥物關聯的研究相對較少。本文利用提出的多重網絡算法處理了難以整合的異構數據，研究基因和藥物之間的關聯，取得了一定的效果，后續可以在這方面做進一步的研究。由于肺癌的患病基數大，數據來源廣泛，因此本文選擇了肺癌進行研究。當其他疾病與肺癌有類似的數據集時，該研究方法也同樣適用于其他疾病。本文選擇的數據集中基因的數量遠大于藥物的數量，導致基因表達相似矩陣和藥物響應相似矩陣維數差異較大，對實驗結果的準確性會產生一定的影響，后續會加以改進。
利益沖突聲明：本文全體作者均聲明不存在利益沖突。
