The probe intensities stored in an AffyBatch object (i.e. mixture.batch) can be accessed by the function intensity(). Scatterplots of probe intensities represent an interesting view on your data, as shown below. However, modern Affymetrix chips contain several millions of probes. This makes traditional scatterplots very heavy to display or print and hard to overview. To overcome the problem, we can use a hexagonal binning plot. In brief, the plotting plane is broken into hexagons and the number of 'dots' that fall in each hexagon is counted. The hexagons are then colored according to the number of points they represent. The darker the hexagon, the more points it represents.

  h1 <- hexbin(intensity(mixture.batch)[, 1], intensity(mixture.batch)[,2])
  plot(h1, main = "Raw intensities")
  dev.new()
  h2 <- hexbin(log2(intensity(mixture.batch)[, 1]), log2(intensity(mixture.batch)[,2]))
  plot(h2, main = "Log transformed intensities")

A superposition of the density estimates for the intensities on each chip is also a very helpful plot:

  dev.new()
  par(mfrow=c(2,1))
  hist(mixture.batch,main="Raw intensities",log=FALSE)
  hist(mixture.batch,main="Log2 transformed intensities",log=TRUE)

2.6 Normalization (Scaling)

The general purpose of normalization is to adjust (or correct) a signal in order to make the comparison with other signals more meaningful. In the case of microarrays, we need to make the results from each of our arrays comparable with the rest. The normalization step is done at the probe level.

Qspline normalization

The normalization method Qspline is one of the normalization methods available in the affy package.

  mixture.qsp <- normalize(mixture.batch, method="qspline")

A plot of the distribution of the probe intensities is a useful diagnostic plot:

  dev.new()
  par(mfrow=c(2,1))
  hist(mixture.qsp, main = "Normalized by qspline")
  hist(mixture.batch, main="Un-normalized intensities")

2.7 Expression Index Calculation

Data analysts are usually interested in the expression values rather than individual probe values. The complete preprocessing of raw Affymetrix data into expression values can be done using background correction, normalization, perfect match correction and summary value computation. Several options exist for each step. To simplify things, we will skip the background correction here, and just focus on normalization and expression index calculation. One way to preprocess affymetrix chips is through the expresso() function, which performs both background correction and normalization as well as expression index calculation. It is slow, however, which is why we only use a small data set in this exercise.

  mixture.ei.es <- expresso(mixture.batch, bg.correct = F, normalize = F, pmcorrect.method = "pmonly", summary.method = "liwong")
  mixture.ei <- exprs(mixture.ei.es)

The expresso() function takes an object of class 'AffyBatch', but returns an object of class 'ExpressionSet'. Objects of class 'ExpressionSet' are easily converted to matrices using the function exprs().

You can easily see the difference after expression index calculation by comparing the number of data points for each chip before and after:

  length(intensity(mixture.batch[,1]))
  length(mixture.ei[,1])

See how the number of data points has decreased.

3 Effect of Normalization in Detail [Optional Exercise]

Because the GeneLogic data set is constructed around a dilution experiment (same samples used, but diluted to different concentrations), the expression levels should correlate with the degree of dilution as given by this vector:

  dilution <- c(1,1,0.75,0.75,0.5,0.5,0.25,0.25,0,0)

We can use this information to verify that normalization indeed does remove noise and not biologically relevant information. But first we need to calculate some expression indexes for qsplined data. For speed, the normalization was omitted before when you calculated the expression indexes above, but that is okay, for we need both normalized and un-normalized data here.

  mixture.qsp.es <- expresso(mixture.batch, bg.correct = F, normalize.method = "qspline", pmcorrect.method = "pmonly", summary.method = "liwong")
  mixture.qsp <- exprs(mixture.qsp.es)

Now calculate the correlations for the two data sets with the vector vec:

  unnorm.cor <- apply(mixture.ei,1, cor, y=dilution)
  norm.cor <- apply(mixture.qsp, 1, cor, y=dilution) 

Visualize the difference in correlation for the normalized and un-normalized data with vec in the following 2 histograms:

  dev.new()
  par(mfrow=c(1,2))
  hist(unnorm.cor, xlim=range(-1,1))
  hist(norm.cor, xlim=range(-1,1))

4 cDNA arrays [Optional Exercise]

cDNA microarrays differ much from Affymetrix arrays, and the strategies to normalize data are different. For this reason Bioconductor contains a package called 'marray' addressing cDNA arrays. As stated previously, however, using the 'marray' package to process cDNA arrays can be very cumbersome. Because of the lack of standardization, you would have to manually input information on slide layouts, choice of probes, type of equipment used, etc. (all the stuff given in the CDF files from Affymetrix, see 2.2 above). Indeed, the marray package is mostly useful for laboratories who produce thousands of cDNA arrays with the same setup. If, however, you want to use the qspline algorithm to normalize a cDNA array data set, you can do this by loading the cDNA data set into R as a matrix, and use the function normalize.qspline(). The only required argument would be the data matrix. Nevertheless, here is a small optional exercise showing how using marray would look

4.1 Load the package and data

The package marray addresses normalization for cDNA arrays. We will use the dataset swirl. More information about it can be found using the help system.

  library(marray)
  data(swirl)

4.2 Observing your arrays

Once again, observing the spatial distribution of your intensities can show interesting features (mainly artefacts)

  maImage(swirl)

The function maImage() has several options. You can see them using the help system, and eventually want to explore some of them.

4.3 Normalization

For cDNA array data the boxplot is often used.

  maBoxplot(swirl)

You can normalize cDNA array data like this:

  swirl.norm <- maNormMain(swirl)

You can observe the effect of normalization:

  maBoxplot(swirl.norm)

Here again, there are different strategies to normalize, some being very computer intensive. Help files for maNormMain() can tell you more.

5 Bibliography and links

1. Laurent Gautier, Leslie Cope, Benjamin M. Bolstad and Rafael A. Irizarry: affy -analysis of Affymetrix data at the probe level. Bioinformatics, 2003
2. Robert Gentleman and Ross Ihaka. R: A language for data analysis and graphics. Journal of Computational and Graphical Statistics, 5(3), 1996.
3. Steen Knudsen: A Biologistis Guide to Analysis of DNA Microarray Data. Wiley, New York, 2002.
4. C Workman, LJ Jensen, H Jarmer, R Berka, L Gautier, HB Nielsen, HH Saxild, C Nielsen, S Brunak, and S Knudsen. A new non-linear normalization method to reduce variability in dna microarray experiments. Genome Biology, 2002.
5. http://www.bioconductor.org/