Difference between revisions of "AffyVsInHouse.R"

From Organic Design wiki
 
({{R}})
 
(3 intermediate revisions by the same user not shown)
Line 1: Line 1:
 +
# {{R}}
 
library(limma)
 
library(limma)
 
library(affy)
 
library(affy)
Line 9: Line 10:
 
# -------------------------------- Affymetrix --------------------------------- #
 
# -------------------------------- Affymetrix --------------------------------- #
  
−
if(0) { # Change for HortResearch
+
if(1) { # Change for HortResearch
 
   dataDir <- "/Volumes/HD2/Max Planck/Data/Affy/DayNight/Celfiles"
 
   dataDir <- "/Volumes/HD2/Max Planck/Data/Affy/DayNight/Celfiles"
 
} else {
 
} else {
Line 31: Line 32:
 
library(limma)
 
library(limma)
  
−
if(0) { # Change for HortResearch
+
if(1) { # Change for HortResearch
−
    
+
   dataDir <- "/Volumes/HD2/Max\ Planck/HortResearch/VariabilityStudy/Data"
 
}else {
 
}else {
 
   dataDir <- "/Users/admin/Desktop/Directories/VariabilityStudy/Data"
 
   dataDir <- "/Users/admin/Desktop/Directories/VariabilityStudy/Data"
Line 44: Line 45:
 
pairs(log2(RG$R), pch=".")
 
pairs(log2(RG$R), pch=".")
 
pairs(log2(RG$G), pch=".")
 
pairs(log2(RG$G), pch=".")
−
 
+
dev.off()
 
RG <- RG[,c("AC3","AC4")]
 
RG <- RG[,c("AC3","AC4")]
  
Line 56: Line 57:
 
# Normalization (loess, printtiploess)
 
# Normalization (loess, printtiploess)
 
MA      <- normalizeWithinArrays(RG, method="loess", bc.method="none")
 
MA      <- normalizeWithinArrays(RG, method="loess", bc.method="none")
 +
MA      <- normalizeBetweenArrays(MA, method="scale")
 
RG.norm <- RG.MA(MA)
 
RG.norm <- RG.MA(MA)
 +
 +
# summary stats
 +
iqr <- function(x, qlims=c(0.25, 0.75)) {
 +
  IQR <- quantile(x, qlims[2]) - quantile(x, qlims[1])
 +
  return(IQR)
 +
}
 +
 +
# Use median and mad for skewed chisq distributions
 +
calcStats <- function(x, type="median", format = "%2.3f", addText = TRUE) {
 +
  xStat <- c()
 +
  if(type=="median") {
 +
    medx  <- sprintf(format, median(x))
 +
    if(addText){
 +
      xStat <- paste(expression(median(x)), "=", medx, sep=" ")
 +
    } else {
 +
      xStat <- medx
 +
    }
 +
  } else {
 +
    madx  <- sprintf(format, mad(x))
 +
    if(addText){
 +
      xStat <- paste(expression(mad(x)), "=", madx, sep=" ")
 +
    } else {
 +
      xStat <- madx
 +
    }
 +
  }
 +
  return(xStat)
 +
}
 +
# Test
 +
calcStats(1:10, type="median")
 +
median(1:10)
 +
calcStats(1:10, type="mad")
 +
mad(1:10)
  
 
# --------------------------- Comparison plots -------------------------------- #
 
# --------------------------- Comparison plots -------------------------------- #
Line 68: Line 102:
 
X11(xpos=0, ypos=0, width=size, height=size)
 
X11(xpos=0, ypos=0, width=size, height=size)
 
X11(xpos=600, ypos=0, width=size, height=size)
 
X11(xpos=600, ypos=0, width=size, height=size)
 +
X11(xpos=0, ypos=600+45, width=size/2, height=size/2)
 +
 
dev.list()
 
dev.list()
  
Line 84: Line 120:
 
Rmean <- apply(log2(RG.norm$R), 1, mean)
 
Rmean <- apply(log2(RG.norm$R), 1, mean)
 
Gmean <- apply(log2(RG.norm$G), 1, mean)
 
Gmean <- apply(log2(RG.norm$G), 1, mean)
 +
 +
RG1mean <- apply(log2(cbind(RG.norm$R[,1], RG.norm$G[,2])),1, mean)
 +
RG2mean <- apply(log2(cbind(RG.norm$R[,2], RG.norm$G[,1])),1, mean)
  
 
dev.set(2)
 
dev.set(2)
 
hist(emean, main=AffyMain, xlim=c(4,16), xlab="Average signal")
 
hist(emean, main=AffyMain, xlim=c(4,16), xlab="Average signal")
 +
text(x=15, y=1700, label=calcStats(emean, type="median"), adj=c(1,0))
 +
text(x=15, y=1550, label=calcStats(emean, type="mad"),  adj=c(1,0))
 +
 
dev.set(3)
 
dev.set(3)
 
hist(Rmean, main=InHouseMain, xlim=c(4,16), xlab="Cy5 Average signal")
 
hist(Rmean, main=InHouseMain, xlim=c(4,16), xlab="Cy5 Average signal")
 +
text(x=8, y=1700, label=calcStats(Rmean, type="median",
 +
format="%2.1f"), adj=c(1,0))
 +
text(x=8, y=1550, label=calcStats(Rmean, type="mad", format="%2.1f"),
 +
adj=c(1,0))
 +
 
Sys.sleep(1)
 
Sys.sleep(1)
 
hist(Gmean, main=InHouseMain, xlim=c(4,16), xlab="Cy3 Average signal")
 
hist(Gmean, main=InHouseMain, xlim=c(4,16), xlab="Cy3 Average signal")
 +
text(x=8, y=1700, label=calcStats(Gmean, type="median",
 +
format="%2.1f"), adj=c(1,0))
 +
text(x=8, y=1550, label=calcStats(Gmean, type="mad", format="%2.1f"),
 +
adj=c(1,0))
 +
 +
dev.set(4)
 +
# Dye pair not much different in mean signal
 +
hist(RG1mean, main=InHouseMain, xlim=c(4,16), xlab="Cy3/Cy5 Average signal")
 +
text(x=11, y=1700, label=calcStats(RG1mean, type="median"), adj=c(1,0))
 +
text(x=11, y=1550, label=calcStats(RG1mean, type="mad"),  adj=c(1,0))
 +
 +
hist(RG2mean, main=InHouseMain, xlim=c(4,16), xlab="Cy3/Cy5 Average signal")
 +
text(x=11, y=1700, label=calcStats(RG2mean, type="median"), adj=c(1,0))
 +
text(x=11, y=1550, label=calcStats(RG2mean, type="mad"),  adj=c(1,0))
  
 
# SD distributions
 
# SD distributions
Line 96: Line 157:
 
Rsd <- apply(log2(RG.norm$R), 1, sd)
 
Rsd <- apply(log2(RG.norm$R), 1, sd)
 
Gsd <- apply(log2(RG.norm$G), 1, sd)
 
Gsd <- apply(log2(RG.norm$G), 1, sd)
 +
 +
RG1sd <- apply(log2(cbind(RG.norm$R[,1], RG.norm$G[,2])), 1, sd)
 +
RG2sd <- apply(log2(cbind(RG.norm$R[,2], RG.norm$G[,1])), 1, sd)
 +
  
 
dev.set(2)
 
dev.set(2)
 
hist(esd, main=AffyMain, breaks = 12, xlim=c(0,1.5), xlab="Standard deviation")
 
hist(esd, main=AffyMain, breaks = 12, xlim=c(0,1.5), xlab="Standard deviation")
 +
text(x=1, y=8000, label=calcStats(esd, type="median"), adj=c(1,0))
 +
text(x=1, y=7300, label=calcStats(esd, type="mad"),  adj=c(1,0))
 +
 
dev.set(3)
 
dev.set(3)
−
hist(Rsd, main=InHouseMain, breaks = 100, xlim=c(0,1.5), xlab="Cy3 interspot standard deviation")  
+
hist(Rsd, main=InHouseMain, breaks = 100, xlim=c(0,1.5), xlab="Cy3 interspot standard deviation")
 +
text(x=1, y=3000, label=calcStats(Rsd, type="median"), adj=c(1,0))
 +
text(x=1, y=2700, label=calcStats(Rsd, type="mad"),  adj=c(1,0))
 
Sys.sleep(1)
 
Sys.sleep(1)
 
hist(Gsd, main=InHouseMain, breaks = 100, xlim=c(0,1.5), xlab="Cy5 interspot standard deviation")
 
hist(Gsd, main=InHouseMain, breaks = 100, xlim=c(0,1.5), xlab="Cy5 interspot standard deviation")
 +
text(x=1, y=3000, label=calcStats(Gsd, type="median"), adj=c(1,0))
 +
text(x=1, y=2700, label=calcStats(Gsd, type="mad"),  adj=c(1,0))
  
−
# summary stats
+
dev.set(4)
−
iqr <- function(x, qlims=c(0.25, 0.75)) {
+
# Dye pair seems to add about another 10% noise (0.107/0.113, 0.107/0.114)
−
  IQR <- quantile(x, qlims[2]) - quantile(x, qlims[1])
+
hist(RG1sd, main=InHouseMain, breaks = 100, xlim=c(0,1.5), xlab="Cy3/Cy5 interspot standard deviation")
−
  return(IQR)
+
text(x=1, y=3000, label=calcStats(RG1sd, type="median"), adj=c(1,0))
−
}
+
text(x=1, y=2700, label=calcStats(RG1sd, type="mad"),  adj=c(1,0))
−
 
 
−
mean(esd)
 
−
mean(Rsd)
 
−
mean(Gsd)
 
−
 
 
−
sqrt(mean(esd^2))
 
−
sqrt(mean(Rsd^2))
 
−
sqrt(mean(Gsd^2))
 
−
 
 
−
iqr(esd)
 
−
iqr(Rsd)
 
−
iqr(Gsd)
 
  
 +
hist(RG2sd, main=InHouseMain, breaks = 100, xlim=c(0,1.5), xlab="Cy3/Cy5 interspot standard deviation")
 +
text(x=1, y=3000, label=calcStats(RG2sd, type="median"), adj=c(1,0))
 +
text(x=1, y=2700, label=calcStats(RG2sd, type="mad"),  adj=c(1,0))
  
 
# CV
 
# CV
 
dev.set(2)
 
dev.set(2)
−
hist(esd/emean, main=AffyMain, breaks = 12, xlim=c(0,0.2), xlab="Coefficient of variation")
+
hist(esd/emean, main=AffyMain, breaks = 12, xlim=c(0,0.2),
 +
xlab="Coefficient of variation")
 +
text(x=0.15, y=6000, label=calcStats(esd/emean, type="median"), adj=c(1,0))
 +
text(x=0.15, y=5500, label=calcStats(esd/emean, type="mad"),  adj=c(1,0))
 +
 
 
dev.set(3)
 
dev.set(3)
 
hist(Rsd/Rmean, main=InHouseMain, breaks=80, xlim=c(0,0.2), xlab="Cy5 interspot coefficient of variation")
 
hist(Rsd/Rmean, main=InHouseMain, breaks=80, xlim=c(0,0.2), xlab="Cy5 interspot coefficient of variation")
 +
text(x=0.15, y=3000, label=calcStats(Rsd/Rmean, type="median"), adj=c(1,0))
 +
text(x=0.15, y=2700, label=calcStats(Rsd/Rmean, type="mad"),  adj=c(1,0))
 
Sys.sleep(1)
 
Sys.sleep(1)
 
hist(Gsd/Gmean, main=InHouseMain, breaks=80, xlim=c(0,0.2), xlab="Cy3 interspot coefficient of variation")
 
hist(Gsd/Gmean, main=InHouseMain, breaks=80, xlim=c(0,0.2), xlab="Cy3 interspot coefficient of variation")
 +
text(x=0.15, y=3000, label=calcStats(Gsd/Gmean, type="median"), adj=c(1,0))
 +
text(x=0.15, y=2700, label=calcStats(Gsd/Gmean, type="mad"),  adj=c(1,0))
 +
 +
dev.set(4)
 +
# Checking a Dye swap, also about a 10% error added dye to dye swap (0.009/0.01)
 +
hist(RG1sd/RG1mean, main=InHouseMain, breaks=80, xlim=c(0,0.2), xlab="Cy3/Cy5 interspot coefficient of variation")
 +
text(x=0.15, y=3000, label=calcStats(RG1sd/RG1mean, type="median"), adj=c(1,0))
 +
text(x=0.15, y=2700, label=calcStats(RG1sd/RG1mean, type="mad"),  adj=c(1,0))
 +
 +
hist(RG2sd/RG2mean, main=InHouseMain, breaks=80, xlim=c(0,0.2), xlab="Cy3/Cy5 interspot coefficient of variation")
 +
text(x=0.15, y=3000, label=calcStats(RG2sd/RG1mean, type="median"), adj=c(1,0))
 +
text(x=0.15, y=2600, label=calcStats(RG2sd/RG1mean, type="mad"),  adj=c(1,0))
 +
 +
 +
# ------------------------- Summary diagnostics ------------------------------- #
 +
# Means
 +
signal <-
 +
  c(
 +
    calcStats(emean, type="median", addText=FALSE),
 +
    calcStats(Rmean, type="median", addText=FALSE),
 +
    calcStats(Gmean, type="median", addText=FALSE),
 +
    calcStats(RG1mean, type="median", addText=FALSE),
 +
    calcStats(RG2mean, type="median", addText=FALSE),
 +
    )
 +
noise <-
 +
  c(
 +
    calcStats(emean, type="mad", addText=FALSE),
 +
    calcStats(Rmean, type="mad", addText=FALSE),
 +
    calcStats(Gmean, type="mad", addText=FALSE),
 +
    calcStats(RG1mean, type="mad", addText=FALSE),
 +
    calcStats(RG2mean, type="mad", addText=FALSE)
 +
    )
 +
 +
xmean <- cbind(signal, noise)
 +
rownames(xmean) <- c("Affy", "Cy5", "Cy3","Cy3/Cy5 1", "Cy3/Cy5 2")
 +
 +
#SD's
 +
signal <-
 +
  c(
 +
    calcStats(esd, type="median", addText=FALSE),
 +
    calcStats(Rsd, type="median", addText=FALSE),
 +
    calcStats(Gsd, type="median", addText=FALSE),
 +
    calcStats(RG1sd, type="median", addText=FALSE),
 +
    calcStats(RG2sd, type="median", addText=FALSE),
 +
    )
 +
noise <-
 +
  c(
 +
    calcStats(esd, type="mad", addText=FALSE),
 +
    calcStats(Rsd, type="mad", addText=FALSE),
 +
    calcStats(Gsd, type="mad", addText=FALSE),
 +
    calcStats(RG1sd, type="mad", addText=FALSE),
 +
    calcStats(RG2sd, type="mad", addText=FALSE)
 +
    )
 +
 +
xsd <- cbind(signal, noise)
 +
rownames(xsd) <- c("Affy", "Cy5", "Cy3","Cy3/Cy5 1", "Cy3/Cy5 2")
 +
 +
# CV's
 +
signal <-
 +
  c(
 +
    calcStats(esd/emean, type="median", addText=FALSE),
 +
    calcStats(Rsd/Rmean, type="median", addText=FALSE),
 +
    calcStats(Gsd/Gmean, type="median", addText=FALSE),
 +
    calcStats(RG1sd/RG1mean, type="median", addText=FALSE),
 +
    calcStats(RG2sd/RG2mean, type="median", addText=FALSE),
 +
    )
 +
noise <-
 +
  c(
 +
    calcStats(esd/emean, type="mad", addText=FALSE),
 +
    calcStats(Rsd/Rmean, type="mad", addText=FALSE),
 +
    calcStats(Gsd/Gmean, type="mad", addText=FALSE),
 +
    calcStats(RG1sd/RG1mean, type="mad", addText=FALSE),
 +
    calcStats(RG2sd/RG2mean, type="mad", addText=FALSE)
 +
    )
 +
 +
xcv <- cbind(signal, noise)
 +
rownames(xcv) <- c("Affy", "Cy5", "Cy3","Cy3/Cy5 1", "Cy3/Cy5 2")
 +
 +
print(xmean)
 +
print(xsd)
 +
print(xcv)

Latest revision as of 02:16, 29 May 2007

Code snipits and programs written in R, S or S-PLUS library(limma) library(affy)

packageDescription("limma", field="Version") packageDescription("affy", field="Version")

vignette("affy")

  1. -------------------------------- Affymetrix --------------------------------- #

if(1) { # Change for HortResearch 

 dataDir <- "/Volumes/HD2/Max Planck/Data/Affy/DayNight/Celfiles"

} else { 

 dataDir  <- "/Users/admin/Desktop/DayNight/Celfiles/"
 Sys.putenv("DISPLAY"=":0")

}

dset<- ReadAffy(filenames=file.path(dataDir, dir(dataDir, pattern=".CEL")), widget = F) # loads CEL files into an affybatch object

un <- ".CEL" # remove extra names sampleNames(dset) <- gsub(un, "", sampleNames(dset))

  1. Obtaining indexes of sampleNames (affy slides) of interest

technicalreps <- grep("00 G048", sampleNames(dset)) techset <- dset[,technicalreps]

  1. Normalization

erma <- rma(techset)

  1. --------------------------------- In house ---------------------------------- #

library(limma)

if(1) { # Change for HortResearch 

 dataDir <- "/Volumes/HD2/Max\ Planck/HortResearch/VariabilityStudy/Data"

}else { 

 dataDir <- "/Users/admin/Desktop/Directories/VariabilityStudy/Data"

} files <- dir(dataDir, pattern="gpr")

  1. Examine genepix, genepix.median

RG <- read.maimages(files, path=dataDir, source="genepix", wt.fun=wtflags(0))

  1. Visually AC3 and AC4 most similar

pairs(log2(RG$R), pch=".") pairs(log2(RG$G), pch=".") dev.off() RG <- RG[,c("AC3","AC4")]

  1. All spots

nrow(RG)

  1. Good spots

apply(RG$weights, 2, sum)

  1. Bad spots

nrow(RG) - apply(RG$weights, 2, sum)

  1. Normalization (loess, printtiploess)

MA <- normalizeWithinArrays(RG, method="loess", bc.method="none") MA <- normalizeBetweenArrays(MA, method="scale") RG.norm <- RG.MA(MA)

  1. summary stats

iqr <- function(x, qlims=c(0.25, 0.75)) { 

 IQR <- quantile(x, qlims[2]) - quantile(x, qlims[1])
 return(IQR)

}

  1. Use median and mad for skewed chisq distributions

calcStats <- function(x, type="median", format = "%2.3f", addText = TRUE) { 

 xStat <- c()
 if(type=="median") {
   medx  <- sprintf(format, median(x))
   if(addText){
     xStat <- paste(expression(median(x)), "=", medx, sep=" ")
   } else {
     xStat <- medx
   }
 } else {
   madx  <- sprintf(format, mad(x))
   if(addText){
     xStat <- paste(expression(mad(x)), "=", madx, sep=" ")
   } else {
     xStat <- madx
   }
 }
 return(xStat)

}

  1. Test

calcStats(1:10, type="median") median(1:10) calcStats(1:10, type="mad") mad(1:10)

  1. --------------------------- Comparison plots -------------------------------- #
  1. Setup

size <- 8 AffyMain <- "Affymetrix technical replication" InHouseMain <- "In house technical replication"

graphics.off() X11(xpos=0, ypos=0, width=size, height=size) X11(xpos=600, ypos=0, width=size, height=size) X11(xpos=0, ypos=600+45, width=size/2, height=size/2)

dev.list()

  1. Pairs plots

pairs(log2(RG.norm$R), pch=".") pairs(log2(RG.norm$G), pch=".")

dev.set(2) plot(exprs(erma), main = AffyMain, pch=".") dev.set(3) plot(log2(RG.norm$R), main = InHouseMain, pch=".")

  1. Histograms
  2. Mean distributions

emean <- apply(exprs(erma), 1, mean) Rmean <- apply(log2(RG.norm$R), 1, mean) Gmean <- apply(log2(RG.norm$G), 1, mean)

RG1mean <- apply(log2(cbind(RG.norm$R[,1], RG.norm$G[,2])),1, mean) RG2mean <- apply(log2(cbind(RG.norm$R[,2], RG.norm$G[,1])),1, mean)

dev.set(2) hist(emean, main=AffyMain, xlim=c(4,16), xlab="Average signal") text(x=15, y=1700, label=calcStats(emean, type="median"), adj=c(1,0)) text(x=15, y=1550, label=calcStats(emean, type="mad"), adj=c(1,0))

dev.set(3) hist(Rmean, main=InHouseMain, xlim=c(4,16), xlab="Cy5 Average signal") text(x=8, y=1700, label=calcStats(Rmean, type="median", format="%2.1f"), adj=c(1,0)) text(x=8, y=1550, label=calcStats(Rmean, type="mad", format="%2.1f"), adj=c(1,0))

Sys.sleep(1) hist(Gmean, main=InHouseMain, xlim=c(4,16), xlab="Cy3 Average signal") text(x=8, y=1700, label=calcStats(Gmean, type="median", format="%2.1f"), adj=c(1,0)) text(x=8, y=1550, label=calcStats(Gmean, type="mad", format="%2.1f"), adj=c(1,0))

dev.set(4)

  1. Dye pair not much different in mean signal

hist(RG1mean, main=InHouseMain, xlim=c(4,16), xlab="Cy3/Cy5 Average signal") text(x=11, y=1700, label=calcStats(RG1mean, type="median"), adj=c(1,0)) text(x=11, y=1550, label=calcStats(RG1mean, type="mad"), adj=c(1,0))

hist(RG2mean, main=InHouseMain, xlim=c(4,16), xlab="Cy3/Cy5 Average signal") text(x=11, y=1700, label=calcStats(RG2mean, type="median"), adj=c(1,0)) text(x=11, y=1550, label=calcStats(RG2mean, type="mad"), adj=c(1,0))

  1. SD distributions

esd <- apply(exprs(erma), 1, sd) Rsd <- apply(log2(RG.norm$R), 1, sd) Gsd <- apply(log2(RG.norm$G), 1, sd)

RG1sd <- apply(log2(cbind(RG.norm$R[,1], RG.norm$G[,2])), 1, sd) RG2sd <- apply(log2(cbind(RG.norm$R[,2], RG.norm$G[,1])), 1, sd)


dev.set(2) hist(esd, main=AffyMain, breaks = 12, xlim=c(0,1.5), xlab="Standard deviation") text(x=1, y=8000, label=calcStats(esd, type="median"), adj=c(1,0)) text(x=1, y=7300, label=calcStats(esd, type="mad"), adj=c(1,0))

dev.set(3) hist(Rsd, main=InHouseMain, breaks = 100, xlim=c(0,1.5), xlab="Cy3 interspot standard deviation") text(x=1, y=3000, label=calcStats(Rsd, type="median"), adj=c(1,0)) text(x=1, y=2700, label=calcStats(Rsd, type="mad"), adj=c(1,0)) Sys.sleep(1) hist(Gsd, main=InHouseMain, breaks = 100, xlim=c(0,1.5), xlab="Cy5 interspot standard deviation") text(x=1, y=3000, label=calcStats(Gsd, type="median"), adj=c(1,0)) text(x=1, y=2700, label=calcStats(Gsd, type="mad"), adj=c(1,0))

dev.set(4)

  1. Dye pair seems to add about another 10% noise (0.107/0.113, 0.107/0.114)

hist(RG1sd, main=InHouseMain, breaks = 100, xlim=c(0,1.5), xlab="Cy3/Cy5 interspot standard deviation") text(x=1, y=3000, label=calcStats(RG1sd, type="median"), adj=c(1,0)) text(x=1, y=2700, label=calcStats(RG1sd, type="mad"), adj=c(1,0))

hist(RG2sd, main=InHouseMain, breaks = 100, xlim=c(0,1.5), xlab="Cy3/Cy5 interspot standard deviation") text(x=1, y=3000, label=calcStats(RG2sd, type="median"), adj=c(1,0)) text(x=1, y=2700, label=calcStats(RG2sd, type="mad"), adj=c(1,0))

  1. CV

dev.set(2) hist(esd/emean, main=AffyMain, breaks = 12, xlim=c(0,0.2), xlab="Coefficient of variation") text(x=0.15, y=6000, label=calcStats(esd/emean, type="median"), adj=c(1,0)) text(x=0.15, y=5500, label=calcStats(esd/emean, type="mad"), adj=c(1,0))

dev.set(3) hist(Rsd/Rmean, main=InHouseMain, breaks=80, xlim=c(0,0.2), xlab="Cy5 interspot coefficient of variation") text(x=0.15, y=3000, label=calcStats(Rsd/Rmean, type="median"), adj=c(1,0)) text(x=0.15, y=2700, label=calcStats(Rsd/Rmean, type="mad"), adj=c(1,0)) Sys.sleep(1) hist(Gsd/Gmean, main=InHouseMain, breaks=80, xlim=c(0,0.2), xlab="Cy3 interspot coefficient of variation") text(x=0.15, y=3000, label=calcStats(Gsd/Gmean, type="median"), adj=c(1,0)) text(x=0.15, y=2700, label=calcStats(Gsd/Gmean, type="mad"), adj=c(1,0))

dev.set(4)

  1. Checking a Dye swap, also about a 10% error added dye to dye swap (0.009/0.01)

hist(RG1sd/RG1mean, main=InHouseMain, breaks=80, xlim=c(0,0.2), xlab="Cy3/Cy5 interspot coefficient of variation") text(x=0.15, y=3000, label=calcStats(RG1sd/RG1mean, type="median"), adj=c(1,0)) text(x=0.15, y=2700, label=calcStats(RG1sd/RG1mean, type="mad"), adj=c(1,0))

hist(RG2sd/RG2mean, main=InHouseMain, breaks=80, xlim=c(0,0.2), xlab="Cy3/Cy5 interspot coefficient of variation") text(x=0.15, y=3000, label=calcStats(RG2sd/RG1mean, type="median"), adj=c(1,0)) text(x=0.15, y=2600, label=calcStats(RG2sd/RG1mean, type="mad"), adj=c(1,0))


  1. ------------------------- Summary diagnostics ------------------------------- #
  2. Means

signal <- 

 c(
   calcStats(emean, type="median", addText=FALSE),
   calcStats(Rmean, type="median", addText=FALSE),
   calcStats(Gmean, type="median", addText=FALSE),
   calcStats(RG1mean, type="median", addText=FALSE),
   calcStats(RG2mean, type="median", addText=FALSE),
   )

noise <- 

 c(
   calcStats(emean, type="mad", addText=FALSE),
   calcStats(Rmean, type="mad", addText=FALSE),
   calcStats(Gmean, type="mad", addText=FALSE),
   calcStats(RG1mean, type="mad", addText=FALSE),
   calcStats(RG2mean, type="mad", addText=FALSE)
   )

xmean <- cbind(signal, noise) rownames(xmean) <- c("Affy", "Cy5", "Cy3","Cy3/Cy5 1", "Cy3/Cy5 2")

  1. SD's

signal <- 

 c(
   calcStats(esd, type="median", addText=FALSE),
   calcStats(Rsd, type="median", addText=FALSE),
   calcStats(Gsd, type="median", addText=FALSE),
   calcStats(RG1sd, type="median", addText=FALSE),
   calcStats(RG2sd, type="median", addText=FALSE),
   )

noise <- 

 c(
   calcStats(esd, type="mad", addText=FALSE),
   calcStats(Rsd, type="mad", addText=FALSE),
   calcStats(Gsd, type="mad", addText=FALSE),
   calcStats(RG1sd, type="mad", addText=FALSE),
   calcStats(RG2sd, type="mad", addText=FALSE)
   )

xsd <- cbind(signal, noise) rownames(xsd) <- c("Affy", "Cy5", "Cy3","Cy3/Cy5 1", "Cy3/Cy5 2")

  1. CV's

signal <- 

 c(
   calcStats(esd/emean, type="median", addText=FALSE),
   calcStats(Rsd/Rmean, type="median", addText=FALSE),
   calcStats(Gsd/Gmean, type="median", addText=FALSE),
   calcStats(RG1sd/RG1mean, type="median", addText=FALSE),
   calcStats(RG2sd/RG2mean, type="median", addText=FALSE),
   )

noise <- 

 c(
   calcStats(esd/emean, type="mad", addText=FALSE),
   calcStats(Rsd/Rmean, type="mad", addText=FALSE),
   calcStats(Gsd/Gmean, type="mad", addText=FALSE),
   calcStats(RG1sd/RG1mean, type="mad", addText=FALSE),
   calcStats(RG2sd/RG2mean, type="mad", addText=FALSE)
   )

xcv <- cbind(signal, noise) rownames(xcv) <- c("Affy", "Cy5", "Cy3","Cy3/Cy5 1", "Cy3/Cy5 2")

print(xmean) print(xsd) print(xcv)

Retrieved from "https://organicdesign.nz/wiki/index.php?title=AffyVsInHouse.R&oldid=70839"