R包biomaRt: 转换ID、注释基因、GO通路
2020-04-21 本文已影响0人
没有猫但是有猫饼
从我写的RNA-seq摸索:4. edgeR/limma/DESeq2差异基因分析→ggplot2作火山图→biomaRt转换ID并注释中提取出来这部分,方便以后修改补充
1 我们利用useMart()
函数选择“ENSEMBL_MART_ENSEMBL”,并将其赋值给my_mart
对象
library('biomaRt')
library("curl")
my_mart <-useMart("ensembl")
在ensembl数据库中包含了77个数据集,可用下面这样的方式查看
datasets <- listDatasets(my_mart)
View(datasets)
datasets
2 选择一个数据集datasset
,这里选人类的
my_dataset <- useDataset("hsapiens_gene_ensembl",
mart = my_mart)
3 💥根据ensembl ID
获取基因名、描述或染色体信息
💥💥💥这里前半部分有误!请一定往下看解决办法
my_newid <- getBM(attributes = c("ensembl_gene_id","external_gene_name","description","chromosome_name"),
filters = "ensembl_gene_id",
values = newinput,
mart = my_dataset)
image.png
💥这里一直报错,并且输出的为内容为0行
💥找到原因是:EBI数据库💥没有小数点💥,所以需要进一步替换为整数的形式。需要把小数点去掉!!这个很重要,所以需要加一个步骤
①还是将差异文件的行名提取出来
inputdata <- as.data.frame(row.names(deseq_res))
②这里将匹配到的.
以及后面的数字连续匹配并替换为空,并重新赋值,一定要是data.frame
格式
newinput <- as.data.frame(gsub("\\.\\d*", "", inputdata[,1]))
③getBM()
转换ID
1)
attributes
参数:用来指定输出的数据类型,就是你要什么,比如entrezgene,hgnc_id。忘记的话可以用listAttributes(你的dataset名字)
查看
2)filters
参数:用来指定数据的输入类型,比如你的原始信息是基因的ensembl ID
,并且有这些基因的染色体位置
信息,那么此处的filter就是ensembl_gene_id
和chromosome_name
等。
3)values
参数:就是你待转换ID
的数据
4)mart
参数:此前定义的数据库,此处就是my_dataset
那么在我这里:
attributes
:我想要输出"ensembl_gene_id",转换后的"external_gene_name",转换后的"description"
filters
:我的原始信息"ensembl_gene_id"
mart
:之前建立的数据库
listAttributes(你的dataset)
可以查看可供选择的attributes
listAttributes(my_dataset)
my_result <- getBM(attributes = c("ensembl_gene_id","external_gene_name","description"),
filters = "ensembl_gene_id",
values = newinput,
mart = my_dataset)
ID转换成功后
这样就完成了对
ensembl_id
的转化和注释
4 最后需要把结果文件deseq_res
和注释文件my_result
两者merge
起来
merge
需要有相同的gene_id
💥但是一定要看看自己文件里的gene_id
是不是一致,如果有一个为小数,就要再添加一列取整后的gene_id
① deseq_res
中gene_id
有小数点 所以再加一列变成new_deseq_res
,新增加的列名为gene_new_id
new_deseq_res <- as.data.frame(deseq_res)
new_deseq_res$gene_new_id <- gsub("\\.\\d*", "", deseq_res$gene_id)
② 修改一下列名,把含有小数点的列命名为gene_all_id
,取整后的为gene_id
,这一步是为了方便merge
colnames(new_deseq_res) <- c('baseMean', 'log2FoldChange','lfcSE','stat','pvalue','padj','gene_all_id','gene_id')
new_deseq_res
③ merge
两个文件,即new_deseq_res
和my_resullt
,生成final_res
文件
by = intersect(names(x), names(y))
为取两个文件所有列名中列名相同的那列!
final_res <- merge(my_result, new_deseq_res, by = intersect(names(my_result), names(new_deseq_res)))
write.table(final_res, 'C:/Users/wang/Desktop/final_result.txt',row.names = FALSE, sep = '\t', quote = FALSE)
结果文件
5 还可以找到某个基因所在的通路GO号
① 选出要查找的基因
#举个例子
entrez = c("673", "837")
② 利用ensembl
构建my_mart
和my_dataset
my_mart <-useMart("ensembl")
#`listDatasets()`可以查看可用的`datasets`
datasets <- listDatasets(my_mart)
View(datasets)
#构建`my_dataset`
my_dataset <- useDataset("hsapiens_gene_ensembl",
mart = my_mart)
③ 查看可输出的attributes
listAttributes(my_dataset)
④ 查找GOid
GOid <- getBM(attributes = c('entrezgene_id', 'go_id'),
filters = 'entrezgene_id',
values = entrez,
mart = my_dataset)
结果
6 与5
相反,可以通过所在的通路GO号
找到某个基因
getBM(attributes = c('entrezgene_id', 'ensembl_gene_id'),
filters = 'go',
values = 'GO:0005524',
mart = my_dataset)