RNA-seqGTF

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_idchromosome_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_resgene_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_resmy_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_martmy_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)
上一篇下一篇

猜你喜欢

热点阅读