11、权重残差图、RLE和NUSE

　　affyPLM包可以对芯片原始数据进行拟合回归，最后得到芯片权重（Weights）残差（Residuals）图、相对对数表达（RLE，Relative log expression）箱线图、相对标准差（NUSE，Normalized unscaled standard errors）箱线图

 

以样品GSM286757.CEL、GSM286756.CEL、GSM286763.CEL、GSM286762.CEL、GSM286759.CEL、GSM286760.CEL、GSM286765.CEL、GSM286766.CEL为例：

 

library(affyPLM) 

rawData<-ReadAffy("GSM286757.CEL","GSM286756.CEL",

"GSM286763.CEL","GSM286762.CEL",

"GSM286759.CEL","GSM286760.CEL",

"GSM286765.CEL","GSM286766.CEL")

Pset <- fitPLM(rawData)

 

boxplot(Pset,col=c(1:8),main="NUSE")    ## NUSE图

Mbox(Pset,col=c(1:8),main="RLE")      ## RLE图

 

image(Pset,type="weights",which=1,main="Weights") ## 权重图

image(Pset, type="resids", which=2, main="Residuals") ## 残差图

image(Pset, type="sign.resids", which=2, main="Residuals sign") ## 符号残差图

 

 

从上面的代码可以看出，经过了fitPLM的计算得到了权重参差、相对对数表达、相对标准差的数据，但是这些数据在Pset中是怎样存储的呢？

 

 

 

 

运行下面的代码可以看清楚：

model=PM ~ -1 + probes + samples

model.param=verify.model.param(rawData,model)

variable.type <- verify.variable.types (model,c(default="factor"))

constraint.type <- verify.constraint.types(model,c(default="contr.treatment"))

 

 

n.probesets <- length(geneNames(rawData))

R.model <- PLM.designmatrix3(rawData,model,variable.type=variable.type,constraint.type=constraint.type)

output <- verify.output.param()

modelparam <- verify.model.param(rawData,model,model.param=model.param)

background.param <- verify.bg.param(R.model, "RMA.2",background.param = list())

normalize.param <- verify.norm.param(R.model, "quantile",normalize.param=list())

 

Fitresults <- .Call("R_rlm_PLMset_c",pm(rawData),mm(rawData),

                                     probeNames(rawData),

                                     n.probesets,

                                     R.model,

                                     output,

                                     modelparam,

                                     TRUE, 

                                     "RMA.2",

                                     background.param, 

                                     TRUE, 

                                     "quantile", 

                                     normalize.param,

                                     0,

                                     PACKAGE="affyPLM")

 

 

其中：

 

一、Fitresults[[4]]和NUSE有关，它是一个阵列，行数是探针组数目，列数是样品数，长度是探针组数目*样品数

 

grp.rma.se1.median <- apply(Fitresults[[4]], 1,median,na.rm=TRUE)

grp.rma.rel.se1.mtx <- sweep(Fitresults[[4]],1,grp.rma.se1.median,FUN='/')

## 以上2步操作是让每一行都减去该行的中位数

 

boxplot(grp.rma.rel.se1.mtx,col=c(1,2,3,4,5,6,7,8),main="NUSE")

 

二、Fitresults[[1]]和RLE有关，它是一个阵列，行数是探针组数目，列数是样品数，长度是探针组数目*样品数

 

medianchip <- apply(Fitresults[[1]], 1, median)

M <- sweep(Fitresults[[1]],1,medianchip,FUN='-')

## 以上2步操作是让每一行都减去该行的中位数

 

boxplot(M,col=c(1,2,3,4,5,6,7,8),main="RLE")

 

三、Fitresults[[3]][[1]]和权重图有关，它是一个阵列，行数是PM探针数目，列数是样品数，长度是探针数目*样品数。在本例中，PM探针数目是604258，样品数是8，那么第一个样品的权重值是Fitresults[[3]][[1]][,1]，长度为探针数目604258

 

## 查看第一个样品的前50个权重值

> Fitresults[[3]][[1]][,1][1:50]

1007_s_at 1007_s_at 1007_s_at 1007_s_at 1007_s_at 1007_s_at 1007_s_at 

1.0000000 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000 

