GEO数据挖掘—处理Affymetrix的表达量芯片的cel文件

rm(list = ls())
#打破下载时间的限制,改前60秒，改后10w秒
options(timeout = 100000) 
#options(scipen = 20)#不要以科学计数法表示

#传统下载方式
library(GEOquery)
eSet = getGEO("GSE28413", destdir = '.', getGPL = F)
#网速太慢，下不下来怎么办
#1.从网页上下载/发链接让别人帮忙下，放在工作目录里
#2.试试geoChina,只能下载2019年前的表达芯片数据
#library(AnnoProbe)
#eSet = geoChina("GSE7305") #选择性代替第8行
#研究一下这个eSet
class(eSet)
length(eSet)
eSet = eSet[[1]] 
class(eSet)
#(1)提取表达矩阵exp
exp <- exprs(eSet)
#⭐第一个要检查的地方👇，表达矩阵行列数，正常是几万行，列数=样本数，
#如果0行说明不是表达芯片或者是遇到特殊情况，不能用此流程分析
dim(exp)
#二个要检查的地方👇
range(exp)#看数据范围决定是否需要log，是否有负值，异常值，如有负值，结合箱线图进一步判断
[1]  5.000 13.449
#可能要修改的地方👇
#exp = log2(exp+1) #需要log才log，不需要log要注释掉这一句
#第三个要检查的地方👇
boxplot(exp,las = 2) #看是否有异常样本

> summary(as.vector(exp))
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  5.000   5.000   5.000   5.776   5.987  13.449 

library(oligo)
celFiles<-list.celfiles('GSE28413_RAW/',listGzipped=T,
                          full.name=TRUE)
celFiles
options(BioC_mirror="https://mirrors.westlake.edu.cn/bioconductor")
exon_data<-oligo::read.celfiles(celFiles)
dim(exprs(exon_data))
exprs(exon_data)[1:4,1:4]
###标准化(一步完成背景校正、分位数校正标准化和RMA(robustmultichipaverage)表达整合)
exon_data_rma<-oligo::rma(exon_data)
exp_probe<-Biobase::exprs(exon_data_rma)
exp_probe[1:4,1:4]

dat=exp_probe
dim(dat)#看一下dat这个矩阵的维度
dat[1:4,1:4]#查看dat这个矩阵的1至4行和1至4列，逗号前为行，逗号后为列

exp <- dat
range(exp)
boxplot(exp,las = 2)

summary(as.vector(exp))

> summary(as.vector(exp))
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
 0.9982  2.6379  3.9583  4.5262  5.9866 13.4490