Workshop.R

Code snipits and programs written in R, S or S-PLUS

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

setwd("C:/DATA/Microarray")

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 False discovery rate control (FDR) p < 0.05
2. This is wrong, it should be the adusted p value list which comes out of topTable
3. 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)