1007_s_at 1007_s_at 1007_s_at 1007_s_at 1007_s_at 1007_s_at 1007_s_at 

1.0000000 0.8790510 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000 

1007_s_at 1007_s_at   1053_at   1053_at   1053_at   1053_at   1053_at 

1.0000000 1.0000000 1.0000000 1.0000000 1.0000000 0.2046946 1.0000000 

  1053_at   1053_at   1053_at   1053_at   1053_at   1053_at   1053_at 

1.0000000 1.0000000 1.0000000 0.4951793 1.0000000 1.0000000 1.0000000 

  1053_at   1053_at   1053_at   1053_at    117_at    117_at    117_at 

1.0000000 0.4462245 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000 

   117_at    117_at    117_at    117_at    117_at    117_at    117_at 

1.0000000 1.0000000 0.6728794 1.0000000 0.3876992 0.8266238 0.7217806 

   117_at    117_at    117_at    117_at    117_at    117_at    121_at 

1.0000000 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000 

   121_at 

1.0000000 

 

 

## 绘制权重图的颜色

col.weights <- terrain.colors(25) 

> col.weights

 [1] "#00A600FF" "#10AC00FF" "#20B100FF" "#32B700FF" "#45BD00FF" "#59C300FF"

 [7] "#6DC900FF" "#83CE00FF" "#9AD400FF" "#B2DA00FF" "#CBE000FF" "#E6E600FF"

[13] "#E6D612FF" "#E7C924FF" "#E8BF36FF" "#E9B848FF" "#EAB35AFF" "#EBB16DFF"

[19] "#ECB27FFF" "#EDB592FF" "#EEBCA5FF" "#EFC5B8FF" "#F0D1CBFF" "#F1E0DFFF"

[25] "#F2F2F2FF"

 

 

这里有25个颜色，[1]~[25]，从草绿色渐变到橘黄色，再渐变到接近白色。越小的权重值分配到的颜色越接近草绿色，越大的权重值越接近白色。这样，每个PM探针都有了对应的权重值和颜色，绘制成图像就是权重图了。

 

四、Fitresults[[8]] [[1]]和残差图有关系，原理和权重图是一样的。

 

符号残差图根据Fitresults[[8]] [[1]]的数据的正负号来确定颜色，正数红色，0白色，负数蓝色。对Fitresults[[8]] [[1]]进行sign(Fitresults[[8]] [[1]])*(log2(abs(Fitresults[[8]] [[1]])+1))计算后得到的数据用于残差图，残差图的颜色是：

 

col.resids <- pseudoPalette(low="blue",high="red",mid="white")

> col.resids

 [1] "#0000FF" "#0B0BFF" "#1515FF" "#2020FF" "#2A2AFF" "#3535FF" "#4040FF"

 [8] "#4A4AFF" "#5555FF" "#6060FF" "#6A6AFF" "#7575FF" "#8080FF" "#8A8AFF"

[15] "#9595FF" "#9F9FFF" "#AAAAFF" "#B5B5FF" "#BFBFFF" "#CACAFF" "#D4D4FF"

[22] "#DFDFFF" "#EAEAFF" "#F4F4FF" "#FFFFFF" "#FFFFFF" "#FFF4F4" "#FFEAEA"

[29] "#FFDFDF" "#FFD5D5" "#FFCACA" "#FFBFBF" "#FFB5B5" "#FFAAAA" "#FF9F9F"

[36] "#FF9595" "#FF8A8A" "#FF8080" "#FF7575" "#FF6A6A" "#FF6060" "#FF5555"

[43] "#FF4A4A" "#FF4040" "#FF3535" "#FF2B2B" "#FF2020" "#FF1515" "#FF0B0B"

[50] "#FF0000"

 

颜色从蓝色渐变到红色，再渐变到白色。越小的残差值分配到的颜色越接近蓝色，

 

这里有50个颜色，[1]~[50]，从蓝色渐变到红色，再渐变到接近白色。越小的残差值分配到的颜色越接近蓝色，越大的权重值越接近白色。这样，每个PM探针都有了对应的残差值和颜色，绘制成图像就是残差图了。

 

转载于:https://www.cnblogs.com/xianwen/p/4522359.html