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°1$logFC>0.5~"up",
deg1$P.Value<0.05°1$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包。
参考链接：
