搜索
查看: 2439|回复: 3

[basic] 关于clusterprofiler的一些疑问(结果和教程不一致)

[复制链接]

1

主题

3

帖子

32

积分

新手上路

Rank: 1

积分
32
发表于 2018-5-12 07:06:16 | 显示全部楼层 |阅读模式
本帖最后由 tuanzi 于 2018-5-12 07:09 编辑

首先我读取了基因symbol
[AppleScript] 纯文本查看 复制代码
library('clusterProfiler')
library('org.Hs.eg.db')
setwd("C:/Users/lenovo/Desktop/Diabetic children/DEG(logFC=0.5,P=0.01)/DEG")
rawgene=read.csv("Geneup.csv",header=TRUE,row.names=1,check.names = FALSE)

之后用Jimmy大神写的转换代码转化了ID
[AppleScript] 纯文本查看 复制代码
geneIDannotation <- function(geneLists,name=T,map=T,ensemble=F,accnum=F){
  ## input ID type : So far I just accept entrezID or symbol
  ## default, we will annotate the entrezID and symbol, chromosone location and gene name 

  suppressMessages(library("org.Hs.eg.db"))
  all_EG=mappedkeys(org.Hs.egSYMBOL) 
  EG2Symbol=toTable(org.Hs.egSYMBOL)
  if( all(! geneLists %in% all_EG) ){
    inputType='symbol'
    geneLists=data.frame(symbol=geneLists)
    results=merge(geneLists,EG2Symbol,by='symbol',all.x=T)
  }else{
    inputType='entrezID'
    geneLists=data.frame(gene_id=geneLists)
    results=merge(geneLists,EG2Symbol,by='gene_id',all.x=T)
  }

  if ( name ){
    EG2name=toTable(org.Hs.egGENENAME)
    results=merge(results,EG2name,by='gene_id',all.x=T)
  }
  if(map){
    EG2MAP=toTable(org.Hs.egMAP)
    results=merge(results,EG2MAP,by='gene_id',all.x=T)
  }
  if(ensemble){
    EG2ENSEMBL=toTable(org.Hs.egENSEMBL)
    results=merge(results,EG2ENSEMBL,by='gene_id',all.x=T)
  }
  if(accnum){
    EG2accnum=toTable(org.Hs.egREFSEQ) 
    results=merge(results,EG2MAP,by='gene_id',all.x=T)
  }
  return(results)
}
gene=na.omit(geneIDannotation(as.matrix(rawgene)[,1]))

结果如下
[AppleScript] 纯文本查看 复制代码
> gene
      gene_id         symbol                                                                   gene_name cytogenetic_location
1       10004       NAALADL1                         N-acetylated alpha-linked acidic dipeptidase like 1              11q13.1
2   100129195    ZSCAN16-AS1                                                     ZSCAN16 antisense RNA 1               6p22.1
3   100129387     GABPB1-AS1                                                      GABPB1 antisense RNA 1              15q21.2
4   100130452   LOC100130452                                                uncharacterized LOC100130452               2q33.1
5   100131211          NEMP2                                nuclear envelope integral membrane protein 2               2q32.2
6   100287569      LINC00173                                  long intergenic non-protein coding RNA 173             12q24.22
7   100499489   LOC100499489                                                uncharacterized LOC100499489              10p12.2
8   100506639   LOC100506639                                                uncharacterized LOC100506639                 5p12
9   100506755       MIR497HG                                               mir-497-195 cluster host gene              17p13.1
10  100507217      LINC01578                                 long intergenic non-protein coding RNA 1578              15q26.1
11  100507334   LOC100507334                                               two pore channel 3 pseudogene                 2q13
12  100507362      LINC01015                                 long intergenic non-protein coding RNA 1015               6p22.1
13  100526772 TMEM110-MUSTN1                                                  TMEM110-MUSTN1 readthrough               3p21.1
14  100528020        FAM187A                                family with sequence similarity 187 member A             17q21.31
15  100532735    INO80B-WBP1                                     INO80B-WBP1 readthrough (NMD candidate)               2p13.1
16  100616371        MIR4746                                                               microRNA 4746              19p13.3

接下来是富集分析
[AppleScript] 纯文本查看 复制代码
ego_cc <- enrichGO(gene[,1], OrgDb=org.Hs.eg.db, ont = "CC", pAdjustMethod = "BH", minGSSize = 5, pvalueCutoff = 0.05, qvalueCutoff = 0.05, readable = TRUE)

结果如下
[AppleScript] 纯文本查看 复制代码
                   ID                  Description GeneRatio   BgRatio      pvalue  p.adjust    qvalue
GO:0016607 GO:0016607                nuclear speck    10/184 381/18698 0.004577086 0.4266889 0.4097028
GO:0005681 GO:0005681         spliceosomal complex     6/184 190/18698 0.011403463 0.4266889 0.4097028
GO:0035371 GO:0035371         microtubule plus-end     2/184  18/18698 0.013286690 0.4266889 0.4097028
GO:0098562 GO:0098562 cytoplasmic side of membrane     6/184 199/18698 0.014064395 0.4266889 0.4097028
GO:0071013 GO:0071013 catalytic step 2 spliceosome     4/184 101/18698 0.017696205 0.4266889 0.4097028
GO:0030863 GO:0030863        cortical cytoskeleton     4/184 105/18698 0.020106456 0.4266889 0.4097028


其中pvalue,qvalue,不符合p/qcutoff=0.01的要求,而且padj和pvalue都是相同的。我在用enrichGO的结果也不能直接用
[AppleScript] 纯文本查看 复制代码
barplot(ego_cc, showCategory=15,title="EnrichmentGO_CC")
我是按照http://www.biotrainee.com/thread-1084-1-1.html进行操作的。
请教各位我的问题出在哪里?



回复

使用道具 举报

4

主题

51

帖子

482

积分

管理员

Rank: 9Rank: 9Rank: 9

积分
482
发表于 2018-5-12 11:02:03 | 显示全部楼层
你的基因都有问题呀!!!
基因的ID不能太大,一百万以上的都是人类研究太少的基因, 你看名字就应该明白
一般来说要20000以下的,这些基因ID是有顺序的。
回复 支持 反对

使用道具 举报

1

主题

3

帖子

32

积分

新手上路

Rank: 1

积分
32
 楼主| 发表于 2018-5-12 18:37:03 | 显示全部楼层
bioinfotrainee 发表于 2018-5-12 11:02
你的基因都有问题呀!!!
基因的ID不能太大,一百万以上的都是人类研究太少的基因, 你看名字就应该明白
...

感谢您的回复,我之前也在奇怪为什么会这样。但是还有一个疑问,就是为什么qcutoff已经设置了但是在结果中没有过滤掉?
回复 支持 反对

使用道具 举报

0

主题

2

帖子

117

积分

注册会员

Rank: 2

积分
117
发表于 2019-3-23 15:48:48 | 显示全部楼层
tuanzi 发表于 2018-5-12 18:37
感谢您的回复,我之前也在奇怪为什么会这样。但是还有一个疑问,就是为什么qcutoff已经设置了但是在结果 ...

我遇到和楼主同一样的问题,请问楼主解决了吗
回复 支持 反对

使用道具 举报

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

本版积分规则

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

GMT+8, 2019-6-26 17:13 , Processed in 0.034912 second(s), 25 queries .

Powered by Discuz! X3.2

© 2001-2013 Comsenz Inc.