探针注释文件中没有基因名字怎么办？（二）

前面给大家介绍了，如果从GEO数据库下载的芯片数据，配套的☞探针注释文件中没有基因名字怎么办？生信交流群有小伙伴用上次讲到的AnnoProbe这个R包去注释GPL16686这个芯片平台就报错了，看来这个包也不是万能的。我们先来看看GPL16686的注释文件长什么样

从这个表中我们可以看到有，染色体上的坐标，正负链，还有GB_ACC。其实这个GB_ACC就是基因的转录本ID号了，它的全称是gene bank accession。如果你点击这里的GB_ACC，例如（NR_046018）就能跳转到相应的NCBI网页

你就可以看到这个基因的名字叫DDX11L1，有人说那我就一个一个去点不就可以得到每一探针对应的基因名字了吗？这也是一种方法，但是应该是最后不得已才会这么做。因为这里的探针有50000+个。了解小编的人都知道，小编一直奉行，“能有代码解决的问题，绝对不手动去做”。

接下来，小编就来用R解决这个问题。其实本质上就是一个ID转换问题，我们知道ID转换需要有一个map文件，里面包含不同ID之间的对应关系。这个map文件我们可以从UCSC这个网站上下载得到。

下面是下载地址

http://hgdownload.soe.ucsc.edu/goldenPath/hg38/database/

这里提供各种各样的文件，我们需要下载的是refGene.txt.gz这个文件。

下面是refGene.txt.gz这个文件的结构，里面包含了我们需要的两种ID号，分别在第2列和第13列

有了这个文件，接下来我们就可以用R代码来做ID转换了，给我们的探针注释文件加上一列gene symbol

#读入refGene.txt文件
anno=read.table("refGene.txt",sep="\t")
#提取第2列和第13列，分别为GB_ACC和gene symbol
anno=anno[,c(2,13)]
#给这两列命名
names(anno)=c("acc","symbol")
#去重复
anno=unique(anno)
#将GB_ACC设置成行名
rownames(anno)=anno$acc

#读取探针注释文件内容
probe=read.table("GPL16686.txt",header=T,sep="\t",comment.char = "#",stringsAsFactors = F)
#通过GB_ACC获取对应的gene symbol
symbol=as.character(anno[probe$GB_ACC,"symbol"])
#将NA转换成空
symbol[is.na(symbol)]=""
#将gene symbol这一列加入到原来的探针注释文件中
probe$symbol=symbol
#保存包含gene symbol的探针的注释文件
write.table(file="GPL16686_with_symbol.txt",probe,quote=F,sep="\t",row.names=F)

得到的结果如下，最后一列就是我们用R加上去的gene名字

参考下文获取完整数据和代码

探针注释文件中没有基因名字怎么办？（二）

欢迎关注公众号 “生信交流平台”，我们创建了交流群，供大家交流讨论生信问题。