GPL平台的soft文件提供的注释信息到底准确吗

这个月初,我推出3个R包,

  • 第一个是整合全部的bioconductor里面的芯片探针注释包。

  • 第二个是整合全部GPL的soft文件里面的芯片探针注释包。

  • 第三个是下载全部的GPL的soft文件里面的探针碱基序列比对后注释包。

配合着详细的介绍:

因为这些包暂时托管在GitHub平台,但是非常多的朋友访问GitHub困难,尤其是我打包了好几百个GPL平台的注释信息后, 我的GitHub包变得非常臃肿,大家下载安装困难,所以我重新写一个精简包。也在:芯片探针ID的基因注释以前很麻烦 和 :芯片探针序列的基因注释已经无需你自己亲自做了, 里面详细介绍了。

最重要的是idmap函数

安装方法说到过:芯片探针序列的基因注释已经无需你自己亲自做了,  使用起来也非常简单:

library(AnnoProbe)
ids=idmap('GPL570',type = 'soft')
head(ids)

仅仅是一句话,就拿到了这个平台的探针的注释信息。需要注意的是,这个函数的type参数,其实是有3个选择,这里我演示的是选择soft这个来源的基因注释信息。

然后就有粉丝反馈说她自己感兴趣的朋友,不同来源出现了很诡异的结果,让她无法理解,希望我解读一下:

  • https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GPL6947

在GEO数据库里面访问该平台的主页,可以看到下面的注释信息

这个信息就是前面我们使用的idmap函数的type参数选择了soft这个选项后的结果。

看看3个参数的差异

那就使用 idmap函数吧,这里的type参数选择3次,返回不同的结果:

library(AnnoProbe)
ids1=idmap('GPL6947',type = 'bioc')
head(ids1)
ids2=idmap('GPL6947',type = 'soft')
head(ids2)
ids3=idmap('GPL6947',type = 'pipe')
head(ids3)

length(unique(ids1[,1]))
length(unique(ids2[,1]))
length(unique(ids3[,1]))

简单从数量上,就能看到显著的差异了

而且可以看到,其中第三种注释结果里面的探针是冗余的,也就是说,一个探针会对应到多个基因

> length(unique(ids1[,1]))
[1] 29457
> length(unique(ids2[,1]))
[1] 49576
> length(unique(ids3[,1]))
[1] 39892

其实这个是合理的状态,因为很多基因的坐标其实是有overlap的,比如miRNA在某个基因内部,某个基因的假基因或者lncRNA也可能会有overlap的。

先比较bioc和soft的注释差异

其中bioc的来源就是该平台对应的bioconductor里面的芯片探针注释包的信息的提取,而soft就是我们前面说的在GEO数据库里面访问该平台的主页看到的注释信息的提取。

ids2_filter=ids2[nchar(ids2[,2])>1,]
m1=merge(ids1,ids2_filter,by.x = 'probe_id',by.y = 'ID')
table(m1[,2]==m1[,3]) 
#FALSE  TRUE 
# 6124 22967 
m1_not=m1[m1[,2] !=m1[,3],]
tmp1=annoGene(m1_not[,2],'SYMBOL')
tmp2=annoGene(m1_not[,3],'SYMBOL')

可以看到,两个注释结果冲突的比例还好,毕竟2.2万多探针都是一致的,然后我们仔细检查了那6千多个冲突了的 探针,发现很诡异的现象。右边,已也就是来源于soft的注释信息显示的,都是一些奇奇怪怪的基因名字。

所以我对它们进行了gencode数据库的注释,发现来源于soft的注释信息的那些奇奇怪怪的基因名字,都是无法去找到记录的,手动搜索一下,发现都是一些被废弃掉的基因名字。

所以,我们的结论是,soft就是我们前面说的在GEO数据库里面访问该平台的主页看到的注释信息的提取,应该是非常的过时了。选择这个方法,是下下策。

