GEO代码

rm(list = ls())#清空环境变量
options(stringsAsFactors = F)##字符不作为因子读入

####getGPL获得平台的注释信息，但下载速度会慢很多
####而且注释文件格式大多不如bioconductor包好用
#目前暂且把分析一个GSE数据集称为一个项目
#在分析一个项目时总不可能一次性做完，更多情况是经常会中途暂停，所以如果下次再打开该project时，为避免再次下载GSE数据集，所以用了以下代码
library(GEOquery)
g <- 'GSE42872'
G <- paste(g,"Rdata",sep = ".")
if(!file.exists(G)){#如果该GSE文件不存在，就下载，并保存
  gset<-getGEO(g,destdir='.',
               AnnotGPL=F,
               getGPL=F)
  save(gset,file=G)
} else load(G)#如果该文件存在，就加载
#这里认识一下对象Large ExpressionSet，即gset[[1]]
#gset是包含了一个Large ExpressionSet对象的列表，而同时Large ExpressionSet对象也是包含很多数据的列表
#所以gset是一个包含了一个列表的列表
lareset <- gset[[1]]#获取对象Large ExpressionSet
lareset@experimentData
#了解该实验数据类型（芯片或测序）、实验平台、样本名称
#了解实验背景知识（abstract(lareset)）、实验设计
annotation(lareset)#可获得实验平台
lareset@annotation#可获得实验平台
lareset@assayData#含有表达矩阵
lareset@assayData[["exprs"]]
exprs(lareset)#两种方式提取表达矩阵
lareset@phenoData#含有分组信息
pData(lareset)#直接提取分组信息
#lareset$XXX等价于lareset[[]]
express<- exprs(gset[[1]])#获取表达矩阵
pdata <- pData(gset[[1]])#获取分组信息

ex <- express
qx <- as.numeric(quantile(ex, c(0., 0.25, 0.5, 0.75, 0.99, 1.0), na.rm=T))
LogC <- (qx[5] > 100) ||
  (qx[6]-qx[1] > 50 && qx[2] > 0) ||
  (qx[2] > 0 && qx[2] < 1 && qx[4] > 1 && qx[4] < 2)
#LogC为FALSE时，代表表达矩阵不需要log2处理
#LogC为TRUE时，代表表达矩阵需要log2处理
#当数据需要log2处理时，因为数据中某个样本的某个基因表达量可能为0，那么log2将会得出负无穷大，针对这种情况有以下两种处理
#第一种：把表达量为0的赋值成NaN，经过log2处理，该值仍为NaN
if (LogC) { ex[which(ex <= 0)] <- NaN
express_log2 <- log2(ex)
print("log2 transform finished")}else{print("log2 transform not needed")}
#第二种：把表达量加1
if (LogC) { 
express_log2 <- log2(ex+1)
print("log2 transform finished")
}
else{
print("log2 transform not needed")
}
#根据pdata分组信息进行分组
#归一化处理
#