Array analysis for Erin Dunn ([email protected]) at MGH. Contact Meeta Mistry ([email protected]) for additional details. Request from client was:
It doesn't look like BrainCloud is an ideal dataset because of the issues related to the PMI. Given this, I think we should probably focus on BrainSpan and run a parallel set of analyses to see if things look cleaner in that new dataset.
- grab the BrainSpan data set and metadata
- re-run the QC of the metadata and basic expression analysis
- isolate one brain region only: Orbitofrontal Cortex and Amygdala
library(ggplot2)
library(gtable)
library(scales)
library(RColorBrewer)
library(GEOquery)
library(affy)
library(arrayQualityMetrics)
library(reshape)
library(xtable)
library(ruv)
library(limma)
library(Biobase)
library(gridExtra)
library(stringr)
library(knitr)
library(png)
library(sva)
source("http://dl.dropboxusercontent.com/u/4253254/Resources/functions.r")- get base directory for analyses
- specify data and results directories
- specify column headers used in metadata file
# Setup directory variables
baseDir <- "."
dataDir <- file.path(baseDir, "data")
metaDir <- file.path(dataDir, "meta")
resultsDir <- file.path(baseDir, "results")RMA background corrected data, with quantile normalization and log2-transformation. The median of all probe sets within one gene (transcript cluster) was used as the estimate of gene expression. In the supplement '...a total of 17,565 mainly protein coding genes were surveyed'
# Load GEO data
gse_gene <- getGEO(filename = file.path(dataDir, "geo/GSE25219-GPL5175_series_matrix.txt.gz"))
# gse_probe <- getGEO(filename=file.path(dataDir,
# 'geo/GSE25219-GPL5188_series_matrix.txt.gz'))pheno_gene <- read.delim(file.path(metaDir, "brainspan_samples_metadata.txt"),
row.names = 1)
meta_donor <- read.delim(file.path(metaDir, "brainspan_donor_metadata.txt"),
row.names = 1)
# Subset data by region and hemisphere
meta_ofc <- pheno_gene[which(pheno_gene$BrainRegion == "OFC" & pheno_gene$Hemisphere ==
"R"), ]
meta_ofc$PMI <- as.numeric(as.character(meta_ofc$PMI))
meta_ofc$Stage <- factor(meta_ofc$Stage)
meta_amy <- pheno_gene[which(pheno_gene$BrainRegion == "AMY" & pheno_gene$Hemisphere ==
"R"), ]Phenotype data is loaded in from which we can isolate our two brain regions of interest. In total there are 57 donors, and from each multiple samples were taken from various brain regions.
Exploring the metadata for consistency, and generating a quick overview:
# Gender distribution
gender <- rbind(table(meta_ofc$Sex), table(meta_amy$Sex))
row.names(gender) <- c("Orbitofrontal cortex", "Amygdala")
kable(gender, format = "html")| F | M | |
|---|---|---|
| Orbitofrontal cortex | 17 | 19 |
| Amygdala | 18 | 19 |
# Age table
age.table <- as.character(unique(meta_donor$Period))
age.table <- strsplit(age.table, ",")
age.table <- do.call("rbind", age.table)
row.names(age.table) <- rep("", nrow(age.table))
colnames(age.table) <- c("Stage", "Description")
kable(data.frame(age.table), format = "html", row.names = FALSE)| Stage | Description |
|---|---|
| 1 | Embryonic |
| 2 | Early fetal |
| 3 | Early fetal |
| 4 | Early mid-fetal |
| 5 | Early mid-fetal |
| 6 | Late mid-fetal |
| 7 | Late fetal |
| 8 | Neonatal and early infancy |
| 9 | Late infancy |
| 10 | Early childhood |
| 11 | Middle and late childhood |
| 12 | Adolescence |
| 13 | Young adulthood |
| 14 | Middle adulthood |
| 15 | Late adulthood |
# Age distribution
ggplot(meta_ofc, aes(Stage)) + geom_bar() + ggtitle("Orbitofrontal Cortex: Age distribution") +
xlab("Developmental Stage") + ylab("Number of Samples") + theme(axis.text.x = element_text(colour = "grey20",
size = 15, angle = 0, hjust = 0.5, vjust = 0.5, face = "plain"), plot.title = element_text(size = rel(2)),
axis.title = element_text(size = rel(1.25)))
# pH measurements
ggplot(na.omit(meta_ofc), aes(x = Stage, y = pH, fill = Stage)) + geom_boxplot() +
ggtitle("Orbitofrontal Cortex: pH levels") + xlab("Stage") + guides(fill = FALSE) +
theme(plot.title = element_text(size = rel(2)), axis.title = element_text(size = rel(1.5)),
axis.text = element_text(size = rel(1.25)))
# Time between death and sample collection; remove non-numeric values
ggplot(na.omit(meta_ofc), aes(x = Stage, y = PMI, fill = Stage)) + geom_boxplot() +
ggtitle("Orbitofrontal Cortex: Postmortem Intervals") + xlab("Stage") +
ylab("Postmortem Interval") + guides(fill = FALSE) + theme(legend.position = "none",
plot.title = element_text(size = rel(2)), axis.title = element_text(size = rel(1.5)),
axis.text = element_text(size = rel(1.25)))
# RNA Integrity
ggplot(na.omit(meta_ofc), aes(x = Stage, y = RIN, fill = Stage)) + geom_boxplot() +
ggtitle("Orbitofrontal Cortex: RIN") + xlab("Stage") + ylab("Postmortem Interval") +
guides(fill = FALSE) + theme(legend.position = "none", plot.title = element_text(size = rel(2)),
axis.title = element_text(size = rel(1.5)), axis.text = element_text(size = rel(1.25)))
As we have seen with the Braincloud data, there is a positive correlation observed with devleopemental stage and PMI and a negative correlations with RIN. Not surprising, fetal samples are more likely to have been obtained faster. The authors also acknolwedge this in the manuscript and account for it by incorporating both as covariates in the model.
ArrayQualityMetrics QC report for GSE25219
arrayQualityMetrics(expressionset = eset.ofc, intgroup = c("fetal"), outdir = "./results/report_OFC",
force = TRUE, do.logtransform = FALSE)
A false color heatmap of the distance between arrays demonstrates a high degree of similarity among fetal samples and likewise with non-fetal samples. Two fetal samples clustering with postnatal, simialr to dendorgram above. Remove these samples.
Density plots (smoothed histograms) for all arrays follow a similar distribution shape and range. The shape of distribution is questionable with a fairly heavy right tail.
Try checking the quality on a handful of .CEL files and see if we see the same dsitribution. If not, it might be better to work directly from .CEL files.
Reading in : ./data/geo/CEL/GSM703924_HSB92-OFC-R.CEL.gz Reading in : ./data/geo/CEL/GSM703972_HSB97-OFC-R.CEL.gz Reading in : ./data/geo/CEL/GSM704048_HSB100-OFC-R.CEL.gz Reading in : ./data/geo/CEL/GSM704080_HSB102-OFC-R.CEL.gz Reading in : ./data/geo/CEL/GSM704110_HSB103-OFC-R.CEL.gz
| SampleName | BrainCode | BrainRegion | Hemisphere | Sex | Age | Stage | PMI | pH | RIN | |
|---|---|---|---|---|---|---|---|---|---|---|
| GSM703924 | HSB92_OFC_R | HSB92 | OFC | R | M | 21 PCW | 6 | 4 | 6.65 | 9.1 |
| GSM703972 | HSB97_OFC_R | HSB97 | OFC | R | F | 17 PCW | 5 | 1 | NA | 9.9 |
| GSM704048 | HSB100_OFC_R | HSB100 | OFC | R | F | 19 PCW | 6 | 4 | 6.56 | 9.5 |
| GSM704080 | HSB102_OFC_R | HSB102 | OFC | R | F | 21 PCW | 6 | >13 | 5.89 | 9.4 |
| GSM704110 | HSB103_OFC_R | HSB103 | OFC | R | M | 13 PCW | 4 | 1.5 | NA | 9.3 |
The raw data seems better than what we obtained from GEO, with higher signal intensities. Two samples show a slightly wider distribution and another is skewed to the left. Will check how this changes after normalization.
The data was normalized for differential gene expression analysis using Robust Multichip Average (RMA) in the oligo BioConductor package. Here, RMA normalizes the intensity values at the probe level, and collapses probes into "core" transcripts based on annotations provided by Affymetrix.
geneSummaries <- rma(affyRaw, target = "core", background = T, normalize = T)Repeat the previous QC using the normalized data, and we see that the distributions are similar to what we had found originally pulled from GEO. One option is removing those wide distribution samples.
In the Kang et al study, age was evaluated using ANOVA and both RIN and PMI were included as covariates. To stay consistent, we performed an ANCOVA the same on Orbitofrontal cortex data, modeling Age/Developemental Stage as a continuous variable.
# Remove outlier samples
remove <- c(which(label(ddata)$x == 1), which(label(ddata)$x == 2))
# Remove NA values
meta_new <- meta_new[-remove, ]
meta_new <- meta_new[which(!is.na(meta_new$PMI)), ]
data_new <- exprs(eset.ofc)[, rownames(meta_new)]
# Update expression set
exprs(eset.ofc) <- data_new
pData(eset.ofc) <- meta_new
# Model fit
mod <- model.matrix(~Stage + PMI + RIN, pData(eset.ofc))
fit <- lmFit(eset.ofc, mod)
fit <- eBayes(fit)
topStage <- topTable(fit, coef = 2, number = nrow(exprs(eset.ofc)), adjust.method = "BH")
hist(topStage$P.Value, col = "grey", border = F, main = "P-value distribution: Age",
xlab = "P-value")
topPMI <- topTable(fit, coef = 3, number = nrow(exprs(eset.ofc)), adjust.method = "BH")
hist(topPMI$P.Value, col = "grey", border = F, main = "P-value distribution: PMI",
xlab = "P-value")
topRIN <- topTable(fit, coef = 4, number = nrow(exprs(eset.ofc)), adjust.method = "BH")
hist(topRIN$P.Value, col = "grey", border = F, main = "P-value distribution: RIN",
xlab = "P-value")
Alot of differentially expressed genes with Age/Developemntal Stage (3020) even at a quite stringent threshold (padj < 0.001). Suprisingly few significant changes associated with RIN (880) and none associated with postmortem interval.
Take the top 6 probes that are affected by age. Even though the changes are not identical, we still see a similarity in the trend as we did with the Braincloud data.
