Affymetrix

##Introduction and Instructions

Welcome to the Affymetrix Tutorial. This tutorial will guide you through a basic Affymetrix experiment using freely available Affymetrix microarray data from Gene Expression Omnibus (GEO). Below is the recommended progression for advancing through this tutorial:

(1) Go through the Shared Features page to familiarize yourself with the methods shared across all microarray experiments
(2) Return to this page to go through the annotated code along with accompanying tutorial slides to become acquainted with a typical Affymetrix microarray experiment
(3) Download Affymetrix.R to work through the tutorial directly on your own computer/server
(4) Use Affymetrix.R as a template for your own Affymetrix microarray experiment, and refer to the annotated code and tutorial slides as needed

---

Tutorial Slides

In the future, this page will hold links to the more useful slides that accompany this presentation

---

Annotated Code

####Experiment Progression

Here are the steps for a basic Affymetrix microarray experiment:

1. Get Data
2. QC on Non-Normalized Data
3. Normalization
4. Batch Correction
5. Outlier Removal
6. QC on Normalized Data
7. Covariate Analysis
8. Annotating Probes
9. Collapse Rows
10. Differential Expression Analysis

1. Get Data

For this tutorial we will be using data from GEO, experiment GSE20295. This is a Parkinson's disease multi-region brain study done with human post-mortem brain tissue. In the tutorial, we will only be looking at the prefrontal cortex.

It is key to download the affy R library for your Affymetrix experiment, since it contains functions for reading the raw data into your R session (ReadAffy) and normalization of Affymetrix arrays (rma).

For your Affymetrix experiment, you have two choices for getting your data into your R session:

(a) You load your data directly from a source (such as GEO or ArrayExpress) into your R session
(b) You download your data to your local computer or server before loading it into your R session

Option (a) can fully utilize tools such as the GEOquery R library. Specific sources often have R libraries for downloading microarray data directly into your R session. In the code tutorial, we use GEOquery for option (a).

Option (b) is more suitable if you are using data from your own experiment, OR you need to select specific parts of an entire dataset to analyze. Sometimes with large datasets (> 500 Mb) it is better to 'pre-filter', or only extract/download data that you need - in this experiment, we only want prefrontal cortex.

You only need to use option (a) OR option (b) to get data into your R session.

#####(a) Load your data directly from a source into your R session (using GEOquery R library)

Use getGEO() to get phenotype data from the entire experiment, and create your datMeta.

gse <- getGEO("GSE20295", GSEMatrix =TRUE,getGPL=FALSE) datMeta = pData(gse[[1]])
rownames(datMeta) = datMeta[,2]
datMeta$title = gsub(" ","_",datMeta$title)
idx = which(datMeta$source_name_ch1 == "Postmortem brain prefrontal cortex")

Create an index idx to later use in keeping only prefrontal cortex samples.
Next, use getGEOSuppFiles() to directly download .CEL.gz files onto your computer from the whole experiment.

getGEOSuppFiles("GSE20295")

Use software such as 7-zip to extract your new "GSE20295_RAW" folder. Also, take care to manually change GSM506039_1_1364_BA9_Pm.CEL.gz to GSM506039_1364_BA9_Pm.CEL.gz before moving on. This is an upload error in GEO - this kind of mistake definitely happens all the time.

filesPFC = paste(datMeta$geo_accession,"_",datMeta$title,".CEL.gz",sep="")[idx]
data.affy = ReadAffy(celfile.path = "./GSE20295/GSE20295_RAW", filenames = filesPFC)
datExpr = exprs(data.affy)
datMeta = datMeta[idx,]

ReadAffy() created an affy object data.affy that only contains the prefrontal cortex samples that we are interested in. It is important for your affy object to only contain the samples you are investigating, since the normalization function rma() needs this affy object as input.

Finally, check datMeta/datExpr ordering and reformat as necessary.

GSM = rownames(pData(data.affy))
GSM = substr(GSM,1,9)
idx = match(GSM, datMeta$geo_accession)
datMeta = datMeta[idx,]

datMeta = datMeta[,-c(3:7,14:36)]
datMeta$characteristics_ch1 = gsub("disease state: control","CTL",datMeta$characteristics_ch1)
datMeta$characteristics_ch1 = gsub("disease state: Parkinson's disease","PKD",datMeta$characteristics_ch1)
datMeta$characteristics_ch1.1 = gsub("gender: male","M",datMeta$characteristics_ch1.1)
datMeta$characteristics_ch1.1 = gsub("gender: female","F",datMeta$characteristics_ch1.1)
datMeta$characteristics_ch1.2 = gsub("age: ","",datMeta$characteristics_ch1.2)
datMeta$characteristics_ch1.3 = gsub("brain region: ","",datMeta$characteristics_ch1.3)
colnames(datMeta)[5:8] = c("Dx","Sex","Age","Region")
datMeta$Dx = as.factor(datMeta$Dx)

Your datMeta and datExpr should look like this:[datMeta and datExpr Example](datExpr/datMeta example slide link)

#####(b) Download your data to your local computer or server before loading it into your R session

Use getGEO() to get phenotype data from the entire experiment, and create your datMeta.

gse <- getGEO("GSE20295", GSEMatrix =TRUE,getGPL=FALSE) datMeta = pData(gse[[1]])
rownames(datMeta) = datMeta[,2]
datMeta$title = gsub(" ","_",datMeta$title)
idx = which(datMeta$source_name_ch1 == "Postmortem brain prefrontal cortex")

Create an index idx to later use in keeping only prefrontal cortex samples.

