Affymetrix芯片储存着大量的生物信息学数据，因此有必要从实战出发的角度，汇总下Affymetrix芯片处理的流程。下面以GSE1438为例

1.数据下载及读取

rm(list = ls())
library(GEOquery)
library(affyPLM)
library(affy)
gse="GSE1428"
#rawdata <- getGEOSuppFiles(gse,fetch_files = F) #下载原始数据
gset <- getGEO(gse,getGPL = F,destdir = ".",)
#setwd(gse)
#filenames <- list.files(pattern = "*.gz$") #批量查找并列出后缀为.gz的文件
#sapply(filenames, gunzip) #批量解压

data.raw <- ReadAffy() #去读当前目录下所有的CEL文件
sampleNames(data.raw)

##  [1] "GSM23735.CEL" "GSM23736.CEL" "GSM23737.CEL" "GSM23738.CEL" "GSM23740.CEL"
##  [6] "GSM23741.CEL" "GSM23742.CEL" "GSM23743.CEL" "GSM23744.CEL" "GSM23745.CEL"
## [11] "GSM23746.CEL" "GSM23747.CEL" "GSM23748.CEL" "GSM23749.CEL" "GSM23750.CEL"
## [16] "GSM23751.CEL" "GSM23752.CEL" "GSM23753.CEL" "GSM23754.CEL" "GSM23755.CEL"
## [21] "GSM23756.CEL" "GSM23757.CEL"

sampleNames(data.raw) <- sub(pattern = "\\.CEL",replacement = "",sampleNames(data.raw))

2.数据预处理

常用的质量控制的指标：平均数法、RLE、NUSE和RNA降解曲线 根据以上指标综合决定实验是否合格，并提出质量不合格的样品。

library(simpleaffy)
data.qc <- qc(data.raw) #平均值法
plot(data.qc)

simpleaffy-qc

Pset <- fitPLM(data.raw)
library(RColorBrewer)
colors <- brewer.pal(12,"Set3")
Mbox(Pset,ylim=c(-1,1),main="RLE",las=3,col=colors)

RLE

boxplot(Pset,ylim=c(0.95,1.22),col=colors,mian="NUSE",las=3)

NUSE

data.deg <- AffyRNAdeg(data.raw)
plotAffyRNAdeg(data.deg,col=colors)
legend("topleft",rownames(pData(data.raw)),col = colors,lwd=1,inset = 0.05,cex = 0.5)

RNA降解曲线

可以看出，这个芯片的整体检查率并不太高，且GSE23740、GSM23745、GSM23746、GSM23750、GSM2375和GSM23757的RLE和NUSE偏离中心太多，整体RNA降解斜率偏低。在实际科研中，我们最好寻找高质量的芯片。

考虑到整体芯片质量不佳，过滤后剩余的样本数会比较少，下面就假装质量还可以进行下游分析（请大家谅解！）

3.RMA标准化

data.rma <- rma(data.raw)

## Background correcting
## Normalizing
## Calculating Expression

data.eset <- exprs(data.rma) #表达矩阵

#整理分组信息
pd <- pData(gset[[1]])
pd <- pd %>% 
  select(geo_accession,description)
colnames(pd) <- c("id","type")
pd$type <- sapply(strsplit(pd$type,"(",fixed = T),"[",1)
pd$type <- trimws(pd$type)
pd$type <- factor(pd$type,levels = c("Young","Older"))
pd <- pd[order(pd$type),]
data.eset <- data.eset[,pd$id]
all(colnames(data.eset)==pd$id)

## [1] TRUE

boxplot(data.eset,col=colors,las=3,mian="after-RMA") #经过RMA标准化后

RMA标准后箱式图

4.差异分析

group_list <- factor(pd$type,levels = c("Young","Older"))
names(group_list) <- pd$id
design <- model.matrix(~group_list)

library(limma)
fit <- lmFit(data.eset,design)
fit1 <- eBayes(fit)
options(options=4)
deg <- topTable(fit1,coef=2,n=Inf) %>% 
  rownames_to_column("probe_id")

5.探针ID注释

gpl <- gset[[1]]@annotation
library(AnnoProbe)
probe2id <- idmap(gpl)
deg1 <- probe2id %>% 
  inner_join(deg,by="probe_id") %>% 
  select(-probe_id) %>% 
  arrange(desc(logFC)) %>% 
  distinct(symbol,.keep_all=T)
deg1$type <- case_when(deg1$P.Value<0.05&deg1$logFC>0.5~"up",
                       deg1$P.Value<0.05&deg1$logFC< -0.5~"down",
                       T~"stable")

6.火山图

ggplot(deg1,aes(x=logFC,y=-log10(P.Value),color=type))+
  geom_point(alpha=0.4,size=3.5)+
  scale_color_manual(values=c("#546de5", "#d2dae2","#ff4757"))+
  geom_vline(xintercept=c(-1,1),lty=4,col="black",lwd=0.8) +
  geom_hline(yintercept = -log10(0.05),lty=4,col="black",lwd=0.8) +
  labs(x="log2(fold change)",
       y="-log10 (p-value)")+
  theme_bw()+
  theme(plot.title = element_text(hjust = 0.5), 
        legend.position="right", 
        legend.title = element_blank())

火山图

当然affy包主要针对的是旧版的Affymetrix芯片，如hgu95/95和hgu133系列。下一篇我们来看看oligo包。

参考链接：

R语言_Affymetrix芯片数据处理

用affy包读取affymetix的基因表达芯片数据-CEL格式数据