前面已經(jīng)講述了R包clusterProfiler的GO富集分析方法动雹，本篇繼續(xù)演示R包goseq的GO富集分析款筑。
相比clusterProfiler中的GO富集分析方法疏橄，goseq的特別之處在于，不再使用超幾何分布（Hyper-geometric distribution）檢驗，而是使用了Wallenius non-central hyper-geometric分布。除了背景基因的數(shù)量外鱼填，它還同時考慮了基因長度信息，認為從某個類別中抽取個體的概率與從某個類別之外抽取一個個體的概率是不同的毅戈，這種概率的不同是通過對基因長度的偏好性進行估計得到的苹丸，從而其認為能更為準確地計算出 GO term被差異基因富集的概率。
goseq包的安裝
對于goseq的安裝也很簡單苇经，一般情況下谈跛，直接通過Bioconductor安裝goseq就可以了。
#bioconductor 安裝
#install.packages('BiocManager') #需要首先安裝 BiocManager塑陵，如果尚未安裝請先執(zhí)行該步
BiocManager::install('goseq')
goseq的GO富集分析（有參向）
這里均對于有參考基因組的情況而言的，無參分析暫不涉及蜡励。
就以人類轉(zhuǎn)錄組數(shù)據(jù)為例展示GO分析的過程吧令花。
1、準備輸入數(shù)據(jù)
需要準備兩類數(shù)據(jù)凉倚。
（1）待分析的基因名稱兼都，例如這里以人類參考基因組hg38版本的ensembl id為例。把基因名稱以一列的形式排開稽寒，放在一個文本文件中（例如命名“gene_select.txt”）扮碧。Excel中查看，就是如下示例這種樣式。
（2）所有基因長度信息蚓土，既包含待分析的基因，也要包含背景基因赖淤，根據(jù)所選擇的參考基因組定義蜀漆。例如這里統(tǒng)計了人類參考基因組hg38版本的所有基因長度，第一列是基因名稱咱旱，第二列是基因長度确丢，放在一個文本文件中（例如命名“gene_length.all.txt”）。Excel中查看吐限，就是如下示例這種樣式鲜侥。
2诸典、goseq的GO富集分析
首先讀取基因列表和長度文件描函，預處理數(shù)據(jù)。
#讀取待GO分析的基因名稱搂赋，該示例來自人參考基因組 hg38 的基因名稱
de_gene <- as.vector(read.delim('gene_select.txt', sep = '\t')[[1]])
#讀取人類 hg38 的所有基因長度信息
gene_list <- read.delim('gene_length.all.txt', sep = '\t', row.names = 1)
#設置待富集的基因為 1赘阀，背景基因為 0
genes <- rep(0, nrow(gene_list))
names(genes) <- rownames(gene_list)
genes[de_gene] <- 1
head(genes) #處理后的待分析數(shù)據(jù)
結(jié)果獲得一組向量結(jié)構(gòu)在R中存儲。向量中記錄了所有基因的名稱脑奠，以及其是作為背景基因存在（0）還是作為富集基因存在（1）基公。
隨后加載包goseq，并使用goseq的內(nèi)部函數(shù)enrichGO()即可完成GO富集分析宋欺。
library(goseq)
#根據(jù)基因長度加權(quán)
pwf <- nullp(DEgenes = genes, bias.data = gene_list$gene_length, plot.fit = FALSE)
head(pwf)
#GO 富集分析
#hg38 是內(nèi)置的人類 hg38 基因組參考庫
#ensGene 意為基因名稱使用 ensembl id 作匹配進行富集
pvals <- goseq(pwf, 'hg38', 'ensGene')
head(pvals)
#手動對 P 值作個校正轰豆，例如 FDR 校正
pvals$FDR <- p.adjust(pvals$over_represented_pvalue, method = 'fdr')
pvals <- pvals[c(1,2,8,4,5,6,7)] #第三列的p值就不要了
head(pvals)
對于各列內(nèi)容：
category，富集到的GO id齿诞；
over_represented_pvalue酸休，富集的p值；
FDR祷杈，p調(diào)整值斑司；
numDEInCat和numInCat，分別為富集到該GO條目中的基因數(shù)目但汞，以及該條目中背景基因總數(shù)目宿刮；
term，富集到的GO的描述信息私蕾；
ontology僵缺，GO分類BP（生物學過程）、CC（細胞組分）或MF（分子功能）踩叭。
3磕潮、手動添加基因名稱
我們看到翠胰，上述goseq的默認輸出中，只統(tǒng)計了富集到該GO條目中的基因數(shù)量自脯，并未把具體的基因名稱展示出來之景。如果需要，基因名稱我們來手動匹配下冤今。
#添加基因闺兢，參考方案
#https://www.biostars.org/p/355247/
getGeneLists <- function(pwf, goterms, genome, ids){
gene2cat <- getgo(rownames(pwf), genome, ids)
cat2gene <- split(rep(names(gene2cat), sapply(gene2cat, length)),
unlist(gene2cat, use.names = FALSE))
out <- list()
for(term in goterms){
tmp <- pwf[cat2gene[[term]],]
tmp <- rownames(tmp[tmp$DEgenes > 0, ])
out[[term]] <- tmp
}
out
}
goterms <- pvals$category
goList <- getGeneLists(pwf, goterms, 'hg38', 'ensGene')
#默認將富集到GO的基因名稱添加在最后一列，然后輸出到本地
pvals$Ensembl_ID <- sapply(pvals$category, function(x) paste(goList[[x]], collapse = '; '))
write.table(pvals, 'goseq.txt', sep = '\t', row.names = FALSE, quote = FALSE)
內(nèi)容格式見上文戏罢，最后一列添加了富集到該條目的基因名稱信息屋谭。
由于輸入文件使用的ensembl id，因此最后展示的也是ensembl id龟糕。如果期望使用其它的基因名稱桐磁，如通俗的symbol id等類型，除了更改輸入文件為使用symbol id的基因名稱做分析外讲岁，還可以通過基因名稱轉(zhuǎn)換的方式對ensembl id和symbol id作個匹配轉(zhuǎn)換我擂。
