KEGG的翻车之旅
2019-12-01 本文已影响0人
佳名
7,744的基因集竟然一个通路都没有富集出来。
kegg <- enrichKEGG(gene= testdata$entrezgene_id,
organism = "hsa",
keyType = "kegg")
于是使用在线富集
kegg <- enrichKEGG(gene= testdata$entrezgene_id,
organism = "hsa",
keyType = "kegg",
use_internal_data = TRUE)
结果提示没有安装KEGG.db，于是安装。
BiocManager::install('KEGG.db')
重新富集
kegg <- enrichKEGG(gene= testdata$entrezgene_id,
organism = "hsa",
keyType = "kegg",
use_internal_data = TRUE)
富集结果为13个通路，修改参数，重新富集
kegg <- enrichKEGG(gene= testdata$entrezgene_id,
organism = "hsa",
keyType = "kegg", maxGSSize =100000,
pvalueCutoff = 0.5,
use_internal_data = TRUE)
富集结果为49个通路
testdata1 <- data.frame(testdata$log2FoldChange)
testdata1
row.names(testdata1) <- testdata$entrezgene
结果报错：
Error in `.rowNamesDF<-`(x, value = value) : 不允许有重复的'row.names'
In addition: Warning message:
non-unique values when setting 'row.names':
重新ID转换
human <- useMart("ensembl","hsapiens_gene_ensembl")
geneID <- getBM(attributes=c("external_gene_name","entrezgene_description","entrezgene_id"),
filters = "external_gene_name",values = resdata[,1], mart = human)
行名一直有重复项，换clusterProfiler
ENTREZID<- bitr(resdata[,1], fromType = "SYMBOL", toType=c("ENTREZID","GENENAME"),
OrgDb = org.Hs.eg.db, drop = TRUE)
成功了，但是ID转换没有getBM多，仔细比对，原来像MT-ND这种线粒体上的基因，bitr函数不识别，bitr识别的是ND。biomaRt比较坑的是，许多ID转换出来是错的，如“SMN1”和“SMN2"这两个基因得到的entrezgene_id都是6606，甚至毫无相干的两个基因，TRIM74和STAG3L3得到的都是378108，连entrezgene_description都是一样的。
使用HTseq的时候，列名为gene_name，出现27353个统计单位；
列名为gene_id，出现58884个统计单位，也就是说gene_name和gene_id并不是一一对应的关系。使用gene_id遇到的情况是ID转换时出现两个gene_name。
ID转换的时候使用clusterProfiler的bitr函数。
使用两份数据进行测试，一份是小鼠，一份是人
Note：resdata1[,1]为数据框第一列，即我要转换的ID，小鼠20009个ID，人 18732个ID。
1 小鼠
1.1 使用biomaRt包
library('biomaRt')
mart <- useMart("ensembl","mmusculus_gene_ensembl")
ENTREZID <- getBM(attributes=c("external_gene_name",
"entrezgene_description",
"ensembl_gene_id",
"entrezgene_id"),
filters = "ensembl_gene_id",
values = resdata1[,1], mart = mart)
#去除重复的ENSEMBL
ENTREZID=ENTREZID[!duplicated(ENTREZID$ensembl_gene_id),]
#为了和resdata1合并，需要重新组成数据框
data <- data.frame(Row.names=ENTREZID$ensembl_gene_id,
entrezgene_id,
gene_name=ENTREZID$external_gene_name,
description=ENTREZID$entrezgene_description)
#合并
data1 <-merge(data,resdata1,by="Row.names")
#测试
testdata <-subset(data1,entrezgene_id!= 'NA')
nrow(testdata)
#[1] 14519
1.2 使用bitr函数
library('clusterProfiler')
library('org.Mm.eg.db')#version=3.10.0
ENTREZID<- bitr(resdata1[,1], fromType = "ENSEMBL",
toType=c('SYMBOL',"ENTREZID","GENENAME"),
OrgDb = org.Mm.eg.db, drop = FALSE)
#去除重复的ENSEMBL
ENTREZID=ENTREZID[!duplicated(ENTREZID$ENSEMBL),]
data <- data.frame(Row.names=ENTREZID$ENSEMBL,
SYMBOL=ENTREZID$SYMBOL,
entrezgene_id=ENTREZID$ENTREZID,
description=ENTREZID$GENENAME)
#合并
data1 <-merge(data,resdata1,by="Row.names")
#测试
testdata <-subset(data1,entrezgene_id!= 'NA')
nrow(testdata)
#[1] 15554
2 人类
2.1 使用biomaRt包
library('biomaRt')
mart <- useMart("ensembl","hsapiens_gene_ensembl")
ENTREZID <- getBM(attributes=c("external_gene_name",
"entrezgene_description",
"ensembl_gene_id",
"entrezgene_id"),
filters = "ensembl_gene_id",
values = resdata1[,1], mart = mart)
#去除重复的ENSEMBL
ENTREZID=ENTREZID[!duplicated(ENTREZID$ensembl_gene_id),]
data <- data.frame(Row.names=ENTREZID$ensembl_gene_id,
entrezgene_id,
gene_name=ENTREZID$external_gene_name,
description=ENTREZID$entrezgene_description)
#合并
data1 <-merge(data,resdata1,by="Row.names")
#测试
testdata <-subset(data1,entrezgene_id!= 'NA')
nrow(testdata)
#[1] 13472
2.2 使用bitr函数
library('clusterProfiler')
library('org.Hs.eg.db')#version=3.10.0
ENTREZID<- bitr(resdata1[,1], fromType = "ENSEMBL",
toType=c('SYMBOL',"ENTREZID","GENENAME"),
OrgDb = org.Hs.eg.db, drop = FALSE)
#去除重复的ENSEMBL
ENTREZID=ENTREZID[!duplicated(ENTREZID$ENSEMBL),]
data <- data.frame(Row.names=ENTREZID$ENSEMBL,
SYMBOL=ENTREZID$SYMBOL,
entrezgene_id=ENTREZID$ENTREZID,
description=ENTREZID$GENENAME)
#合并
data1 <-merge(data,resdata1,by="Row.names")
#测试
testdata <-subset(data1,entrezgene_id!= 'NA')
nrow(testdata)
#[1] 14896
