目的是为了找出表达矩阵中有哪些基因与目的基因有相关性。
#=======================================================
#=======================================================
library(GEOquery)
rm(list=ls())
library(dplyr)
library(tidyr)
library(Biobase)
library(limma)
setwd('D:\\SCIwork\\F23\\GSE48780')
gsename = "GSE48780"
# 下载基因芯片数据,destdir参数指定下载到本地的地址
gse<- getGEO(gsename, destdir = ".")
##根据GSE号来下载数据,下载_series_matrix.txt.gz
gpl<- getGEO('GPL570', destdir = ".")
##根据GPL号下载的是芯片设计的信息, soft文件
gse <- getGEO(filename = 'GSE48780_series_matrix.txt.gz')
gpl <- getGEO(filename = 'GPL570.soft')
# 查看列名
colnames(Table(gpl))
Table(gpl)[1:10,1:6] # 前10行前6列信息
gpl <- gpl@dataTable@table
colnames(gpl)
gpl <- gpl %>%
dplyr::select(ID, "Gene Symbol")
write.csv(gpl,"GPL.csv", row.names = F)
# gse中的行名ID与gene name的对应关系
genename = read.csv("GPL.csv")
genename <- genename%>%
tidyr::separate(Gene.Symbol,
into = c('Gene', 'Symbol'),
sep='\\///')%>%
dplyr::select(ID,Gene )
##########################################################################################
##
###########################################################################################
setwd('D:\\SCIwork\\F23\\GSE48780')
# 构建表达矩阵
exprSet <- as.data.frame(exprs(gse)) # 得到表达矩阵,行名为ID,需要转换
# 转换ID为gene name
exprSet$ID = rownames(exprSet)
express = merge( x=genename, y=exprSet, by="ID")
express$ID = NULL
express[which(is.na(express),arr.ind = T)]<-0 #结合which进行缺失替代
exprSet <- aggregate(x = express[,2:ncol(express)],
by = list(express$Gene),
FUN = max)
head(exprSet)
exprSet <- as.data.frame(exprSet)
exprSet <-exprSet[-1,]
names(exprSet)[1] <- 'ID'
rownames(exprSet) <- exprSet$ID
exprSet$ID <- NULL
write.csv(exprSet, file = 'exprSet.csv')
save(exprSet, file = 'exprSet.Rdata')