标签归档:Biomart

topGO包进行-GO富集分析-做有向无环图

GOdata = new(“topGOdata”, ontology = “MF”, allGenes = geneList,annot = annFUN.gene2GO, gene2GO = geneID2GO)

利用topGO进行分析,最重要的是构建topGO对象,构建topGO需要两个参数:

1,topGO需要基因和GO号的对应关系

2,基因列表,用来标记背景基因(所有基因)及差异基因

一,获取ensembl和GO号的对应关系: geneID2GO

如果你有现成的gene id和go id的对应关系,文件格式为

gene_ID制表符GO_ID1, GO_ID2, GO_ID3, ...

每行,则可以利用readMappings读取该文件(topGO包里面的函数)。

而我有一个cuffdiff的结果文件gene_exp.diff,用的是ensembl的注释GFF文件。首先需要获得ensembl和GO的对应关系,这里利用biomaRt包(前面的文章中有讲,这里不详细介绍)。

library(biomaRt)
genes = useEnsembl(biomart="ensembl",dataset="hsapiens_gene_ensembl")
# 得到go信息和gene
gene2goInfo <- getBM(attributes=c('ensembl_gene_id','go_id','entrezgene','name_1006','go_linkage_type','namespace_1003'), mart = genes)
# 过滤
gene2goInfo=gene2goInfo[gene2goInfo$go_id != "", ]
# 上述是一个gene对应一个go id,需要合并为一个gene对应多个go id,利用by函数(神器)
geneID2GO = by(gene2goInfo$go_id, gene2goInfo$ensembl_gene_id, function(x) as.character(x))

geneID2GO是符合topGO要求的,gene2goInfo和geneID2GO的格式如下

> head(gene2goInfo)
  ensembl_gene_id      go_id entrezgene                         name_1006
3 ENSG00000281614 GO:0005886       3635                   plasma membrane
4 ENSG00000281614 GO:0005829       3635                           cytosol
5 ENSG00000281614 GO:0005515       3635                   protein binding
6 ENSG00000281614 GO:0007165       3635               signal transduction
7 ENSG00000281614 GO:0005856       3635                      cytoskeleton
8 ENSG00000281614 GO:0050852       3635 T cell receptor signaling pathway
  go_linkage_type     namespace_1003
3             IEA cellular_component
4             TAS cellular_component
5             IPI molecular_function
6             TAS biological_process
7             IEA cellular_component
8             TAS biological_process

> head(geneID2GO)
$ENSG00000000003
[1] "GO:0004871" "GO:0005515" "GO:0005887" "GO:0070062" "GO:0007166"
[6] "GO:0043123" "GO:0039532" "GO:1901223"

$ENSG00000000005
[1] "GO:0016021" "GO:0005737" "GO:0005515" "GO:0005635" "GO:0016525"
[6] "GO:0001886" "GO:0001937" "GO:0071773" "GO:0035990"

二,构建基因列表:geneList

diff=read.table("./gene_exp.diff",sep="t",header=TRUE)
# 差异表达基因
interesting_genes=factor(diff[diff$significant=="yes",]$gene_id)
# 所有基因
all_genes <- diff$gene_id
# 构建基因列表
geneList <- factor(as.integer (all_genes %in% interesting_genes))
names(geneList)=all_genes

三,构建topGO对象

source("https://bioconductor.org/biocLite.R")
chooseBioCmirror()
biocLite("topGO")
biocLite("org.Hs.eg.db")
biocLite("GO.db")
biocLite("Rgraphviz")
library(topGO)

GOdata <- new("topGOdata", ontology = "MF", allGenes = geneList,annot = annFUN.gene2GO, gene2GO = geneID2GO)

Building most specific GOs …..
( 3893 GO terms found. )

Build GO DAG topology ……….
( 4356 GO terms and 5490 relations. )

Annotating nodes ……………
( 16722 genes annotated to the GO terms. )

四,基因富集分析