You can download the prefrontal cortex expression data directly from the GEO page GSE20295. For GSE20295_RAW.tar choose 'custom' download, and only select files with 'BA9' and '.CEL.gz' in the name. Download these files into a new directory 'GSE20295_PFC'. Use software such as 7-zip to extract your new "GSE20295_PFC/GSE20295_RAW" folder. Then, take care to manually change GSM506039_1_1364_BA9_Pm.CEL.gz to GSM506039_1364_BA9_Pm.CEL.gz before moving on. This is an upload error in GEO - this kind of mistake definitely happens all the time.

Next, simply use ReadAffy() to read in your pre-selected prefrontal cortex samples from your directory.

data.affy = ReadAffy(celfile.path = "./GSE20295_PFC/GSE20295_RAW")
datExpr = exprs(data.affy)
datMeta = datMeta[idx,]

ReadAffy() created an affy object data.affy that only contains the prefrontal cortex samples that we are interested in. It is important for your affy object to only contain the samples you are investigating, since the normalization function rma() needs this affy object as input.

Finally, check datMeta/datExpr ordering and reformat as necessary.

GSM = rownames(pData(data.affy))
GSM = substr(GSM,1,9)
idx = match(GSM, datMeta$geo_accession)
datMeta = datMeta[idx,]

datMeta = datMeta[,-c(3:7,14:36)]
datMeta$characteristics_ch1 = gsub("disease state: control","CTL",datMeta$characteristics_ch1)
datMeta$characteristics_ch1 = gsub("disease state: Parkinson's disease","PKD",datMeta$characteristics_ch1)
datMeta$characteristics_ch1.1 = gsub("gender: male","M",datMeta$characteristics_ch1.1)
datMeta$characteristics_ch1.1 = gsub("gender: female","F",datMeta$characteristics_ch1.1)
datMeta$characteristics_ch1.2 = gsub("age: ","",datMeta$characteristics_ch1.2)
datMeta$characteristics_ch1.3 = gsub("brain region: ","",datMeta$characteristics_ch1.3)
colnames(datMeta)[5:8] = c("Dx","Sex","Age","Region")
datMeta$Dx = as.factor(datMeta$Dx)

Your datMeta and datExpr should look like this:[datMeta and datExpr Example](datExpr/datMeta example slide link)

2. QC on Non-Normalized Data

Before normalizing our data, we want to take a look at the raw data to make sure that it looks OK. We do Quality Control (QC) tests to do this. First, we log2 normalize our data just for these tests (and check out the dimensions):

datExpr = log2(datExpr) dim(datExpr)

There are three simple and effective tests to do:

(a) Boxplot
(b) Histogram
(c) MDS Plot

#####(a) Boxplot

A boxplot gives us an idea of how much gene expression each sample has as a whole. We want to make sure that the means of the samples are generally in a line.

boxplot(datExpr,range=0, col = as.numeric(datMeta$Dx), xaxt='n', xlab = "Array", main = "Boxplot Pre-Normalization", ylab = "Intensity")
legend("topright",legend = levels(datMeta$Dx),fill = as.numeric(as.factor(levels(datMeta$Dx))))

[Boxplot Example Pre-Normalization](link to boxplot slide)

#####(b) Histogram

A histogram shows us the distribution of amount of expression for each probe within a sample. The x-axis is the amount of expression, and the y-axis is what percentage of probes had this x-value of expression. Each line is a different sample. We want to make sure that the histograms of the samples are generally overlapping.

i=1
plot(density((datExpr[,i]),na.rm=T),col = as.numeric(datMeta$Dx[i]),main = "Hist of Log2 Exp", xlab="log2exp",xlim=c(4,16),ylim=c(0,0.5))
for(i in 2:29){lines(density((datExpr[,i]),na.rm=T), col = as.numeric(datMeta$Dx)[i],)}
legend("topright",legend = levels(datMeta$Dx),fill = as.numeric(as.factor(levels(datMeta$Dx))))

[Histogram Example Pre-Normalization](link to histogram slide)

#####(c) MDS Plot

An MDS (Multi-Dimensional Scaling) plot uses principal coordinate analysis to find the two vectors that capture the most amount of variance in a dataset. It is a visual representation of how similar samples are in their expression. We just want to make sure that there isn't one sample way off by itself away from all of the other samples. Usually samples will cluster based on batch at this stage, which we verify later on in the analysis.

mds = cmdscale(dist(t(datExpr)),eig=TRUE)
plot(mds$points,col=as.numeric(datMeta$Dx),pch=19)
legend("topright",legend = levels(datMeta$Dx),fill = as.numeric(as.factor(levels(datMeta$Dx))))

[MDS Example Pre-Normalization](link to mds slide)

3. Normalization

Normalization helps to remove experimental error from microarray experiments. RMA (robust multichip average) normalization is a tool designed specifically for Affymetrix arrays. It performs background correction, probe summarization, quantile normalization and log2 transformation. Background correction eliminates microarray 'noise', ie. removes the probes that bind to nothing but emit signal. Probe summarization takes the expression from each probe and summarizes it to its 'probe set' using median polishing.

datExpr = rma(data.affy, background=T, normalize=T, verbose=T)
datExpr = exprs(datExpr)

[Boxplot Post-Normalization Example](link to slide) [Histogram Post-Normalization Example](link to slide) [MDS Plot Post-Normalization Example](link to slide)

4. Batch Correction

5. Outlier Removal

6. QC on Normalized Data

7. Covariate Analysis

8. Annotating Probes

9. Collapse Rows

10. Differential Expression Analysis