Repository navigation
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
In the future, this page will hold links to the more useful slides that accompany this presentation
####Experiment Progression
Here are the steps for a basic Affymetrix microarray experiment:
1. Get Data
2. QC on Raw 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
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)
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)
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. Each row in datExpr will be a probe set after normalization.
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)
Sequencing run batches usually contribute the most to variance within an Affymetrix microarray dataset. It is important to remove the 'batch effects' before analyzing the dataset for phenotypic differences. First, we extract batch information (the sequencing run date) from our affy object.
batch = protocolData(data.affy)$ScanDate
batch = substr(batch,1,8)
batch = as.factor(batch)
table(batch)
datMeta$Batch = batch
plot(mds$points,col=as.numeric(datMeta$Batch),pch=19)
legend("topright",legend = levels(datMeta$Batch),fill = as.numeric(as.factor(levels(datMeta$Batch))))
[MDS Plot Post-Normalization by Batch Example](link to slide)
As can be seen in the MDS plot, the samples are clearly separating by batch at this stage.
We use ComBat from the sva R library to remove batch effects from our dataset. It uses either parametric or non-parametric empirical Bayes frameworks for adjusting data for batch effects. We need to remove all singular batches before moving on, since batch effects cannot be properly measured (and consequently removed) from a single batch. From our table(batch) command, we can see that we have one singular batch, "09/04/03".
to_remove = (datMeta$Batch == "09/04/03")
datExpr = datExpr[,!to_remove]
datMeta = datMeta[!to_remove,]
datMeta$Batch = droplevels(datMeta$Batch)
We can now remove batch contributions from our expression data.
mod = model.matrix(~datMeta$Dx)
batch = as.factor(datMeta$Batch)
datExpr.combat = ComBat(dat = datExpr,batch = batch,mod = mod)
datExpr = datExpr.combat
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))))
plot(mds$points,col=as.numeric(datMeta$Batch),pch=19)
legend("topright",legend = levels(datMeta$Batch),fill = as.numeric(as.factor(levels(datMeta$Batch))))
[MDS Plot Post-ComBat Dx Example](link to slide)
[MDS Plot Post-ComBat Batch Example](link to slide)
As can be seen in the MDS plots, ComBat effectively removes the variance due to batch from our dataset, which allows the variance due to phenotype (Control v. Parkinson's Disease) to come to the forefront.
Outliers can detract from 'true' phenotypic signal in a dataset. An outlier will commonly be a sample that doesn't have much variance in gene expression at all. We identify these outliers with network connectivity based statistics, and then remove these outliers from the data. This approach is especially effective when preparing data for a gene co-expression network experiment, which utilize tools such as Weighted Gene Co-Expression Network Analysis (WGCNA).
First, we plot do hierarchical clustering on our samples to visually identify potential outliers.
colnames(datExpr)=datMeta$geo_accession
tree = hclust(dist(t(datExpr)), method="average")
plot(tree)
[Sample Tree Example](link to slide)
Evidently, sample GSM506064 does not cluster with the other samples. Let's see if network statistics can confirm that this sample is an outlier.
normadj = (0.5 + 0.5*bicor(datExpr))^2
netsummary = fundamentalNetworkConcepts(normadj)
C = netsummary$Connectivity
Z.C = (C-mean(C))/sqrt(var(C))
to_keep = abs(Z.C) < 2
table(to_keep)
colnames(datExpr)[!to_keep]
Network analysis indicates that sample GSM506044 is an outlier in our dataset, which disagrees with the hierarchical clustering tree. This is OK though, since the connectivity measure is much more indicative of a sample that does not contribute to the group variance. The hierarchical clustering measure is more indicative of similarity in expression, which is not as important to a microarray experiment as variance across phenotypes.
So we choose to remove sample GSM506044 from our dataset.
datExpr = datExpr[,to_keep]
datMeta = datMeta[to_keep,]
We now perform the same three QC tests as before on our normalized, batch corrected, and outlier removed dataset.
dim(datExpr)
dim(datMeta)
#####(a) Boxplot
boxplot(datExpr,range=0, col = as.numeric(datMeta$Dx), xaxt='n', xlab = "Array", main = "Boxplot Normalized", ylab = "Intensity")
legend("topright",legend = levels(datMeta$Dx),fill = as.numeric(as.factor(levels(datMeta$Dx))))
[Boxplot Finalized Dataset Example](link to slide)
The medians/means of the samples should be about the same.
#####(b) Histogram
i=1
plot(density((datExpr[,i]),na.rm=T),col = as.numeric(datMeta$Dx[i]),main = "Hist of Log2 Exp", xlab="log2 exp")
for(i in 2:27){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 Finalized Dataset Example](link to slide)
The different sample lines should overlap much more strongly than before.
#####(c) MDS Plot
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 Plot Finalized Dataset Example](link to slide)
The control and Parkinson's disease samples should clearly separate from each other.
It’s important that all biological and technical covariates are not confounded by group/disease state. We don’t want the other factors (ie. sex, age) to be correlated with disease state (Dx). Thus, we want the p-value of the t-test (or F-test, if you have multiple disease states) of the coefficients predicting disease state from factors such as age, gender, batch, etc. to be greater than 0.05. If this is not the case for one of your factors, then you will need to take steps to remove confounders. This step of the analysis is discussed in more detail in Shared Features.
par(mfrow=c(2,2))
par(mar=c(3,2,3,2))
plot(datMeta$Dx, ylab="Number", main="Subjects")
for(i in c(6,7,9)){
if( i == 6 || i == 9 ){
print(paste(i,"Character Graph",sep=" "))
A = anova(lm(as.numeric(as.factor(datMeta[,i])) ~ datMeta$Dx)); p = A$"Pr(>F)"[1]
plot(as.factor(datMeta[,i]) ~ datMeta$Dx, main=paste(colnames(datMeta)[i]," p=", signif(p,2)), ylab="", xlab="")
}
else{
print(paste(i,"Number Graph",sep=" "))
A = anova(lm(as.numeric(datMeta[,i]) ~ datMeta$Dx)); p = A$"Pr(>F)"[1]
plot(as.numeric(as.character(datMeta[,i])) ~ datMeta$Dx, main=paste(colnames(datMeta)[i]," p=", signif(p,2)), ylab="", xlab="")
}
}
[Covariate Plot Example](link to slide)
Lucky for us, we do not have any obvious confounds in our dataset.