【原】学徒笔记——芯片数据的注释文件获取

目前芯片数据的分析流程都可以通过AnnoProbe包简单的完成上游分析，包括表达矩阵获取、分组方案的构建和数据注释，但是也存在一些平台的数据无法被该包直接获取。AnnoProbe获取注释信息的方式是通过对信息文件中的GPL字段信息，直接从数据库下载相关编号，但是意外总会发生。以下是我在几个GSE数据分析中遇到的情况总结：

一、文章出错

写文章确实是个严谨的事，但是万一呢，有时候做个脑瘤的分析整个糖尿病的编号在里面，也是大受震撼，一般来说起码都是一个物种的，平台一不一致问题不大的样子。通篇检查一下，可能就是差那么一位数，但是一定有写对的地方。

二、手动下载文件

直接在 GEO 平台搜索对应的 GPL 编号，通过 AnnoProbe 包的 checkGPL 函数检查一下，返回的 FALSE，即R包的数据库里找不到这个平台的注释文件，所以要去手动下载然后读取，但是这种野生的soft文件都比较放纵：

1、symble编号淹没在各种信息里

举个栗子，GSE90604，通过GEO平台的查询可以得知这组数据集的平台是GPL17692

library(AnnoProbe)
checkGPL("GPL17692")
[1] FALSE


然后通过GEO平台去看看到底怎么回事，

是的，并不是下面的这个soft_formatted_family_file(s)文件，而是这个button，实际上这个Data_table就是下载文件的预览，可以看到gene_assignment一列里用 //符号隔开了很多的信息，选择SYMBLE的那一段，提取出来

library(data.table)
b=fread('GPL17586-45144.txt',data.table = F)[,c(2,8)] #D读取下载到本地的soft文件
head(b)
library(stringr)
b$gene=str_split(b$gene_assignment,' // ',simplify = T)[,2]将基因信息列按照 // 符号分离，抽取symble编号
ids=b[c(1,3)]


这样子就可以提取到一个正常的soft文件了，后续按照既有流程走下去就好。

2、就是不给你常规symble号，就是玩儿

事情还没结束，同一篇文章里还有另外一个坑，GSE49810，常规检查一下

library(AnnoProbe)
checkGPL("GPL15648")
[1] FALSE


问题不大，再来，下载soft文件之后，看文件预览

好像没什么乱七八糟的符号，好像问题不大，开始常规操作

a = getGEO(gpl_number,destdir = ".") #读取本地文件到一个列表
b = a@dataTable@table         #提取需要的元素
colnames(b)
ids = b[,c("ID","Description")]
colnames(ids) = c("probe_id","symbol")#指定列名
ids = ids[ids$symbol!="" & !str_detect(ids$symbol,"///"),] #去除空的信息行


注释好了，做差异分析，然后转换ID准备做富集，但是注释失败，全部转换失败！意外总是有的，关键在于问题在哪儿。回去看这个预览文件，基因信息部分并不是SYMBLE编号，symble编号都是几个大写字母，都是一些缩写，然而预览文件里完全不是，都是基因全名。大佬在这时发出了光辉，clusterprofiler预判了这种情况，bitr函数中fromtype参数允许多种输入参数，先查看数据库中有哪些数据

keytypes(org.Hs.eg.db)
 [1] "ACCNUM"       "ALIAS"        "ENSEMBL"      "ENSEMBLPROT" 
 [5] "ENSEMBLTRANS" "ENTREZID"     "ENZYME"       "EVIDENCE"    
 [9] "EVIDENCEALL"  "GENENAME"     "GENETYPE"     "GO"          
[13] "GOALL"        "IPI"          "MAP"          "OMIM"        
[17] "ONTOLOGY"     "ONTOLOGYALL"  "PATH"         "PFAM"        
[21] "PMID"         "PROSITE"      "REFSEQ"       "SYMBOL"      
[25] "UCSCKG"       "UNIPROT"     


其中就有GENENAME的选项

library(clusterProfiler)
library(org.Hs.eg.db)
df <- bitr(unique(deg$symbol), fromType = "GENENAME",
           toType = c( "ENTREZID"),
           OrgDb = org.Hs.eg.db)
df2 <- bitr(unique(deg$symbol), fromType = "GENENAME",
            toType = c( "SYMBOL"),
            OrgDb = org.Hs.eg.db)


然后查看，信息出来了，但是仍然提示有64%左右的数据没有注释上，所以还是有大量的数据损失了，因此严格来说只成功了一半，还有一种比较将就的通关方式，即使用bitr函数的 drop = F 参数，这个参数的功能是保留所有没注释上的基因，但是会造成后续的分析结果没有基因名，操作空间有限。

（其实学徒的知识背景仍然是不足，这个芯片里面的探针本身就是基因ID，就是基因的entrez ID ，完全不需要通过 GENENAME来进行桥接 ）

3、其实是有的，但是又不完全有

这次有问题的是GPL14550这个平台，按照常规的流程先用geochina函数去统一下载所有的表达矩阵和相关信息，再根据idmap函数自动下载注释文件，显示下载失败，找不到网址，那就上GEO数据库里去找吧，然后就看到了

好家伙，数据库里提供的文件直接什么都没有，整个 GENE_SYMBOL 列都是空的，不要慌，遇到什么奇葩都先打开 google 看看，然后看到了技能树论坛里曾经有人做过这个平台的注释（遇到arguments must have same length-QA-生信技能树 (biotrainee.com)）大大的不对劲，仔细看这个帖子里描述，他并没有遇到下载的问题，所以试试老一点命令，因为流程里的是操作是通过 idmap() 函数来做的，这个函数隶属于 AnnoProbe ，然后尝试如下命令：

library(GEOquery)
> checkGPL("GPL14550")   ##通过命令检查这个平台是否存在下载列表里
[1] TRUE     ##TRUE表示存在，对应FALSE显示数据库里没有收录

> getGEO("GPL14550")  ##然后下载成功了
File stored at: 
C:\Users\Dreamon\AppData\Local\Temp\Rtmpa61tEz/GPL14550.soft


然后去下载地址查看一下下载文件，是正常的！然后通过读取本地注释文件的方式导入流程

b = getGEO("GPL14550",destdir = "./")
d = b@dataTable@table


所以这个注释实际在数据库里是收纳了的，只是这个数据库不是 GEO 是 GEOquery 的，所以姜还是老的辣。这边建议不是很常见的GPL编号如果 idmap() 找不到的话，请给 checkGPL() 一个机会，也给 getGEO() 一个面子。

4、Google

实际上在最近的学习中，同学间这类的情况也偶有出现，不过解决办法没有这个折腾，很直接，谷歌一下，相信自己，坑里一定有垫底的，也一定会有出坑的梯子！

后续还有什么奇葩注释的操作再更新。

文末友情推荐

与十万人一起学生信，你值得拥有下面的学习班：

• 数据挖掘（GEO,TCGA,单细胞）2021第5期
• 生信爆款入门-2021第5期