关于clusterprofiler的一些疑问（结果和教程不一致）

tuanzi

首先我读取了基因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进行操作的。
请教各位我的问题出在哪里？

bioinfotrainee

你的基因都有问题呀！！！
基因的ID不能太大，一百万以上的都是人类研究太少的基因， 你看名字就应该明白
一般来说要20000以下的，这些基因ID是有顺序的。

tuanzi

bioinfotrainee 发表于 2018-5-12 11:02
你的基因都有问题呀！！！
基因的ID不能太大，一百万以上的都是人类研究太少的基因， 你看名字就应该明白
...

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

沉默的胡杨

tuanzi 发表于 2018-5-12 18:37
感谢您的回复，我之前也在奇怪为什么会这样。但是还有一个疑问，就是为什么qcutoff已经设置了但是在结果 ...

我遇到和楼主同一样的问题，请问楼主解决了吗