共表達網絡分析識別心肌分化過程中的關鍵因子

文/張藝

心臟疾病伴隨著心肌細胞的壞死，已經壞死的心肌細胞會完全失去收縮功能，無法起死回生。心肌細胞再生能夠使心臟得到修復，使受損的心臟得到功能的完善,直接重編程能夠為這一目標的實現提供幫助。與干細胞技術相比，直接重編程把一類終末分化細胞直接轉化為其他類終末細胞，而不用經歷干細胞階段。在很短的時間內就可以獲得功能更佳的心肌細胞，誘導的時間也會很短。采用自體細胞不會受到倫理因素的限制，而且可能進行在體內直接轉分化，不必進行細胞移植，能夠避免形成腫瘤的風險、存活率低以及免疫排斥反應的種種問題。直接重編程誘導方向較為精確，也不會產生靶細胞之外的細胞類型。不過，心肌細胞直接重編程技術還不成熟，無法得到到成熟的心肌細胞，僅僅有一小部分細胞擁有自發電活動，自發收縮的功能；因為心肌細胞自身很難進行增殖，所以要獲取大量數目的細胞用于臨床是很困難的。

圖1：WGCNA算法基本步驟

圖2：共表達網絡

Davis等[5]在1987年的發現奠定了細胞重編程的基礎。2010年，Ieda等[6]受到多種轉錄因子組合誘導多能干細胞技術的啟發，首次報道了通過過表達三個心臟發育的核心轉錄因子Gata4、Mef2C、Tbx5(GMT)，可在體外成功誘導成纖維出現心肌樣細胞表型。之后，Song等[7]在經典三因子的基礎上加入Hand2（GHMT），可以使成年心臟成纖維細胞及鼠尾成纖維細胞，更有效地轉化為心肌樣細胞；Protze等[8]用胚胎成纖維細胞，篩選重編程因子，發現用Mef2C、Myocd、Tbx5(MMT)組合能夠上調心肌細胞的表達；Wada等[9]在GMT的基礎上，加入Mesp1和Myocd（GMTMM），能使人成纖維細胞表現出更多心肌細胞特點。整體來看，核心的重編程因子，能在一定程度上，把成纖維細胞轉為心肌細胞的狀態，但是效率并不高,如何高效地獲取更多功能完備的心肌細胞成為一個難題。因此提高成纖維細胞直接重編程為心肌細胞的效率，需要更深入研究轉錄機制，尋找更高效率的轉錄因子組合。

單細胞基因組測序技術在近些年來發展很快。傳統RNA測序方法，可單次處理大量的細胞，但得到的是平均表達水平。不存在兩個完全相同的細胞。單細胞RNA測序（scRNA-seq）可以挖掘出單個細胞獨具的基因表達變化。利用scRNA-seq技術，在單細胞的水平，分析心肌細胞發育的機制，會更有利于挖掘出重要的轉錄因子，重新對轉錄因子進行組合，從而提高重編程為心肌細胞的效率。

圖3：心肌相關聚類模塊

1 數據和方法

1.1 數 據

數據來自NCBI，數據庫編號為GSE76118，是2016年Sean等[11]人的工作，采集了胚胎發育第8.5天（e8.5）, e9.5和e10.5天的心臟單細胞轉錄組數據。

1.2 方法

1.2.1 數據處理

選取233個組織和細胞類型相關的候選核心轉錄因子，從下載好的RNA-Seq數據中篩選出這些候選轉錄因子的所有單細胞數據。過濾掉零表達和低表達的轉錄因子，最終剩余1258個轉錄因子作為研究對象。原始數據的基因表達值單位是CPM(count-per-million),在接下來的數據分析過程中需要進行對數化，取log2(CPM+1)。

1.2.2 共表達網絡構建方法

構建共表達網絡的步驟為，根據基因與基因之間的相關性關系判定兩基因之間的權重大小。這個關系描述了基因表達模式之間的相似性。相關性的計算方法有斯皮爾曼相關系數、皮爾遜相關系數等等。每個節點代表一個基因，每條邊代表共表達關系的強弱，利用相關性來構建出網絡。相關系數值分布于-1到1之間。在無向網絡中，取相關系數的絕對值。

大量的文獻選擇利用皮爾森相關系數來計算基因之間的權重，這樣的做法存在一個明顯的缺陷，閾值的取值通常是人為定義，取值過小容易造成網絡體積過大，冗余信息干擾到關鍵基因的發掘，取值過大則容易丟失有用信息。為避免這樣的處境，權重基因共表達網絡分 析（Weighted Gene Co-Expression Network Analysis，WGCNA)，利用軟閾值的判定方法，能夠對這個問題進行化解。

作為一種系統生物學分析方法，WGCNA構建了一個冪指數鄰接函數，對相關系數取n次冪，來進行加權。用鄰接系數αij=power(Sij,β)=|Sij|β，在加權系數β的作用下，相關系數的分布，趨向于無尺度分布。利用WGCNA分析方法構建加權共表達網絡的步驟如圖1所示。

