Q：Microarry数据标准化流程？
本文为记录阅读王诗翔老师博客《Bioconductor分析基因芯片数据》的读书笔记。
早期的生物信息学数据主要是由芯片（microarray）产生的，尽管其逐渐被测序所取代， 但是在价格方面还是有非常大的优势，目前仍然是许多科研项目的主力军。聊起科研就不得不提 芯片制造商的三大巨头：Affymetrix、Agilent和Illumina。
最初的芯片主要是为了检测mRNA的表达谱而开发的，在当时几乎就是Affymetrix的天下。所以 目前有许多数据是Affymentrix时代的产物，因此研究其表达谱芯片的分析和处理非常有必要。
用于研究基因谱的芯片有两类：cDNA芯片和寡核苷酸芯片。
Affymetrix表达谱芯片大致可分为旧版和新版两类，旧版一般指Human Genome U133(hgu)系列， 更新版如人类的HST2.0，常的处理R包：affy、simpleaffy和oligo。
1.快速入门
affy包是Bioconductor平台开发的，主要用来读取Affymetrix芯片CEL原始文件，预处理，获得表达矩阵。
学习affy前首先要明确下列术语：
-
探针（probe）：包含25个碱基的寡核苷酸；
-
由于探针长度较短，为保证杂交的特异性，affy公司设计了两种探针，一类完全匹配（perfect match），另一类不匹配（mismatch），其信号值称为PM和MM。二者成对出现， 且有共同的affyID。 其实PM和MM除第13个碱基外，其余完全一致，只是在MM中把PM的第13个碱基替换成了互补碱基而已。
-
CEL文件：储存每个探针的信号值和定位信息
-
CDF文件：储存每个探针对的定位信息
下面以CCL包为例，完成预处理，最终获得基因表达矩阵。
rm(list = ls())
library(pacman)
p_load(CLL) #载入CLL包时，会自动加载affy包
data("CLLbatch") #CLL自带示例CEL格式数据
class(CLLbatch) #数据格式为affyBatch
## [1] "AffyBatch"
## attr(,"package")
## [1] "affy"
CLLrma <- rma(CLLbatch) # rma预处理
## Warning: replacing previous import 'AnnotationDbi::tail' by 'utils::tail' when
## loading 'hgu95av2cdf'
## Warning: replacing previous import 'AnnotationDbi::head' by 'utils::head' when
## loading 'hgu95av2cdf'
##
## Background correcting
## Normalizing
## Calculating Expression
eset <- exprs(CLLrma) #取出表达矩阵
eset[1:6,1:6]
## CLL10.CEL CLL11.CEL CLL12.CEL CLL13.CEL CLL14.CEL CLL15.CEL
## 100_g_at 7.495697 7.945159 7.861043 7.990252 7.889718 7.788036
## 1000_at 7.251120 8.298723 8.473531 8.131041 8.051400 7.874776
## 1001_at 4.457362 4.517529 4.357791 4.651236 4.473598 4.767024
## 1002_f_at 3.984693 3.981675 4.065171 4.132480 4.065022 3.821567
## 1003_s_at 6.437046 6.200826 6.412218 6.313572 6.105120 6.085442
## 1004_at 7.574146 8.054949 8.276920 8.105801 7.990090 7.596703
2.数据预处理
预处理的目的是将探针水平的数据转换成基因表达数据，数据结构包含AffyBatch类和ExpressionSet类：
-
AffyBatch →存储探针水平数据（相当于读取CEL文件）；
-
ExpressionSet →存储表达水平数据（基因表达矩阵）。
预处理通过质量控制（剔除不合格的芯片数据）→ 标准化（将所有芯片中的基因表达值转化到一个可以比较的水平）→ 后续分析
2.1 示例数据
rm(list = ls())
library(pacman)
p_load(CLL)
data("CLLbatch")
data.class(CLLbatch)
## [1] "AffyBatch"
data("disease") #CLL对应的临床信息
head(disease)
## SampleID Disease
## 1 CLL10 <NA>
## 2 CLL11 progres.
## 3 CLL12 stable
## 4 CLL13 progres.
## 5 CLL14 progres.
## 6 CLL15 progres.
phenoData(CLLbatch) #表型信息
## An object of class 'AnnotatedDataFrame'
## sampleNames: CLL10.CEL CLL11.CEL ... CLL9.CEL (24 total)
## varLabels: sample
## varMetadata: labelDescription
featureData(CLLbatch)
## An object of class 'AnnotatedDataFrame': none
annotation(CLLbatch) #注释信息
## [1] "hgu95av2"
protocolData(CLLbatch)
## An object of class 'AnnotatedDataFrame': none
exprs_matrix <- assayData(CLLbatch)[[1]] #表达矩阵
exprs_matrix[1:6,1:6]
## CLL10.CEL CLL11.CEL CLL12.CEL CLL13.CEL CLL14.CEL CLL15.CEL
## 1 183.0 113.0 119.0 130.0 133.0 141.0
## 2 11524.8 6879.8 7891.3 8627.5 8205.3 5354.3
## 3 301.0 146.0 133.0 160.0 153.0 130.0
## 4 11317.0 6747.0 7916.0 8616.0 7865.0 5382.0
## 5 115.0 82.0 78.8 101.0 65.3 69.3
## 6 139.0 107.0 121.0 145.0 126.0 98.0
exprs(CLLbatch)[1:6,1:6]
## CLL10.CEL CLL11.CEL CLL12.CEL CLL13.CEL CLL14.CEL CLL15.CEL
## 1 183.0 113.0 119.0 130.0 133.0 141.0
## 2 11524.8 6879.8 7891.3 8627.5 8205.3 5354.3
## 3 301.0 146.0 133.0 160.0 153.0 130.0
## 4 11317.0 6747.0 7916.0 8616.0 7865.0 5382.0
## 5 115.0 82.0 78.8 101.0 65.3 69.3
## 6 139.0 107.0 121.0 145.0 126.0 98.0
可以看出assayData和exprs方法调取的结果是一致的。
-
assayData槽是AffyBatch类必不可少的，其返回结果的第一个元素是矩阵类型，用于保存 基因表达矩阵。行对应不同的探针组，用无重复的索引值表示；列对应样本。用exprs()调取的 正是此基因表达矩阵。
2.2 质量控制
质量控制主要集中在CEL文件的处理，可以分为3个不同层次：
-
直接观察：image()函数
-
平均值法：simpleaffy包
-
数据拟合法：affyPLM包
image():直接看芯片上所有位点的灰度值
image(CLLbatch[,1])
image.png
这种方法只是对芯片信号强度有一个总体认识：图像特别黑，信号强度低；图像特别亮，信号 可能过饱和。很明显，这种方法主观性太强，不能量化，不建议用。比较简单的方法是基于各种 平均值。
qc():平均值法
类似平均值法有一个共同特点：
假设一组实验中每个芯片数据对于某个平均值指标都相差不大。
p_load(simpleaffy)
library(CLL)
data("CLLbatch")
data.qc <- qc(CLLbatch) #质量控制
plot(data.qc)
image.png
上图显示的是全部质量控制总览图： * 第1列是样本名称； * 第2列是两个数字，上面是检出率，下面是平均背景噪声； * 第3列（“QC Stats”）:最下面的横轴是尺度因子等对应的坐标，-3~3。共用到3项指标，尺度因子、 GAPDH 3’/5’比值和actin 3’/5’比值。出现红色的“bioB”，说明该样品中未检测到BioB。
简单来说，所有指标出现蓝色正常，红色可能存在质量问题。一般来说，一个芯片的各项指标都不太正常，尤其是BioB无法检测到，建议判定该芯片失败。比如CLL15.CEL数据检出率明显低于其他样品，actin3/actin5远大于3，且未检测到BioB，可判定其为无效。
affyPLM：数据拟合
library(CLL)
library(affyPLM)
## Loading required package: preprocessCore
data("CLLbatch")
Pset <- fitPLM(CLLbatch) #对数据集做回归分析，结果为PLMset类
image(CLLbatch[,1])
image.png
image(Pset,type="weights",which=1,main="Weights") #根据计算结果画权重图
image.png
image(Pset,type="resids",which=1,main="Residuals") #残差图
image.png
image(Pset, type="sign.resids",which=1, main="Residuals.sign")
image.png
affyPLM包数据拟合时引入了加权最小二乘法来进行回归，理论上，权重和残差的分布是随机的， 应该看到绿色均匀分布的权重图和红蓝均匀分布的残差图。
相对对数表达箱线图（RLE）反应了基因表达量的一致性趋势，定义为一个探针组在某个样品的表达值除以该探针组在所有样品中表达值的中位数后取对数。(即log10(x/median(x)))
如果用RLE来控制数据质量，每个样品的中心应该非常接近纵坐标0的位置。
library(RColorBrewer)
library(affyPLM)
colors <- brewer.pal(12,"Set3")
Mbox(Pset,ylim=c(-1,1),col=colors,main="RLE",las=3)
image.png
NUSE是一种比RLE更为灵敏的质量检测手段。NUSE定义是探针组在某个样品的PM值的标准差除以该探针值在 各样品中的PM标准差的中位值。(即sd(PM)/median(sd(PM)))，如果芯片质量都非常可靠，他们的sd应该非常接近，NUSE值都在1附近。
boxplot(Pset,ylim=c(0.95,1.22),col=colors,main="NUSE",las=3)
image.png
可以看出CLL1和CLL10质量有显著问题，需要舍弃。
RNA降解是影响芯片质量的很重要因素。因为RNA是从5‘端开始降解的，理论上5’端的荧光强度应该低于 3‘端的荧光强度。如果RNA降解曲线的斜率太小，甚至接近0，很可能RNA降解太严重，需要作为坏数去除（下图CLL13）。
library(affy)
data.deg <- AffyRNAdeg(CLLbatch)
plotAffyRNAdeg(data.deg,col=colors)
legend("topleft",rownames(pData(CLLbatch)),col = colors,lwd=1,inset = 0.05,cex = 0.5)
image.png
CLLbatch <- CLLbatch[,-match(c("CLL10.CEL","CLL1.CEL","CLL13.CEL"),sampleNames(CLLbatch))]
去掉3个质量差的。
RLE、NUSE都是基于平均值的思想，此外还可以从芯片间的相互关系来检测芯片的质量（物以类聚，人以群分）。即组内芯片应该聚在一起，组间明显分开。Pearson线性相关系数就是常用指标，但是在实际科研中，往往不是直接查看相关系数矩阵，而是根据相关系数导出的距离矩阵，进行聚类分析或主成分分析对样品归类并可视化。
library(pacman)
p_load(gcrma)
p_load(graph)
p_load(affycoretools)
p_load(CLL)
data(disease)
CLLgcrma <- gcrma(CLLbatch) #用gcrma法预处理
## Adjusting for optical effect.....................Done.
## Computing affinities
## Loading required package: AnnotationDbi
## Loading required package: stats4
## Loading required package: IRanges
## Loading required package: S4Vectors
##
## Attaching package: 'S4Vectors'
## The following objects are masked from 'package:dplyr':
##
## first, rename
## The following object is masked from 'package:tidyr':
##
## expand
## The following object is masked from 'package:base':
##
## expand.grid
##
## Attaching package: 'IRanges'
## The following object is masked from 'package:simpleaffy':
##
## members
## The following object is masked from 'package:grDevices':
##
## windows
## The following objects are masked from 'package:dplyr':
##
## collapse, desc, slice
## The following object is masked from 'package:purrr':
##
## reduce
##
## Attaching package: 'AnnotationDbi'
## The following object is masked from 'package:dplyr':
##
## select
## .Done.
## Adjusting for non-specific binding.....................Done.
## Normalizing
## Calculating Expression
eset <- exprs(CLLgcrma) #提取表达矩阵
pearson_cor <- cor(eset)
dist.lower <- as.dist(1-pearson_cor)
hc <- hclust(dist.lower,"ave") #聚类
plot(hc)
image.png
# PCA主成分分析
colnames(eset);
## [1] "CLL11.CEL" "CLL12.CEL" "CLL14.CEL" "CLL15.CEL" "CLL16.CEL" "CLL17.CEL"
## [7] "CLL18.CEL" "CLL19.CEL" "CLL20.CEL" "CLL21.CEL" "CLL22.CEL" "CLL23.CEL"
## [13] "CLL24.CEL" "CLL2.CEL" "CLL3.CEL" "CLL4.CEL" "CLL5.CEL" "CLL6.CEL"
## [19] "CLL7.CEL" "CLL8.CEL" "CLL9.CEL"
samplenames <- sub(pattern = "\\.CEL",replacement = "",colnames(eset))
groups <- factor(disease[,2])
plotPCA(eset,addtext=samplenames,groups=groups,groupnames=levels(groups))
image.png
从聚类结果看，progress和stable组根本不能很好区分，但还不能判定实验失败。理论上讲，如果总体上两组数据是分开的，那么说明我们关心的导致癌症从progress到stable的因素起到了主要因素，如果不是，需要具体问题具体分析。CLL数据来自不同的个体，很可能个体差异起到主导作用，才会出现未达到预期聚类。
因此，只有当聚类图中有明显的类别差异时，才适合考虑去除个别不合群的样品；如果整体分类被打乱，则不能简单判定所有样品都有问题。
注意：使用PCA分析时，必须考虑前两个主成分是否具有代表性，即查看前2个主成分的累计贡献率，如果低于60%，考虑使用多维尺度分析来构建分类图。
2.3 背景校正、标准化和汇总
芯片经过质量控制，去除不合格的样品后，还需要经过北京校正、标准化和汇总3个步骤。
-
背景校正：实际科研中，经常发现大量的MM（Mismatch）信号比PM（Perfect Match）信号还有强。这时候就需要用统计模型来去除背景噪声，即背景校正（Background correction）。
-
标准化：消除非实验误差，让组间测量具有可比性。根据假设：bulk normalization，即仅有一小部分基因表达值在不同条件下有差异，因此使用所有基因表达值作为参考进行标准化。
最后用一定的统计方法，将前面获得的荧光信号值从探针水平，汇总到探针组水平，此过程成为汇总（Summariztion）。
MAS5预处理不指定任何芯片，而是设定一个假想的芯片最为参考，这样每个芯片都可以单独计算，这样有新的数据加入时可不必重新计算已经标准化过的数据。
RMA预处理默认算法时分位数方法。
下面用MAS方法进行背景校正和标准化：
CLLmas5 <- bg.correct(CLLbatch,method = "mas") #使用mas方法背景校正
data_mas5 <- normalize(CLLmas5,method="constant") #constant方法标准化
head(pm(data_mas5)/pm(data_mas5),5) #查看每个样品的缩放系数
## CLL11.CEL CLL12.CEL CLL14.CEL CLL15.CEL CLL16.CEL CLL17.CEL CLL18.CEL
## 175218 1 1 1 1 1 1 1
## 356689 1 1 1 1 1 1 1
## 227696 1 1 1 1 1 1 1
## 237919 1 1 1 1 1 1 1
## 275173 1 1 1 1 1 1 1
## CLL19.CEL CLL20.CEL CLL21.CEL CLL22.CEL CLL23.CEL CLL24.CEL CLL2.CEL
## 175218 1 1 1 1 1 1 1
## 356689 1 1 1 1 1 1 1
## 227696 1 1 1 1 1 1 1
## 237919 1 1 1 1 1 1 1
## 275173 1 1 1 1 1 1 1
## CLL3.CEL CLL4.CEL CLL5.CEL CLL6.CEL CLL7.CEL CLL8.CEL CLL9.CEL
## 175218 1 1 1 1 1 1 1
## 356689 1 1 1 1 1 1 1
## 227696 1 1 1 1 1 1 1
## 237919 1 1 1 1 1 1 1
## 275173 1 1 1 1 1 1 1
芯片内标准化方法针对双通道数据，又可分为全局化方法（Global normaliztion）和荧光强度依赖的方法（Intensity-dependent normalzation）。
全局化方法：假设红光信号强度与绿光信号强度呈正比例关系，即R=kG（R:红色信号强度；G:绿色信号强度）。
下面来比较不同算法的处理效果：
rm(list = ls())
library(gcrma)
library(affyPLM)
library(RColorBrewer)
library(CLL)
data("CLLbatch")
colors <- brewer.pal(12,"Set3")
CLLmas5 <- mas5(CLLbatch)
## background correction: mas
## PM/MM correction : mas
## expression values: mas
## background correcting...done.
## 12625 ids to be processed
## | |
## |####################|
CLLrma <- rma(CLLbatch)
## Background correcting
## Normalizing
## Calculating Expression
CLLgcrma <- gcrma(CLLbatch)
## Adjusting for optical effect........................Done.
## Computing affinities.Done.
## Adjusting for non-specific binding........................Done.
## Normalizing
## Calculating Expression
## hist
hist(CLLbatch,main="orignal",col=colors)
legend("topright",rownames(pData(CLLbatch)),col=colors,
lwd = 1,inset = 0.05,cex = 0.5,ncol = 3)
image.png
hist(CLLmas5,main="MAS 5.0",xlim=c(-150,2^10),col=colors)
image.png
hist(CLLrma,main="RMA",col=colors)
image.png
hist(CLLgcrma,main="gcRMA",col=colors)
image.png
简要说明：RMA算法将多条曲线重合在一起，有利于进行下一步差异分析，但有双峰出现，不符合高斯正态分布。gcRMA算法稍微好些，但不意味着gcRMA总是优于RMA算法。
## boxplot
boxplot(CLLbatch, col=colors, las=3, main="original")
image.png
boxplot(CLLmas5, col=colors, las=3, ylim=c(0,1000), main="MAS 5.0")
image.png
boxplot(CLLrma, col=colors, las=3, main="RMA")
image.png
boxplot(CLLgcrma, col=colors, las=3, main="gcRMA")
image.png
简要说明： 3种算法各样本中值比较接近,当MAS5、gcRMA都有一定拖尾现象，总体来说RMA效果比较令人满意。
2.4 预处理的一体化计算
在实际应用中，一般都是预处理一条龙服务，RMA和MAS5比较常用。
image.png
预处理方法
MAS5和RMA的异同点：
-
MAS5对每个芯片单独进行标准化；RMA采用多芯片模型，对所有芯片一起标准化。
-
MAS5利用MM探针去除背景噪声，而RMA不使用MM信息，而是基于PM信号分布采用随机模型估算表达量。
注意：RMA处理后的数据是经过log2为底的对数转换的，而MAS5不是。
3.完整流程
rm(list = ls())
library(limma)
library(gcrma)
library(CLL)
library(RColorBrewer)
library(affyPLM)
data("CLLbatch")
data("disease")
#RLE和NUSE质控
Pset <- fitPLM(CLLbatch) #对数据进行回归分析
colors <- brewer.pal(12,"Set3")
Mbox(Pset,ylim=c(-1,1),col=colors,main="RLE",las=3) #RLE法质量控制
boxplot(Pset,ylim=c(0.95,1.22),col=colors,main="NUSE",las=3) #NUSE质量控制
#RNA降解情况
library(affy)
data.deg <- AffyRNAdeg(CLLbatch)
plotAffyRNAdeg(data.deg,col=colors)
legend("topleft",rownames(pData(CLLbatch)),col = colors,lwd=1,inset = 0.05,cex = 0.5)
CLLbatch <- CLLbatch[,-match(c("CLL10.CEL","CLL1.CEL","CLL13.CEL"),sampleNames(CLLbatch))]
CLLrma <- rma(CLLbatch)
sampleNames(CLLrma) <- gsub(".CEL$","",sampleNames(CLLrma))
disease <- disease[match(sampleNames(CLLrma),disease[,"SampleID"]),]
eset <- exprs(CLLrma) #表达矩阵
boxplot(eset,col=colors)
#分组信息
disease <- factor(disease[,"Disease"])
design <- model.matrix(~-1+disease)
con.matrix <- makeContrasts(contrasts = "diseaseprogres.-diseasestable",
levels = design)
#差异分析
fit <- lmFit(eset,design)
fit1 <- contrasts.fit(fit,con.matrix)
fit2 <- eBayes(fit1)
deg <- topTable(fit2, coef=1,
n=nrow(fit2), lfc=log2(1.5))
diff <- deg[deg[,"P.Value"]<0.01,]
# ID注释
library(annotate)
affydb <- annPkgName(CLLbatch@annotation,type = "db")
library(affydb,character.only = T)
diff$symbols <- getSYMBOL(rownames(diff),affydb) #获取Gene Symbol
diff$ENtrezID <- getEG(rownames(diff),affydb) #获取Entrez ID
参考链接：
Bioconductor分析基因芯片数据
