Microarray Data Analysis

Load Packages

suppressPackageStartupMessages({
    library(here)
    library(oligo)
})

Load the dataset

• GSM651311 Keratinocytes, untreated 24h, rep1
• GSM651316 Keratinocytes, DEX-treated 24h, rep1
• GSM651321 Keratinocytes, untreated 24h, rep2
• GSM651326 Keratinocytes, DEX-treated 24h, rep2

dataset <- read.celfiles(list.files(here("data/microarray"), full.names = TRUE))

## Loading required package: pd.hg.u95av2

## Loading required package: RSQLite

## Loading required package: DBI

## Platform design info loaded.

## Reading in : /home/jialin/Courses/code6150/data/microarray/GSM651309.CEL.gz
## Reading in : /home/jialin/Courses/code6150/data/microarray/GSM651310.CEL.gz
## Reading in : /home/jialin/Courses/code6150/data/microarray/GSM651311.CEL.gz
## Reading in : /home/jialin/Courses/code6150/data/microarray/GSM651312.CEL.gz
## Reading in : /home/jialin/Courses/code6150/data/microarray/GSM651313.CEL.gz
## Reading in : /home/jialin/Courses/code6150/data/microarray/GSM651314.CEL.gz
## Reading in : /home/jialin/Courses/code6150/data/microarray/GSM651315.CEL.gz
## Reading in : /home/jialin/Courses/code6150/data/microarray/GSM651316.CEL.gz
## Reading in : /home/jialin/Courses/code6150/data/microarray/GSM651317.CEL.gz
## Reading in : /home/jialin/Courses/code6150/data/microarray/GSM651318.CEL.gz
## Reading in : /home/jialin/Courses/code6150/data/microarray/GSM651319.CEL.gz
## Reading in : /home/jialin/Courses/code6150/data/microarray/GSM651320.CEL.gz
## Reading in : /home/jialin/Courses/code6150/data/microarray/GSM651321.CEL.gz
## Reading in : /home/jialin/Courses/code6150/data/microarray/GSM651322.CEL.gz
## Reading in : /home/jialin/Courses/code6150/data/microarray/GSM651323.CEL.gz
## Reading in : /home/jialin/Courses/code6150/data/microarray/GSM651324.CEL.gz
## Reading in : /home/jialin/Courses/code6150/data/microarray/GSM651325.CEL.gz
## Reading in : /home/jialin/Courses/code6150/data/microarray/GSM651326.CEL.gz
## Reading in : /home/jialin/Courses/code6150/data/microarray/GSM651327.CEL.gz
## Reading in : /home/jialin/Courses/code6150/data/microarray/GSM651328.CEL.gz

#dataset <- dataset[, c(paste0(c("GSM651310", "GSM651320", "GSM651315", "GSM651325"), ".CEL.gz"))]
dataset <- dataset[, c(paste0(c("GSM651311", "GSM651321", "GSM651316", "GSM651326"), ".CEL.gz"))]
sampleNames(dataset) <- c("control_rep1", "control_rep2", "treat_rep1", "treat_rep2")
pData(dataset)$group <- c("control", "control", "treat", "treat")
dataset

## ExpressionFeatureSet (storageMode: lockedEnvironment)
## assayData: 409600 features, 4 samples 
##   element names: exprs 
## protocolData
##   rowNames: control_rep1 control_rep2 treat_rep1 treat_rep2
##   varLabels: exprs dates
##   varMetadata: labelDescription channel
## phenoData
##   rowNames: control_rep1 control_rep2 treat_rep1 treat_rep2
##   varLabels: index group
##   varMetadata: labelDescription channel
## featureData: none
## experimentData: use 'experimentData(object)'
## Annotation: pd.hg.u95av2

Quality control

boxplot(dataset, target = "core")

oligo::MAplot(dataset, pairs = TRUE, ylim = c(-1, 1))

library(arrayQualityMetrics)
arrayQualityMetrics(expressionset = dataset,
    outdir = "docs/microarray_qualitymetrics",
    force = TRUE, do.logtransform = TRUE,
    intgroup = c("group"))

RMA

The RMA method proceeds with background subtraction, normalization and summarization using a deconvolution method for background correction, quantile normalization and the RMA (robust multichip average) algorithm for summarization.

edata <- oligo::rma(dataset)

## Background correcting
## Normalizing
## Calculating Expression

edata

## ExpressionSet (storageMode: lockedEnvironment)
## assayData: 12625 features, 4 samples 
##   element names: exprs 
## protocolData
##   rowNames: control_rep1 control_rep2 treat_rep1 treat_rep2
##   varLabels: exprs dates
##   varMetadata: labelDescription channel
## phenoData
##   rowNames: control_rep1 control_rep2 treat_rep1 treat_rep2
##   varLabels: index group
##   varMetadata: labelDescription channel
## featureData: none
## experimentData: use 'experimentData(object)'
## Annotation: pd.hg.u95av2

Access the quality after normalization

library(ggplot2)
library(cowplot)

## 
## Attaching package: 'cowplot'

