GO/KEGG富集分析enrichment analysis 92

本教程使用GOstats来对给定的基因symbols做富集分析。

> library("org.Hs.eg.db")
> library("GSEABase")
> library("GOstats")
> genes <- c("AREG", "FKBP5", "CXCL13", "KLF9", "ZC3H12A", "P4HA1", "TLE1", "CREB3L2", "TXNIP", "PBX1", "GJA1", "ITGB8", "CCL3", "CCND2", "KCNJ15", "CFLAR", "CXCL10", "CYSLTR1", "IGFBP7", "RHOB", "MAP3K5", "CAV2", "CAPN2", "AKAP13", "RND3", "IL6ST", "RGS1", "IRF4", "G3BP1", "SEL1L", "VEGFA", "SMAD1", "CCND1", "CLEC3B", "NEB", "AMD1", "PDCD4", "SCD", "TM2D3", "BACH2", "LDLR", "BMPR1B", "RFXAP", "ASPH", "PTK2B", "SLC1A5", "ENO2", "TRPM8", "SATB1", "MIER1", "SRSF1", "ATF3", "CCL5", "MCM6", "GCH1", "CAV1", "SLC20A1")
> ##GO BP enrichment analysis
> goAnn <- get("org.Hs.egGO")
> universe <- Lkeys(goAnn)
> head(universe)
[1] "1"         "10"        "100"       "1000"      "10000"     "100008586"
> entrezIDs <- mget(genes, org.Hs.egSYMBOL2EG, ifnotfound=NA)
> head(entrezIDs)
$AREG
[1] "374"
 
$FKBP5
[1] "2289"
 
$CXCL13
[1] "10563"
 
$KLF9
[1] "687"
 
$ZC3H12A
[1] "80149"
 
$P4HA1
[1] "5033"
 
> entrezIDs <- as.character(entrezIDs)
> head(entrezIDs)
[1] "374"   "2289"  "10563" "687"   "80149" "5033" 
> params <- new("GOHyperGParams",
+                  geneIds=entrezIDs,
+                  universeGeneIds=universe,
+                  annotation="org.Hs.eg.db",
+                  ontology="BP",
+                  pvalueCutoff=0.01,
+                  conditional=FALSE,
+                  testDirection="over")
> over <- hyperGTest(params)
> library(Category)
> glist <- geneIdsByCategory(over)
> glist <- sapply(glist, function(.ids) {
+ 	.sym <- mget(.ids, envir=org.Hs.egSYMBOL, ifnotfound=NA)
+ 	.sym[is.na(.sym)] <- .ids[is.na(.sym)]
+ 	paste(.sym, collapse=";")
+ 	})
> head(glist)
                                                                                                                                                                                         GO:0002690 
                                                                                                                                                              "PTK2B;CXCL10;CCL3;CCL5;VEGFA;CXCL13" 
                                                                                                                                                                                         GO:0002688 
                                                                                                                                                              "PTK2B;CXCL10;CCL3;CCL5;VEGFA;CXCL13" 
                                                                                                                                                                                         GO:0048518 
"RHOB;ASPH;ATF3;CCND1;BMPR1B;CAV1;CAV2;CCND2;PTK2B;GJA1;IL6ST;CXCL10;IRF4;LDLR;SMAD1;MAP3K5;PBX1;RFXAP;CCL3;CCL5;SLC20A1;TLE1;CLEC3B;VEGFA;CFLAR;CXCL13;TXNIP;CYSLTR1;AKAP13;MIER1;CREB3L2;ZC3H12A" 
                                                                                                                                                                                         GO:1901623 
                                                                                                                                                                           "PTK2B;CCL3;CCL5;CXCL13" 
                                                                                                                                                                                         GO:0002687 
                                                                                                                                                              "PTK2B;CXCL10;CCL3;CCL5;VEGFA;CXCL13" 
                                                                                                                                                                                         GO:2000403 
                                                                                                                                                                           "PTK2B;CCL3;CCL5;CXCL13" 
> bp <- summary(over)
> bp$Symbols <- glist[as.character(bp$GOBPID)]
> head(bp)
      GOBPID       Pvalue OddsRatio    ExpCount Count Size                                        Term
1 GO:0002690 3.490023e-08  39.24490  0.19267044     6   52 positive regulation of leukocyte chemotaxis
2 GO:0002688 9.278865e-08  32.80297  0.22601725     6   61          regulation of leukocyte chemotaxis
3 GO:0048518 1.170957e-07   4.28014 13.56103476    32 3660   positive regulation of biological process
4 GO:1901623 1.176400e-07 128.80174  0.04816761     4   13         regulation of lymphocyte chemotaxis
5 GO:0002687 1.496978e-07  30.05918  0.24454325     6   66  positive regulation of leukocyte migration
6 GO:2000403 2.969857e-07  96.58170  0.05928321     4   16 positive regulation of lymphocyte migration
                                                                                                                                                                                            Symbols
