KEGG通路是在功能分析中非常有效的富集方法，通過富集分析，我們可以有效識別出和目的基因（interest gene）顯著相關的生物學功能。常用的富集手段有超幾何分佈，fisher精確檢驗，二項分佈等等。例如DAVID軟體就是採用校正的Fisher exact test演算法。然而這種通路富集分析的方法存在一定的侷限性，那就是獲得的結果只能是一條完整的大通路，如代謝通路或其他調控通路，但是通過觀察，我們發現在每條通路中仍可進一步分為多個小的通路，並參與調控某些特定的功能，我們稱之為——子通路。
軟體工具#KEGG子通路拆分流程#
分析流程
1. 匯入資料
匯入資料為基因列表形式，即由使用者提供需要進行富集的基因列表，基因為gene symbol形式，格式如圖
所有的通路XML格式分析檔案已從KEGG通路資料庫下載，作為底層分析資料進行呼叫。命令列操作介面如下
使用者輸入包含基因列表的檔名，篩選顯著子通路的閾值，以及計算顯著性的演算法。
軟體工具#KEGG子通路拆分流程#
軟體工具#KEGG子通路拆分流程#
2. ID轉換
由於KEGG資料庫使用的是自己的ID，通路中包括多種元素，如‘hsa’開頭的代表蛋白質，‘cpd’開頭的代表化合物等等。在通路XML格式檔案中，所有的通路元素都被標記為KEGG的ID，如下所示
第一列數字代表元素ID，第二列為對應的基因或化合物資訊，有些ID可能對應多個基因。在富集過程中我們主要關注的是蛋白和基因，所以在拆分子通路的過程中，我們過濾掉化合物等非基因元素，最後將所有的元素ID轉換成gene symbol。所有的通路內基因編號，ID，symbol等資訊均來自KEGG資料庫主頁下載獲得。
軟體工具#KEGG子通路拆分流程#
3. 拆分子通路
拆分過程包括以下兩部分：
第一步，我們通過正則表示式，從XML通路檔案中提取所有的基因互作對，即對應著通路中通過邊連線的兩個基因，他們之間的邊可能代表啟用，抑制，磷酸化等，這些具體的互作型別不在我們考慮範圍內，我們只想知道兩個基因是否有互作關係。並且最重要的一點，從XML檔案獲取的互作對資訊是有向互作，即它們之間的邊代表了方向性，這保證了我們在後面從根節點向葉子節點grow的過程可以有向的進行。
第二步，收集所有存在互作的非孤立基因，獲得通路內的基因列表，對每一個基因進行迭代過程。在每一次迭代中，我們都假設這個基因為根節點，將所有與之互作的基因作為他的葉子節點，以此類推，直到獲得最末端的葉子，迭代終止。這樣我們就獲得了一條子通路。
4. 合併子通路
在拆分子通路過程中，由於每個基因節點都被作為根節點進行了一次迭代過程，所以會出現一種情況，那就是一條子通路的根節點很可能是另一條子通路的葉子節點，兩條子通路具有重複性。舉例來說就是一條子通路是從樹幹到葉子，而另一條子通路是這條樹幹上的一條樹枝。用圖來說明就是下面的情況
如圖，在由ABC構成的通路中，我們可以拆分得到下面三個子通路，分別是ABC, AC, BC。但是BC這條子通路其實是ABC的一部分，所以我們要將這樣的子通路進行合併，最後我們期望得到的結果是1. ABC 2. AC。
軟體工具#KEGG子通路拆分流程#
5. 顯著性計算
我們從一條KEGG通路中獲得了X條子通路，但是並不是每一條子通路都和我們的interest gene顯著相關，有些子通路中甚至不包含任何interest gene。為了識別出被interest gene顯著富集到的重要子通路，這裡我們採用超幾何分佈的演算法。對於一個KEGG通路，我們計算通路中完整基因作為背景資料，對於M個interest gene，我們計算在一條由N個基因構成的子通路中，至少有K個基因屬於interest gene的概率。輸出結果如下
使用者可根據右側顯著性P值篩選重要的子通路，若interest gene顯著富集到某一大通路，則這些顯著的子通路則提示了interest gene具體的生物活動和作用機制。
軟體工具#KEGG子通路拆分流程#
6. FDR校正
所有經過超幾何分佈或Fisher精確檢驗計算的子通路顯著性P值均經過BH多重檢驗法進行校正，校正前和校正後的P值資訊均儲存在結果檔案‘pathID_subpath_sig.txt’檔案。
由於對所有KEGG通路進行拆分和子通路顯著性富集後由於獲得資料和檔案較多，所以這裡增加了一個顯著性閾值的篩選，例如使用者在命令列設定閾值0.05，則只輸出顯著的主通路和子通路資訊。輸出檔案以每個主通路命名，包括每條子通路資訊，富集到子通路內的基因個數，基因列表，顯著性P值，FDR值，格式如下
完整程式碼下載地址
軟體工具#KEGG子通路拆分流程#