用R进行microarray分析（2）

归一化探针水平的数据并计算基因表达数据

x<-expresso(a,bg.correct = TRUE,bgcorrect.method = "mas",
normalize.method = "constant",pmcorrect.method = "pmonly",
summary.method = "medianpolish")

函数 expresso 封装了计算基因表达数据过程中的归一化、汇总探针集、背景校正、PM 探针校正等 4 个过程。
本例中 bg.correct = TRUE 表示进行背景校正，采用 mas 方法,
normalize.method = "constant", 归一化采用 constant 方法;
pmcorrect.method = "pmonly",PM 校正采用 pmonly 方法;
summary.method = "medianpolish", 汇总探针集采用 medianpolish 方法。

我们看一看其他还存在什么方法：
背景校正：

bgcorrect.methods()

bgcor.png

归一化方法：

normalize.method(a)

normalize_method.png

PM校正方法：

 pmcorrect.methods()

pmcor.png

汇总探针集方法：

express.summary.stat.methods()

exp_method.png

将分析结果保存到文件中

write.exprs(x,file="expr.txt")
#也可以保存为csv文件：write.exprs(x,file="expr.csv")

我们用excel打开：第一列为该芯片的探针号，我们可以根据基因芯片平台的信息对应至相应的基因。后面6列为不同样本的表达量。

exp.png

绘图了解芯片数据概况

#microarray fig
par(mfrow=c(2,3))
image(a[,c(1,2,3,4,5,6)])

由于数据较多，生成图像可能比较缓慢。

microarray.png

par(mfrow=c(2,3))
hist(log2(intensity(a[,1])),breaks=100,col="blue",main="GSM712682")
hist(log2(intensity(a[,2])),breaks=100,col="blue",main="GSM712683")
hist(log2(intensity(a[,3])),breaks=100,col="blue",main="GSM712684")
hist(log2(intensity(a[,4])),breaks=100,col="blue",main="GSM712688")
hist(log2(intensity(a[,5])),breaks=100,col="blue",main="GSM712689")
hist(log2(intensity(a[,6])),breaks=100,col="blue",main="GSM712690")

由于图片保存的较小，且数据较为密集，所以此处蓝色看起来如同黑色。

microarray_2.png

hist((as.numeric(exprs(x),col="blue")),main="Gene expression")

microarray_3.png

探针水平数据

boxplot(a,col="red",main="Probe Level Data",names=c("82", "83", "84", "88", "89", "90"))
#82，83...表示GSM号后两位

Rplot.png

基因表达数据

boxplot(data.frame(exprs(x)),col="blue",main="Gene Expression Data",names=c("82", "83", "84", "88", "89", "90"))
#进行log2变换
boxplot(log2(data.frame(exprs(x))),col="blue",main="Gene Expression Data (log2)",names=c("82", "83", "84", "88", "89", "90"))

Rplot01.png
Rplot02.png

我们还可以分析两个样本之间的探针情况。但由于本例中样本数量较多不进行详细说明。

#[]中的参数是选择具体的样本
plot(exprs(a)[,c(1,4)],log="xy",pch=".",main=" All probes")#散点图显示所有探针
plot(pm(a)[,c(1,4)],log="xy",pch=".",main=" Matched Probes")# 图中仅显示匹配探针
plot(mm(a)[,c(1,4)],log="xy",pch=".",main=" Unmatched Probes")#图中仅显示失配探针
plot(exprs(x)[,1:2],pch=".",main="Gene expression data")#基因表达情况

probe.png
image.png

MVA图

上三角对M-A作图，M值为两块芯片对应探针的比率值的对数；A值是对应探针数据的对数均值
下三角显示对应M-A图中M值的中值和四分位距
对角线为芯片名称

par(mfrow=c(1,1))
pms<-pm(a[,c(1,4)])
mva.pairs((pms))

MVA.png