趙燕偉，王振興，王熙梓，賈曉東，李明君

(聊城市人民醫(yī)院,山東 聊城 252000)

食管癌是最常見、最具侵襲性的惡性腫瘤之一。每年全世界大約有40萬名新確診患者并且約有30萬患者死亡[1]。據美國國家癌癥研究院報道的數據顯示食管癌患者的五年生存率僅為19.9%。近年來食管癌的治療效果得到了很大改善，但食管癌的預后仍然不能令人滿意。化療和基于放療的治療方法已應用于許多患者，根據目前醫(yī)學研究的統(tǒng)計分析，食管癌手術前的新輔助化療或聯合放化療的患者比單獨手術的患者的生存期會有所增加[2]。但是，對患者治療的方法是因人而異的，有些患者的化療效果很好，有些患者可能效果較差。所以確定哪些人適合使用哪種治療手段就顯得十分重要，患者也能避免嘗試不同無效的治療手段時帶來的副作用的傷害。所以，可以預測化療反應的分子標記物需要被更多的研究，找到適當的分子標記物，更好的識別哪些藥物對患者有效，哪些藥物患者會產生耐藥性，從而在克服潛在的化療耐藥性的情況下提供新的治療手段。

微小核糖核酸(miRNA)是一種非編碼RNA，和轉錄因子功能相似，可以調節(jié)信使RNA(mRNA)的表達，miRNA起到的是抑制基因表達的功能，從而實現對各種生物學功能的調控，目前很多研究已經發(fā)現miRNA和轉錄因子在癌癥的發(fā)生發(fā)展中起到關鍵的作用[3-4]。miRNA和轉錄因子作為生物標記物的研究也很多，研究人員發(fā)現在食管癌患者和正常樣本中，多個miRNA表達異常，這些表達的差異直接影響了腫瘤的生長、增值、侵襲，miR-485-5p就可以使乙酰氨基葡糖轉移酶水平下降從而抑制食管癌的增殖和侵襲，miR-502通過促進AKT磷酸化介導食管癌細胞增殖[5-6]。此外還有研究發(fā)現miRNA和轉錄因子的異常表達還影響著化學療法的療效，miR-205納米制劑可以使得前列腺癌對化療敏感，miR-146a通過下調CHOP的表達誘導機體耐藥[7-9]。在多種癌癥的研究中有轉錄因子和miRNA作為標志物的相關報道眾多，miR-590-3p通過抑制Wnt/beta-catenin通路內的WIF1和DKK1基因促進了結直腸癌的增殖，轉錄因子ZNF217是乳腺癌進展過程中的預后生物標志物和治療靶點[10-12]，但是到目前為止miRNA是如何影響著食管癌對藥物的耐藥性的研究還很少。

本研究中我們整合識別了對不同藥物敏感及耐藥的細胞系，找到了在敏感和耐藥細胞系中起到關鍵作用的基因集合，之后我們整合了多個轉錄因子 及miRNA對靶基因調控關系數據庫，構建了包含轉錄因子及miRNA的轉錄調控網絡，這些調控關系是實驗證實的，之后為了確定這些調控關系在食管癌中是被真正的激活，我們使用了食管癌患者的數據對網絡內每一對關系進行相關性計算，得到在食管癌中顯著相關的關系對，篩選后的網絡進一步進行分析。使用不同藥物作用的下耐藥和敏感細胞系進行差異分析得到的差異表達基因，篩選出藥物特異的調控模塊，從根本上找到那些影響細胞系耐藥相關的基因是被哪些轉錄因子或者miRNA調控的，這些轉錄因子及miRNA可能成為新的藥物耐藥性預測的標記物，為后續(xù)的藥物研發(fā)及患者的治療提供幫助。

1 材料方法

1.1 敏感及耐藥細胞系識別

