映射codelink芯片ID
GE公司codelink芯片(这个破产品已经停产了)上的探针ID,需要转成别的ID,看了一些在线的转换ID工具,都不支持,探针ID基本上都只支持affy的。
想起了bioconductor里的biomaRt,这个包可以检索BioMart数据库,这个数据库里有N多种ID。试了一下,果然没问题。
#加载biomaRt包
library(biomaRt)
#选取数据库
ensembl = useMart("ensembl", dataset = "hsapiens_gene_ensembl")
#获取可供检索的attribute
attributes=listAttributes(mart=ensembl)
#看一下attributes里是否有codelink,得到以下的结果。
attributes[grep("codelink",attributes[,1]), ]
name description
21 codelink Codelink ID
#再搜一下其它我想要的ID
#读进含有codelink ID的文件。
mrna_id <- read.table("mrna_id.txt")
#进行ID映射。
idmap <- getBM(attributes=c("codelink","refseq_dna","external_gene_id",
"embl","hgnc_id","hgnc_symbol"),
filters ="codelink", values=mrna_id[,1],
mart=ensembl,output="list")
乱78糟的ID都可以通过biomart进行ID映射。。生物上各种各样的ID实在是太乱了。。。
没想到程序跑完竟然报错,囧~~~
找到了一个codelink的注释包,hwgcod:
GE CodeLink Human Whole Genome Bioarray (~55 000 human genes) Annotation Data (hwgcod)
写了下面段小代码,算是解决了。5万多个点,有些没对上,有些却是一对多。这种结果很正常。。囧~~~
library(hwgcod)
mrna_id <- read.table("mrna_id.txt",header=FALSE, stringsAsFactors=FALSE)
accnum <- as.list(hwgcodACCNUM)
accnum <- accnum[!is.na(accnum)]
map <- function(name){
accnum[[name]]
}
mrna_acc <- lapply(mrna_id[,1],map)
writetofile <- function(idx){
line = ""
for(n in mrna_acc[[idx]])
line=paste(line,n)
write(line,"mrna_acc_id.txt",append=TRUE)
}
lapply(1:length(mrna_id[[1]]),writetofile)