从ensembl的ID到其转录本坐标

假设通过了某种分析(差异分析，peaks富集等等)得到了指定的基因集，但是是以ensembl数据库ID命名的，比如：

[1] "ENSMUSG00000000031" "ENSMUSG00000000792" "ENSMUSG00000001027" "ENSMUSG00000001508" "ENSMUSG00000002007"
[6] "ENSMUSG00000002985"


想得到它们这些ID对应的转录本的起始终止坐标。

第一步：把ensembl的ID对应到entrez ID

需要指定ensembl数据库ID来源的物种的org包

library(clusterProfiler)
library(org.Mm.eg.db)
geneList <- sample(mappedRkeys(org.Mm.egENSEMBL),100)
g=bitr( geneList , fromType = "ENSEMBL",
         toType = c("ENTREZID", "SYMBOL"),
         OrgDb = org.Mm.eg.db)[,2]


第二步：根据entrez ID来获取转录本信息

这个需要借助 TXDB 系列包

library("TxDb.Mmusculus.UCSC.mm10.knownGene")
txdb=TxDb.Mmusculus.UCSC.mm10.knownGene
columns(txdb)
select(txdb, keys = g, columns=c( "TXCHROM","TXSTART","TXEND","TXNAME", "TXSTRAND"), keytype="GENEID")


当然了，这个时候的转录本ID是UCSC数据库的，不过一般来说，坐标就够了，转录本ID其实是用不上的。

第4步：包装成函数

get_bed <- function(file){
 library(clusterProfiler)
 library(org.Mm.eg.db)
 library("TxDb.Mmusculus.UCSC.mm10.knownGene")
 txdb=TxDb.Mmusculus.UCSC.mm10.knownGene
 columns(txdb)
 g=bitr( geneList , fromType = "ENSEMBL",
              toType = c("ENTREZID", "SYMBOL"),
              OrgDb = org.Mm.eg.db)[,2]

bed=select(txdb, keys = g, columns=cl, keytype="GENEID")[,c(3,5,6,1,2,4)]
 write.table(bed,file = paste0(file,'.bed'),quote = F,col.names = F,row.names = F,sep = '\t')
}


画外音：坐标需要排序

当我把制作好的bed文件拿去做下游分析的时候，发现，报错，原因是bed文件并没有排序，O(∩_∩)O哈哈~

 keepChr=paste0('chr',c(1:22,"X","Y","M"))
 bed=bed[bed$TXCHROM %in% keepChr,]
 bed$TXCHROM=as.factor(bed$TXCHROM)
 levels(bed$TXCHROM)=keepChr
 sort_bed=bed[order(bed$TXCHROM,bed$TXSTART),]


所以又加上了几行代码，搞定！