1.2.3 從網絡中尋找功能聚類模塊

利用cytoScape中的插件MCODE（Molecular Complex Detection），在龐大的網絡中進行聚類，構建功能模塊。MCODE計算網絡中各個節點信息。以該點為種子節點，能擴展出的最大k值k-core，此K-core的密度以及該節點的分數值。節點的分數反映出，這個節點及其周邊的節點的密集度。之后，從分數值最大的節點開始，調用函數getClusterCore( )，利用該節點為種子節點，進行擴展，然后加入符合參數條件的鄰接節點。最后，據參數的要求進行后續的處理，得到功能模塊。

1.2.4 篩選網絡中樞紐基因

雖然多細胞生物中的所有細胞基本共享相同的DNA序列，但在每一種細胞類型中，僅有一組特定的基因被主動表達。不同的細胞類型被不同的活性基因組合和它們之間的相互作用所管控。在每個細胞內，一組靶基因和調節基因，即轉錄因子（TF），彼此相互作用并形成基因調控網絡。基因調控網絡在拓撲上包括高度連接的組件和幾個具有低連接性的節點。識別網絡中的關鍵基因，能夠更好的理解細胞命運的轉變挑戰在于如何找到掌控細胞命運轉變的基因。樞紐基因根據網絡拓撲性結合生物學意義來進行挑選。

對網絡中的節點使用cytoHubba進行網絡拓撲參數計算。cytoHubba包括節點連接度、邊緣滲出組件、最大鄰居組件、最大鄰居組件的密度、最大團中心性、瓶頸值、偏心度、緊密度、發散性、中介性、應力11個節點參數的計算。結合生物學意義，選用以下幾個參數來度量網絡中轉錄因子的重要程度。

假定G=(V,E)是個無向圖，V代表網絡中的點集，E是邊集。對于一個集合S，用|S|來表征它的基數（集合中的元素個數）。

a. 度（Degree，Deg）:

Deg(v)=|N(v)|，節點的度越高，該節點相連的節點越多，代表其重要性越高。

b. 緊密度（Closeness, Clo）：

c. 中介性（Betweenness，BC）

d. 瓶頸值（Bottleneck, BN）

Ts表示以s為根節點的最短路徑樹。BN(v)=∑s∈Vps(v)，當多于|V(Ts)|/4條從s出發的路徑在v點匯合，ps(v)=1；否則ps(v)=0。經過所有的瓶頸節點有一個共同特點，就是它們的鄰居節點可分為多個不相交的節點集合，這些節點集合需要互相通信交流的時候，只能經過該瓶頸節點。故瓶頸值越大，重要性也越高。

2 結果

表1：7個樞紐基因

2.1 共表達網絡構建

從預處理后的數據中篩選出1258個轉錄因子，并提取這些轉錄因子的表達矩陣，通過上述方法對網絡進行構建，如圖2。

2.2 鑒別共表達網絡中心肌相關模塊

通過cytoScape中MCODE工具的應用，在共表達網絡中識別到了8個共表達模塊，通過分析網絡模塊，圖3所示的模塊與心臟發育的過程緊密相關。

2.3 篩選樞紐基因

從圖3模塊中挑選樞紐基因，得到影響心肌發育的樞紐基因。

首先通過cytoHubba計算每個節點的Degree、Closeness、Betweenness、Bottleneck等參數值，評定出候選的18個樞紐基因。經過R語言的GOstats包對18個樞紐基因進行GO富集分析，最終綜合篩選出Ctnnb1、Nme2、Nf ia、Gata6、Hmgb1、Hmga1、Notch1（如表1）這7個樞紐因子在心肌細胞生長過程中發揮著極其重要的作用。

3 討論

通過對心肌發育的單細胞轉錄組數據進行預處理，利用加權基因共表達網絡分析算法，對1258個轉錄因子構建共表達網絡，通過MCODE聚類方法識別到心肌相關的模塊，通過網絡拓撲性的特征參數計算和GO富集分析最終篩選出重要的候選轉錄因子。其中，Notch1轉錄因子已被研究證明，對胚胎期心臟的發育有影響，并且和成年期動物心血管系統疾病的關系十分密切，Gata6轉錄因子對心肌能量代謝及細胞增殖具有重要作用。深入研究這些轉錄因子將是下一步工作的重要研究內容，應用生物學實驗方法確認這些因子在重編程過程中的作用也是今后工作的重點。

總體來說，直接重編程和干細胞技術相較，跳過了干細胞的階段，在短短的時間內可以獲得功能更佳的心肌細胞，誘導時間大大減少。不過心臟直接重編程技術還面臨著許多挑戰，對心臟直接重編程分子機制的了解還不夠深入。研究表明小鼠和人類的重編程過程也有所不同，如轉錄因子不同，需要的時間不同等。從目前的技術到轉化為臨床應用還是有很大距離，仍存在很多問題，需要利用更多新穎的計算機技術方法來進行持續深入的挖掘與探索。