2022-05-06

# 首先加载依赖的R包
library(annaffy)
library(affy)
library(annotate)
library(GEOquery)
library(hgu133plus2.db)
library(org.Hs.eg.db)

# 设置RAW文件所在文件夹的路径
setwd("~\\Desktop\\data")
path <- getwd()
# gse给出所有的待处理GEO文件名称
gse <- c("GSE1234", "GSE345", "GSE4456")
for (n in gse) {
  setwd(path)
  # dir.create函数创建与GSE文件同名的文件夹，以便压缩包的解压和文件提取
  dir.create(n)
  # 设置默认路径
  setwd(paste(path, n, sep = "/"))
  # 对压缩文件解压
  untar(paste(paste(path, n, sep = "/"), "RAW.tar", sep = "_"))
  # 提取所有的命名中带有gz的文件
  files <- dir(pattern="gz$") ##加载文件
  # 将files中的文件进行合并
  sapply(files, gunzip) ##合并文件
  # 将所有名字带有CEL的文件存放到filelist变量中，每个CEL文件就是一个样本的芯片数据
  filelist <- dir(pattern="CEL$")
  # 使用ReadAffy函数读取每个CEL文件，并对芯片数据进行初步处理
  data <- ReadAffy(filenames=filelist)
  # 使用rma函数对芯片数据进行标化
  eset <- rma(data)
  # expres函数用于提取eset对象中的表达数据
  eset.e <- data.frame(exprs(eset))
  # 提取表达数据的基因信息
  affydb <- annPkgName(data@annotation,type="db")
  require(affydb, character.only=TRUE)
  genes <- as.character(aafSymbol(as.character(rownames(eset)),affydb))
  eset.e$genes <- genes
  # 至此实现了对GEO原始数据的RMA标准化处理，处理后的表达谱数据为eset.e
  write.csv(eset.e, "rma_exp.csv", row.names = F) # 导出RMA标化后的数据

  # 有些研究发现可以将RMA标化方法与MAS5标化方法进行结合，具体的实现方式如下

  # 使用mas5calls函数对data进行处理，得到mas5标化后的标化数据
  calls <- mas5calls(data) # get PMA calls
  # exprs函数用于提取表达谱
  calls <- exprs(calls)
  # 对数据calls进行筛选，最终得到数据即为在RMA标化的基础上，进一步考虑了mas5
  absent <- rowSums(calls == 'A') # how may samples are each gene 'absent' in all
  absent <- which (absent == ncol(calls)) # which genes are 'absent' in all samples
  rmaFiltered <- eset[-absent,] # filters out the genes 'absent' in all samples 
  filtered <- data.frame(exprs(rmaFiltered))
  symbols<-as.character(aafSymbol(as.character(rownames(rmaFiltered)),affydb))
  filtered$genes <- symbols
  write.csv(filtered, "rma_mas5_exp.csv", row.names = F)
}
rm(list = ls())