EDA for Microarray data

Exploratory data analysis (EDA) can be used to get an idea of what your data looks like. Here we are analyzing microarray data from 8 samples: two groups of four. A first step in any analysis of genomics data is to learn its general properties and search for problematic samples.

# library(devtools)
# install_github("dagdata","genomicsclass")
library(dagdata)
#load data
data(SpikeInEDA)
#histogramm of the first array: not very informative in linear scale
hist(int[,1])
#Better picture using log2 scale
hist(log2(int[,1]))

for(i in 1:ncol(int))
  if(i==1) plot(density(log2(int[,i])),col=(i==4)+1, main=NULL) else lines(density(log2(int[,i])),col=(i==4)+1)

Note that one histogram (we highlighted it by making it red) looks different: it has a different shape from the rest. Is this sample different from the rest in any significant way? If we compute the correlation between this sample and the rest it is not very different and all very high.

signif(cor(int),2)

##                 2353m99hpp_av08 2353n99hpp_av08 2353o99hpp_av08
## 2353m99hpp_av08            1.00            0.99            0.99
## 2353n99hpp_av08            0.99            1.00            0.99
## 2353o99hpp_av08            0.99            0.99            1.00
## 2353p99hpp_av08            0.99            0.99            0.99
## 2353q99hpp_av08            0.98            0.98            0.98
## 2353r99hpp_av08            0.98            0.98            0.98
## 2353s99hpp_av08            0.98            0.98            0.98
## 2353t99hpp_av08            0.98            0.98            0.98
##                 2353p99hpp_av08 2353q99hpp_av08 2353r99hpp_av08
## 2353m99hpp_av08            0.99            0.98            0.98
## 2353n99hpp_av08            0.99            0.98            0.98
## 2353o99hpp_av08            0.99            0.98            0.98
## 2353p99hpp_av08            1.00            0.98            0.98
## 2353q99hpp_av08            0.98            1.00            0.99
## 2353r99hpp_av08            0.98            0.99            1.00
## 2353s99hpp_av08            0.98            0.99            0.99
## 2353t99hpp_av08            0.98            0.99            0.99
##                 2353s99hpp_av08 2353t99hpp_av08
## 2353m99hpp_av08            0.98            0.98
## 2353n99hpp_av08            0.98            0.98
## 2353o99hpp_av08            0.98            0.98
## 2353p99hpp_av08            0.98            0.98
## 2353q99hpp_av08            0.99            0.99
## 2353r99hpp_av08            0.99            0.99
## 2353s99hpp_av08            1.00            0.99
## 2353t99hpp_av08            0.99            1.00

##we don't need to show all the points so we take random samples
#of size 10000
library(rafalib)
splot<-function(x,y,...){
  ind<-sample(length(x),10000)
  x=x[ind];y=y[ind] 
  plot(x,y,...)
}  
#scatterplot between sample 1 and 2
splot(log2(int[,1]),log2(int[,2]), 
      main="sample 1 and 2") 
#Scatterplot between sample 1 and 4
splot(log2(int[,1]),log2(int[,4]),
      main="sample 1 and 4")

Note that samples 1 through 4 are replicates and should produce the same values up to measurement error. Scatterplots and correlation are not the best tools to detect problems. For example 1,2,3,4 and 100,200,300,400 two lists with very different values have perfect correlation. A better measure is the differences between the values and therefore a better plot is a rotation of the scatter plot containing the differences (log ratios) on the y-axis and the averages (in the log scale) on the x-axis. This plot is a refereed to as an MA-plot (plot of difference versus the average). 

maplot<- function(x,y,...) splot((x+y)/2,y-x,...)

maplot(log2(int[,1]),log2(int[,2]),xlab="A",ylab="M", main="Sple 1 and 2", ylim=c(-2,2))
maplot(log2(int[,1]),log2(int[,3]),xlab="A",ylab="M",main="Sple 1 and 3",ylim=c(-2,2))
maplot(log2(int[,1]),log2(int[,4]),xlab="A",ylab="M", main="Sple 1 and 4", ylim=c(-2,2))

Now the problem is obvious. It turns out this samples comes from an array for which a spatial problem can be detected at the original image level. We actually have the grid locations for these measurements and can recreate the image.

##we are doing this for two arrays 1 and 4
library(matrixStats) ##need rowMedians
library(RColorBrewer)
for(i in c(1,4)){
  r=log2(int[,i])-rowMedians(log2(int)) 
  ## r are residuals from median array
  ## to avoind outliers taking over colors of image
  ### define a MAX
  MAX<-1
  r[r>MAX]<-MAX
  r[r< -MAX] <- -MAX
  ##we now that every other column is skipped
  mat=matrix(NA,max(locations[,1]),max(locations[,2]+1)/2)
  for(j in 1:nrow(locations)){
     mat[locations[j,1],(locations[j,2]+1)/2]<-r[j]
  }
  image(mat,col=brewer.pal(11,"RdBu"))
}

On the second image we can clearly see the spatial pattern (blue are positive residuals, red are negative)

Licence

Licence

References

https://github.com/genomicsclass

Download raw data

The raw data files for this lab are in the rawdata repository, available here:
https://github.com/genomicsclass/rawdata.
Click Download ZIP in order to download all the files, unzip this file which should result in a rawdata-master folder. Rename this folder to rawdata.