其次比较bioc和pipe的注释差异

其中bioc的来源就是该平台对应的bioconductor里面的芯片探针注释包的信息的提取,而pipe是我们自己下载全部的GPL的soft文件里面的探针碱基序列比对后注释结果,理论上这个应该是最新的,但是需要实际数据来证明!

table( paste(ids1[,1],ids1[,2],sep = '_') %in% paste(ids3[,1],ids3[,2],sep = '_') )
# FALSE  TRUE 
# 2826 26631 
ids1_not_ids3=ids1[!paste(ids1[,1],ids1[,2],sep = '_') %in% paste(ids3[,1],ids3[,2],sep = '_'),]
m2=merge(ids1_not_ids3,ids3,by='probe_id')
tmp1=annoGene(m2[,2],'SYMBOL')
tmp2=annoGene(m2[,3],'SYMBOL')

m2_not=m2[!m2[,2] %in% tmp1[,1],]
m2_ok=m2[m2[,2] %in% tmp1[,1],]

可以看到,在bioc来源的注释结果里面只有不到3千的探针是冲突了,同样的道理我们进行gencode数据库的注释这样的检验,发现是前者有38%的探针注释到的基因名字是无法在gencode数据库里面进行映射。

但是其余的可以映射的那些探针,就说明了bioc和pipe的注释是真的冲突了

探索冲突的那些探针注释

因为冲突的探针有八百多个,虽然说在3万探针里面,占比不大,其实是可以忽略的,但是本着探索的精神,我们还是仔细看看:

那就拿这个ILMN_1652209探针吧,可以看到它是可以比对到参考基因组的3个位置:

ILMN_1652209    16  chr14   19083887    1   50M TGAATCACAGCACCTTCAAATATGGCCACTTCCAAGGCAGAATTAAACAT
ILMN_1652209    256 chr14   19316112    1   50M ATGTTTAATTCTGCCTTGGAAGTGGCCATATTTGAAGGTGCTGTGATTCA
ILMN_1652209    272 chr22   15806451    1   50M TGAATCACAGCACCTTCAAATATGGCCACTTCCAAGGCAGAATTAAACAT

在blat工具也可以看到同样的结果:

也就是说,这个探针居然是有方向的,至少我的注释里面并没有去考虑这一点。

chr14    19062316    19131167    DUXAP9  transcribed_processed_pseudogene
chr14    19268853    19337730    DUXAP10 transcribed_processed_pseudogene
chr14    19301704    19316914    BMS1P17 transcribed_unprocessed_pseudogene
chr22    15784959    15829984    DUXAP8  lncRNA
chr22    15805263    15820884    BMS1P22 transcribed_unprocessed_pseudogene

那么,为什么bioc来源会给出一个BMS1的注释呢?连染色体都不一样啊!

chr10    42782795    42834937    BMS1    protein_coding

这个实在是很诡异哦,然后我谷歌了一下,哈哈,大家都使用了这个错误的注释!

但是我们看soft来源,这个探针是LOC441969,简单搜索可以知道,它与BMS1并不是一回事,就是在14q11.1位置,并不是在chr10上面。

LOC441969 valid_symbol LOC441969 similar to BMS1-like, ribosome assembly protein; ribosome biogenesis protein BMS1 homolog. 

实际上这个LOC441969的 中英文全称:类似核糖体生物发生蛋白BMS1 同源物 similar to Ribosome biogenesis protein BMS1 homolog,OMIM/位置:14q11.1

但是冲突并不可怕,因为大部分的基因其实并不是大家感兴趣的,研究下去也走不通的。所以重点是看它们到底是能不能注释到基因ID甚至GO/KEGG数据库功能。

了解Illumina HumanHT-12 平台

还是蛮出名的表达芯片平台,尤其是其 Illumina HumanHT-12 V4.0 expression beadchip ,数据集之多可以在表达芯片领域排到前10 。

(0)

相关推荐