skr！GEO芯片数据的探针ID转换

今天这个帖子，对我来说很有意义，我在最后要介绍一位心中的男神。

我每天努力学习就是为了消除心中的疑惑。 几年前，当我跟着别人的代码跑了个流程后，突破不了的是两个事情

• 1.如何不借助 excel 的拖拽，对芯片进行分组，

• 2.如何方便优雅全能地进行探针转换。

那个时候，R语言基础很差，处理不了数据，很有挫败感，所以就停止了R语言的学习。 直到我碰到了那个男人。

现在我解决了这个问题，分享给大家。

第一种，我们可以直接用平台的数据

进入官网 https://www.ncbi.nlm.nih.gov/geo/ 

知道平台是GPL6244 

这时候我们进入R语言，用GEOquery中的getGEO可以获得探针和基因名的信息 网络不好的可以直接略过

1. library(GEOquery)

2. GPL6244 <>'GPL6244',destdir ='.')

转换成数据框形式,有3万行，12列

GPL6244_anno <->Table(GPL6244)

查看内容，我们发现基因名称藏在了gene_assignment这一列的中间 

所以我们要把他和第一列id提取出来

1. library(dplyr)

2. library(tidyr)

3. probe2symbol_df <- gpl6244_anno="" %="">%

4.  select(ID,gene_assignment) %>%

5.  filter(gene_assignment != '---') %>%

6.  separate(gene_assignment,c('drop','symbol'),sep='//') %>%

7.  select(-drop)

看一下，数据已经被提取出来了。 

假如getGEO这一步网络不好呢

library(GEOquery)
GPL6244 <>'GPL6244',destdir ='.')

我们在这个一开始的这个页面下载平台的soft文件 

点击soft文件 

下载解压 然后用data.table这个包中的fread即可阅读进来，注意，skip这个参数十分重要！！

1. GPL6244_anno <>'./GSE42872_family.soft/GSE42872_family.soft',skip ='ID')

得到了GPL6244_anno,我们又可以运行下面的代码提出探针和基因名称对应的关系了。

library(dplyr)
library(tidyr)
probe2symbol_df <- gpl6244_anno="" %="">% 
  select(ID,gene_assignment) %>% 
  filter(gene_assignment != '---') %>% 
  separate(gene_assignment,c('drop','symbol'),sep='//') %>% 
  select(-drop)

第二种，使用R包

R有很多注释包，首先我们要知道什么平台用什么R包。 我这里整理了一个对应关系 

读入R语言

1. platformMap <->'platformMap.txt')

为了方便查询，我写了以下的代码，每次只要修改平台的名称即可

index <->'GPL6244'
paste0(platformMap$bioc_package[grep(index,platformMap$gpl)],'.db')

输出的是，

'hugene10sttranscriptcluster.db'

安装并加载这个包

1. BiocInstaller::biocLite('hugene10sttranscriptcluster.db')

2. library(hugene10sttranscriptcluster.db)

使用toTable函数找到对应关系

probe2symbol_df <->get('hugene10sttranscriptclusterSYMBOL'))

那探针ID转换的两种方法就讲完了

正式进行探针ID之间的转换

我们现在用的是矩阵文件exprSet, 

我们要达到的效果是这个样的 

其实很简单，就是把两个数据框按照探针名称合并就可以了， 同时在这里，我们还需要把重复的探针去掉，方法就是计算每一个探针的平均值，留下最大的那个。

先调整两个数据框的探针名字和类型，方便合并

1. names(exprSet)[1] <->1]

2. exprSet$probe_id <->as.character(exprSet$probe_id)

然后，

重点就来了，下面的代码恢弘大气，一次性完成这两个事情。

library(dplyr)
exprSet <- exprset="" %="">% 
  inner_join(probe2symbol_df,by='probe_id') %>% #合并探针的信息
  select(-probe_id) %>% #去掉多余信息
  select(symbol, everything()) %>% #重新排列，
  mutate(rowMean =rowMeans(.[grep('GSM', names(.))])) %>% #求出平均数(这边的.真的是画龙点睛)
  arrange(desc(rowMean)) %>% #把表达量的平均值按从大到小排序
  distinct(symbol,.keep_all = T) %>% # symbol留下第一个
  select(-rowMean) %>% #反向选择去除rowMean这一列
  tibble::column_to_rownames(colnames(.)[1]) # 把第一列变成行名并删除

真的很skr！

这个代码的产生过程，以后再说。

这里的%>% 相当于Linux中的xargs，

Hadley Wickham真的是天才！

如果不是他的这个神包tidyverse以及健明给的鼓励 我可能到现在还入不了R的门！