|
前面給大家介紹了，如果從GEO數(shù)據(jù)庫下載的芯片數(shù)據(jù)，配套的?探針注釋文件中沒有基因名字怎么辦？有小伙伴用上次講到的AnnoProbe這個(gè)R包去注釋GPL16686這個(gè)芯片平臺(tái)就報(bào)錯(cuò)了，看來這個(gè)包也不是萬能的。我們先來看看GPL16686的注釋文件長(zhǎng)什么樣 從這個(gè)表中我們可以看到有，染色體上的坐標(biāo)，正負(fù)鏈，還有GB_ACC。其實(shí)這個(gè)GB_ACC就是基因的轉(zhuǎn)錄本ID號(hào)了，它的全稱是gene bank accession。如果你點(diǎn)擊這里的GB_ACC，例如（NR_046018）就能跳轉(zhuǎn)到相應(yīng)的NCBI網(wǎng)頁 你就可以看到這個(gè)基因的名字叫DDX11L1，有人說那我就一個(gè)一個(gè)去點(diǎn)不就可以得到每一探針對(duì)應(yīng)的基因名字了嗎？這也是一種方法，但是應(yīng)該是最后不得已才會(huì)這么做。因?yàn)檫@里的探針有50000+個(gè)。了解小編的人都知道，小編一直奉行，“能有代碼解決的問題，絕對(duì)不手動(dòng)去做”。 接下來，小編就來用R解決這個(gè)問題。其實(shí)本質(zhì)上就是一個(gè)ID轉(zhuǎn)換問題，我們知道ID轉(zhuǎn)換需要有一個(gè)map文件，里面包含不同ID之間的對(duì)應(yīng)關(guān)系。這個(gè)map文件我們可以從UCSC這個(gè)網(wǎng)站上下載得到。 下面是下載地址 http://hgdownload.soe./goldenPath/hg38/database/ 這里提供各種各樣的文件，我們需要下載的是refGene.txt.gz這個(gè)文件。 下面是refGene.txt.gz這個(gè)文件的結(jié)構(gòu)，里面包含了我們需要的兩種ID號(hào)，分別在第2列和第13列 有了這個(gè)文件，接下來我們就可以用R代碼來做ID轉(zhuǎn)換了，給我們的探針注釋文件加上一列g(shù)ene 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")#去重復(fù)anno=unique(anno)#將GB_ACC設(shè)置成行名rownames(anno)=anno$acc#讀取探針注釋文件內(nèi)容probe=read.table("GPL16686.txt",header=T,sep="\t",comment.char = "#",stringsAsFactors = F)#通過GB_ACC獲取對(duì)應(yīng)的gene symbolsymbol=as.character(anno[probe$GB_ACC,"symbol"])#將NA轉(zhuǎn)換成空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)得到的結(jié)果如下，最后一列就是我們用R加上去的gene名字 |
|
|
