|
前面给大家介绍了，如果从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./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 symbolanno=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 symbolsymbol=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名字 |
|
|