Fisher检验,通过2*2的列联表的进行计算。当然还有其他检验,比如KS,KS.elim。

—————-非差异表达基因—差异表达基因
注释到A通路——– 20 ———– 50
没有注释到A通路—- 1870 ——— 80

resultFisher <- runTest(GOdata, algorithm = "classic", statistic = "fisher")
allRes <- GenTable(GOdata,classicFisher = resultFisher, orderBy = "classicFisher", topNodes = 10)

根据orderBy的参数对象进行排序,得到结果如下

GO.ID Term Annotated Significant
1 GO:0032794 GTPase activating protein binding 14 2
2 GO:0015278 calcium-release channel activity 16 2
3 GO:0099604 ligand-gated calcium channel activity 16 2
4 GO:0005527 macrolide binding 19 2
5 GO:0005528 FK506 binding 19 2
6 GO:0001225 RNA polymerase II transcription coactiva… 1 1

Expected classicFisher
1 0.05 0.0012
2 0.06 0.0016
3 0.06 0.0016
4 0.07 0.0023
5 0.07 0.0023
6 0.00 0.0038

五,生成有向无环图(directed acycline praph,DAG)

svg("go.dag.svg")
showSigOfNodes(GOdata, score(resultFisher), firstSigNodes = 5, useInfo = 'all')
dev.off()

继续阅读

利用biomaRt包下载HGMD公开版的突变位点

前两面文章介绍了Ensembl的biomart,相信你对biomart应该有所了解了,在此再介绍一种方法,即通过R语言包biomaRt下载HGMD的数据。

HGMD的最新数据是需要购买授权才行,公开版信息不仅滞后,而且不能下载,不能得到基因组位置,在biostart上看到有人说Ensembl整合了HGMD的公开版,心想能获得公开版的数据也不错,于是采用biomaRt包下载。

各位不要高兴,最终的结果是,只得到了所有突变的基因组位置,未能下到具体的突变类型,以及与表型的关系。不过能下载基因组的位置,也算不错,结合对这些位置的注释,能获取不少信息。如果您对这些位置的利用有更多或者更好的想法,欢迎与我讨论。

1,安装biomaRt包

source("http://bioconductor.org/biocLite.R")
chooseCRANmirror()
chooseBioCmirror()
biocLite("biomaRt")

2,显示ensembl的biomart

library(biomaRt)
listEnsembl()
#如果要显示特定版本的,添加version参数
#listEnsembl(version=78)


继续阅读

安装BioMart Perl及利用BioMart Perl API下载数据

上篇介绍了如何利用ensembl的biomart服务下载ensembl gene id与NCBI entrez gene id的对应关系时,最后一步是保存result。biomart也提供通过biomart-perl,在本地通过perl脚本下载,并通过标准输出到终端上。biomart提供生成好的perl脚本,只需在选择好相关attribute和filter之后,点击中间上方的perl即可。
biomart-perl-api
继续阅读

如何获取Ensembl gene id和NCBI的gene id及与HGNC的对应关系

Ensembl和NCBI都是盛名的基因组研究机构,提供相关的基因组结构注释文件,比如gtf或者gff,但注释的id却不是统一的。比如基因ID,Ensembl有Ensembl gene id,NCBI有entrez gene id。不同的人用的基因注释文件来源不同,就需要进行转换。本文主要讲如何利用Ensembl的Biomart,下载对应关系。

Biomart整合了各种生物学注释数据,提供了易于操作的界面,在线提供批量下载,以加速科学研究。Ensembl已应用biomart提供相关服务。

The BioMart project provides free software and data services to the international scientific community in order to foster scientific collaboration and facilitate the scientific discovery process. The project adheres to the open source philosophy that promotes collaboration and code reuse.

Ensembl的biomart网址为http://asia.ensembl.org/index.html

第一步,选择相应的数据库

选择ensemble gene 83

step1-select-database

选择homo sapiens gene

step1-select-database

继续阅读