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 分析
- library(goseq)
- data(genes)
- head(genes)
- 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 参数的输入值
举例:
- genes <- c("ENSG00000124208", "ENSG00000182463", "ENSG00000124201", "ENSG00000124205", "ENSG00000124207")
- getgo(genes,'hg19','ensGene')
这里显示了每一个 gene 参与到的所有 GO ID
3. 通过 getlength 函数检索 gene 的长度
getlength 用法:
1 getlength(genes, genome, id)
结果: 结果是一个向量, 包含所有基因的长度, 如果某个基因的长度无法检索到, 用 NA 代替. 这个向量是 nullp 函数中 bias.data 的输入值.
举例:
- genes <- c("ENSG00000124208", "ENSG00000182463", "ENSG00000124201", "ENSG00000124205", "ENSG00000124207")
- 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 的输入, 用来进一步作图.
举例:
- data(genes)
- 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 这三个参数就可以
举例:
- data(genes)
- pwf <- nullp(genes,'hg19','ensGene')
- pvals <- goseq(pwf,'hg19','ensGene')
- head(pvals)
这里的选择 over_represented_pvalues<0.05 就是具有统计学意义的 GO ID 了
1 enriched.GO<-pvals[pvals$over_represented_pvalue<0.05,]
来源: http://www.bubuko.com/infodetail-3358319.html