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

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

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

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

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

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

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

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

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

知道平台是GPL6244 

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

1. library(GEOquery)

2. GPL6244 <-getGEO('GPL6244',destdir =".")

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

1. 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这一步网络不好呢

1. library(GEOquery)

2. GPL6244 <-getGEO('GPL6244',destdir =".")

我们在这个一开始的这个页面下载平台的soft文件 点击soft文件 下载解压 然后用data.table这个包中的fread即可阅读进来，注意，skip这个参数十分重要！！

1. GPL6244_anno <-data.table::fread("./GSE42872_family.soft/GSE42872_family.soft",skip ="ID")

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

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)

第二种，使用R包

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

1. platformMap <- data.table::fread("platformMap.txt")

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

1. index <- "GPL6244"

2. paste0(platformMap$bioc_package[grep(index,platformMap$gpl)],".db")

输出的是，

"hugene10sttranscriptcluster.db"

安装并加载这个包

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

2. library(hugene10sttranscriptcluster.db)

使用toTable函数找到对应关系

1. probe2symbol_df <- toTable(get("hugene10sttranscriptclusterSYMBOL"))

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

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

我们现在用的是矩阵文件exprSet, 我们要达到的效果是这个样的 

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

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

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

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

然后，

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

1. library(dplyr)

2. exprSet <- exprSet %>%

3.  inner_join(probe2symbol_df,by="probe_id") %>% #合并探针的信息

4.  select(-probe_id) %>% #去掉多余信息

5.  select(symbol, everything()) %>% #重新排列，

6.  mutate(rowMean =rowMeans(.[grep("GSM", names(.))])) %>% #求出平均数(这边的.真的是画龙点睛)

7.  arrange(desc(rowMean)) %>% #把表达量的平均值按从大到小排序

8.  distinct(symbol,.keep_all = T) %>% # symbol留下第一个

9.  select(-rowMean) %>% #反向选择去除rowMean这一列

10.  tibble::column_to_rownames(colnames(.)[1]) # 把第一列变成行名并删除

真的很skr！

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

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

Hadley Wickham真的是天才！

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