從癌癥藥物敏感性基因組學數據庫獲(GDSC)[11]得了全部的食管癌細胞系表達數據及藥物作用后的IC50值數據。然后分別對每一種藥物治療的所有細胞系進行劃分，分為敏感和耐藥細胞系。一個藥物作用于多個細胞系，每個細胞系的IC50值有所差異，通過計算同種藥物作用下不同細胞系IC50的均值及標準差，確定了分類標準，細胞系的IC50值大于均值加上0.8倍的標準差為耐藥細胞系，IC50值小于均值減去0.8倍的標準差為敏感細胞系。為了保證后續(xù)計算的統(tǒng)計學準確性，要求耐藥和敏感細胞系的數量均不小于5個，根據兩個測度識別了食管癌敏感及耐藥細胞系。

1.2 耐藥性顯著相關基因集合

通過識別得到的多組藥物相關的敏感和耐藥細胞系進行差異表達計算，使用的方法為t檢驗，選取顯著性結果的閾值為P< 0.01，得到每組藥物的耐藥性相關基因集合。對全部13種藥物的分類完成的細胞系進行了相關計算，得到13組差異表達基因。

1.3 食管癌轉錄調控網絡構建

為了識別出哪些miRNA和轉錄因子調控藥物耐藥及功能相關的基因，從多個數據庫整理了實驗證實的轉錄調控關系，miRNA調控靶基因的關系對來自miRTarBase[13]、TarBase[14]，轉錄因子對靶基因及miRNA的調控關系來自于TRANSFAC[15]，還整合了單獨的轉錄因子對miRNA調控數據庫TransmiR[16]，所有的調控關系都是經過實驗證實的，初始的轉錄調控網絡內一共得到5 265條調控關系對，2 772個節(jié)點。

從TCGA數據庫下載了食管癌的mRNA和miRNA表達數據[17]，首先對下載的數據的miRNA名稱注釋為成熟miRNA名稱，基因的ENSG名稱注釋成標準的基因縮寫名稱。使用得到的食管癌數據對網絡內的每一對關系進行相關性計算，計算方法選擇的是皮爾森相關，相關閾值選取P< 0.05。進行篩選后得到的網絡稱之為食管癌特異的轉錄調控網絡。網絡內剩余2 445對調控關系，1 612個節(jié)點。

1.4 耐藥相關模塊挖掘

為了得到各種藥物耐藥相關的轉錄調控模塊，將敏感和耐藥細胞系計算得到的差異表達基因投入網絡，根據關聯有罪理論，以及盡可能的找到網絡內調控種子基因的關鍵基因，本研究以差異基因作為種子基因，進行擴一步搜索其相關調控基因，獲取了藥物特異的調控網絡。網絡并非為全聯通網絡，由多個分離的組分構成，由于網絡內節(jié)點度比較高的基因起到重要作用，這里選取模塊內的最大組分，最大組分的模塊內包含最核心的基因以及最相比其他組分最多的節(jié)點數量，包含更為豐富的信息，最后剔除其余散點，這一模塊稱為相關藥物的耐藥模塊。

1.5 功能富集分析及藥物作用相關通路檢驗分析

功能富集分析使用的是在線的網站工具DAVID(Database for Annotation, Visualization, and Integrated Discovery)[5]，將模塊內的miRNA剔除，使用模塊內全部的轉錄因子及靶基因進行富集分析，分析GO_BP功能及KEGG通路的富集結果。

2 結果分析

2.1 食管癌特異轉錄調控網絡

通過整合多個數據庫內的miRNA及轉錄因子對靶基因的調控關系，構建了所有關系均為實驗證實的轉錄調控網。因為很多基因在不同組織內是特異性表達的，所以想通過真實的食管癌數據對調控關系進行驗證，下載了TCGA的173例食管癌患者腫瘤數據，對所有的調控關系對進行相關性計算，使用的方法為Pearson相關，經過計算選取P值顯著的關系對作為食管癌中能行使特定調控功能的關系對，由這些關系所組成的網絡稱之為食管癌特異的轉錄調控網絡。我們可以基于這一網絡去探索網絡內某個重要基因的調控關系，去進一步研究該基因之所以能行駛功能是不是調控了其他基因導致的，或者從另一角度找到該基因的表達值改變的原因，是否是由于受其它基因的調控。構建了這樣的食管癌特異的網絡對后續(xù)的關鍵調控子相關調控機制的研究提供了便利(見圖1)。

