功能富集分析流程（一）

功能富集分析流程是怎樣的？

接下來，我們以KEGG數(shù)據(jù)庫為例，帶大家梳理下功能富集分析流程：

1、教你怎么將KEGG數(shù)據(jù)庫本地化。

2、如何自己定義功能富集的背景基因集（以物種人類為例，選取一個平臺檢測的所有基因和全部的人類已知基因做背景基因有什么區(qū)別么？差異在哪里？）

3、如何自己寫代碼實現(xiàn)功能富集核心原理。

小編有話說

首先做一個基因集合的功能富集，必須搞清楚一件事情，基因功能富集和基因功能注釋的差別：前者是一個 gene set 是否共同參與完成了一個生物學(xué)過程，而功能注釋只是針對single gene來說的，這個基因可能參與的生物學(xué)過程，有一些什么功能，只需要看這個基因是否在某個通路中出現(xiàn)，一個簡單的map過程而已，沒有用到統(tǒng)計學(xué)模型，不需要顯著性計算。

而評價一個gene set是否共同參與完成了某個生物學(xué)過程，這有很多統(tǒng)計學(xué)模型和算法來評估這個事情，比如最常見的超幾何模型和GSEA功能富集的排秩原理。這里小編會以超幾何原理為例，告訴你怎么自己寫代碼實現(xiàn)這個功能富集過程，加深對功能富集的理解。

第一步 KEGG數(shù)據(jù)庫本地化

》獲取KEGG通路中所有通路列表

1.1?KEGG官網(wǎng)https://www.genome.jp/kegg/，首先我們看一下KEGG數(shù)據(jù)庫中的Pathway。

KEGG將pathway分為了七大類，每一大類下面有小類，大家可以自己去詳細(xì)看看是什么樣的。

1.?Metabolism?代謝	代謝
2.?Genetic?Information?Processing	遺傳信息處理
3.?Environmental?Information?Processing	環(huán)境信息處理
4.?Cellular?Processes	細(xì)胞過程
5.?Organismal?Systems	組織系統(tǒng)
6.?Human?Diseases	疾病
7.?Drug?Development	藥物

1.2?進(jìn)入如下這個網(wǎng)址：

https://www.kegg.jp/kegg/rest/keggapi.html#list，點擊如圖所示鏈接得到KEGG所有人類通路的列表，截至20190418，收錄了330條通路。

通路列表如下，第一列是KEGG通路ID號，第二列是通路名稱。

path:hsa00010	Glycolysis / Gluconeogenesis – Homo sapiens (human)
path:hsa00020	Citrate cycle (TCA cycle) – Homo sapiens (human)
path:hsa00030	Pentose phosphate pathway – Homo sapiens (human)
path:hsa00040	Pentose and glucuronate interconversions – Homo sapiens (human)
path:hsa00051	Fructose and mannose metabolism – Homo sapiens (human)
path:hsa00052	Galactose metabolism – Homo sapiens (human)
path:hsa00053	Ascorbate and aldarate metabolism – Homo sapiens (human)
path:hsa00061	Fatty acid biosynthesis – Homo sapiens (human)
path:hsa00062	Fatty acid elongation – Homo sapiens (human)

1.3?獲取任意一個pathway下載網(wǎng)址

進(jìn)入網(wǎng)址https://www.kegg.jp/kegg/kegg3a.html，點擊圖片中框框位置，得到通路下載鏈接：

https://www.kegg.jp/kegg-bin/download?entry=hsa00010&format=kgml

1.4? 根據(jù)上面下載連接的特點和1.2通路列表id生成所有的KEGG通路下載地址

https://www.kegg.jp/kegg-bin/download?entry=hsa00010&format=kgml

即替換這個網(wǎng)址中entry=hsa00010即可，R代碼如下：

#?生成某一通路的鏈接,最后用迅雷批量下載xml文件
rm(list=ls())
setwd(“/path”)

#?導(dǎo)入kegg文件（通路ID,通路名字）
kegg?<-?read.table(“KEGG_pathwayID.txt”,header=F,sep=“\t”,quote=“\””,stringsAsFactors=F)
dim(kegg)?

#?生成某一通路的鏈接(手動網(wǎng)頁復(fù)制)
original_URL?<-?“http://www.kegg.jp/kegg-bin/download?entry=hsa00010&format=kgml”??
replaced_pathwayId?<-?gsub(“path:”,“”,kegg[,1])

replaced_URL?<-?sapply(replaced_pathwayId,function(x){gsub(“hsa00010”,x,original_URL)})

hsa_keggURLs?<-?cbind(pathwayId=replaced_pathwayId,pathwayName=kegg[,2],pathwayURL=replaced_URL)

write.table(hsa_keggURLs,file=“hsa_keggURLs.txt”,sep=“\t”,row.names=F,quote=F)

然后使用wget 命令下載hsa_keggURLs.txt文件中所有通路的xml文件，或者使用迅雷下載也是可以噠！

本期功能富集分析流程先講到這里，下期我們講解xml文件的特點和如何根據(jù)通路的xml文件提取每個通路中注釋的geneid，并獲得對應(yīng)關(guān)系。