转自“果子学生信”
一、首先获取用于转换的probe2symbol文件:
第1种方法 直接用平台数据(网络不好的可以直接略过)
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)
第2种方法 使用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之间的转换
1.获取矩阵文件exprSet
2.调整两个数据框的探针名字和类型,方便合并
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以及健明!!
网友评论