从我写的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)
