Repository navigation
Illumina
Welcome to the Illumina Tutorial. This tutorial will guide you through a basic Illumina experiment using freely available Illumina 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 Illumina microarray experiment
(3) Download Illumina.R to work through the tutorial directly on your own computer/server
(4) Use Illumina.R as a template for your own Illumina microarray experiment, and refer to the annotated code and tutorial slides as needed
This tutorial is a great place to start learning about Illumina microarray experiments. However, it was designed as a 'starting point' and not as an encyclopedia of all of the problems you might encounter in your experiment. After becoming comfortable with the basics discussed here, you should have the knowledge and skills to overcome your own roadblocks in experimentation. Please refer to the Shared Features page for a discussion of how to approach unforeseen problems that can (and will!) come up in your experiments.
This tutorial assumes you have a basic understand of programming and the use of the R language, including how to install packages from CRAN and bioconductor, and how to set up an appropriate development environment for R. This tutorial will make use of the latest packages available. Some features will briefly explained, but links will be provided to further documentation which you are encouraged to look over.
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. Collapse Rows
9. Differential Expression Analysis
This tutorial will use a dataset from GEO (GSE29378). This dataset is a from post-mortem hippocampal brain tissue collected from Alzheimer's patients and controls. Two hippocampal regions, CA1 and CA3, were collected. For the sake of simplicity, only CA3 will be examined in this tutorial.
We will begin by downloading the dataset using the GEOQuery library. The object returned by the getGEO() command from GEOQuery is an ExpressionSet object provided by the Biobase package, which is automatically loaded when the GEOQuery library is loaded. ExpressionSet objects are a very useful abstraction to store both expression data and metadata, and we will make extensive use of them throughout the tutorial. However, the ExpressionSet object returned by getGEO() contains the "final" expression vales used by the original study authors for their analysis. This means these values have already been processed through a previous pipeline. Please see this slide for an explanation of why you should rerun raw data in your new pipeline.
Fortunately GEO normally requires study authors to upload the original, unprocessed data as supplemental files. We will use the getGEOSuppFiles() function also provided by GEOQuery. This function will download the supplemental files as a .gz file, which will need to be extracted.
Before we run any commands, we need to load the libraries we will use for this analysis in a specific order, due to conflicting names between some packages.
library(limma)
library(lumiHumanIDMapping)
library(lumiHumanAll.db)
library(annotate)
library(lumi)
library(GEOquery)
library(ggplot2)
library(Cairo)
library(sva)
library(WGCNA)
library(reshape2)
library(plyr)
library(dplyr)
library(broom)
library(magrittr)
library(purrr)
library(stringr)
library(readr)We will learn more about many of these packages as we advance, but the packages in the bottom half must be loaded after the ones in the first half or you will encounter errors.
geo.object <- getGEO('GSE29378', destdir = '.')[[1]]
geo.supplemental <- getGEOSuppFiles('GSE29378', baseDir = './')When the downloads have completed (this can take some time depending the size of the file and the speed of your Internet connection), you need to extract the raw data from the supplemental files. You can use a GUI based program such as 7-zip, or you can use the following command in any Unix-based environment
gunzip GSE29378/GSE29378_non-normalized.txt.gz
tar xf GSE29378/GSE29378_RAW.tarThe raw expression data is located in GSE29378_non-normalized.txt.gz. The confusingly named GSE29378_RAW.tar does not infact contain any expression data but instead contains platform annotation information. This can be useful in the event that there are difficulties getting the correct annotations for the specific platform used in the study.
Now that we have the data downloaded, we need to extract the relevent subject metadata (such as diagnosis, age, sex, brain region, etc.) from the ExpressionSet object returned by getGEO(). We can get this information using the pData() function from Biobase, which will return a data frame where the rows correspond to each of the samples. If you examine this data frame, you will there are many columns, the large majority of which are irrelevant for our analysis. The title column actually contains most of the information we need, but it is joined together in one large string. We will split it up and then format each column appropriately, with the help of several packages. Note that the batch number does need to be exracted from a different column.
pdata.geo <- pData(geo.object) %>% select(title) %>% map(str_split_fixed, " - ", 5) %>% reduce(rbind) %>% data.frame
colnames(pdata.geo) <- c("Diagnosis", "Subject.ID", "Region", "Sex", "Age")
pdata.geo$Sample.ID <- sampleNames(geo.object)
pdata.geo$Age %<>% str_replace("yr", "") %>% as.character %>% as.numeric
pdata.geo$Subject.ID %<>% str_replace("#", "")
pdata.geo$Subject.ID <- str_replace(pdata.geo$Subject.ID, "#", "")
pdata.geo$Region %<>% str_replace("hippocampus ", "")
pdata.geo$Batch <- factor(geo.object$characteristics_ch1.10)We now have a data frame pdata.geo which contains the necessary metadata about our samples. Now we just need to get the raw data imported and we are ready to go! However, this will actually prove challenging for two reasons
- The column names of the raw data are not formatted correctly
- There are extra samples included in the raw data which are not in our metadata, and must be removed
The improperly formatted column names will cause the lumiR() function we would normally use to import the data to fail. Specifically, the detection p-value columns all have the same name, and the average signal intensity columns are missing 'AVG_SIGNAL' from their name. We will fix these problems by reading in the file as tab-separated values with read_tsv, reformatting the column names, and then writing it back as tab-separated values. Please see this slide for more information Illumina microarray columns.
geo.reformat <- read_tsv('./save/GSE29378/GSE29378_non-normalized.txt')
pval.columns <- str_detect(colnames(geo.reformat), "Detection Pval") %>% which
colnames(geo.reformat)[pval.columns] <- paste(colnames(geo.reformat)[pval.columns - 1], "Detection Pval", sep = ".")
colnames(geo.reformat)[pval.columns - 1] %<>% paste("AVG_SIGNAL", sep = ".")
write.table(geo.reformat, "./save/geo_reformat.tsv", row.names = FALSE, sep = '\t')We can now read the raw expression data in using lumiR(), which import them as an ExpressionSet.
lumi.raw <- lumiR('./save/geo_reformat.tsv', lib.mapping = 'lumiHumanIDMapping') The lumiR() function has read in our raw data, and by providing a mapping library, it will also convert the Illumina probe IDs to nuIDs, which are a more consistent way of determing which gene is actually measured by a given probe. Please see this slide for more information about Illumina probe annotation.
Unfortunately, this raw data contains more samples than the metadata we have in pdata.geo. We will need to remove these extra samples. The only way we can identify these samples is by parsing their column names, and matching them to corresponding metadata in pdata.geo. We will also add this metadata to our new ExpressionSet object.
pdata.lumi <- map(sampleNames(lumi.raw), str_split_fixed, "\\.", 7) %>% reduce(rbind) %>% data.frame
colnames(pdata.lumi) <- c("Diagnosis", "Subject.ID", "Region", "Sex", "Age", "X1", "X2")
pdata.lumi$sampleID <- lumi.raw$sampleID
pdata.lumi$Region %<>% str_replace("b", "")
pdata.lumi$Age %<>% as.character %>% as.integer
pdata.lumi$Diagnosis %<>% revalue(c(A = "Alzheimers", C = "Control"))
pData(lumi.raw) <- select(pdata.lumi, -(X1:X2))
sampleNames(lumi.raw) <- pdata.lumi$sampleIDWe are able to extract batch number, which will be enough for us to identify the samples to keep. In this case, we will make a unique identifier for each sample by combining the subject ID number with the brain region, separated by an _. We must also remove two technical replicates.
lumi.names <- paste(pdata.lumi$Subject.ID, pdata.lumi$Region, sep = "_")
geo.names <- paste(pdata.geo$Subject.ID, pdata.geo$Region, sep = "_")
known.names <- lumi.names %in% geo.names
rep.names <- str_detect(pdata.lumi$sampleID, "rep")
lumi.known <- lumi.raw[,known.names & !rep.names]
lumi.known$Batch <- str_replace(pdata.geo$Batch, "^.*: ", "") %>% factor %>% as.integer
sampleNames(lumi.known) <- paste(lumi.known$Subject.ID, lumi.known$Region, sep = "_")Now that we have the raw data for the correct set of samples, we will select only the ones from the CA3 region.
lumi.ca3 <- lumi.known[,lumi.known$Region == "CA3"]We will normalize our the intensity levels in our expression to account the large variability normally seen in gene expression. In this case, we can only use log2 normalization.
lumi.log2 <- lumiT(lumi.ca3, "log2")Raw Illumina microarray data which has been properly exported through GenomeStudio will contain more quality control columns which allows for the use of the more robust variance-stabilized transformation (VST) normalization method, which would be run as with "vst" substited for "log2" in the above command. Please see this slide for more information about normalization methods in Illumina arrays.
We will use three simple plots to examine the quality of our data:
- Boxplot
- Histogram
- MDS Plot
First we will prepare a data frame which will can be used in ggplot for making several of our figures.
expr.log2 <- exprs(lumi.log2) %>% t %>% data.frame(Sample.Name = sampleNames(lumi.log2), Batch = lumi.log2$Batch, Diagnosis = lumi.log2$Diagnosis)
expr.log2.melt <- melt(expr.log2, id = c("Sample.Name", "Batch", "Diagnosis"))
colnames(expr.log2.melt)[4:5] <- c("Symbol", "Intensity")Then we will make a boxplot of all of the probe intensities in every sample. The sample means should be roughly in line. This can be very slow and computationally intensive if you have a large number of samples and/or probes. Furthermore, we intensionally want to save this as a JPG (an image format with lossy compression) because a higher quality format will be very large and slow to load on most computers.
p <- ggplot(expr.log2.melt, aes(x = Sample.Name, y = Intensity, fill = factor(Batch))) + geom_boxplot() + theme_bw()
p <- p + theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust = 1.0))
p <- p + ggtitle("Log2 normalized signal intensity") #+ ylab("Intensity") + xlab("Sample")
p <- p + theme(legend.position = "none", panel.grid.major = element_blank(), panel.grid.minor = element_blank())
ggsave(filename = "boxplot_normalized.jpg", plot = p, width = 15 , height = 8)
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. We also expect the histogram to be right skewed, with most probes have fairly low expression and a small number with very high expression.
p <- ggplot(expr.log2.melt, aes(Intensity, group = Sample.Name, col = Diagnosis)) + geom_density() + theme_bw()
p <- p + theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())
p <- p + ggtitle("Histogram of Log2 Expression") + ylab("Density") + xlab("Log2 Expression")
CairoPDF("histogram_log2", height = 5, width = 9)
print(p)
dev.off()
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 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. However, in large samples clustering of any kind will be difficult to see and does not indicate any problem with the data.
mds.log2 <- exprs(lumi.log2) %>% t %>% dist %>% cmdscale(eig = TRUE)
mds.log2.plot <- data.frame(Sample.Name = rownames(mds.log2$points), Diagnosis = lumi.log2$Diagnosis, Batch = factor(lumi.log2$Batch), Component.1 = mds.log2$points[,1], Component.2 = mds.log2$points[,2])
p <- ggplot(mds.log2.plot, aes(x = Component.1, y = Component.2, col = Diagnosis)) + geom_point()
p <- p + theme_bw() + theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())
p <- p + xlab("Component 1") + ylab("Component 2") + ggtitle("MDS of Diagnosis")
CairoPDF(file = "mds_diagnosis", height = 6, width = 7)
print(p)
dev.off()
Since we have not found any problems with our data yet, we are ready to move on to the next steps. First we need to normalize between arrays to account for individual inter-array variability.
lumi.norm <- lumiN(lumi.log2, "rsn")The robust spline normalization (RSN) method has been used successfully with many Illumina microarray datasets. Please see this slide for more information about normalization methods.
Illumina microarrays provide an additional unique piece of information in the form of probe detection scores, which measure the confidence with which a given intensity was measured. This in the form of a p-value, provided for each probe in each sample. Probes in which no samples meet the the detection score cutoff are dropped because they could not be reliably detected in any sample. Please see this slide for more information about detection scores.
The recommended detection score cutoff is generally 0.01, the default setting used by the lumi function detectionCall(). This function returns the number of probes which pass the detection score threshold. We can then drop those probes from our ExpressionSet object.
lumi.cutoff <- detectionCall(lumi.norm) #Get the count of probes which passed the detection threshold per sample
lumi.expr <- lumi.norm[which(lumi.cutoff > 0),] #Drop any probe where none of the samples passed detection thresholdWe will also drop any probes which are unannotated. The gene symbols for each probe are retrieved using getSYMBOL from the annotate package.
symbols.lumi <- getSYMBOL(rownames(lumi.expr), 'lumiHumanAll.db') %>% is.na
lumi.expr.annot <- lumi.expr[!symbols.lumi,] Running microarray batches contribute substantial variance to many Illumina microarray experiments. We can examine an MDS plot of the batches to see if there is any clustering.
p <- ggplot(mds.log2.plot, aes(x = Component.1, y = Component.2, col = Batch)) + geom_point()
p <- p + theme_bw() + theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())
p <- p + xlab("Component 1") + ylab("Component 2") + ggtitle("MDS of Batch")
CairoPDF(file = "mds_batch", height = 6, width = 7)
print(p)
dev.off()
Although it is difficult to observe any clear clustering, we still need to remove any batch effects that are present. We will do this using ComBat() from the sva package. Note that if your data has any single batch samples, these must be dropped because the batch effects cannot be distinguished from other differences.
model.combat <- model.matrix(~ Diagnosis + Sex + Age, data = pData(lumi.expr.annot)) %>% data.frame
expr.combat <- ComBat(dat = exprs(lumi.expr.annot), batch = factor(lumi.expr.annot$Batch), mod = model.combat)
lumi.combat <- lumi.expr.annot
exprs(lumi.combat) <- expr.combat We can remake our MDS plots with the batch-effect corrected expression values.
mds.combat <- exprs(lumi.combat) %>% t %>% dist %>% cmdscale(eig = TRUE)
mds.combat.plot <- data.frame(Sample.Name = rownames(mds.combat$points), Diagnosis = lumi.combat$Diagnosis, Batch = factor(lumi.combat$Batch), Component.1 = mds.combat$points[,1], Component.2 = mds.combat$points[,2])
p <- ggplot(mds.combat.plot, aes(x = Component.1, y = Component.2, col = Diagnosis)) + geom_point()
p <- p + theme_bw() + theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())
p <- p + xlab("Component 1") + ylab("Component 2") + ggtitle("MDS of Diagnosis")
CairoPDF(file = "mds_diagnosis_combat", height = 7, width = 7)
print(p)
dev.off()
p <- ggplot(mds.combat.plot, aes(x = Component.1, y = Component.2, col = Batch)) + geom_point()
p <- p + theme_bw() + theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())
p <- p + xlab("Component 1") + ylab("Component 2") + ggtitle("MDS of Batch")
CairoPDF(file = "mds_batch_combat", height = 6, width = 7)
print(p)
dev.off()
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 hierarchical clustering and 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). In fact, we use the WGCNA R library to extract connectivity data and determine our outliers.
Hierarchical clustering can be used to visually identify outliers.
tree.combat <- exprs(lumi.combat) %>% t %>% dist %>% hclust(method = "average")
CairoPDF("clustering_combat", width = 13, height = 10)
plot(tree.combat, main = "Hierarchical Clustering Sammples")
dev.off()
We can see that 1294_CA3 is a clear outlier. We can now check if connectivity-based statistics also indicate this is an outlier, and if there are any other outliers.
normalized.adjacency <- (0.5 + 0.5 * bicor(exprs(lumi.combat))) ^ 2
network.summary <- fundamentalNetworkConcepts(normalized.adjacency)
connectivity <- network.summary$Connectivity
connectivity.zscore <- (connectivity - mean(connectivity)) / sqrt(var(connectivity))
connectivity.plot <- data.frame(Sample.Name = names(connectivity.zscore), Z.score = connectivity.zscore, Sample.Num = 1:length(connectivity.zscore))
p <- ggplot(connectivity.plot, aes(x = Sample.Num, y = Z.score, label = Sample.Name)) + geom_text(size = 4, colour = "red")
p <- p + geom_hline(aes(yintercept = -2))
p <- p + theme_bw() + theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())
p <- p + xlab("Sample Number") + ylab("Z score") + ggtitle("Sample Connectivity")
CairoPDF("connectivity_combat", width = 10, height = 10)
print(p)
dev.off()
We can then remove any samples in which the absolute value of the Z-score of the connectivity is greater than 2, and then remake our QC plots.
lumi.rmout <- lumi.combat[,abs(connectivity.zscore) < 2]expr.rmout <- exprs(lumi.rmout) %>% t %>% data.frame(Sample.Name = sampleNames(lumi.rmout), Batch = lumi.rmout$Batch, Diagnosis = lumi.rmout$Diagnosis)
expr.rmout.melt <- melt(expr.rmout, id = c("Sample.Name", "Batch", "Diagnosis"))
colnames(expr.rmout.melt)[4:5] <- c("Symbol", "Intensity")
p <- ggplot(expr.rmout.melt, aes(x = Sample.Name, y = Intensity, fill = factor(Batch))) + geom_boxplot() + theme_bw()
p <- p + theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust = 1.0))
p <- p + ggtitle("Log2 normalized signal intensity") + ylab("Intensity") + xlab("Sample")
p <- p + theme(legend.position = "none", panel.grid.major = element_blank(), panel.grid.minor = element_blank())
ggsave(filename = "boxplot_rmout.jpg", plot = p, family = "Oxygen", width = 15 , height = 8)
p <- ggplot(expr.rmout.melt, aes(Intensity, group = Sample.Name, col = Diagnosis)) + geom_density() + theme_bw()
p <- p + theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())
p <- p + ggtitle("Histogram of Log2 Expression") + ylab("Density") + xlab("Log2 Expression")
CairoPDF("histogram_rmout", height = 5, width = 9)
print(p)
dev.off()
mds.rmout <- exprs(lumi.rmout) %>% t %>% dist %>% cmdscale(eig = TRUE)
mds.rmout.plot <- data.frame(Sample.Name = rownames(mds.rmout$points), Diagnosis = lumi.rmout$Diagnosis, Batch = factor(lumi.rmout$Batch), Component.1 = mds.rmout$points[,1], Component.2 = mds.rmout$points[,2])
p <- ggplot(mds.rmout.plot, aes(x = Component.1, y = Component.2, col = Diagnosis)) + geom_point()
p <- p + theme_bw() + theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())
p <- p + xlab("Component 1") + ylab("Component 2") + ggtitle("MDS of Diagnosis")
CairoPDF(file = "mds_diagnosis_rmout", height = 6, width = 7)
print(p)
dev.off()
We don’t want the other factors (ie. sex, age) to be correlated with disease state (Diagnosis), since this will confound our analysis (we won't be able to tell the difference between covariate effects). Thus, we want the p-value of the t-test (or ANOVA F-test, if you have multiple disease states) of the coefficients from a linear model predicting factors such as age, gender, batch, etc. from disease state 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 .
lm.age <- lm(Age ~ Diagnosis, data = pData(lumi.rmout)) %>% anova %>% tidy
lm.sex <- lm(as.numeric(Sex) ~ Diagnosis, data = pData(lumi.rmout)) %>% anova %>% tidy
lm.batch <- lm(Batch ~ Diagnosis, data = pData(lumi.rmout)) %>% anova %>% tidy
p <- ggplot(pData(lumi.rmout), aes(x = Diagnosis, y = Age)) + geom_boxplot()
p <- p + theme_bw() + theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())
p <- p + theme(axis.title = element_blank()) + ggtitle(paste("Age (p < ", round(lm.age$p.value[1], 3), ")", sep = ""))
CairoPDF("age_boxplot", height = 6, width = 6)
plot(p)
dev.off()
p <- ggplot(pData(lumi.rmout), aes(x = Diagnosis, fill = Sex)) + geom_bar()
p <- p + theme_bw() + theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())
p <- p + theme(axis.title = element_blank(), legend.title = element_blank())
p <- p + ggtitle(paste("Sex (p < ", round(lm.sex$p.value[1], 3), ")", sep = ""))
CairoPDF("sex_barplot", height = 6, width = 6)
plot(p)
dev.off()
p <- ggplot(pData(lumi.rmout), aes(x = Diagnosis, fill = factor(Batch))) + geom_bar()
p <- p + theme_bw() + theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())
p <- p + theme(axis.title = element_blank(), legend.title = element_blank())
p <- p + ggtitle(paste("Batch (p < ", round(lm.batch$p.value[1], 3), ")", sep = ""))
CairoPDF("batch_barplot", height = 6, width = 6)
plot(p)
dev.off()
Many probes are duplicates for the same genes. We want to collapse our dataset down so that it we only have one probe per gene. This will simplify downstream analyses where duplicate genes may cause problems. The collapseRows() function from the WGCNA package can be used to do this very efficiently. You can choose how exactly to summarize probes, but the default setting works well in most cases. The default setting is to take the probe set with the maximum mean (across samples) as the final ensembl gene row for your datExpr.
expr.symbols <- featureNames(lumi.rmout) %>% getSYMBOL('lumiHumanAll.db') %>% factor
expr.collapse <- collapseRows(exprs(lumi.rmout), rowGroup = expr.symbols, rowID = rownames(exprs(lumi.rmout)))
lumi.collapse <- ExpressionSet(assayData = expr.collapse$datETcollapsed, phenoData = phenoData(lumi.rmout))