搜索
查看: 638|回复: 8

clusterprofiler>gseGO>genelist 怎么搞 多谢大家了【R语言】

[复制链接]

2

主题

13

帖子

139

积分

注册会员

Rank: 2

积分
139
发表于 2017-9-30 10:41:42 | 显示全部楼层 |阅读模式
目的:我现在手上有一堆geneID(symbol格式),想用clusterprofiler package做kegg 和 GO 的 GSEA enrichment
打算:用geo数据库得到表达量差异的fold change的排序的genelist,之后带入package跑。不知这种想法可以不

大家好,我最近在用Y叔的神包cluster profiler想做个gsea的enrichment,其他的都好说,但是,那个genelist我搞不清楚啊

在这个里面我用的gene列表是Y叔在clusterprofiler中id转换的示例的基因,用的代码大多数也是基于也是y叔讲的例子
[Plain Text] 纯文本查看 复制代码
library(clusterProfiler)
x <- c("GPX3",  "GLRX",   "LBP",   "CRYAB", "DEFB1", "HCLS1",   "SOD2",   "HSPA2",
       "ORM1",  "IGFBP1", "PTHLH", "GPC3",  "IGFBP3","TOB1",    "MITF",   "NDRG1",
       "NR1H4", "FGFR3",  "PVR",   "IL6",   "PTPRM", "ERBB2",   "NID2",   "LAMB1",
       "COMP",  "PLS3",   "MCAM",  "SPP1",  "LAMC1", "COL4A2",  "COL4A1", "MYOC",
       "ANXA4", "TFPI2",  "CST6",  "SLPI",  "TIMP2", "CPM",     "GGT1",   "NNMT",
       "MAL",   "EEF1A2", "HGD",   "TCN2",  "CDA",   "PCCA",    "CRYM",   "PDXK",
       "STC1",  "WARS",  "HMOX1", "FXYD2", "RBP4",   "SLC6A12", "KDELR3", "ITM2B")
eg = bitr(x, fromType="SYMBOL", toType=c("ENTREZID","GO"), OrgDb="org.Hs.eg.db")
head(eg)

eg_new = eg[c("ENTREZID","GO")]


结果> head(eg_new)
  ENTREZID         GO
1     2878 GO:0000302
2     2878 GO:0004602
3     2878 GO:0004602
4     2878 GO:0005515
5     2878 GO:0005576
6     2878 GO:0005615
>


If user have GO annotation data (in data.frame format with first column of gene ID and second column of GO ID), they can use enricher() and gseGO() functions to perform over-representation test and gene set enrichment analysis.

If genes are annotated by direction annotation, it should also annotated by its ancestor GO nodes (indirect annation). If user only has direct annotation, they can pass their annotation to buildGOmap function, which will infer indirection annotation and generate a data.frame that suitable for both enricher() and  gseGO().

这里面应该讲的是如果有GO注释的数据(数据框形式,第一列是gene ID,第二列是GO ID),因为我不清楚这个geneID到底是symbol还是entreID,于是两个都试了,都失败了。而且根据文档的上下文我猜这里的geneID应该是entreID。

之后跑gseGO

[Plain Text] 纯文本查看 复制代码
ego3 <- gseGO(geneList     = eg_new,
              OrgDb        = org.Hs.eg.db,
              ont          = "CC",
              nPerm        = 1000,
              minGSSize    = 100,
              maxGSSize    = 500,
              pvalueCutoff = 0.05,
              verbose      = TRUE)


得到报错
> ego3 <- gseGO(geneList     = eg_new,
+               OrgDb        = org.Hs.eg.db,
+               ont          = "CC",
+               nPerm        = 1000,
+               minGSSize    = 100,
+               maxGSSize    = 500,
+               pvalueCutoff = 0.05,
+               verbose      = TRUE)
Error in `[.data.frame`(x, order(x, na.last = na.last, decreasing = decreasing)) :
  undefined columns selected
>


所以  大神们啊  我该咋办啊 我这卡了一周多了   
谢谢大家



上一篇:16S扩增子数据拆分问题
下一篇:细胞死亡
回复

使用道具 举报

29

主题

120

帖子

914

积分

高级会员

Rank: 4

积分
914
发表于 2017-9-30 11:12:36 | 显示全部楼层
你物种是人,有org.Hs.eg.db包,就不要用自行新建geneid2goid的矩阵了(这个一般适用于非模式物种),直接用你最后一个代码就行了,修改下geneList,这个只是一个向量,包含你的gene symbol
回复 支持 反对

使用道具 举报

4

主题

22

帖子

491

积分

中级会员

Rank: 3Rank: 3

积分
491
发表于 2017-10-3 17:13:12 | 显示全部楼层
Y叔公众号里面这样写的:
GSEA输入的geneList,有三个特性:

数值型向量,可以是fold change或者其它的数值型变量都可以
命名,每一个数字都有一个对应的名字,就是相应的基因ID了
排序,数字是从高到低排序的
假设你的数据是一个csv文件,这个文件至少应该有两个columns,一个是基因ID(不允许有duplication),第二个是相应的表达量或fold change之类的数字型变量。

那么你的geneList可以这样搞出来:

d = read.csv(your_csv_file)
## assume 1st column is ID
## 2nd column is FC

