oligo包处理原始芯片数据

# 1.载入包

rm(list=ls()) #   清空变量
library(oligo)
library(GEOquery)

# if (!requireNamespace("BiocManager", quietly = TRUE))
#   install.packages("BiocManager")
# 
# BiocManager::install("oligo")
# BiocManager::install("GEOquery")

# 2.下载原始文件

# 以GSE156199为例
gse_id <- 'GSE156199'
setwd("/Users/zhengxueming/scripts/R_scripts/test/")
baseDir <- "/Users/zhengxueming/scripts/R_scripts/test/"
getGEOSuppFiles(gse_id, makeDirectory = TRUE, baseDir = baseDir,
                fetch_files = TRUE, filter_regex = NULL)

如果是压缩文件，需要先解压到当前文件夹。

# 3. 得到原始数据

workDir <- file.path(baseDir, gse_id)
celfiles <- list.files(workDir, "\\.gz$")  # 匹配以.gz 结尾的文件。
data.raw <- read.celfiles(filenames = file.path(workDir, celfiles))
# read.xysfiles (NimbleGen芯片)
class(data.raw) # "GeneFeatureSet"

# 4. 芯片质量评估

fit <- fitProbeLevelModel(data.raw)
# 按样本名画图，可以给不同分组标记不同颜色
groups <- c("group1","group1","group2","group2")
boxplot(fit, names = NA, col = as.factor(groups))
legend("topright", legend = unique(groups), fill = as.factor(unique(groups)),
       box.col = NA, bg = "white", inset = 0.01)

par(mfrow = c(2, 2))
MAplot(data.raw, pairs = FALSE)

# 5.表达量计算

data.eset <- rma(data.raw) # rma算法，返回值经过log2变换
## Background correcting
## Normalizing
## Calculating Expression
class(data.eset)  #ExpressionSet 对象

data.exprs <- exprs(data.eset)  #得到表达矩阵
class(data.exprs)  # matrix" "array"

sample.names <- colnames(data.exprs)  # sample name
probe_ids <- rownames(data.exprs) # probe id

# 6.和GEO下载的处理后表达矩阵比较

# 下载得到表达矩阵
library(GEOquery)
gse <- getGEO('GSE156199',getGPL = FALSE) # 同时下载了GPL文件
exprSet <- exprs(gse[[1]]) # 基因表达矩阵
# 比较
dim(exprSet)
dim(data.exprs)

exprSet[1:5,1:4]
data.exprs[1:5,1:4]

summary(exprSet)
summary(data.exprs)

# sessionInfo()