1                                                                                                                                                               PTK2B;CXCL10;CCL3;CCL5;VEGFA;CXCL13
2                                                                                                                                                               PTK2B;CXCL10;CCL3;CCL5;VEGFA;CXCL13
3 RHOB;ASPH;ATF3;CCND1;BMPR1B;CAV1;CAV2;CCND2;PTK2B;GJA1;IL6ST;CXCL10;IRF4;LDLR;SMAD1;MAP3K5;PBX1;RFXAP;CCL3;CCL5;SLC20A1;TLE1;CLEC3B;VEGFA;CFLAR;CXCL13;TXNIP;CYSLTR1;AKAP13;MIER1;CREB3L2;ZC3H12A
4                                                                                                                                                                            PTK2B;CCL3;CCL5;CXCL13
5                                                                                                                                                               PTK2B;CXCL10;CCL3;CCL5;VEGFA;CXCL13
6                                                                                                                                                                            PTK2B;CCL3;CCL5;CXCL13
> dim(bp)
[1] 525   8
> ##KEGG enrichment analysis
> keggAnn <- get("org.Hs.egPATH")
> universe <- Lkeys(keggAnn)
> params <- new("KEGGHyperGParams", 
+                     geneIds=entrezIDs, 
+                     universeGeneIds=universe, 
+                     annotation="org.Hs.eg.db", 
+                     categoryName="KEGG", 
+                     pvalueCutoff=0.01,
+                     testDirection="over")
> over <- hyperGTest(params)
> kegg <- summary(over)
> library(Category)
> glist <- geneIdsByCategory(over)
> glist <- sapply(glist, function(.ids) {
+ 	.sym <- mget(.ids, envir=org.Hs.egSYMBOL, ifnotfound=NA)
+ 	.sym[is.na(.sym)] <- .ids[is.na(.sym)]
+ 	paste(.sym, collapse=";")
+ 	})
> kegg$Symbols <- glist[as.character(kegg$KEGGID)]
> kegg
  KEGGID       Pvalue OddsRatio ExpCount Count Size                                   Term                                    Symbols
1  04510 9.643801e-05  7.873256 1.124361     7  200                         Focal adhesion    CCND1;CAPN2;CAV1;CAV2;CCND2;ITGB8;VEGFA
2  04060 5.487123e-04  5.821855 1.489779     7  265 Cytokine-cytokine receptor interaction BMPR1B;IL6ST;CXCL10;CCL3;CCL5;VEGFA;CXCL13
3  04062 3.739699e-03  5.486219 1.062521     5  189            Chemokine signaling pathway              PTK2B;CXCL10;CCL3;CCL5;CXCL13
> library("pathview")
> gIds <- mget(genes, org.Hs.egSYMBOL2EG, ifnotfound=NA)
> gEns <- unlist(gIds)
> gene.data <- rep(1, length(gEns))
> names(gene.data) <- gEns
> for(i in 1:3){pv.out <- pathview(gene.data, pathway.id=as.character(kegg$KEGGID)[i], species="hsa", out.suffix="pathview", kegg.native=T)}

写成函数

enrich <- function(entrezIDs, orgDbName="org.Hs.eg.db", pvalueCutoff=.01){
	require(orgDbName, character.only=TRUE)
	require("GSEABase")
	require("GOstats")
	require("Category")
	goAnn <- get(gsub(".db", "GO", orgDbName))
	universe <- Lkeys(goAnn)
	onto <- c("BP", "MF", "CC")
	res <- lapply(onto, function(.onto){
		param <- new('GOHyperGParams',
					 geneIds= entrezIDs,
					 universeGeneIds=universe,
					 annotation=orgDbName,
					 ontology=.onto,
					 pvalueCutoff=pvalueCutoff,
					 conditional=FALSE,
					 testDirection="over")
		over <- hyperGTest(param)
		glist <- geneIdsByCategory(over)
		glist <- sapply(glist, function(.ids) {
			.sym <- mget(.ids, envir=get(gsub(".db", "SYMBOL", orgDbName)), ifnotfound=NA)
			.sym[is.na(.sym)] <- .ids[is.na(.sym)]
			paste(.sym, collapse=";")
		})
		summary <- summary(over)
		if(nrow(summary)>1) summary$Symbols <- glist[as.character(summary[, 1])]
		summary
	})
	names(res) <- onto
	keggAnn <- get(gsub(".db", "PATH", orgDbName))
	universe <- Lkeys(keggAnn)
	param <- new("KEGGHyperGParams",
				 geneIds=entrezIDs,
				 universeGeneIds=universe,
				 annotation=orgDbName,
				 categoryName="KEGG",
				 pvalueCutoff=pvalueCutoff,
				 testDirection="over")
	over <- hyperGTest(param)
	kegg <- summary(over)
	glist <- geneIdsByCategory(over)
	glist <- sapply(glist, function(.ids) {
		.sym <- mget(.ids, envir=get(gsub(".db", "SYMBOL", orgDbName)), ifnotfound=NA)
		.sym[is.na(.sym)] <- .ids[is.na(.sym)]
		paste(.sym, collapse=";")
	})
	kegg$Symbols <- glist[as.character(kegg$KEGGID)]
	res[["kegg"]] <- kegg
	res
}
 