## The following object is masked from 'package:ggplot2':
## 
##     ggsave

plot_grid(
    local({
        PCA_raw <- prcomp(t(log2(exprs(dataset))), scale = FALSE)
    
        dataGG <- data.frame(PC1 = PCA_raw$x[,1], PC2 = PCA_raw$x[,2],
                         Group = pData(dataset)$group)
            
        qplot(PC1, PC2, data = dataGG, color = Group,
              main = "PCA plot of the raw data (log-transformed)", asp = 1.0, geom = "text",
        label = sampleNames(dataset)) + scale_colour_brewer(palette = "Set2")
    }),
    local({
        PCA <- prcomp(t(exprs(edata)), scale = FALSE)
        dataGG <- data.frame(PC1 = PCA$x[,1], PC2 = PCA$x[,2],
                  Group = pData(dataset)$group)
        qplot(PC1, PC2, data = dataGG, color =  Group,
            main = "PCA plot of the normalized data", asp = 1.0, geom = "text",
            label = sampleNames(edata)) +
            scale_colour_brewer(palette = "Set2")
    })
)

par(mfrow = c(1,2))
boxplot(dataset)
boxplot(edata)

par(mfrow = c(1,1))

oligo::MAplot(edata, pairs = TRUE, ylim = c(-1, 1))

library(RColorBrewer)
library(pheatmap)
dists <- as.matrix(dist(t(exprs(edata)), method = "manhattan"))
diag(dists) <- NA
hmcol <- colorRampPalette(rev(brewer.pal(9, "PuOr")))(255)
pheatmap(dists, col = rev(hmcol), clustering_distance_rows = "manhattan",
                    clustering_distance_cols = "manhattan")

Filter based on intensity

Microarray data commonly show a large number of probes in the background intensity range. They also do not change much across arrays. Hence they combine a low variance with a low intensity. We want to filter these results as they may contribute to false positive results in the differential expression analysis. The bars represent the distribution of median intensities, the red vertical line represents the threshold value. If more than two samples of a gene have intensities larger than the threshold value, the gene will be kept.

edata_medians <- rowMedians(exprs(edata))
hist_res <- hist(edata_medians, 100, col="#e7efd8", freq = FALSE,
                 main = "Histogram of the median intensities",
                 xlab = "Median intensities")
emp_mu <- hist_res$breaks[which.max(hist_res$density)]
emp_sd <- mad(edata_medians)/2
prop_cental <- 0.50
lines(sort(edata_medians),
      prop_cental*dnorm(sort(edata_medians), mean = emp_mu, sd = emp_sd),
    col = "grey10", lwd = 3)

#cut_val <- 0.05 / prop_cental
thresh_median <- qnorm(0.05 / prop_cental, emp_mu, emp_sd)
abline(v = thresh_median, lwd = 4, col = "red")

samples_cutoff <- 2
idx_thresh_median <- apply(exprs(edata), 1, function(x){
                   sum(x > thresh_median) >= samples_cutoff})
table(idx_thresh_median)

## idx_thresh_median
## FALSE  TRUE 
##  4060  8565

edata <- subset(edata, idx_thresh_median)

Identification of differentially expressed genes

Create a design matrix. We will also consider the batch effects between replicates in order to remove them.

library(limma)

## 
## Attaching package: 'limma'

## The following object is masked from 'package:oligo':
## 
##     backgroundCorrect

## The following object is masked from 'package:BiocGenerics':
## 
##     plotMA

f <- factor(c("control", "control", "treat", "treat"))
batch <- factor(c("rep1", "rep2", "rep1", "rep2"))
design <- model.matrix(~ 0 + f + batch)
colnames(design)

## [1] "fcontrol"  "ftreat"    "batchrep2"

colnames(design) <- c("control", "treat", "batch")
design

##   control treat batch
## 1       1     0     0
## 2       1     0     1
## 3       0     1     0
## 4       0     1     1
## attr(,"assign")
## [1] 1 1 2
## attr(,"contrasts")
## attr(,"contrasts")$f
## [1] "contr.treatment"
## 
## attr(,"contrasts")$batch
## [1] "contr.treatment"

We can fit the linear model, define appropriate contrast to test the hypothesis on treatment effect and compute the moderated t–statistics by calling the eBayes function.

data.fit <- lmFit(exprs(edata), design)
head(data.fit$coefficients)

##            control    treat      batch
## 100_g_at  7.797102 8.130695 -0.2555991
## 1000_at   7.593515 7.814241  0.1087141
## 1002_f_at 4.612296 4.664025 -0.1417638
## 1003_s_at 6.091525 6.084982 -0.2208167
## 1004_at   5.848525 5.979026 -0.1795733
## 1005_at   8.166263 9.442745 -0.3113774

contrast.matrix <- makeContrasts(treat-control,levels=design)
data.fit.con <- contrasts.fit(data.fit,contrast.matrix)
data.fit.eb <- eBayes(data.fit.con)

Volcano plot to show the distribution of fold change and p value.

volcanoplot(data.fit.eb,highlight=10)

