功能富集分析流程（一）

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

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

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

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

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

小編有話說(shuō)

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

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

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

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

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

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

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

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

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

通路列表如下，第一列是KEGG通路ID號(hào)，第二列是通路名稱。

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?獲取任意一個(gè)pathway下載網(wǎng)址

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

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

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

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

即替換這個(gè)網(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)?

#?生成某一通路的鏈接(手動(dòng)網(wǎng)頁(yè)復(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文件的特點(diǎn)和如何根據(jù)通路的xml文件提取每個(gè)通路中注釋的geneid，并獲得對(duì)應(yīng)關(guān)系。