注：(a)中網絡內橙色的點為miRNA，藍色的點為轉錄因子，綠色點為靶基因.

2.2 食管癌耐藥模塊

通過對藥物作用的食管癌細胞系進行區(qū)分，將同一種藥物作用下不同耐藥能力的細胞系進行區(qū)分，計算同種藥物作用下耐藥和敏感細胞系的差異表達基因，這些基因在藥物作用的時候使細胞系產生耐藥反應，從而使得藥效降低。共篩選了滿足分析條件的食管癌藥物13種(見表1)。

表1 13種藥物中敏感耐藥細胞系數量Table 1 The number of sensitive resistant cell lines in 13 drugs

其中包含Trametinib、Topotecan、AZD7762和Afatinib等在內的10種藥物均顯示已經應用于食管癌患者的治療[18-26]，其余3種藥物Teniposide、Dasatinib和Dactolisib也有研究發(fā)現其潛在的食管癌治療的價值，而且發(fā)現了不同患者使用該藥物時敏感性有一定差異[27-29]。對每一種藥物進行差異表達分析，獲得了13組差異表達基因集合。我們發(fā)現在不同藥物的作用下，部分基因在兩組樣本中差異表達頻率很高，這些基因可能影響多種藥物在食管癌中的耐藥性。尤其是FANCD2OS、GUCY2F、HES7三個基因，分析得知其在不同的6種藥物中均顯著的差異表達，通過文獻查找也發(fā)現在以往的研究及試驗中證實，這些基因發(fā)生突變和多個癌癥的發(fā)生發(fā)展顯著相關，而且影響患者的預后[30-33]。繼續(xù)對這些耐藥相關的基因的研究，從而發(fā)現是否是因為這些耐藥相關基因發(fā)生了突變的患者預后較差的原因是否為對治療藥物產生了耐藥性導致的。

這些耐藥相關的基因是否受其他基因調控，或者是否依靠調控其他基因的表達從而真正實現患者耐藥性的產生，將差異表達基因映射到食管癌特異的調控網絡內，從網絡內提取出和投入的種子基因最臨近的基因集合，將篩選后的最大聯調網絡提取出來，作為食管癌各種藥物特異的調控模塊，選擇網絡內的最大連通組分，既保留了網絡內的關鍵調控因子，也盡可能的保留了食管癌耐藥相關的差異表達基因，從而能繼續(xù)分析模塊的功能，以及提取每一步的調控關系，找到關鍵的調控因子。根據上文體積的方法，最終挖掘出13種食管癌耐藥模塊(見圖2)。

圖2 13種藥物耐藥特異模塊Fig.2 Resistance-related characteristic modules of 13 drugs

因為每個藥物耐藥基因不同，模塊的大小也有所差異。發(fā)現大多數模塊內，耐藥相關基因都位于模塊的下游，也就是這些耐藥相關基因的表達受到了一系列轉錄因子及miRNA的調控，上游的一個或幾個關鍵的調控子發(fā)生改變，可能導致耐藥相關基因的異常表達的出現。

挖掘后的模塊進行功能分析，分析每個模塊所行使的功能，我們使用是在線的功能注釋工具DAVID，通過分析結果我們發(fā)現，很多模塊的功能主要富集與細胞增殖及血管生成上，這些功能與腫瘤的發(fā)生發(fā)展相關，有趣的是我們也發(fā)現了吉西他濱模塊內的CRHBP和TP53均是細胞對藥物應答功能相關基因[34]，該藥物相關模塊內基因的富集結果(見圖3)。

圖3 吉西他濱模塊內基因富集分析結果Fig.3 Results of gene enrichment analysis in Gemcitabine module

2.3 關鍵調控基因對藥物耐藥性的影響

