Following the limma protocol

For today’s class we will follow the protocol for two-colors microarrays

• Define the experiment
• Read the hybridization data
• Normalize each slide internally
• Normalize between arrays
• Fit a linear model
• Evaluate statistical significance

Targets

It is recommended to start with a description of the experiments

We write a tab-separated file, indicating at least

• Filename
• Cy3
• Cy5

                                 Filename    Cy3    Cy5
1       GSM3303967_Dcg2699_vs_WT_I.gpr.gz mutant     wt
2      GSM3303968_Dcg2699_vs_WT_II.gpr.gz mutant     wt
3 GSM3303969_Dcg2699_vs_WT_III_csw.gpr.gz     wt mutant

Reading the targets file

library(limma)
targets <- readTargets("targets.txt")
kable(targets)

This example was taken from GSE117566 in NCBI GEO

Reading the arrays

RG <- read.maimages(targets$Filename, source = "genepix")

Read GSM3303967_Dcg2699_vs_WT_I.gpr.gz 
Read GSM3303968_Dcg2699_vs_WT_II.gpr.gz 
Read GSM3303969_Dcg2699_vs_WT_III_csw.gpr.gz 

class(RG)

[1] "RGList"
attr(,"package")
[1] "limma"

names(RG)

[1] "R"       "G"       "Rb"      "Gb"      "targets" "genes"   "source" 
[8] "printer"

dim(RG$R)

[1] 45220     3

Genes

head(RG$genes, n=15)

   Block Row Column              ID            Name
1      1   1      1 GE_BrightCorner GE_BrightCorner
2      1   1      2 GE_BrightCorner GE_BrightCorner
3      1   1      3      DarkCorner      DarkCorner
4      1   1      4      DarkCorner      DarkCorner
5      1   1      5      DarkCorner      DarkCorner
6      1   1      6      DarkCorner      DarkCorner
7      1   1      7      DarkCorner      DarkCorner
8      1   1      8      DarkCorner      DarkCorner
9      1   1      9      DarkCorner      DarkCorner
10     1   1     10      DarkCorner      DarkCorner
11     1   1     11      DarkCorner      DarkCorner
12     1   1     12      DarkCorner      DarkCorner
13     1   1     13      DarkCorner      DarkCorner
14     1   1     14      DarkCorner      DarkCorner
15     1   1     15      DarkCorner      DarkCorner

Comparing R and G channels

for(i in 1:3) plot(RG$R[,i], RG$G[,i])

Using Logarithms

for(i in 1:3) plot(log(RG$R[,i]), log(RG$G[,i]))

Looking at the slides

for(i in 1:3) imageplot(RG$R[,i], RG$printer)
for(i in 1:3) imageplot(RG$G[,i], RG$printer)

Using log

for(i in 1:3) imageplot(log(RG$R[,i]), RG$printer)
for(i in 1:3) imageplot(log(RG$G[,i]), RG$printer)

Background

for(i in 1:3) imageplot(log(RG$Rb[,i]), RG$printer)
for(i in 1:3) imageplot(log(RG$Gb[,i]), RG$printer)

Background correction

The easiest way to normalize is to subtract the background from the signal \[R_\text{corrected} = R - Rb\]

RG_corr <- backgroundCorrect(RG, method = "subtract")
class(RG_corr)

[1] "RGList"
attr(,"package")
[1] "limma"

names(RG_corr)

[1] "R"       "G"       "targets" "genes"   "source"  "printer"

Histogram

There are negative numbers. Not good for logarithms

Histogram of log

Warning in log(RG_corr$R[, 1]): NaNs produced

Avoiding negative values

A better correction requires a better model

RG_corr <- backgroundCorrect(RG, method = "normexp")

Array 1 corrected
Array 2 corrected
Array 3 corrected
Array 1 corrected
Array 2 corrected
Array 3 corrected

Normexp uses a model of background

Background corrected images

for(i in 1:3) imageplot(log(RG$R[,i]), RG$printer)
for(i in 1:3) imageplot(log(RG$G[,i]), RG$printer)

Comparing corrected

for(i in 1:3) plot(log(RG_corr$R[,i]), log(RG_corr$G[,i]))

Changing perspective

We care about the difference of expression \[M = \frac{\log(R)-\log(G)}{2}\]

For symmetry, we also keep the average expression \[ A = \frac{\log(R)+\log(G)}{2}\]

MA plot

for(i in 1:3) plotMA(RG_corr, array=i)

Correcting banana shape

MA <- normalizeWithinArrays(RG_corr, method="loess")
for(i in 1:3) plotMA(MA, array=i)

Comparing densities

Think of these as softened histograms

Between array normalization

MA.q <- normalizeBetweenArrays(MA, method = "quantile")
plotDensities(MA.q)

Let’s drop control spots

gal <- readGAL("GPL16989_Genepix_GAL_4plex_Design_041487.gal.gz")
# num_replicas <- table(MA$genes$ID)
# controls <- names(num_replicas)[num_replicas>=10]
MA <- MA[gal$ControlType == "false", ]
MA <- MA[MA$genes$ID != "EMPTY", ]
MA.q <- normalizeBetweenArrays(MA, method = "quantile")