GOseq的介紹
GOseq是一個R包，用於尋找GO terms，即基因富集分析。此方法基於 Wallenius non-central hyper-geometric distribution。相對於普通的超幾何分布(Hyper-geometric distribution)，此分布的特點是從某個類別中抽取個體的概率與從某個類別之外抽取一個個體的概率是不同的，這種概率的不同是通過對基因長度的偏好性進行估計得到的，從而能更為准確地計算出 GO term 被差異基因富集的概率。
1.GOseq的安裝
>BiocManager::install("goseq")
2.參考數據集
這里我們采用GOseq包里的內置數據集genes來做GO分析
1 library(goseq) 2 data(genes) 3 head(genes) 4 str(genes)
這里genes數據集是EMSEMBL gene的向量集合，其中1代表差異表達
3.通過getgo函數獲得GO terms
getgo的用法：
1 getgo(genes, genome, id,fetch.cats=c("GO:CC","GO:BP","GO:MF"))
genes:genes是輸入的gene向量或列表
genome:參考基因組，比如hg38,hg19
id:輸入基因的類型,比如ensGene
fetch,cats:fetch.cats是"GO:CC", "GO:BP", "GO:MF" & "KEGG"的一系列組合
這里用supportedOrganisms()函數來查看支持的genome和id
結果：結果是一個列表，包含每一個gene對應的所有的GO ID,這個值是goseq函數中gene2cat參數的輸入值
舉例：
1 genes <- c("ENSG00000124208", "ENSG00000182463", "ENSG00000124201", "ENSG00000124205", "ENSG00000124207") 2 getgo(genes,'hg19','ensGene')
這里顯示了每一個gene參與到的所有GO ID
3.通過getlength函數檢索gene的長度
getlength用法：
1 getlength(genes, genome, id)
結果：結果是一個向量，包含所有基因的長度，如果某個基因的長度無法檢索到，用NA代替。這個向量是nullp函數中bias.data的輸入值。
舉例：
1 genes <- c("ENSG00000124208", "ENSG00000182463", "ENSG00000124201", "ENSG00000124205", "ENSG00000124207") 2 getlength(genes,'hg19','ensGene')
這里基因長度出現了3036.5，是因為這里基因長度取得是轉錄本長度的中位數。
4.使用nullp函數(Probability Weighting Function)計算概率加權函數
nullp函數介紹：
Calculates a Probability Weighting Function for a set of genes based on a given set of biased data (usually gene length) and each genes status as differentially expressed or not.
nullp函數用法：
1 nullp(DEgenes, genome, id, bias.data=NULL,plot.fit=TRUE)
DEgenes：DEgenes的格式是一個二元向量，其中1代表差異表達，0代表非差異表達，還有包括gene id，格式與內置數據集genes一樣
bias.data:bias.data是一個數值向量，通常是基因轉錄本長度的中位數，單位是bp.如果設置bias.data=NULL，nullp函數將通過getlength函數來獲取gene的長度。所以這里默認設置為bias.data=NULL
plot.fit:plot.fit這里將pwf作圖，默認設置為plot,fit=TRUE
一般nullp函數后面的參數都選擇默認，只用選擇設置DEgenes, genome和 id即可
結果：結果是一個數據框，行名為gene id,列名為"DEgenes", "bias.data" 和 "pwf"，這個數據框對象是goseq函數的輸入，用來計算富集的GO terms，也可以作為plotPWF的輸入，用來進一步作圖。
舉例：
1 data(genes) 2 pwf <- nullp(genes, 'hg19', 'ensGene')
5.使用goseq函數進行GO富集分析
goseq函數介紹：
Does selection-unbiased testing for category enrichment amongst differentially expressed (DE) genes for RNA-seq data. By default, tests gene ontology (GO) categories, but any categories may be tested.
goseq函數用法：
1 goseq(pwf, genome, id, gene2cat = NULL,test.cats=c("GO:CC", "GO:BP", "GO:MF"),method = "Wallenius", repcnt = 2000, use_genes_without_cat=FALSE)
pwf：這里的pwf是由nullp函數得到的結果，為一個數據框
gene2cat:這里的gene2cat是由getgo函數得到的結果，如果設置gene2cat=NULL，goseq函數將會自動地用getgo函數來獲得GO ID，默認設置是gene2cat=NULL
method:這里method有三種選擇，"Wallenius", "Sampling" 和 "Hypergeometric".這里"Sampling" 和 "Hypergeometric"方法幾乎從沒被使用過
一般goseq函數后面的參數都可以選擇默認，只用選擇pwf,genome和id這三個參數就可以
舉例：
1 data(genes) 2 pwf <- nullp(genes,'hg19','ensGene') 3 pvals <- goseq(pwf,'hg19','ensGene') 4 head(pvals)
這里的選擇over_represented_pvalues<0.05就是具有統計學意義的GO ID了
1 enriched.GO<-pvals[pvals$over_represented_pvalue<0.05,]