Then we sort the results by their absolute t-statistics.

top <- topTable(data.fit.eb, number = Inf)
head(top)

##              logFC  AveExpr         t      P.Value    adj.P.Val        B
## 37989_at -4.187867 5.058882 -21.61960 3.275846e-17 2.805762e-13 28.37765
## 38428_at -4.238118 5.624551 -17.23860 5.383859e-15 2.020411e-11 23.91131
## 40522_at  3.109750 8.149514  17.02754 7.076745e-15 2.020411e-11 23.66414
## 36629_at  2.923410 7.890901  16.27198 1.929437e-14 4.131408e-11 22.75146
## 41193_at -2.827488 6.189725 -15.87952 3.299633e-14 5.652271e-11 22.25952
## 1081_at  -3.472805 7.263600 -14.50306 2.373227e-13 3.387781e-10 20.43040

Check how many results can we get if we use a p value cutoff by 0.001 or an adjusted p value cutoff by 0.05.

table(top$adj.P.Val < 0.05)

## 
## FALSE  TRUE 
##  8117   448

table(top$P.Value < 0.001)

## 
## FALSE  TRUE 
##  8205   360

We would like to visualize the distribution of p value with a histogram.

hist(top$P.Value, col = brewer.pal(3, name = "Set2")[1], breaks = 100,
      main = "Control vs Treatment of DEX", xlab = "p-values")

Annotating genes

We need to annotate the gene names by the probe IDs.

suppressPackageStartupMessages({
    library(hgu95av2.db)
})
get_symbol <- function(probeid) {
    ans <- mapIds(hgu95av2.db::hgu95av2.db,
                  keys = probeid, keytype = "PROBEID", column = "SYMBOL", multiVals = "first")
    unname(ans)
}
get_genename <- function(probeid) {
    ans <- mapIds(hgu95av2.db::hgu95av2.db,
                  keys = probeid, keytype = "PROBEID", column = "GENENAME", multiVals = "first")
    unname(ans)
}
top$symbol <- get_symbol(rownames(top))

## 'select()' returned 1:many mapping between keys and columns

top$gene_name <- get_genename(rownames(top))

## 'select()' returned 1:many mapping between keys and columns

Generate a table of differentially expressed genes

We will use a cutoff of adjusted p value by 0.05.

selected <- top[, c("symbol", "gene_name", "logFC", "P.Value", "adj.P.Val")]
#top <- cbind(data.frame(probeID = rownames(top), stringsAsFactors = FALSE), top)
#rownames(top) <- NULL
selected <- selected[selected$adj.P.Val < 0.05,]
nrow(selected)

## [1] 448

We will only show the first 100 genes.

knitr::kable(selected[1:100,])

Gene clustering

The following is a heatmap of gene expression values of the differentially expressed genes.

library(pheatmap)
mat <- exprs(edata)[rownames(selected), ]
rownames(mat) <- selected$symbol
pData(edata)$sample_name <- rownames(pData(edata))
anno <- as.data.frame(pData(edata))[, c("group", "sample_name")]
pheatmap(mat, annotation_col = anno, show_rownames = FALSE,
         main = "Heatmap of gene expression levels")

We can roughly divide the genes (probes) into three clusters, and the control/treat samples are also clearly separated.

PCA plot

library(ggfortify)
autoplot(prcomp(exprs(edata)[rownames(selected), ]))

The PCA plot did not provide much useful information.

Pathway analysis

library(ReactomePA)

## 

## ReactomePA v1.24.0  For help: https://guangchuangyu.github.io/ReactomePA
## 
## If you use ReactomePA in published research, please cite:
## Guangchuang Yu, Qing-Yu He. ReactomePA: an R/Bioconductor package for reactome pathway analysis and visualization. Molecular BioSystems 2016, 12(2):477-479

entrez_ids <- mapIds(hgu95av2.db::hgu95av2.db, keys = rownames(selected),
                     keytype = "PROBEID", column = "ENTREZID", multiVals = "first")

## 'select()' returned 1:many mapping between keys and columns

entrez_ids <- entrez_ids[!is.na(entrez_ids)]
reactome_enrich <- enrichPathway(gene = entrez_ids, organism = "human")

barplot(reactome_enrich, showCategory = Inf)

dotplot(reactome_enrich, showCategory = Inf)

ReactomePA::emapplot(reactome_enrich, showCategory = Inf)

up_regulated <- selected[selected$logFC > 0, ]
up_regulated <- up_regulated[rownames(up_regulated) %in% names(entrez_ids), ]
ReactomePA::emapplot(
    enrichPathway(gene = entrez_ids[rownames(up_regulated)], organism = "human"),
    showCategory = Inf)

## Using `nicely` as default layout

down_regulated <- selected[selected$logFC < 0, ]
down_regulated <- down_regulated[rownames(down_regulated) %in% names(entrez_ids), ]
ReactomePA::emapplot(
    enrichPathway(gene = entrez_ids[rownames(down_regulated)], organism = "human"),
    showCategory = Inf)