source-code.Rmd

---
title: Differential Gene Expression Data from the Human Central Nervous System across Alzheimer’s Disease, Lewy Body Diseases, and the Amyotrophic Lateral Sclerosis and Frontotemporal Dementia Spectrum
author: Ayush Noori, Aziz M. Mezlini, Bradley T. Hyman, Alberto Serrano-Pozo, Sudeshna Das
bibliography: references.bib
csl: neurobiology-of-disease.csl
output:
  prettydoc::html_pretty:
    theme: cayman
    highlight: github
    toc: yes
editor_options: 
  chunk_output_type: inline
---

```{r setup, include = FALSE}
knitr::opts_chunk$set(message=FALSE, warning=FALSE)
```

# Overview

>In [Noori et al. (2020)](https://doi.org/10.1016/j.nbd.2020.105225), we performed a systematic review and meta-analysis of human CNS transcriptomic datasets in the public [Gene Expression Omnibus (GEO)](https://www.ncbi.nlm.nih.gov/geo/) and [ArrayExpress](https://www.ebi.ac.uk/arrayexpress/) repositories across Alzheimer's disease (AD), Lewy body diseases (LBD), and the amyotrophic lateral sclerosis and frontotemporal dementia (ALS-FTD) disease spectrum. Here, we provide the source code to generate the data quality reports and perform differential expression analyses for these datasets. Please see our accompanying publication in *Data In Brief* for additional information.


# Setup

Install [R](https://www.r-project.org/) to run our code. Requisite packages are loaded.

```{r load-packages}

# base functions for R/Bioconductor
library(Biobase)

# query Gene Expression Omnibus or ArrayExpress
library(GEOquery)
library(ArrayExpress)

# data pre-processing
library(oligo)
library(arrayQualityMetrics)
library(sva)

# differential expression analysis
library(limma)

# probe annotation
library(org.Hs.eg.db)
library(AnnotationDbi)

# data visualization
library(ggplot2)
library(EnhancedVolcano)

```


# Query Repository

Specified transcriptomics datasets are retrieved by querying either the Gene Expression Omnibus or ArrayExpress repositories. The following script performs differential expression analysis on [GSE1297](https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE1297), however, this analysis is consistent across all included datasets.

```{r load-data}

# load series and platform data from GEO
gset = getGEO("GSE1297", GSEMatrix =TRUE, AnnotGPL=TRUE)
if (length(gset) > 1) idx = grep("GPL96", attr(gset, "names")) else idx = 1
gset = gset[[idx]]

# clean column names
fvarLabels(gset) = make.names(fvarLabels(gset))

# sample groups generated by GEO2R utility
gsms = "1XX111XXX11X0X00000XXXX0XX1X00X"
sml = c()
for (i in 1:nchar(gsms)) { sml[i] = substr(gsms,i,i) }

# eliminate samples marked as "X"
sel = which(sml != "X")
sml = sml[sel]
gset = gset[ ,sel]

```


# RMA Normalization

Missing values are removed.

```{r remove-missing}

sumNA = apply(gset, 1, function(x) sum(is.na(x)))
idxNA = which(sumNA > 0)
if(length(idxNA) > 0) {gset = gset[-idxNA, ]}

```

Robust Multichip Average (RMA) normalization is performed using the `oligo` package [@carvalho_framework_2010; @irizarry_exploration_2003].

``` {r rma-normalization}

# subset expression matrix from ExpressionSet
ex = exprs(gset)

# detect if normalization is required
qx = as.numeric(quantile(ex, c(0., 0.25, 0.5, 0.75, 0.99, 1.0), na.rm=T))
LogC = (qx[5] > 100) ||
  (qx[6]-qx[1] > 50 && qx[2] > 0) ||
  (qx[2] > 0 && qx[2] < 1 && qx[4] > 1 && qx[4] < 2)

# perform RMA normalization
if (LogC) {
  ex = apply(ex, 2, as.numeric)
  exRMA = oligo::basicRMA(ex, rownames(gset))
  exprs(gset) = exRMA
}

```

Disease labels are generated and the `ExpressionSet` is reordered.

``` {r disease-labels}

# map group label to AD vs. CTRL
dis = sml; dis[dis == 1] = "AD"; dis[dis == 0] = "CTRL"
dis = factor(dis, levels=c("CTRL", "AD"))
rel = order(dis)

# reorder ExpressionSet and label
gset = gset[, rel]; sml = sml[rel]; dis = dis[rel]
pData(gset)$Condition = dis

```


# Quality Control

The `arrayQualityMetrics` package is used to generate data quality reports. Outliers are detected via boxplots, MA plots, and inter-array distance comparison [@kauffmann_microarray_2010; @kauffmann_arrayqualitymetrics_2009]. **The full data quality report for GSE1297 is available [here](https://ayushnoori.github.io/nd-diff-expr/QC).**

``` {r qc-report}

# generate quality control report
qc = arrayQualityMetrics(gset, outdir = "QC", reporttitle ="GSE1297 QC - CA1, GPL96", force = TRUE)

# identify outliers
array = qc$arrayTable[, 3:5]
keep = apply(array == "", 1, all)
remove = rownames(pData(gset))[!keep]

```

Following the above normalization, the expression matrix is visualized in a boxplot which indicates aberrant samples [@wickham_ggplot2_2016].

``` {r boxplot}

# parse data for boxplot
ex = exprs(gset)
samples = rep(rownames(pData(gset)), each = nrow(gset))
samples = factor(samples, levels = unique(samples))
pheno = rep(as.character(dis), each = nrow(gset))
bp = data.frame(samples, c(ex), pheno)
bp[["outlier"]] = bp$samples %in% remove

# create boxplot
box = ggplot(bp, aes(x = samples, y = ex, fill = pheno, linetype = outlier)) +
  geom_boxplot(outlier.alpha = 0.1, na.rm = TRUE) + 
  scale_fill_manual(values =  c("#FE7B72", "#84B2D7")) +
  scale_linetype_manual(values = c("solid", "dashed")) + 
  labs(title = "GSE1297 Box Plot - CA1, GPL96",
       x = "Samples", y = "Probe Intensity", fill = "Condition", linetype = "Outlier?") +
  theme(plot.title = element_text(hjust = 0.5, size = 16, face="bold"),
        axis.title.x = element_text(size=12, face="bold"),
        axis.title.y = element_text(size=12, face="bold"),
        legend.title = element_text(size=10, face="bold"))

```

Samples which fail to pass one or more of the three outlier detection steps in `arrayQualityMetrics` are omitted.

``` {r remove-outliers}

gset = gset[, keep]
sml = sml[keep]
dis = dis[keep]

```

Probes with low expression are capped at the 20th percentile. This threshold is represented on the final boxplot.

``` {r filter-low, fig.width=17, fig.height=8, dpi=600}

# cap probes with low-expression
filter = exprs(gset) # with outliers removed
thresh = quantile(c(filter), 0.20)
filter[which(filter < thresh)] = thresh
exprs(gset) = filter

# boxplot with low-expression threshold
box + geom_hline(yintercept=thresh, linetype="dashed", color = "red", size = 0.5)

```


# Surrogate Variable Analysis

Model matrices are created.

``` {r model-matrices}

# set group names
sml = paste("G", sml, sep="")
fl = as.factor(sml)
gset$description = fl

# create full model matrix
mod = model.matrix(~ description, gset)

# create null model matrix
mod0 = model.matrix(~1, gset)

```

The number of significant surrogate variables ($\le$ 4) is estimated using the asymptotic approach proposed by @leek_asymptotic_2011, then SVA is performed [@leek_sva_2012; @leek_capturing_2007].

``` {r sva}

# estimate # of surrogate variables
n.sv = num.sv(exprs(gset), mod, method="leek")
if(n.sv > 4) { n.sv = 4 }

# perform SVA
svobj = sva(exprs(gset), mod, mod0, n.sv=n.sv)
modSv = cbind(mod, svobj$sv)

```


# Differential Expression Analysis

Differential expression analysis is performed using the `limma` package [@phipson_robust_2016; @ritchie_limma_2015].

``` {r limma-analysis}

# generate linear model fit
fit = lmFit(gset, modSv)

# compute contrasts
cont.matrix = cbind("C1"=c(0,1,rep(0,svobj$n.sv)))
fit2 = contrasts.fit(fit, cont.matrix)

# empirical Bayes statistics for DE
fit2 = eBayes(fit2, 0.01)

```

The ranked table of all differentially-expressed genes (DEGs) with confidence intervals included is created. Then, the interquartile range (IQR) is calculated from the raw expression matrix.

``` {r calculate-iqr}

# get DE genes
tT = topTable(fit2, adjust="fdr", sort.by="B", number=nrow(fit2), confint=TRUE)

# calculate IQR
ex = exprs(gset)[rownames(tT), ]
IQR = apply(ex, 1, function(x) IQR(x, na.rm = TRUE))
tT$IQR = IQR

```


# Probe Annotation

The annotation package specified below corresponds to Affymetrix Human Genome U133A Array (i.e., GPL96) for GSE1297. Annotation packages for other microarrays are available from [R/Bioconductor](https://www.bioconductor.org/packages/release/data/annotation/). Probes wihtout Entrez ID mapping are removed [@maglott_entrez_2007].

``` {r probe-annotation}

# load annotation package
library(hgu133a.db)

# perform mapping and remove duplicate probes
annot = AnnotationDbi::select(hgu133a.db, tT$ID, c("SYMBOL", "ENTREZID", "GENENAME"))
annot = annot[!duplicated(annot$PROBEID), ]

# join DE genes with annotation
dat = cbind(annot[, c("PROBEID", "ENTREZID", "SYMBOL")], tT[, c("P.Value", "logFC", "CI.L", "CI.R", "AveExpr", "IQR")])
colnames(dat) = c("Probe", "Entrez", "Symbol", "pVal", "logFC", "CI.L", "CI.R", "AveExpr", "IQR")

# remove probes without Entrez ID mapping
dat = dat[!is.na(dat$Entrez), ]

```

In the event that multiple probes map to the same gene, the single probe with the greatest interquartile range (IQR) is retained [@walsh_microarray_2015; @wang_detecting_2012].

``` {r filter-duplicate}

# group by ENTREZ ID, then order by IQR (greatest to least)
dat = dat[order(dat$Entrez, -dat$IQR), ]

# retain probe with the greatest IQR only
dat = dat[!duplicated(dat$Entrez), ]

```

For each gene, the z-score is computed as $\frac{x - \mu}{\sigma}$. Genes are subsequently ordered by absolute z-score.

``` {r z-score}

# compute z-scores
SE = (dat$CI.R - dat$CI.L)/(1.96*2)
dat$Z = dat$logFC/SE

# order by absolute z-score
dat = dat[order(-abs(dat$Z)), ]

```

The final table of DEGs is excerpted below.

``` {r final-table, echo = FALSE}
library(knitr)
library(kableExtra)
kable(dat[1:15, ], row.names = FALSE, align = "c") %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive", "bordered"), font_size=12)
```

Top DEGs are visualized via a volcano plot [@blighe_enhancedvolcano_2019].

``` {r volcano-plot, fig.width=16.7, fig.height=10, dpi=600}

EnhancedVolcano(dat,
                lab = dat$Symbol,
                x = "logFC",
                y = "pVal",
                ylim = c(0, max(-log10(dat[, "pVal"]), na.rm=TRUE)),
                title = "GSE1297 Volcano Plot - CA1, GPL96",
                subtitle = "p-value Threshold = 0.05, FC Threshold = 2",
                caption = NULL,
                pCutoff = 0.05,
                FCcutoff = log2(2))

```


# Citation

**Research Paper:** Noori, A.,<sup>1,2,3,4</sup>, Mezlini, A.M.,<sup>2,3,4,5</sup>, Hyman, B.T.,<sup>2,4,5</sup>, Serrano-Pozo, A.,<sup>2,4,5*</sup>, Das, S.,<sup>2,3,4,5</sup>. 2020. Systematic review and meta-analysis of human transcriptomics reveals neuroinflammation, deficient energy metabolism, and proteostasis failure across neurodegeneration. Neurobiology of Disease  149, 105225. https://doi.org/10.1016/j.nbd.2020.105225

**Methods Paper:** Noori, A.,<sup>1,2,3,4</sup>, Mezlini, A.M.,<sup>2,3,4,5</sup>, Hyman, B.T.,<sup>2,4,5</sup>, Serrano-Pozo, A.,<sup>2,4,5*</sup>, Das, S.,<sup>2,3,4,5</sup>. 2021. Differential gene expression data from human central nervous system across Alzheimer’s disease, Lewy body diseases, and amyotrophic lateral sclerosis and frontotemporal dementia spectrum. In press.

### Affiliations
1.	Harvard College, Cambridge, MA 02138, United States of America
2.	Department of Neurology, Massachusetts General Hospital, Boston, MA 02114, United States of America
3.	MIND Data Science Lab, Cambridge, MA 02139, United States of America
4.	MassGeneral Institute for Neurodegenerative Disease, Charlestown, MA 02129, United States of America
5.	Harvard Medical School, Boston, MA 02115, United States of America

*\* Co-Corresponding Authors*


# References