一文解决大批量基因相关性分析

#=======================================================


#=======================================================

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')