本文采用GSE7014数据集作为示例数据集
删除了35-44样本，仅保留26个样本进行计算
分为T2DM和Control组

首先下载CEl文件和准备Targets文件

截屏2021-06-03 下午9.16.06.png

第一步提取

### 读入Targets
setwd("/Users/apple/Desktop/练手T2DM/T2DM/GSE7014/")
library(limma)
Targets <- readTargets("Targets_GSE7104.txt")  # 不指定行名所在的列

### 读入CEL，使用readaffy
library(affy)
GSE7014_cel <- ReadAffy(filenames = Targets$FileName,
                         celfile.path = "GSE7014_RAW/",
                         phenoData = Targets)
cel <- GSE7014_cel

# 表达矩阵 exprs查看
head(exprs(cel))[, 1:5]
# 样本信息 pData
GSE7014_targets <- pData(cel)
# 样本名称 sampeNames
sampleNames(cel)
# ScanDate 查看样本扫描日期
cel@protocolData@data$ScanDate
# expression matrix 表达矩阵
GSE7014_expr_cel <- exprs(cel)

### 保存数据，备下一步分析用
save(GSE7014_cel, GSE7014_expr_cel, GSE7014_targets, file = "./GSE7014_cel.Rdata")

第二步

setwd("/Users/apple/Desktop/练手T2DM/T2DM/GSE7014/芯片质量评估/")
# 加载Affymetrix芯片GSE7014数据(第一步变量)
load_input <- load("/Users/apple/Desktop/练手T2DM/T2DM/GSE7014/GSE7014_cel.Rdata")
load_input

library(affy)
GSE7014_expr_cel <- exprs(GSE7014_cel)#提取表达矩阵
head(GSE7014_expr_cel)[, 1:5]

GSE7014_targets <- pData(GSE7014_cel) #提取targets（样本分组）数据，在phenoData中的data

第三步芯片数据质量评估 （比较简单，得到一份报告）

# 基于arrayQulitMetrics
if (!requireNamespace("BiocManager", quietly = TRUE))
  install.packages("BiocManager")

BiocManager::install("arrayQualityMetrics")

library(arrayQualityMetrics)
arrayQualityMetrics(GSE7014_cel, 
                    outdir = "./芯片数据质量评估 基于arrayQulitMetrics/",    # 通过该文件夹下的index来看结果
                    force = TRUE,   # 如果本身有个文件夹的话就把文件夹刷新了
                    intgroup = "group",    # 直接读取原始数据GSE7014_cel@phenoData@data[["group"]]
                    do.logtransform = TRUE)
dev.off()   # dev.off报错了也无妨，就是为了让它报错

第四步芯片数据标准化（基于RMA和gcRMA两种方法）

RMA

# RMA
bgcorrect.methods()  # 可以查看背景校正方法
GSE7014_bgc <- affy::bg.correct(GSE7014_cel, method = "rma")  # 最常用rma方法

# 背景矫正后作图观察样本
# boxplot
par(mfrow = c(1, 1))
par(mar = c(9, 5, 3, 3))
par(cex = 0.8)
boxplot(GSE7014_cel,
        las = 2, 
        col = rep(c("blue", "red"), each = 20),  # 20个实验组
        ylab = "Log intensity")

# 绘制PCA
# PCA_new函数
PCA_new <- function(expr, ntop = 500, group, show_name = F){
  library(ggplot2)
  library(ggrepel)
  rv <- genefilter::rowVars(expr)
  select <- order(rv, decreasing = TRUE)[seq_len(min(ntop, length(rv)))]
  # 主成分计算函数
  pca <- prcomp(t(expr[select, ])) # 最核心的函数代码
  percentVar <- pca$sdev^2/sum(pca$sdev^2)
  d <- data.frame(PC1 = pca$x[, 1], 
                  PC2 = pca$x[, 2], 
                  group = group, 
                  name = colnames(expr))
  attr(d, "percentVar") <- percentVar[1:2]
  if (show_name) {
    # 作图函数ggplot
    ggplot(data = d, aes_string(x = "PC1", y = "PC2", color = "group")) + 
      geom_point(size = 2) +
      xlab(paste0("PC1: ", round(percentVar[1] * 100), "% variance")) + 
      ylab(paste0("PC2: ", round(percentVar[2] * 100), "% variance")) +
      geom_text_repel(aes(label = name),
                      size = 3,
                      segment.color = "black",
                      show.legend = FALSE )
  } else {
    ggplot(data = d, aes_string(x = "PC1", y = "PC2",color = "group")) + 
      geom_point(size = 2) +
      xlab(paste0("PC1: ", round(percentVar[1] * 100), "% variance")) + 
      ylab(paste0("PC2: ", round(percentVar[2] * 100), "% variance"))
  }
}


PCA_new(log2(exprs(GSE7014_bgc)), 
        nrow(GSE7014_bgc), 
        group = GSE7014_targets$group,
        show_name = T)
#可以看到PCA图还是没什么改善

gcRMA

# gcRMA
library(gcrma)
GSE7014_gcrma <- gcrma(GSE7014_cel)
GSE7014_expr_gcrma <- exprs(GSE7014_gcrma) # 提取原始数据中的表达矩阵


# 运行以后再做箱线图观察
dev.new()
par(mfrow = c(1, 1))
par(mar = c(9, 5, 3, 3))
par(cex = 0.8)
boxplot(GSE7014_gcrma, # 表达矩阵
        las = 2, 
        col = rep(c("blue", "red"), each = 20),  # 20个实验组
        ylab = "Log intensity")
# dev.off()
# 可以看到GSM161970这个样本离群值非常明显
 

PCA_new(GSE7014_expr_gcrma, 
        nrow(GSE7014_gcrma), 
        group = GSE7014_targets$group,
        show_name = T)
# 可以看到970样本多重分析下来，可以考虑删除这个样本

##之后也可以再重复一次上一步arrayQualityMetrics进行质量评估，验证这一步

第五步缺失值的补充（采用KNN法），选取gcRMA标准化进行（rma也可）

#首先看一下有没有缺失值
sum(is.na(exprs(GSE7014_gcrma)))  # 可以使用RMA函数的结果也可以用gcRMA的结果
# 发现没有缺失值

# 如果有缺失值那么运行下面代码
if (!requireNamespace("BiocManager", quietly = TRUE))
  install.packages("BiocManager")
BiocManager::install("impute")
library(impute)
# KNN法impute函数计算
GSE7014_gcrma <- impute.knn(GSE7014_gcrma, maxp = 30000)  # maxp默认的是1500，好像会经常报错，可以尝试调大到30000
sum(is.na(GSE7014_gcrma)) # 再次查看是否有缺失值

第六步标准化结果可视化

# 再次用boxplot和PCA进行可视化查看

# 运行以后再做箱线图观察
dev.new()
par(mfrow = c(1, 1))
par(mar = c(9, 5, 3, 3))
par(cex = 0.8)
boxplot(GSE7014_gcrma, # 表达矩阵
        las = 2, 
        col = rep(c("blue", "red"), each = 20),  # 20个实验组
        ylab = "Log intensity")
# dev.off()

PCA_new(GSE7014_expr_gcrma, 
        nrow(GSE7014_gcrma), 
        group = GSE7014_targets$group,
        show_name = T)

## 保存变量后续使用
save(GSE7014_expr_gcrma, GSE7014_targets, 
     file = "affymetrix/GSE29450/GSE29450_after_gcrma.Rdata")