Read Agilent data

Agilent data is a two color arrays. For two-color arrays it’s slightly more complicated, because you have a pairing of files (red and green channels). Different scanners spit out different formats for this. So the target information will be a little bit different. limma package is used to read the files. The readTargets function tells you what’s in the red and the green channels.

The function read.maimages(microarray images) is used to read the data. You have to tell it what software produced the images. There’s different software for that, although these days, there’s a default that works pretty broadly. So all we’re doing now is telling it where the files are. We tell it what imaging software was used to produce these files (“genepix”). The function read.maimages, stores red and green separately.

library(limma)
library(rafalib)#installed from github
#Define a base directory
basedir <- "~/hubiC/Documents/R/doc/english/genomics/rawdata/agilent"
setwd(basedir)
#Read sample information table : it tells you what&#39;s in the red channel and the green channel
targets <- readTargets("TargetBeta7.txt")
#Read microarray files : Red and green are stored separately
RG <- read.maimages(targets$FileName, source="genepix")

## Read 6Hs.195.1.gpr 
## Read 6Hs.168.gpr 
## Read 6Hs.166.gpr 
## Read 6Hs.187.1.gpr 
## Read 6Hs.194.gpr 
## Read 6Hs.243.1.gpr

Normalization

Data are normalized using MA.RG function. It stores the information as the log ratio– that’s the M– and the average of the logs, which is the A.

#Normalization
MA <- MA.RG(RG,bc.method="none")

MA plot

Just to give you a quick idea of what you can do with this. You can do an MA plot. You just plot the A versus the M.

#Normalization
MA <- MA.RG(RG,bc.method="none")
#MA plot
plot(MA$A[,1],MA$M[,1])#MA plot for the first sample

Array image

Another nice feature is to make images pretty quickly.

#Array image
imageplot(MA$M[,2], RG$printer, zlim=c(-3,3))

This is a nice example where we see an image where there seems to be a little bit of a problem. If you look at the edge, you see that it's quite green, compared to the middle, which is redder. Then you have too much red on the left. So we're detecting that there's somewhat of a problem in this file.

Licence

Licence

References

https://github.com/genomicsclass

Download raw data

The raw data files for this lab are in the rawdata repository, available here: https://github.com/genomicsclass/rawdata Click Download ZIP in order to download all the files, unzip this file which should result in a rawdata-master folder. Rename this folder to rawdata.

Read Affymeterix CEL files

We start by reading in the sample information table. This is usually created by the person who performed the experiment.

#Set working directory to the the celfiles
basedir <- "~/hubiC/Documents/R/doc/english/genomics/rawdata/celfiles"
setwd(basedir)
library(affy)

#Sample information table : it has file names and data from spiking experiment (spiking concentration)
tab <- read.delim("sampleinfo.txt",check.names=FALSE,as.is=TRUE)
rownames(tab) <- tab$filenames
tab[1:6, 1:3]

##                                       filenames 37777_at 684_at
## 1521a99hpp_av06.CEL.gz   1521a99hpp_av06.CEL.gz     0.00   0.25
## 1532a99hpp_av04.CEL.gz   1532a99hpp_av04.CEL.gz     0.00   0.25
## 2353a99hpp_av08.CEL.gz   2353a99hpp_av08.CEL.gz     0.00   0.25
## 1521b99hpp_av06.CEL.gz   1521b99hpp_av06.CEL.gz     0.25   0.50
## 1532b99hpp_av04.CEL.gz   1532b99hpp_av04.CEL.gz     0.25   0.50
## 2353b99hpp_av08r.CEL.gz 2353b99hpp_av08r.CEL.gz     0.25   0.50

#list all the .cel files that are in the current directory.
fns <- list.celfiles()
fns

## [1] "1521a99hpp_av06.CEL.gz"  "1521b99hpp_av06.CEL.gz" 
## [3] "1532a99hpp_av04.CEL.gz"  "1532b99hpp_av04.CEL.gz" 
## [5] "2353a99hpp_av08.CEL.gz"  "2353b99hpp_av08r.CEL.gz"

#Check whether the filenames are the same in the directory and in the sample info tab
fns %in% tab[,1] ##check

## [1] TRUE TRUE TRUE TRUE TRUE TRUE

#Read cel files in the current directory
ab <- ReadAffy(phenoData=tab)

ReadAffy function creates an AffyBatch object which object contains the information you need.

#Extract the perfect match probe-level intensities
dim(pm(ab))

## [1] 201807      6

#Phenotypic data : sample information 6X17
dim(pData(ab))

## [1]  6 17

#Plateform used for gene information
annotation(ab)

## [1] "hgu95a"

Normalization

The last thing to do here is to turn probe-level information into gene-level information. You can preprocess this probe-level information in many ways. One way you can do it is using this algorithm called rma. It’s going to turn probe-level data into gene-level data, quantile normalization and also background correction.

e <- rma(ab)

## Background correcting
## Normalizing
## Calculating Expression

dim(e)

## Features  Samples 
##    12626        6

You notice that this information is smaller, it's only 12,000 by 6. It's been summarized now.

If you are not interested in probe level data you could use this function :

setwd(basedir)
ejust <- justRMA(filenames=tab[,1],phenoData=tab)
dim(ejust)

Licence

Licence

References

https://github.com/genomicsclass