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

Posted on by

这个月初，我推出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这个来源的基因注释信息。

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

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

image-20191215192209639

这个信息就是前面我们使用的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]))

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

image-20191215192610061

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

> 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的注释信息显示的，都是一些奇奇怪怪的基因名字。

image-20191215194240964

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

image-20191215194334380

所以，我们的结论是，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万探针里面，占比不大，其实是可以忽略的，但是本着探索的精神，我们还是仔细看看：

image-20191215200559677

那就拿这个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工具也可以看到同样的结果：

image-20191215201305157

也就是说，这个探针居然是有方向的，

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

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

image-20191215202315897

但是我们看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 。

实际上，我还是没有解答，到底哪个注释方法是值得我们采纳的，后续推文，精彩继续吧。

GPL6947 Illumina HumanHT-12 V3.0 expression beadchip oligonucleotide beads Homo sapiens 49576 25185 523 Illumina Inc.
Accession Title Technology Taxonomy Data Rows Samples Count Series Count Contact
GPL10558 Illumina HumanHT-12 V4.0 expression beadchip oligonucleotide beads Homo sapiens 48107 83328 2449 Illumina Inc.
​

Comments are closed.