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