go <- enrich(entrezIDs, "org.Hs.eg.db", .05)

92 thoughts on “GO/KEGG富集分析enrichment analysis

  1. Reply 逝水 3月 26,2014 8:20 上午

    glist <- sapply(glist, function(.ids)
    这一步会报错,没有glist对象,是不是这一步有问题?

  2. Reply super 4月 6,2014 2:11 下午

    你好,请问这个GO/KEGG分析的包和DAVID哪个更powerful?我只用DAVID做过,从没用过这个~~

  3. Reply super 4月 6,2014 2:14 下午

    补上条,我还要请教另外还有很多收费软件,如metacore,ipa等 它们比起免费软件的优势在哪里?劣势就不用说了……要花钱

    • Reply admin 4月 7,2014 7:50 上午

      我觉得本质上没有区别,它们都用的是相同的算法。我几乎没有用过收费软件来进行enrichment analysis。我想主要的区别在于人性化方面,比如UI设计,比如图型化输出等等。还有就是,如果你用于商业用途的话,开源软件也很有可能是收费的。

  4. Reply super 4月 8,2014 5:58 上午

    谢谢。我的数据是一组ipi ID的基因,必须先用DAVID转成gene symbol才可以用GOstats,所以我就直接用DAVID了。但我试验了您的数据觉得也很赞毕竟全在R环境下操作。但好像必须限制要gene symbol。
    此外如果我的基因不是人的 比如是老鼠的牛的 找哪个数据库?

    • Reply admin 4月 8,2014 3:52 下午

      其实重要的是得到entrez_id, 示例中使用gene symbol主要是因为这个最常见。如果基因不是人的,而是其它物种的,基本大同小异。把org.Hs.eg.db换成相应的数据库即可,如果是小鼠的就是org.Mm.eg.db,牛的是org.Bt.eg.db.相应的,org.Hs.egGO, org.Hs.egPATH等都换成对应的物种的就可以了,比如org.Mm.egGO…

  5. Reply super 4月 8,2014 7:25 下午

    你好 我问的就是你说的这些我怎么知道 什么牛是Bt 小鼠是Mm这些 哪里有索引?
    此外不用gene_symbol如果用其它表示形式怎么弄?直接写ensembl ID 或ipi 的id就行了么?谢谢

  6. Reply Colin 5月 13,2014 3:12 上午

    在您的博客中受益匪浅,我是一名在读的硕士研究生,想请教您怎么总结这么多R方面的知识,平时的工作中都是完全要用到才去学习吗?而且给我的感觉要熟悉这些不同的包就要看很多文档,做很多调试啊,真实超人,能不能分享些经验!谢谢您

    • Reply admin 5月 13,2014 10:30 上午

      我学习使用R已经3年了。没啥经验,就是每个工作日都用。用得多了,就一通百通了。最重要的是,不断向有他人学习。很多人都有自己的经验,一但他共享给我,我就用心体会就好了。很多知识都是从初学者那里得到的呢。

  7. Reply ? 5月 16,2014 1:38 上午

    GO term enrichment之后,外文文献中一般要做张图(条形图是横着的,对应的纵坐标的GO术语也是横着的),请问下,那张图有什么简单的方法可以做?

  8. Reply ygc 7月 10,2014 1:32 上午

    clusterProfiler包 ygc.name/2012/05/07/clusterprofiler-package-comparing-biological-themes-gene-clusters 会直接返回富集的基因,如果函数调用时参数readable设为TRUE,返回来的就是gene symbol.

    直接,简单,方便。

  9. Reply Mabing 7月 11,2014 3:40 上午

    博主您好,两个问题,请教一下。第一个问题,爬你这段富集分析代码时,这个错误
    > gIds <- mget(genes, org.Hs.egSYMBOL2EG, ifnotfound=NA)
    Error in sqliteExecStatement(con, statement, bind.data) :
    RS-DBI driver: (connection with pending rows, close resultSet before continuing)
    是怎么回事?
    第二问题:博主R生信分析知识非常广泛和扎实,对于像我这种初学者,能给点学习建议么?爬去您的程序运行,排错,换数据,继续跑数据,这样有意义么?谢谢您~

    • Reply admin 7月 15,2014 10:16 下午

      不好意思,没有看到您的留言,所以这么晚才回复。
      1. 数据联接已经打开,关闭再联即可。
      2. 有意义。任何事情都有个过程。一开始不明所以,但是可能增加自己运行的信心与经验。慢慢地,知识丰富起来,就好了。

  10. Reply hiker 7月 15,2014 9:03 下午

    博主好,贴您的代码之后,运行出现这个错误
    > over kegg <- summary(over)
    Error: getAnnMap: package KEGG not available
    博主的代码高深难懂,可否注释再详细些,让俺等小白可以看懂。

  11. Reply boya888 7月 26,2014 8:07 下午

    博主老师好!目前已经得到所有的差异表达的基因,请问GOstats可以作水稻Osativa、玉米maize等作物的转录组(RNA-seq)的富集分析吗?http://bioconductor.org/packages/release/BiocViews.html#___Organism页面只有2个玉米maizecdf、maizeprobe和1个BSgenome.Osativa.MSU.MSU7文件!
    请问怎么办呢?谢谢

    • Reply admin 7月 27,2014 9:14 上午

      是肯定可以的。但是因为这两者都没有org的annotation包,甚至没有db0包,所以需要自己准备一些文件,具体参见:http://www.bioconductor.org/packages/release/bioc/vignettes/GOstats/inst/doc/GOstatsForUnsupportedOrganisms.pdf

  12. Reply boya888 7月 27,2014 8:12 下午

    非常感谢指导!data.frame object格式为
    go_id Evidence gene_id
    1 GO:0008150 ND 1
    2 GO:0002576 TAS 2
    可是植物基因组注释文件GenomeAnnotation.table没有evidence codes (GO的?)这一项
    5: PFAM 6: Panther 7: KOG 8: KEGG ec 9: KEGG Orthology 10: Gene Ontology terms
    GRMZM2G000231_P01 PF00326 PTHR11731,PTHR11731:SF54 GO:0008236,GO:0006508
    GRMZM2G000236_P01 PF00724 PTHR22893,PTHR22893:SF13 KOG0134 GO:0010181,GO:0016491,GO:0055114
    好在keggframeData 只需要KEGG Ids 和 gene IDs两列!
    请问我是不是可以依据GenomeAnnotation.table来分拣出GO 和 KEGG的 GO.db和 KEGG.db呢?
    再次感谢

    • Reply admin 7月 28,2014 2:01 下午

      为什么不去http://www.geneontology.org/page/download-annotations下载Oryza sativa的?这里是http://www.geneontology.org/page/evidence-code-decision-tree,可能有一些代码不在这个tree里面。对于maize的,你可以去http://www.maizegdb.org/documentation/maizemap/ontology.php下载。

  13. Reply boya888 7月 27,2014 8:59 下午

    补充:Gene Ontology terms (NOTE: these are automated results from interpro2go in most genomes, *not* empirically derived)
    难道GO Evidence全是ND ?

  14. Reply Zander 7月 30,2014 10:47 下午

    博主你好,有两个问题请教一下。
    (1)GO enrichment analysis的时候,博主采用org.Hs.egGO,不知道有没有尝试过使用biomaRt包提取与GO-BP有关的genes?
    (2)KEGG enrichment analysis的时候,博主采用org.Hs.egPATH,系统提示KEGG.db的数据是older version,已经过时了,系统建议用KEGGREST或者RECTOME包,好像GOstats包目前还不支持这两个包了,不知道博主有木有解决的方法?

    • Reply admin 7月 31,2014 6:08 上午

      1,可以使用biomaRt提取。完全没有问题。其区别在于使用org.Hs.egGO使用的是本地库,速度快。而使用biomaRt是在线数据库,资料新。
      2,对于KEGG,的确没有什么太好的办法。可以使用KEGGREST来桥接。我本人不推荐使用RECTOME,结果不是很理想。但是不妨试试。但其实,如果是科研用的话,KEGG.db的库也是够的,并不是很旧。

  15. Reply boya888 8月 2,2014 8:30 上午

    请问博主KEGG Orthology(KO)和kegg号的对应关系的问题!http://www.kegg.jp/kegg-bin/show_pathway?@ko00904/reference%3dwhite/default%3d%23bfffbf/K04120/K04125/K14038
    网页上 ko00904对应三个kegg号(可是在线分析下载result.ko文件却没有这些对应信息)。google搜 “kegg id convert”也没有明确的结果!请问(KO)和kegg号的对应关系在哪儿下呢?

    • Reply admin 8月 5,2014 8:22 上午

      这个我没有做过。而且我对KEGG Orthology还真不熟。但是我想,使用http://www.kegg.jp/kegg-bin/view_ortholog_table?orthology=K14038这样的网址,可以得到你想要的任何转换。具体它的对应关系从哪里下载,需要花时间去阅读它的说明文档。

  16. Reply monkeyjun 8月 3,2014 7:00 上午

    谢谢博主老师上次关于mirdeep2的提示, 现在要进一步对miRNA-seq的结果进行分析
    ,就是想问下您您用的R 的版本是多少,因为在查找KEGG包的时候出现“package ‘KEGG’ is not available (for R version 3.0.2) ” 的情况,希望能得到您的指导。
    谢谢。
    一名生物信息本科生

  17. Reply boya888 8月 5,2014 7:47 下午

    博主老师好!我在用gostats作玉米maize的KEGG富集分析时遇到问题!toTable出了问题, frame = toTable(‘/export/home/GOstats/zmRkeggFrame1.txt’)#参照M. Carlson, 2014写的操作手册#
    Error in (function (classes, fdef, mtable) : unable to find an inherited method for function ‘toTable’ for signature ‘”character”’可是能正常读
    > frame head(frame)
    koID GeneID
    1 K00002 GRMZM2G479423_T05
    2 K00004 GRMZM2G154007_T02
    但是frame = toTable就出问题!

    • Reply admin 8月 6,2014 12:38 下午

      你的toTable是从哪里来的?在AnnotationDbi中的toTable只接受Bimap,而你这输入的是个文件名。我不知道你所谓的M. Carlson, 2014写的操作手册是什么,没有链接,我无从考证。

  18. Reply boya888 8月 6,2014 5:56 上午

    补充:我是在服务器管理帐户的R上载入本用户安装的GOstats所需要的包,如category等!
    GOstats能否在作出的pathway中标示上下调的基因(红,上调;绿,下调)呢?谢谢

  19. Reply 天天小白 8月 12,2014 1:18 上午

    请问有没有包 包括KEGG的全部关系 比如以2列的形式给出每个通路的全部基因之间的联系和方向??

  20. Reply 天天小白 8月 15,2014 5:10 上午

    非常感谢~ 已经找到通路和基因的对应关系了~ 请问哪个函数能调出全部相连基因对及其方向?

  21. Reply ff2345 8月 25,2014 2:03 上午

    老师,您好!我做的物种还没有全基因组发表。我用KAAS得到了转录组序列号与KO号的一个列表。请问如何将我的基因序列号与pathway联系起来再做富集分析?

    • Reply admin 8月 25,2014 8:03 上午

      这样的情况我没有遇见过。推荐你两步走:1,阅读文献,看是否有其它人做过类似的事情,如果有,参照前人做法。2,如果没有,将从KAAS得到KO号列表当做全库,将你的基因列表转换成KO号当成抽选库,然后使用http://pgfe.umassmed.edu/ou/archives/3372中的办法计算P值。

  22. Reply ff2345 8月 26,2014 3:05 上午

    KAAS指的是KEGG网站上的自动注释服务器,KO 指的是KEGG pahway上的每个同源蛋白(KEGG Orthology),我的目的是想找到KO和pathway得对应关系。看到您的博客上对KEGG API也有所研究,想问下有没有比较好的方法

    • Reply admin 8月 27,2014 9:51 上午

      这个我没有做过。无能为力了。

    • Reply liyi 5月 25,2017 4:46 上午

      前辈,你好,我关注的物种全基因组去年发表的,但kegg还没收录。我用KAAS得到了整个基因id与KO号的一个列表。然后RNA-seq后得到了一些差异基因,请问接下来具体该如何将我的这些基因序列做kegg富集分析,还有遇到的问题就是kaas注释的KO号对应pathway说明,少于整个基因组id对应的KO号。不知道说清楚没?本人也还没弄大懂,望赐教

      • Reply admin 5月 28,2017 10:01 上午

        这个时候还是使用reactome吧。对于KEGG,你可以利用同源性做一些分析,但是也许不是最佳的方案。关于你说的KO号的问题,我没有能够理解清楚。

  23. Reply ff2345 9月 3,2014 10:00 上午

    上面那个问题我自己解决了,用gage包的kegg.gsets(“ko”)函数,再通过ID的转换就能够计算富集。@boya888
    还有个问题,大部分文章中看得都是FDR,不只是p-value,一直不明白如何从p-value经过Bofforin校正算出FDR值。

  24. Reply boya888 9月 11,2014 7:33 上午

    谢谢博主和ff2345的指点!我们做得是植物的GO和KEGG富集分析!操作手册是How To Use GOstats and Category to do Hypergeometric testing with unsupported model organisms (Carlson, 2014)。
    作植物得自己用AnnotationDbi建annotation 包吧!参照http://grokbase.com/t/r/bioconductor/0943r8h6yq/bioc-hypergtest-different-results-using-different-annotation-packages构建GO.db,可是这个帖子是芯片数据,而我们做得是RNA-seq数据。还有,KEGG.db该怎么构建呢?有没有unsupported model organisms的GO和KEGG富集分析的protocol(类似于博主的帖子或文档)呢?谢谢

  25. Reply Robin Li 9月 18,2014 3:04 上午

    博主你好:
    在您的code中我发现您在前面使用的是envir=org.Hs.egSYMBOL2EG,在后面使用的是envir=org.Hs.egSYMBOL,我help了一下,发现两个envir的help文档是一样的,请问两者之间有区别么?两者之间可以互换?

  26. Reply Robin Li 9月 19,2014 10:01 下午

    这样啊!我 ls(“package:org.Hs.eg.db”)看来一下,是不是所有的后缀加了EG的数据都是和相应的相反,从symbol到id,没有EG的是从id到symble?

  27. Reply huyudong 9月 24,2014 8:37 上午

    老师,你好!请问(1)本例中,为什么需要 gEns <- unlist(gIds), gene.data <- rep(1, length(gEns)), names(gene.data) <- gEns难道不能使用gene.data <-as.character(gIds)?(2)我目前在做酵母的KEGG富集分析时,按照本教程,前面所有的分析都没有问题,但在用pathview做图时,出现问题,提示No of the genes or compounds mapped to the pathway!Argument gene.idtype or cpd.idtype may be wrong. 请问这个可能的原因是什么,关于酵母的gene.idtype我设定为“KEGG”,因为酵母在map时,用的是systematic names(类似:YPR163C),非Entrez ID,请问这个设置正确吗?我使用的注释为org.Sc.sgd.db.

    • Reply admin 9月 24,2014 9:31 下午

      关于问题2,你可以仔细阅读pathview中关于gene.idtype的说明。我不能保证使用KEGG就是对的,但是似乎应该是这样。
      关于问题1,这里的gene.data是需要一系统的gene的表达值,它应该是numeric,它的names是相应的基因名。

  28. Reply huyudong 9月 25,2014 11:16 上午

    嗯,好的,谢谢老师,我再试试!

  29. Reply biogg 11月 26,2014 2:04 上午

    博主,您好!
    在做GO/kegg分析时,Bioconductor包中的函数一般都是调用已有的数据库,对与人或者小鼠常见的物种,做这些分析容易些。但对于一些不常见的物种该怎么做呢?多谢指教!

    • Reply admin 11月 30,2014 3:55 下午

      对于bioconductor没有提供的物种,是有一些困难,但是我的建议是先理解标准的物种怎么做,然后知道为什么了之后,你自然就会做非标准的物种了。如果你已经完全了解了标准物种的做法之后仍然有困难,请再次留言讨论。

      • Reply biogg 12月 4,2014 4:13 上午

        多谢回复,如老师所说了解了标准物种的做法,自然会做非标准物种。我在实际应用中遇到的问题没有解决,应该是我没有真正掌握分析方法,买Bioconductor红皮书学习了一段时间,主要是想做GO和pathway富集分析以及可视化展示,使用GeneAnswers包来显示显著性富集的pathway网络关系。书中敲代码上去,每一步的输入文件详细的打出了解是什么格式,但遇到了一个至今未解决的问题。
        关键的一句代码:y=geneAnswersBuilder(humanGeneInput,”org.Hs.eg.db”,categoryType=”KEGG”,testType=”hyperG”,pvalueT=0.1,geneExpressionProfile=humanExpr,verbose=FALSE)
        getEnrichmentInfo(y)[1:6,]

        输入或者参考的三个文件humanGeneInput和humanExpr,”org.Hs.eg.db”
        > head (humanGeneInput,2)
        EntrezID logFC P.Value
        39400_at 23102 -0.999785 1.482860e-05
        1303_at 6452 -1.343031 1.974284e-05

        humanExpr文件,行名基因名称,列是基因在各个样品中的表达量。

        我使用的是二代测序的数据,物种较偏,根据基因比对nr库的结果提取了每个基因对应的gi号(Entrez ID,humanGeneInput和humanExpr这两个文件顺利的得到了。

        接下来要构造的是非标准物种geneAnswersBuilder需要的annotationLib,说来惭愧,R基础较差,不清楚标准物种如“org.Hs.eg.db”中的文件格式该怎么看,参考过geneAnswersBuilder的一些说明,尝试构造我个人需要的annotationLib,但没有成功。

        我想我遇到的问题应该是比较简单的,但一直没有解决掉,还请指点下,非常感谢。

        • Reply admin 12月 4,2014 9:39 上午

          我建议你先放弃参考红皮书的想法,因为那是为写书而准备的例子,只考虑到了好运行,易上手。你的这个问题比较复杂,你还是自主学习一下GOStats比较容易一些。或者你可以参考ChIPpeakAnno包中相关的部分。你需要了解的关键是超几何分布,理解了它,我想困难就会少很多。

  30. Reply biogg 11月 27,2014 12:57 上午

    看这页的评论学习了!还得再试试!痛苦无助的过程,希望能得到满意的结果啊

  31. Reply yjw 11月 27,2014 3:00 下午

    老师您好,我是您书的读者,按照书上内容做KEGG分析时出现
    yy <- geneAnswersReadable (y,verbose = FALSE)
    Loading required package: KEGG.db

    KEGG.db contains mappings based on older data because the original resource was removed from the the public domain before the most recent update was produced. This package should now be considered deprecated and future versions of Bioconductor may not have it available. Users who want more current data are encouraged to look at the KEGGREST or
    reactome.db packages
    是不是因为KEGG数据库的缘故,现在用geneanswers包做的都不可靠了。
    还有用显著富集的KEGG结果做热图时
    yyy <- geneAnswersSort (yy, sortBy="pvalue")
    看到帮助文档sortby有很多种复习方法,包括genenum, pvalue, foldchange,oddsratio等,可能问题比较傻,不太理解这些方法富集到疾病的原理,为什么选不同方法结果差异很大,那选哪个最好呢?谢谢老师

    • Reply admin 11月 30,2014 4:18 下午

      首先,感谢你认真阅读了书中的内容。
      关于KEGG的问题,可以这么说,你不用太担心富集结果不可靠的问题,最多是这一年内的注释不能及时更新了。现在的确是更鼓励使用KEGGREST这类的包直接从KEGG网站调数据了。但是我也没有摸到最后的解决方案呢。
      geneAnswersSort应该不会影响差异结果,它只是对结果进行排序。但是当你使用热图来显示结果的时候,因为排序不同,你会看到不同的结果,我建议使用pvalue过滤后再使用foldchange排序。

  32. Reply yjw 12月 11,2014 12:24 上午

    老师您好,请教一下您,GOterm中列表只有计数,那在go term上富集的差异表达基因如何查看,或者导出呢?另外是不是count和size差距越大越好呢,还有,ExpCount的大小说明什么呢,我做出来的结果ExpCount很小,接近于0.

    • Reply admin 12月 11,2014 9:04 下午

      你的问题好像不少,主要还是没有理解enrichment 的一些基本概念。最好的办法是阅读一些原始文献。我的博文中有关于如何从go term返回gene的,你搜一下。当然,还有很多办法,回头我写上来。关于结果的解读,请先查阅vignettes

  33. Reply yjw 12月 11,2014 8:40 上午

    另外,老师,请教一下,差异基因的pathway分析能用R做吗,有没有什么好方法,后面的pathview就是实现这一功能的吗

  34. Reply Aure 1月 6,2015 11:13 下午

    博主,发现不同的R版本安装的org.Hs.eg.db的版本就不同,universe的id不同,最终的富集结果也不同,R3.0.2的版本head(universe)的结果和你的相似,R3.1的版本head(universe)的结果是[1] “1” “2” “3” “9” “10” “11”。这个怎么解决呢?另外universe <- Lkeys(org.Hs.egGO)和universe <- Lkeys(org.Hs.egPATH)的结果应该是一致的吗?谢谢

    • Reply admin 1月 7,2015 1:01 下午

      head(universe)的不同很可能只是排序不同。所以不用太介意。universe <- Lkeys(org.Hs.egGO)和universe <- Lkeys(org.Hs.egPATH)的结果应该是不一致的,因为PATH和GO的记录条数可能不一样。

  35. Reply Aure 1月 7,2015 10:47 下午

    谢谢。但是博主,不同的版本的org.Hs.eg.db内包含的基因数目确实存在差异,看了下2.5比2.10少了1000个基因左右,最后的富集结果也有差异,2.10和DAVID的结果更接近些

    • Reply admin 1月 8,2015 9:10 上午

      谢谢分享你的结论。这个具体我没有比较过。这也是生物信息学现在面临的一个问题,不同软件版本之间操作的差异性。所有高质量的文章,一定需要实验结果的验证。对于结果的可重复性问题,只要说明完备,其它的学者可以重复就好。

  36. Reply maria 1月 8,2015 7:45 上午

    楼主,我在运行 genes <- c()时出现错误,
    Error: unexpected symbol in:
    " "TRAPPC2P1", "EDNRA", "EYA4", "CYP3A4", "CCNT1", "KHSRP", "PAX7", "CDK4", "ZAP70", "NOTCH4", "NANOG", "HNF4A", "ULBP2", "MAPT", "MTHFD1L", "KLK10", "KLK6", "MPL", "CYP19A1", "NKX2-1", "SLC5A5
    我输入的基因总共有六百多个,是因为容纳不了这么多基因,还是说它列出来的这些symbol是数据库无法识别的。
    急用啊

  37. Reply maria 1月 8,2015 9:27 上午

    后来我试了,发现
    genes <- c( "CYP1B1", "TRAPPC2P1", "EDNRA", "EYA4", "CYP3A4", "CCNT1", "KHSRP", "PAX7", "CDK4", "ZAP70", "NOTCH4", "NANOG", "HNF4A", "ULBP2", "MAPT", "MTHFD1L", "KLK10", "KLK6", "MPL", "CYP19A1", "NKX2-1", "SLC5A5", "TG", "RAP1A", "STAT3", "BDNF", "BECN1", "TNRC6A", "DTL", "PIK3CD", "SEPT7", "AVEN")
    这样是可以的,也就是说跟字符格式无关,我将六百多个基因分写到两个数据集:genes1,genes2。这样也是行的通的,但是这样两个数据集无法进行下面的步骤。
    我想请问老师您有什么方法可以再让他们合在一起么,我用了merge语句,但是仍旧行不通。或者有没有什么方法可以不用手动输入这么多基因,直接用文件形式导入。
    期待您的回复。

  38. Reply maria 1月 8,2015 9:29 上午

    或者如果您有时间的话,我可以把这些基因发您的邮箱,您帮我看看?

  39. Reply aure 1月 9,2015 7:14 上午

    那么多基因,建议使用scan命令导入文本文件进行识别

  40. Reply Pygathrix 4月 21,2015 5:34 上午

    你好,请问目前的代码是补过上面的bug的了吗?

  41. Reply Pygathrix 4月 21,2015 9:52 上午

    博主好,请指教!

    我按步骤模拟这个脚本的时候碰到一点问题。

    > params over head(bp)
    GOBPID Pvalue OddsRatio ExpCount Count Size
    1 GO:1901623 4.385268e-09 114.089636 0.06502873 5 19
    2 GO:2000403 1.255872e-08 88.714597 0.07871898 5 23
    3 GO:0009893 1.286923e-08 4.817152 10.88717761 30 3181
    4 GO:0042127 1.790927e-08 5.942930 4.83950617 20 1414
    5 GO:1901700 2.842488e-08 6.007788 4.45960152 19 1303
    6 GO:0048247 5.222769e-08 63.847059 0.10267693 5 30
    Term
    1 regulation of lymphocyte chemotaxis
    2 positive regulation of lymphocyte migration
    3 positive regulation of metabolic process
    4 regulation of cell proliferation
    5 response to oxygen-containing compound
    6 lymphocyte chemotaxis
    Symbols
    1 CXCL13;PTK2B;CXCL10;CCL3;CCL5
    2 CXCL13;PTK2B;CXCL10;CCL3;CCL5
    3 CXCL13;AKAP13;PTK2B;GCH1;GJA1;IL6ST;CXCL10;IRF4;ITGB8;AREG;LDLR;SMAD1;MAP3K5;ASPH;ATF3;PBX1;CCND1;RFXAP;RGS1;CCL3;CCL5;CREB3L2;TLE1;CLEC3B;VEGFA;ZC3H12A;CAV1;CAV2;CFLAR;CCND2
    4 CXCL13;TXNIP;PTK2B;GJA1;IGFBP7;IL6ST;CXCL10;AREG;SMAD1;ASPH;ATF3;PBX1;CCND1;CCL5;BMPR1B;KLF9;VEGFA;CAV1;CAV2;CFLAR
    5 CXCL13;TXNIP;PTK2B;GCH1;GJA1;IGFBP7;CXCL10;AREG;RHOB;MAP3K5;ASPH;CCND1;CCL3;CCL5;KLF9;ZC3H12A;CAPN2;CAV1;CAV2
    6 CXCL13;PTK2B;CXCL10;CCL3;CCL5

    这些结果跟帖子里的不一样。

    • Reply admin 4月 21,2015 10:46 上午

      之前讨论过,使用不同版本的org注释,会导致不同的结果,这是可以意料的。但是不知道你的这个结果是不是也出自这个原因。
      我的教程,只是教给你一个思路。明白理解了每一步是什么意思是最重要的。至于结果的可重复性,我想我会在以后的博文中多注意这个问题。

  42. Reply Pygathrix 4月 21,2015 9:54 上午

    博主好,刚才的帖子没贴全,这一步出现了这个问题:

    > params over <- hyperGTest(params)
    Warning messages:
    1: In .local(name, pos, envir, all.names, pattern) :
    ignoring 'pos' argument
    2: In .local(name, pos, envir, all.names, pattern) :
    ignoring 'envir' argument
    3: In .local(name, pos, envir, all.names, pattern) :
    ignoring 'all.names' argument

    但是下面可以继续走下去。不知是缺什么信息?会不会影响结果?

  43. Reply paulchen 11月 7,2016 1:41 上午

    glist <- sapply(glist, function(.ids) {.sym <- mget(.ids, envir=org.Hs.egSYMBOL, ifnotfound=NA).sym[is.na(.sym)] <- .ids[is.na(.sym)]paste(.sym, collapse=";")})这段正确应该怎么写??谢谢老师,上面弄错了

  44. Reply mervin 5月 26,2017 10:56 下午

    Error in getGoToEntrezMap_db(p) :
    The genes you are testing do not have any corresponding GO terms for the ontology you are searching.

    老师您好。我遇到了上面这个错误,遇到没有富集到任何GO的gene list程序会终止,如何才能跳过这种gene list继续运行呢?

    • Reply admin 5月 28,2017 10:06 上午

      你可以试试维护者给出的下面的代码(我抄过来的):
      ontology=”BP”
      Gos <- mget(eset, org.Hs.egGO, ifnotfound=NA) haveGo <- sapply(Gos, function(x) { if (length(x) == 1 && is.na(x)) FALSE else { onts <- Biobase::subListExtract(x, "Ontology", simplify=TRUE) ontology %in% onts }}) eset <- eset[haveGo]

Leave a Reply

  

  

  

%d 博主赞过: