for (i in 1:6)
{
name = paste("histogram",i,".jpg",sep="")
jpeg(name)
hist(data[,i],lwd=2,which='pm',ylab='Density',xlab='Log2 intensities',main=ph@data$sample[i])
dev.off()
}

The which argument allows you to specify if you want to use:

• perfect match probes only: which=’pm’
• mismatch probes only: which=’mm’
• both: which=’both’

Since the use of MM probes is highly controversial, we will work with perfect match probes only.

The lwd argument defines the line width, xlab and ylab define the titles of the axes and main defines the title of the histogram.

We’ll only use the PM intensities

pmexp = pm(data)

We create two empty vectors that will serve as the two columns of the data frame:

• one to store the sample names in, called sampleNames
• one to store the log intensities in , called logs

The initialization of these vectors is essential. If you run into problems when creating these vectors and you have to repeat the plotting, always reinitialize these vectors (recreate them so that they are empty) 

The dim() method returns the number of rows (probes) and columns (samples) of the matrix containing the PM intensities. The number of rows corresponds to the number of PM intensities (probes) that is available for each sample (in this example that’s 251078). The sampleNames vector needs the same number of rows consisting of sample names.

The logs vector will contain the log PM intensities of all six samples stacked into a single column. Up to this point we have extracted raw PM intensities, they are not yet log-transformed! You can do the log transformation using the log2() method.

sampleNames = vector()
logs = vector()
for (i in 1:6)
{
sampleNames = c(sampleNames,rep(ph@data[i,1],dim(pmexp)[1]))
logs = c(logs,log2(pmexp[,i]))
}

If you have 3 groups of 3 replicates, the code is as follows:

sampleNames = vector()
logs = vector()
for (i in 1:9)
{
sampleNames = c(sampleNames,rep(ph@data[i,1],dim(pmexp)[1]))
logs = c(logs,log2(pmexp[,i]))
}

After we have filled the vectors with the data we need, we combine sample names and log intensities into a single data frame:

logData = data.frame(logInt=logs,sampleName=sampleNames)

Now we can create the plot:

dataHist2 = ggplot(logData, aes(logInt, colour = sampleName)) 
dataHist2 + geom_density()

name = "boxplot.jpg"
jpeg(name)
boxplot(data,which='pm',col='red',names=ph@data$sample) 
dev.off()

The which argument allows you to specify if you want to use:

• perfect match probes only: which=’pm’
• mismatch probes only: which=’mm’
• both: which=’both’

Since the use of MM probes is highly controversial, we will work with PM probes only.
The col argument specifies the color of the boxes, the names argument the labels on the X-axis.

We’ll only use the PM intensities

pmexp = pm(data)

We create two empty vectors that will serve as the two columns of the data frame:

• one to store the sample names in, called sampleNames
• one to store the log intensities in , called logs

The initialization of these vectors is essential. If you run into problems when creating these vectors and you have to repeat the plotting, always reinitialize these vectors (recreate them so that they are empty) 

The dim() method returns the number of rows (probes) and columns (samples) of the matrix containing the PM intensities. The number of rows corresponds to the number of PM intensities (probes) that is available for each sample (in this example that’s 251078). The sampleNames vector needs the same number of rows consisting of sample names.

The logs vector will contain the log PM intensities of all six samples stacked into a single column. Up to this point we have extracted raw PM intensities, they are not yet log-transformed! You can do the log transformation using the log2() method.

sampleNames = vector()
logs = vector()
for (i in 1:6)
{
sampleNames = c(sampleNames,rep(ph@data[i,1],dim(pmexp)[1]))
logs = c(logs,log2(pmexp[,i]))
}

If you have 3 groups of 3 replicates, the code is as follows:

sampleNames = vector()
logs = vector()
for (i in 1:9)
{
sampleNames = c(sampleNames,rep(ph@data[i,1],dim(pmexp)[1]))
logs = c(logs,log2(pmexp[,i]))
}

After we have filled the vectors with the data we need, we combine sample names and log intensities into a single data frame:

logData = data.frame(logInt=logs,sampleName=sampleNames)

Now we can create the plot:

dataBox = ggplot(logData,aes(sampleName,logInt))
dataBox + geom_boxplot()

We create two vectors:

• one to store the sample names in, called sampleNames
• one to store the normalized log intensities in , called normlogs

The dim() method returns the number of rows (probe sets) and columns (samples) of the matrix containing the normalized intensities. The number of rows corresponds to the number of normalized intensities (probe sets) in each sample (in this example it’s 22810). The sampleNames vector needs the same number of entries of each sample name stacked into a single column.

The normlogs vector will contain the normalized intensities of all six samples stacked into a single column.

sampleNames = vector()
normlogs = vector()
for (i in 1:6)
{
sampleNames = c(sampleNames,rep(ph@data[i,1],dim(data.matrix)[1]))
normlogs = c(normlogs,data.matrix[,i])
}

If you have 3 groups of 3 samples, the code is as follows:

sampleNames = vector()
normlogs = vector()
for (i in 1:9)
{
sampleNames = c(sampleNames,rep(ph@data[i,1],dim(data.matrix)[1]))
normlogs = c(normlogs,data.matrix[,i])
}

After having filled the vectors with data, we combine sample names and normalized intensities into a single data frame:

normData = data.frame(norm_logInt=normlogs,sampleName=sampleNames)

Now we can create the plot. To compare this plot to the boxplot we made for raw log intensities we have to set them to the same Y-scale. The scale for the normalized data is much lower than these for the raw data so we have to make the box plot setting the Y-scale from 2 to 16. To set the Y-scale we can use the ylim() method.
To add a title to the plot we use the ggtitle() method.

dataBox = ggplot(normData, aes(sampleName,norm_logInt))
dataBox + geom_boxplot() + ylim(2,16) + ggtitle("after normalization")

Additionally, we have to make the box plot of the raw data again this time setting the Y-scale from 2 to 16.

dataBox = ggplot(logData,aes(sampleName,logInt))
dataBox + geom_boxplot() + ylim(2,16) + ggtitle("before normalization")

for (i in 1:6)
{
name = paste("MAplot",i,".jpg",sep="")
jpeg(name)
MAplot(data,which=i)
dev.off()
}

The which argument allows you to specify which array to compare with the median array. So if you have 6 arrays, you have to repeat this command 6 times, varying the value of the which argument from 1 to 6.

for (i in 1:6)
{
name = paste("MAplotnorm",i,".jpg",sep="")
jpeg(name)
MAplot(data.rma,which=i)
dev.off()
}

The which argument allows you to specify which array to compare with the median array. So in our example, we have 6 arrays, so we have to repeat this command 6 times, varying the value of the which argument from 1 to 6.

data.qc = qc(data)

The first quality measure are the average intensities of the background probes on each array.
You can retrieve them by using the avbg() method.
According to Affymetrix guidelines, the average background values of different arrays should be comparable.

avbg(data.qc)

The second quality measure are the scale factors: factors used to equalize the mean intensities of the arrays.
You can retrieve them by using sfs() method.
The scale factors should be within 3-fold of each other.

sfs(data.qc)

The third quality measure are the percent present calls: the percentage of spots that generate a significant signal (significantly higher than background) according to the Affymetrix detection algorithm.
You can retrieve them by using percent.present() method.
The precent present calls should be similar especially among replicates. Very low values (<20%) are a possible indicator of a poor quality array.

percent.present(data.qc)

The last quality measure are the 3’/5′ ratios of the quality control probe sets representing housekeeping genes (genes expressed in any tissue in any organism) like actin and GADPH. For these housekeeping genes 3 probe sets are available: one at the 5′ end of the gene, one at the middle and one at the 3′ end of the gene.
You can retrieve them by using ratios() method.
Affymetrix suggests that 3’/5′ ratios below 3 show acceptable RNA degradation and recommend caution if that value is exceeded for a given array.

ratios(data.qc)

data.rma = rma(data)
data.matrix = exprs(data.rma)

The ExpressionSet class is the superclass of the AffyBatch class meaning that AffyBatches are a special type of ExpressionSets. AffyBatches will therefore have the same characteristics and behaviour as ExpressionSets but AffyBatches will also have a set of specific characteristics and functions that are not shared by ExpressionSets.

library(“gcrma”)

To perform GCRMA normalization:

data.gcrma = gcrma(data)

As said before, GCRMA uses the affinity of each probe during background correction. These probe affinities are stored in an AffyBatch object, called affinity.info.
By default, affinity.info is not specified in the gcrma command: then affinities are computed by gcrma based on the data of the reference experiment.
However, the user can also choose to compute the affinities based on the data of their own experiment and use these affinities during normalization:

my.affinity.info = compute.affinities.local(data)
data.gcrma = gcrma(data,affinity.info=my.affinity.info)

The gcrma command comes with two additional arguments

• fast: the default value is FALSE. If fast is set to TRUE, specific procedures are used to speed up the background correction.
• type: has 3 values
  • fullmodel: uses probe sequences to calculate the affinities and combines the affinities with MM probe signals for background extraction
  • affinities: uses probe sequences to calculate affinities and ignores MM probes
  • MM: background correction based on MM intensities alone

So to speed up calculations use the following command:

data.gcrma = gcrma(data,fadt=TRUE)

bgcorr = pm(bg.correct(data,method="rma"))

First we need to get the data in the correct format for ggplot: a data frame with:

• raw log intensities in column 1
• background corrected log intensities in column 2
• sample names in column 3

We will work with PM intensities only

pmexp = pm(data)

We create three empty vectors that will serve as the three columns of the data frame:

• one to store the sample names in, called sampleNames
• one to store the log intensities in , called logs
• one to store the background corrected log intensities in, called corrlogs

The initialization of these vectors is essential. If you run into problems when creating these vectors and you have to repeat the plotting, always reinitialize these vectors (recreate them so that they are empty). 

The dim() method returns the number of rows (probes) and columns (samples) of the matrix containing the PM intensities. The number of rows corresponds to the number of PM intensities (probes) that is available for each sample (in this example that’s 251078). The sampleNames vector needs the same number of rows consisting of sample names.

The logs and the corrlogs vector will contain the log PM intensities of all six samples stacked into a single column. Up to this point we have extracted raw PM intensities, they are not yet log-transformed! You can do the log transformation using the log2() method.

sampleNames = vector()
logs = vector()
corrlogs = vector()
for (i in 1:6)
{
sampleNames = c(sampleNames,rep(ph@data[i,1],dim(pmexp)[1]))
logs = c(logs,log2(pmexp[,i]))
corrlogs = c(corrlogs,log2(bgcorr[,i]))
}

If you have 3 groups of 3 replicates, the code is as follows:

sampleNames = vector()
logs = vector()
corrlogs = vector()
for (i in 1:9)
{
sampleNames = c(sampleNames,rep(ph@data[i,1],dim(pmexp)[1]))
logs = c(logs,log2(pmexp[,i]))
corrlogs = c(corrlogs,log2(bgcorr[,i]))
}

Then we combine sample names and log intensities into one data frame:

corrData = data.frame(logInt=logs,bgcorr_logInt=corrlogs,sampleName=sampleNames)

Now we can create the plot. We use the geom_abline() method to add a red diagonal to the plot:

dataScatter = ggplot(corrData, aes(logInt,bgcorr_logInt))
dataScatter + geom_point() + geom_abline(intercept=0,slope=1,colour='red') + facet_grid(.~sampleName)

bgcorr = bg.adjust.gcrma(data)

The rest of the plot is created in exactly the same manner as described above for the RMA method.

bgcorr = pm(bg.correct(data,method="rma"))

First we need to get the data in the correct format for ggplot: a data frame with:

• raw log intensities in column 1
• background corrected log intensities in column 2
• sample names in column 3

We will work with PM intensities only

pmexp = pm(data)

We create three empty vectors that will serve as the three columns of the data frame:

• one to store the sample names in, called sampleNames
• one to store the log intensities in , called logs
• one to store the background corrected log intensities in, called corrlogs

The initialization of these vectors is essential. If you run into problems when creating these vectors and you have to repeat the plotting, always reinitialize these vectors (recreate them so that they are empty). 

The dim() method returns the number of rows (probes) and columns (samples) of the matrix containing the PM intensities. The number of rows corresponds to the number of PM intensities (probes) that is available for each sample (in this example that’s 251078). The sampleNames vector needs the same number of rows consisting of sample names.

The logs and the corrlogs vector will contain the log PM intensities of all six samples stacked into a single column. Up to this point we have extracted raw PM intensities, they are not yet log-transformed! You can do the log transformation using the log2() method.

sampleNames = vector()
logs = vector()
corrlogs = vector()
for (i in 1:6)
{
sampleNames = c(sampleNames,rep(ph@data[i,1],dim(pmexp)[1]))
logs = c(logs,log2(pmexp[,i]))
corrlogs = c(corrlogs,log2(bgcorr[,i]))
}

If you have 3 groups of 3 replicates, the code is as follows:

sampleNames = vector()
logs = vector()
corrlogs = vector()
for (i in 1:9)
{
sampleNames = c(sampleNames,rep(ph@data[i,1],dim(pmexp)[1]))
logs = c(logs,log2(pmexp[,i]))
corrlogs = c(corrlogs,log2(bgcorr[,i]))
}

Then we combine sample names and log intensities into one data frame:

corrData = data.frame(logInt=logs,bgcorr_logInt=corrlogs,sampleName=sampleNames)

Now we can create the plot. We use the geom_abline() method to add a red diagonal to the plot:

dataScatter = ggplot(corrData, aes(logInt,bgcorr_logInt))
dataScatter + geom_point() + geom_abline(intercept=0,slope=1,colour='red') + facet_grid(.~sampleName)

bgcorr = bg.adjust.gcrma(data)

The rest of the plot is created in exactly the same manner as described above for the RMA method.

data.matrix["256852_at",]

Raw intensities are stored in data, you can retrieve the raw PM intensities by using the pm() method. By using the probe set ID as a second argument, you can retrieve the PM intensities of the row with this name:

pm(data,"256852_at")

name = "boxplot.jpg"
jpeg(name)
boxplot(data,which='pm',col='red',names=ph@data$sample) 
dev.off()

The which argument allows you to specify if you want to use:

• perfect match probes only: which=’pm’
• mismatch probes only: which=’mm’
• both: which=’both’

Since the use of MM probes is highly controversial, we will work with PM probes only.
The col argument specifies the color of the boxes, the names argument the labels on the X-axis.

We’ll only use the PM intensities

pmexp = pm(data)

We create two empty vectors that will serve as the two columns of the data frame:

• one to store the sample names in, called sampleNames
• one to store the log intensities in , called logs

The initialization of these vectors is essential. If you run into problems when creating these vectors and you have to repeat the plotting, always reinitialize these vectors (recreate them so that they are empty) 

The dim() method returns the number of rows (probes) and columns (samples) of the matrix containing the PM intensities. The number of rows corresponds to the number of PM intensities (probes) that is available for each sample (in this example that’s 251078). The sampleNames vector needs the same number of rows consisting of sample names.

The logs vector will contain the log PM intensities of all six samples stacked into a single column. Up to this point we have extracted raw PM intensities, they are not yet log-transformed! You can do the log transformation using the log2() method.

sampleNames = vector()
logs = vector()
for (i in 1:6)
{
sampleNames = c(sampleNames,rep(ph@data[i,1],dim(pmexp)[1]))
logs = c(logs,log2(pmexp[,i]))
}

If you have 3 groups of 3 replicates, the code is as follows:

sampleNames = vector()
logs = vector()
for (i in 1:9)
{
sampleNames = c(sampleNames,rep(ph@data[i,1],dim(pmexp)[1]))
logs = c(logs,log2(pmexp[,i]))
}

After we have filled the vectors with the data we need, we combine sample names and log intensities into a single data frame:

logData = data.frame(logInt=logs,sampleName=sampleNames)

Now we can create the plot:

dataBox = ggplot(logData,aes(sampleName,logInt))
dataBox + geom_boxplot()

We create two vectors:

• one to store the sample names in, called sampleNames
• one to store the normalized log intensities in , called normlogs

The dim() method returns the number of rows (probe sets) and columns (samples) of the matrix containing the normalized intensities. The number of rows corresponds to the number of normalized intensities (probe sets) in each sample (in this example it’s 22810). The sampleNames vector needs the same number of entries of each sample name stacked into a single column.

The normlogs vector will contain the normalized intensities of all six samples stacked into a single column.

sampleNames = vector()
normlogs = vector()
for (i in 1:6)
{
sampleNames = c(sampleNames,rep(ph@data[i,1],dim(data.matrix)[1]))
normlogs = c(normlogs,data.matrix[,i])
}

If you have 3 groups of 3 samples, the code is as follows:

sampleNames = vector()
normlogs = vector()
for (i in 1:9)
{
sampleNames = c(sampleNames,rep(ph@data[i,1],dim(data.matrix)[1]))
normlogs = c(normlogs,data.matrix[,i])
}

After having filled the vectors with data, we combine sample names and normalized intensities into a single data frame:

normData = data.frame(norm_logInt=normlogs,sampleName=sampleNames)

Now we can create the plot. To compare this plot to the boxplot we made for raw log intensities we have to set them to the same Y-scale. The scale for the normalized data is much lower than these for the raw data so we have to make the box plot setting the Y-scale from 2 to 16. To set the Y-scale we can use the ylim() method.
To add a title to the plot we use the ggtitle() method.

dataBox = ggplot(normData, aes(sampleName,norm_logInt))
dataBox + geom_boxplot() + ylim(2,16) + ggtitle("after normalization")

Additionally, we have to make the box plot of the raw data again this time setting the Y-scale from 2 to 16.

dataBox = ggplot(logData,aes(sampleName,logInt))
dataBox + geom_boxplot() + ylim(2,16) + ggtitle("before normalization")

for (i in 1:6)
{
name = paste("MAplot",i,".jpg",sep="")
jpeg(name)
MAplot(data,which=i)
dev.off()
}

The which argument allows you to specify which array to compare with the median array. So if you have 6 arrays, you have to repeat this command 6 times, varying the value of the which argument from 1 to 6.

for (i in 1:6)
{
name = paste("MAplotnorm",i,".jpg",sep="")
jpeg(name)
MAplot(data.rma,which=i)
dev.off()
}

The which argument allows you to specify which array to compare with the median array. So in our example, we have 6 arrays, so we have to repeat this command 6 times, varying the value of the which argument from 1 to 6.

ph@data[ ,2] = c("control","control","control","mutant","mutant","mutant")
colnames(ph@data)[2]="source"
ph@data

Now you need to tell limma which sample belongs to which group. This info is contained in the second column (the column that we called source) of the PhenoData.

groups = ph@data$source

The names of the groups have to be transformed into factors. Factors in statistics represent the variable that you control: which sample belongs to which group. You can factorize the grouping variable by using the factor() method

f = factor(groups,levels=c("control","mutant"))

Then you need to create a design matrix, a matrix of values of the grouping variable. ANOVA needs such a matrix to know which samples belong to which group. Since limma performs an ANOVA or a t-test (which is just a specific case of ANOVA), it needs such a design matrix. You can create it using the model.matrix() method.
The argument of the model.matrix method is a model formula. The argument of the model.matrix method is a model formula. The tilde (~) in the argument specifies the right hand side of a model equation. If you define the design matrix with ~0, limma will simply calculate the mean expression level in each group.
Provide informative names for each column using the colnames() method.

design = model.matrix(~ 0 + f)
colnames(design) = c("control","mutant")

Essentially, limma will compare the mean expression level in the control samples to the mean expression level in mutant samples and this for each gene. So first of all, limma needs to calculate the mean expression levels using the lmFit() method. This method will fit a linear model (defined in design) to the data to calculate the mean expression level in the control and in the mutant samples:

data.fit = lmFit(data.matrix,design)

You can view the results of the fit. For instance you can show the estimates of the fit for the first 10 probe sets.

data.fit$coefficients[1:10,]

The first column contains the mean log expression in control samples.
The second column contains the mean log expression in mutant samples.

Now you have to tell limma which groups you want to compare. For this you define a contrast matrix defining the contrasts (comparisons) of interest by using the makeContrasts() method. Using the column names of the design matrix you specify the baseline sample (control) you want to compare to and the sample of interest (mutant).

contrast.matrix = makeContrasts(mutant-control,levels=design)
data.fit.con = contrasts.fit(data.fit,contrast.matrix)

Now limma is ready to perform the statistical test to compare the groups. Since we are comparing two groups, it will perform a t-test.
Limma does not perform a standard t-test but it will calculate a moderated t-statistic with shrunken standard deviation for each gene (see slides). An empirical Bayes method is used to shrink the variance of each gene towards a common value for all the genes. This is to lower the influence of very low or very high standard deviations on the t-test. Since the number of replicates is very low, the standard deviations will not be very reliable, ordinary t-statistics are not recommended. However, although there are only a few replications for each gene, the total number of measurements is very large. So information is combined across the genes (i.e., genome-wide shrinkage) to improve performance.
The moderated t-test is performed by using the eBayes() method.

data.fit.eb = eBayes(data.fit.con)

The eBayes() method returns a data frame data.fit.eb containing multiple slots. To get an overview of all the slots, use the names() method:

names(data.fit.eb)

You can retrieve the log fold changes of each gene via the coefficients slot:

data.fit.eb$coefficients[1:10,]

This slot contains the coefficient of the contrast: it’s simply the difference between the mean log expression in mutant samples and the mean log expression in control samples. This difference is usually called a log fold change.

The eBayes() method has performed a moderated t-test on each gene. That means it has calculated a t-statistic and a corresponding p-value for each gene. These values can be found in the t and the p.valuje slot. To view the t-statistics and p-values of the moderated t-test for the first 10 probe sets:

data.fit.eb$t[1:10,]
data.fit.eb$p.value[1:10,]

Remarkably, data.fit.eb also contains F-statistics and p-values for this F-statistic in the F and F.p.value slots while an F-statistic is in an ANOVA. Essentially, a t-test is a special case of an ANOVA used for single comparisons. If you have just a single comparison the F-statistic is the square of the t-statistic. For multiple comparisons, this relation is not true.

ph@data[ ,2] = c("control","control","control","drug","drug","drug","exercise","exercise","exercise")
colnames(ph@data)[2]="source"
ph@data

Since we have 3 groups with 3 replicates each, the factor that determines the grouping will have 3 levels instead of 2, so the code is as follows:

f = factor(groups,levels = c("control","drug","exercise"))

Then you need to create a design matrix, a matrix of values of the grouping variable. ANOVA needs such a matrix to know which samples belong to which group. Since limma performs an ANOVA, it needs such a design matrix. You can create it using the model.matrix() method.
The argument of the model.matrix method is a model formula. The tilde (~) in the argument specifies the right hand side of a model equation. If you define the design matrix with ~0, limma will simply calculate the mean expression level in each group.

design = model.matrix(~ 0 + f) 
colnames(design)=c("control","drug","exercise")
design
data.fit = lmFit(data.rma,design)

Afterwards you can tell limma which groups you want to compare. For this you define a contrast matrix defining the contrasts (comparisons) of interest by using the makeContrasts() method. In this example we make all pairwise comparisons:

contrast.matrix = makeContrasts(exercise-control,drug-exercise,drug-control,levels=design)
data.fit.con = contrasts.fit(data.fit,contrast.matrix)
data.fit.eb = eBayes(data.fit.con)

You can view the results of the ANOVA in the slots of the data.fit.eb object. The statistic that is calculated in ANOVA is the F-statistic, you can find the F-statistic and its corresponding p-value for each gene in the F and F.p.value slots.
ANOVA will just determine if there is a difference between the 3 groups but it will not tell you which groups differ: is it exercise and drug or drug and control ? You have no idea based on the results of the ANOVA. That’s why an ANOVA is always followed by a series of pairwise comparisons. The t-statistics and the resulting p-values of the pairwise comparisons are stored in the t and p.value slots.
Log fold changes can be found in the coefficients slot.

data.fit.eb$coefficients[1:10,]
data.fit.eb$t[1:10,]
data.fit.eb$p.value[1:10,]
data.fit.eb$F[1:10]
data.fit.eb$F.p.value[1:10]

Treatment is the grouping variable dividing the data set into two groups: before and after treatment. To tell limma that your data is paired you just create a second grouping variable called patient:

ph@data[ ,2] = c("before","treated","before","treated","before","treated")
colnames(ph@data)[2]="Treatment"
ph@data[ ,3] = c("patient1","patient1","patient2","patient2","patient3","patient3")
colnames(ph@data)[3]="Patient"
ph@data

Of course, now you need to factorize both grouping variables:

groupsP = ph@data$Patient 
groupsT = ph@data$Treatment
fp = factor(groupsP,levels=c("patient1","patient2","patient3"))
ft = factor(groupsT,levels=c("before","treated"))

Then you need to create a design matrix, a matrix of values of the grouping variable. ANOVA needs such a matrix to know which samples belong to which group. Since limma performs an ANOVA, it needs such a design matrix. You can create it using the model.matrix() method.
The argument of the model.matrix method is a model formula. The tilde (~) in the argument specifies the right hand side of the model equation. The plus (+) is used to combine factors.

paired.design = model.matrix(~ fp + ft)
colnames(paired.design)=c("Intercept","Patient2vsPatient1","Patient3versusPatient1","AftervsBefore")
data.fit = lmFit(data.rma,paired.design)
data.fit$coefficients

Lmfit() will fit a linear model to the data. It will create a data frame called data.fit. The interesting data is in the coefficients table which contains 4 columns:

• the first column contains the intercept of the linear fit, in most cases it has no implicit meaning.
• the second column compares patient 1 and patient 2: it is the difference between the average expression of patient 2 over the two treatments and the average expression level of patient 1 over the two treatments.
• the third column compares patient 1 and patient 3 in the same way.
• the fourth column compares before and after treatment: it calculates the average difference in expression level between after and before treatment. This is the measure that is used in the paired t-test and compared to 0. So these are the values we are interested in. The more this measure differs from 0 the more likely it is that the treatment has an effect.

Performing a moderated paired t-test is now done using eBayes().

data.fit.eb = eBayes(data.fit)
data.fit.eb$p.value[,1:20]

Again, the relevant p-values are in the fourth column, called Treatment. These are the p-values generated by the comparison of after treatment and before treatment.

Treatment is the grouping variable dividing the data set into two groups: before and after treatment. To tell limma that your data is paired you just create a second grouping variable called patient:

ph@data[ ,2] = c("before","during","treated","before","during","treated","before","during","treated")
colnames(ph@data)[2]="Treatment"
ph@data[ ,3] = c("P1","P1","P1","P2","P2","P2","P3","P3","P3")
colnames(ph@data)[3]="Patient"
ph@data

Of course, now you need to factorize both grouping variables:

groupsP = ph@data$Patient 
groupsT = ph@data$Treatment
fp = factor(groupsP,levels=c("P1","P2","P3"))
ft = factor(groupsT,levels=c("before","during","treated"))

Then you need to create a design matrix, a matrix of values of the grouping variable. ANOVA needs such a matrix to know which samples belong to which group. Since limma performs an ANOVA, it needs such a design matrix. You can create it using the model.matrix() method.
The argument of the model.matrix method is a model formula. The tilde (~) in the argument specifies the right hand side of a model equation. The plus (+) is used to combine factors.

paired.design = model.matrix(~ fp + ft)
colnames(paired.design)=c("Intercept","P2vsP1","P3vsP1","DuringvsBefore","AftervsBefore")
data.fit = lmFit(data.rma,paired.design)
data.fit$coefficients

Lmfit() will fit a linear model to the data. It will create a data frame called data.fit. The interesting data is in the coefficients table which contains 5 columns:

• the first column contains the intercept of the linear model. In most cases, it has no intrinsic meaning
• the second column compares patient1 and patient2: it is the difference between the average expression of patient2 over the 3 treatments and the average expression level of patient1 over the 3 treatments.
• the third column compares patient1 and patient3 in the same way.
• the fourth column compares during and before treatment: it calculates the average difference in expression level between during and before treatment. This is the measure that is used in the repeated measures ANOVA. So these are values we are interested in.
• the fifth column compares after and before treatment: it calculates the average difference in expression level between after and before treatment. This is the measure that is used in the repeated measures ANOVA. So these are values we are interested in.

Performing the repeated measures ANOVA is now done using eBayes().

data.fit.eb = eBayes(data.fit)
data.fit.eb$p.value[,1:20]

Again, the relevant p-values are in the fourth and fifth column, called DuringvsBefore and AftervsBefore. These are the p-values generated by the comparison of bfore, during and after treatment.

What if you also want to compare after and during treatment ?
Then you need to create a contrast matrix to tell limma which comparison to make. To do this you use the makeContrasts() method:

cont.matrix = makeContrasts(aftervsduring=after-during,levels=paired.design)
data.contr = contrasts.fit(data.fit,cont.matrix)
data.contr$coefficients
data.fit.eb = eBayes(data.contr)

The coef parameter specifies the column of data.fit.eb that should be used for the plot. In our example in Arabidopsis coef=2 since the second column of data.fit.eb contains the results of the comparison between mutant and control plants.
In our example of paired data data.fit.eb contains four columns:

• the second containing the results of the comparison between patient2 and patient1
• the third containing the results of the comparison between patient3 and patient1
• the fourth containing the results of the comparison between before and after the treatment

So in this example you set coef=2 or coef=3 if you want to find genes that are DE between patients or coef=4 if you want to find genes that are DE as a result of the treatment.

The highlight parameter allows to specify the number of highest scoring genes (on the Y-axis) for which names will be attached on the plot.

name = "Volcano.jpg"
jpeg(name)
volcanoplot(data.fit.eb,coef=2,highlight=10)
dev.off()

• The file argument specifies the file that you want to write to.
• The row.names argument specifies if row names are to be printed (the default is TRUE!).
• The col.names argument specifies if column names are to be printed (the default is TRUE!).
• The quote argument specifies if character or factor columns will be surrounded by double quotes (the default is TRUE!).

write.table(IDs.up,row.names=FALSE,col.names=FALSE,quote=FALSE,file="d:/trainingen/R/upIDs.txt")
write.table(IDs.down,row.names=FALSE,col.names=FALSE,quote=FALSE,file="d:/trainingen/R/downIDs.txt")

• normalized log intensities in the first column
• sample names in the second column
• probe set IDs in the third column

We first have to obtain the normalized intensities of the upregulated genes. The normalized intensities are stored in data.matrix. The downregulated genes are stored in topdowns. We can obtain their probe set IDs via the rownames() method. Then we select their normalized intensites from data.matrix using their probe set IDs:

data.matrix.up = data.matrix[(rownames(topdowns)),]

We create three vectors:

• one to store the sample names in, called sampleNames
• one to store the probe set IDs in, called featureNames
• one to store the normalized log intensities in, called heatlogs

The heatlogs vector will contain the normalized intensities of the upregulated genes in all six samples stacked into a single column.

There are 71 downregulated genes in my data set. I can see that by using the dim() method and looking at the first number it returns (dim(topdowns)[1]). This means that the sampleNames vector needs 254 entries of each sample name stacked into a single column.

The featureNames vector needs the probe set IDs of these genes repeated 6 times and stacked into a single column.

sampleNames = vector()
featureNames = vector()
heatlogs = vector()
for (i in 1:6)
{
sampleNames = c(sampleNames,rep(ph@data[i,1],dim(topdowns)[1]))
featureNames = c(featureNames,rownames(data.matrix.up[1:dim(topdowns)[1],]))
heatlogs = c(heatlogs,data.matrix.up[1:dim(topdowns)[1],i])
}

Then we combine sample names, probe set IDs and normalized intensities into one data frame:

heatData = data.frame(norm_logInt=heatlogs,sampleName=sampleNames,featureName=featureNames)

Now we can create the plot. A heat map can be created using the geom_tile() method. Low normalized intensities will be plotted in green and high normalized intensities will be plotted in red. This can be done by using the scale_fill_gradient() method.

dataHeat = ggplot(heatData, aes(sampleName,featureName))
dataHeat + geom_tile(aes(fill=norm_logInt)) + scale_fill_gradient(low="green", high="red")