GEO芯片数据的探针ID转换

转自“果子学生信”

一、首先获取用于转换的probe2symbol文件：

第１种方法 直接用平台数据(网络不好的可以直接略过)

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

转换成数据框形式

GPL6244_anno <- Table(GPL6244)

基因名隐藏在gene_assignment列的中间，将其与第一列id提取出来

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

如果网络不好，点击下载_family.soft.gz文件，解压后用data.table包中的fread即可读取，注意：：：skip这个参数十分重要！！

GPL6244_anno <-data.table::fread("./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包

platformMap <- data.table::fread("platformMap.txt")
index <- "GPL6244"
paste0(platformMap$bioc_package[grep(index,platformMap$gpl)],".db")

安装并加载输出结果

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

使用toTable函数找到对应关系

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

通过上述两种方式，均可获得用于基因注释的probe2symbol_df，就可进行探针ID之间的转换了

二、进行探针ID之间的转换

１.获取矩阵文件exprSet
２.调整两个数据框的探针名字和类型，方便合并

names(exprSet)[1] <- names(probe2symbol_df)[1]
exprSet$probe_id <- as.character(exprSet$probe_id)

3.①将两个数据框按照探针名称合并；②去除重复探针（计算每一个探针的平均值，留下最大的那个）

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]) # 把第一列变成行名并删除

作者最后感谢了Hadley Wickham的神包tidyverse以及健明！！