本研究對13種藥物進行了具體的分析，研究發(fā)現所有的耐藥相關基因出現在模塊的下游，針對相應的模塊進行了詳細的分析，以拓撲替康模塊為例(見圖4)。

圖4 拓撲替康藥物特異模塊Fig.4 Characteristic module of Topotecan

模塊內黃色的點為差異表達的耐藥相關基因，再對這些差異表達基因研究中我們發(fā)現，以CSF1基因為代表的模塊內葉子節(jié)點上的基因，在很多研究中已經發(fā)現它和藥物的耐藥性顯著相關[35]，而且其行使功能也是由于miR-130 miRNA家族基因調控的[36]，這與我們找到的模塊內的調控關系是一致的。然后我們進一步分析了網絡內最上游的轉錄因子及miRNA，這一模塊內最上游的基因為SP1，其可以通過表達的改變調控相應的通路使得腫瘤患者產生耐藥反應，這樣的通路包括MAPK信號通路[37]，MAPK通路內的RELA基因在模塊內也是由SP1調控在進一步調控其他耐藥相關基因。我們還發(fā)現網絡內有兩個miRNA位于葉子節(jié)點上，hsa-miR-9-5p和hsa-miR-125a-5p這樣的位于葉子節(jié)點上的基因更可能直接行使功能，使機體產生耐藥性。通過對最新的研究的檢索，發(fā)下hsa-miR-9-5p再轉錄后調控中扮演者重要的角色，通過改變DNA拓撲異構酶Ⅱα致使相應的細胞系產生耐藥性[38]。通過研究我們發(fā)現位于上游的基因可以通過一系列的調控關系，逐步的達到調控耐藥相關基因的目的，位于根節(jié)點的基因能直接的反應出機體是否會產生耐藥性，但往往單一的某個基因的異常不足以準確的確定機體是否產生耐藥。所以能通過節(jié)點調控一個或者多個耐藥相關的基因而且行使了重要功能的跟節(jié)點就顯得尤為關鍵，本研究認為這些上游的基因可能成未來耐藥研究的新的生物標記物。通過對模塊的分析，本研究挖掘了不同藥物模塊內的關鍵調控基因，這些基因可能成為潛在的生物標記物(見表2)。

表2 13種藥物模塊內潛在關鍵調控子Table 2 Potential key regulators in 13 drug modules

3 討 論

藥物耐藥性是癌癥治療失敗的一個重要原因，已經有研究發(fā)現了一些轉錄因子和miRNA可以影響藥物的耐藥性的產生，但這些研究大多是針對直接和耐藥相關的基因，對這些基因為什么表達異常研究還有所欠缺，本文通過構建了一個食管癌特異的轉錄調控網絡，基于這一網絡將可以調控耐藥基因的模塊都挖掘出來，得到13種作用于食管癌細胞系的耐藥模塊。而且隨著后續(xù)研究發(fā)現了包括SP1及hsa-miR-21-5p、hsa-miR-1-5p在內的關鍵的轉錄調控通路上游關鍵調控子，這些調控子可能經過一步或是幾步的調控最終改變耐藥基因的表達，從而導致患者耐藥性的產生。各模塊內的關鍵調控子可以作為潛在的生物標志物，用于檢測病人是否適用相應藥物。目前耐藥及敏感樣本的差異分析是根據不同的食管癌細胞系完成的，未來如果數據足夠充足，使用表達模式更復雜的不同癌癥亞型的樣本人群數據進行分析可以進一步提升分析結果的準確性。目前的研究分析的調控關系完全是來源于整理了實驗證實關系的數據庫，隨著研究的不斷進展，以及潛在的調控關系可能還未被研究人員發(fā)現，現在使用的這些調控關系可能不夠完整，會遺漏部分真實存在的調控關系對，導致網絡的連通程度降低，選取網絡最大連通組分的時候會有部分真實的調控關系的遺漏，如果能很好的解決假陽性率的問題，后續(xù)的研究中可以使用軟件預測的調控關系對來代替真實調控關系，將會大大增加轉錄調控網絡的大小以及最終挖掘的知識更加豐富。