hello,hello!各位小伙伴们大家好,我是大家的小编豆豆,最近因为南京疫情,导致很多学校被封了,很多实验样品进不来,所以很多做实验的同学开始学生信。前两天,我妹妹在做GEO数据分析时遇到一点问题,就是将芯片数据的探针ID转化为Gene ID。小编以前也是学数据挖掘出身,知道这个是小伙伴们做GEO数据挖掘的第一道坎,今天小编就来写一个函数帮助小伙伴们快速的解决这个问题。
1.从GEO数据库下载表达矩阵和注释信息(以编号GSE69078为例)
GEO官网:https://www.ncbi.nlm.nih.gov/geo/
2.用R语言获取样本临床信息,并将探针ID转化为转化为Gene symbol
代码语言:javascript复制
# 设置镜像
options(BioC_mirror="https://mirrors.ustc.edu.cn/bioc/")
# 下载R包
if (!require("readr", quietly = TRUE))
install.packages("readr")
if (!require("dplyr", quietly = TRUE))
install.packages("dplyr")
if (!require("tidyr", quietly = TRUE))
install.packages("tidyr")
if (!require("BiocManager", quietly = TRUE))
install.packages("BiocManager")
if (!require("GEOquery", quietly = TRUE))
BiocManager::install("GEOquery")
# 加载R包
library(GEOquery)
# 读取表达矩阵压缩文件
gset = getGEO(filename = './GSE69078_series_matrix.txt.gz', destdir=".",
getGPL = F)
# 获取临床信息
pd_GSE_data = pData(gset)
# 写出样本临床信息
library(readr)
write_tsv(x = data.frame(GSE_id = rownames(pd_GSE_data),pd_GSE_data),file = './GSE69078_pb_info.txt',col_names = T)
# 获取表达矩阵
GSE_data_expr = as.data.frame(exprs(gset))
GSE_data_expr = data.frame(probe_id = rownames(GSE_data_expr),GSE_data_expr)
write_tsv(x = GSE_data_expr,file = './GSE69078_probe_expression.txt',col_names = T)
# 读取探针注释信息
GSE_gpl = read_tsv(file = './GPL570-55999.txt',show_col_types = F,comment = '#')
# 获取探针对应的gene symbol,不同的芯片平台Gene symbol所在的列可能略有不同,大家先看看Gene symbol在那一列,然后在选取探针ID和gene Symbol
GSE_gpl = GSE_gpl[,c(1,11)]
# 去除一个探针对应多个symbol,不同的芯片平台,多个基因分隔符可能不一样,大家根据实际情况修改分隔符
GSE_gpl = GSE_gpl[-which(grepl('///',GSE_gpl$`Gene Symbol`)),]
# 将探针表达矩阵转化为symbol表达矩阵,对于一个symbol对应多个探针,我们去这几个探针的均值
# 将探针表达矩阵转化为gene symbol表达矩阵
probe_annotation = function(matrix,annotate,mathod = c('sum','mean','median','min','max')[2]){
# matrix是一个表达矩阵,第一列为探针ID,其他列为每个探针ID对应样本的表达值
# annotate是探针注释信息,包含两列吗,第一列为探针ID,第二列为探针ID的注释信息
# mathod多个探针ID对应同一个symbol的处理方法,默认为均值
library(dplyr)
library(tidyr)
names(matrix)[1] = 'probe_id'
names(annotate) = c('probe_id','symbol')
if(length(matrix$probe_id)==length(unique(matrix$probe_id))){
if(length(annotate$probe_id)==length(unique(annotate$probe_id))){
dat = merge(x = annotate,y = matrix) %>% dplyr::select(-probe_id)
if(mathod == 'sum'){
dat = aggregate(dat[,-1], by=list(symbol=dat[,1]),sum)
}else if(mathod == 'mean'){
dat = aggregate(dat[,-1], by=list(symbol=dat[,1]),mean)
}else if(mathod == 'min'){
dat = aggregate(dat[,-1], by=list(symbol=dat[,1]),min)
}else if(mathod == 'median'){
dat = aggregate(dat[,-1], by=list(symbol=dat[,1]),median)
}else {
dat = aggregate(dat[,-1], by=list(symbol=dat[,1]),max)
}
dat = as.data.frame(dat)
return(dat)
}else {
print('输入的探针注释的probe ID有重复,请重新输入去重之后的探针注释文件')
}
}else {
print('输入的探针表达矩阵中的probe ID有重复,请重新输入去重之后的探针表达矩阵')
}
}
symbol_exp = probe_annotation(matrix = GSE_data_expr,annotate = GSE_gpl,mathod = 'mean')
write_tsv(x = symbol_exp,file = './GSE69078_symbol_expression.txt',col_names = T)