Workshop.R

1. ----------------------------------------------------------------------------- #
2. Setup
3. ----------------------------------------------------------------------------- #

1. Load add on package 'limma'

library(limma)

1. Help documentation

limmaUsersGuide()

1. Default working directory

getwd()

1. Setting the working directory

1. Relative paths to directory containing GPR files and results

dataDir <- "GPR" resultsDir <- "Results"

1. ----------------------------------------------------------------------------- #
2. Reading in data
3. ----------------------------------------------------------------------------- #

wt.flags <- function(x) {

 as.numeric(x$Flags >= -50)

}

1. Checking that weight function is working

x <- list() x$Flags <- c(0,-50,-75,-100) x wt.flags(x)

1. Creating RGList object for experiment

RG <- read.maimages( dir(dataDir), path=dataDir, source="genepix.median",

                   wt.fun=wt.flags, other.columns="Flags") 

1. Structure of RG

dim(RG) summary(RG)

for(i in names(RG)) {

 pad <- "--------------"
 print(paste(pad, i, pad))
 str(RGi)

}

1. Creating names for accessing each slide

slideNames <- rownames(RG$targets) slideNames

1. ----------------------------------------------------------------------------- #
2. Summary statistics and diagnostics
3. ----------------------------------------------------------------------------- #

1. Range of Intensity channels

summary(RG$R) summary(RG$Rb) summary(RG$G) summary(RG$Gb)

1. Mapping of spot positions to annotation information

RG$genes[1:30,]

1. Flags diagnostics

for( i in slideNames) {

 pad <- "--------------"
 print(paste(pad, i, pad))
 print(table(RG$other$Flags[,i]))

}

1. Weights diagnostics

for( i in slideNames) {

 pad <- "--------------"
 print(paste(pad, i, pad))
 print(table(RG$weights[,i]))

}

1. MA-plots

par(mfrow=c(2,2)) for( i in slideNames ) {

 # the backgroundCorrect function stops automatic background subtraction
 plotMA(backgroundCorrect(RG, method="none"), array=i)

}

par(mfrow=c(2,2)) for( i in slideNames ) {

 imageplot(RG$Rb[,i], layout=RG$printer, low="white", high="red",
           zlim=c(0,500))

} for( i in slideNames ) {

 imageplot(RG$Gb[,i], layout=RG$printer, low="white", high="green",
           zlim=c(0,500))

}

1. Density plot

plotDensities(backgroundCorrect(RG, method="none"))

1. Padding out suspect spots

RG$R[RG$other$Flags == -75] <- NA RG$G[RG$other$Flags == -75] <- NA RG$R[RG$other$Flags == -100] <- NA RG$G[RG$other$Flags == -100] <- NA

1. ----------------------------------------------------------------------------- #
2. Normalization and Diagnostic images
3. ----------------------------------------------------------------------------- #

1. normalizing printtiploess with no background correction

MA <- normalizeWithinArrays(RG, method="printtiploess", bc.method="none") boxplot(MA$M~col(MA$M),names=colnames(MA$M))

1. Checking range of M

summary(MA$M)

1. Diagnostic images after normalization

par(mfrow=c(2,2)) for( i in slideNames ) {

 plotMA(MA, array=i)

}

for( i in slideNames ) {

 imageplot(MA$M[,i], layout=RG$printer)

}

plotDensities(MA)

1. ----------------------------------------------------------------------------- #
2. Analysis
3. ----------------------------------------------------------------------------- #

1. Dye swap design

design <- c(1,-1,-1,1)

fit <- lmFit(MA, design) fit <- eBayes(fit)

1. The top 30 differentially expressed genes

N <- 30 topTable(fit, number=N)

1. Number of genes with FDR < 0.05

sum(fit$p.value<0.05, na.rm=T)

1. Checking results on an MA-plot

ord <- order(fit$lods,decreasing=TRUE) topN <- ord[1:N] plotMA(fit) points(fit$Amean[topN],fit$coef[topN],pch=16,col="red")

1. ----------------------------------------------------------------------------- #
2. Saving topTable results
3. ----------------------------------------------------------------------------- #

1. Create output directory

dir.create(resultsDir) write.table(topTable(fit, number=nrow(MA)), quote=F, sep="\t",

           file=file.path(resultsDir,"genelist.txt"), row.names=F)