## feature 1: numeric vector
geneList = d[,2]

## feature 2: named vector
names(geneList) = as.character(d[,1])

## feature 3: decreasing order
geneList = sort(geneList, decreasing = TRUE)
有了geneList,就可以愉快地用clusterProfiler进行GSEA分析了。
回复 支持 反对

使用道具 举报

2

主题

13

帖子

139

积分

注册会员

Rank: 2

积分
139
 楼主| 发表于 2017-10-10 11:04:07 | 显示全部楼层
xotong 发表于 2017-10-3 17:13
Y叔公众号里面这样写的:
GSEA输入的geneList,有三个特性:

不好意思回复完了,最近给自己挖的坑太多,填的累死了。

言归正传

谢谢你的回复,这个我也有看到,不过我这只有gneelist。想知道这个fold change怎么来。我设想的是fold change是否可以来自于GEO里面的数据,如果可行的话,具体应该怎样做。您要是比较忙,给我说几个key words 我自己先去查。

谢谢,麻烦了
回复 支持 反对

使用道具 举报

2

主题

13

帖子

139

积分

注册会员

Rank: 2

积分
139
 楼主| 发表于 2017-10-10 11:10:18 | 显示全部楼层
anlan 发表于 2017-9-30 11:12
你物种是人,有org.Hs.eg.db包,就不要用自行新建geneid2goid的矩阵了(这个一般适用于非模式物种),直接 ...

你好,谢谢你的回复
我现在的问题就是那个向量哪来,我现在只有一个genelist(里面只有symbol)。那个向量应该是表达量的降序吧,可是我没有啊。是否可以用GEO数据库里的。可是咋用,我不清楚,可以麻烦你给我说一下吗,如果你比较忙,说几个key words我自己去查也可以

谢谢
回复 支持 反对

使用道具 举报

29

主题

120

帖子

914

积分

高级会员

Rank: 4

积分
914
发表于 2017-10-10 19:40:04 | 显示全部楼层
蔡卷子 发表于 2017-10-10 11:10
你好,谢谢你的回复
我现在的问题就是那个向量哪来,我现在只有一个genelist(里面只有symbol)。那个向 ...

你物种是人吧?org.Hs.eg.db会把你的gene symbol转化为合适的gene id的,只要这个gene symbol在org.Hs.eg.db包中有,所以这个gene symbol可以作为你的genelist在gseGO函数中
回复 支持 反对

使用道具 举报

2

主题

13

帖子

139

积分

注册会员

Rank: 2

积分
139
 楼主| 发表于 2017-10-11 10:27:45 | 显示全部楼层
anlan 发表于 2017-10-10 19:40
你物种是人吧?org.Hs.eg.db会把你的gene symbol转化为合适的gene id的,只要这个gene symbol在org.Hs.eg ...

我的物种是人,您是说把这个genelist直接将格式为vetor的gene symbol输入进去吗。可是我这边还是报错。
这个是最后代码
[AppleScript] 纯文本查看 复制代码
kk <- gseGO(geneList = v.gene.list.t,
            OrgDb = org.Hs.eg.db,
            ont = "CC",
            nPerm = 1000,
            minGSSize = 100,
            maxGSSize = 500,
            pvalueCutoff = 0.05,
            verbose = TURE)


这个是genelist的格式和内容
[AppleScript] 纯文本查看 复制代码
> head(v.gene.list.t)
[1] "AARS2"  "AIRE"   "ATM"    "BLM"    "BMP15"  "BMPR1B"
> class(v.gene.list.t)
[1] "character"
回复 支持 反对

使用道具 举报

29

主题

120

帖子

914

积分

高级会员

Rank: 4

积分
914
发表于 2017-10-11 20:32:44 | 显示全部楼层
蔡卷子 发表于 2017-10-11 10:27
我的物种是人,您是说把这个genelist直接将格式为vetor的gene symbol输入进去吗。可是我这边还是报错。
...

不好意思哈,这个是R包的GSEA分析,我之前说错了,不仅需要gene symbol,还需有你在做差异分析的结果中的foldchange,然后将foldchange值赋予genelist向量,并将对应的gene symbol赋予这个genelist向量的名称。如果你没有表达量数据以及foldchange的话,就不要用clusterprofiler包的gesGO函数来做GO富集了,可以选择enrichGO函数,这个就只需要gene symbol就行了。。。之前说错了,对你造成误导了,不好意思。。。
回复 支持 反对

使用道具 举报

2

主题

13

帖子

139

积分

注册会员

Rank: 2

积分
139
 楼主| 发表于 2017-10-17 16:53:15 | 显示全部楼层
anlan 发表于 2017-10-11 20:32
不好意思哈,这个是R包的GSEA分析,我之前说错了,不仅需要gene symbol,还需有你在做差异分析的结果中的 ...

哦哦哦 gseGO的我做出来了  谢谢大神  大神能回答我  我已经很感谢了 谢谢大神
回复 支持 反对

使用道具 举报

您需要登录后才可以回帖 登录 | 立即注册

本版积分规则

QQ|手机版|小黑屋|生信技能树    

GMT+8, 2017-12-18 05:38 , Processed in 0.120890 second(s), 25 queries .

Powered by Discuz! X3.2

© 2001-2013 Comsenz Inc.