Code.Rmd

---
title: "ApoE4 dose effects on serum metabolome in Alzheimer's Disease"
subtitle: "a Data Science approach"
author: "George Miliarakis, MSc Thesis results"
date: "18-03-2024"
output:
    tufte::tufte_html: default
    toc: yes
    toc_depth: 4
    highlight: pygments
editor_options:
  chunk_output_type: inline
---

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

```{r packages, include=FALSE}
require(heatmaply)
require(globaltest)
require(ggplot2)
require(ggthemes)
require(plotly)
require(furrr)
require(caret)
require(dplyr)
require(e1071)
require(FMradio)
require(pROC)
require(xgboost)
require(rpart)
require(DMwR2)
require(nnet)
require(rags2ridges)
```

## Data Exploration

The (interactive) correlation heatmap reveals very high correlation among aminoacids and TG compounds.

```{r heatmap}
load("data/data.Rdata")
X <- df[df$Diagnosis == "Probable AD", 1:230]
C <- round(cor(X), 2)
heatmaply_cor(C, color = "magma", plot_method = "plotly", dendrogram = F, reorderfun = sort.default(d, w), main = "Correlation Heatmap", file = "heatmap.html", colorbar_thickness = 15, colorbar_len = 0.5)
```

## Differential expression of metabolites per ApoE genotype
### Global test in AD

Performing a Global Test of ApoE4 dose-effects, on serum metabolites of AD patients, correcting for sex (Ho: ApoE4 dose has no effect on mean metabolite levels, Ha: it has an effect), yields a significant difference (p = 0.0171). The most significantly affected metabolites are triglycerides and diglycerides (FDR-adjusted p-value <0.05).
```{r gt, fig.height=8}
AD <- subset(df, Diagnosis == "Probable AD")
AD$E4dose <- as.factor(AD$E4dose)

gt <- globaltest::gt(E4dose ~ 1, E4dose ~ . - E4 - Diagnosis - target, data = AD, standardize = TRUE)
gt
covariates <- covariates(gt, sort = T)
res.gt <- extract(covariates) %>%
  sort() %>%
  p.adjust(, method = "FDR") %>%
  globaltest::result() %>%
  mutate_if(is.numeric, round, digits = 3)
res.gt$direction <- as.factor(res.gt$direction)
levels(res.gt$direction) <- c("no ApoE4", "1 ApoE4", "2 ApoE4")
res.gt$alias <- NULL
res.gt <- res.gt[, 1:4]
names(res.gt) <- c("Inheritance", "Assoc. with", "FDR", "p-value")
res.gt <- res.gt[, c(1:(ncol(res.gt) - 2), ncol(res.gt), (ncol(res.gt) - 1))]
filter(res.gt, `p-value` < 0.05)
```

### Nested Linear Models
Metabolite-level models to test for ApoE4 dose effects, correcting for age at diagnosis, sex, smoking status, alcohol consumption, hypertension (and medication), hyperlipidemia (and medication), anticoagulant medication, anti-depressants, mean arterial pressure (MAP) and body mass index (BMI). The F-tests yielded p-values < 0.05, however none of them survived FDR correction. Interestingly, in the SCD group, aminoacids are mostly affected and no TGs or DGs compared to the AD group.
```{r}
nested <- function(Y, x, ...) {
  ### Analysis of Variance (ANOVA)
  x <- as.factor(x)
  df <- cbind(Y, x, ...)
  ncol <- ncol(Y)
  covariates <- names(...)
  F_tests <- furrr::future_map(df[, 1:ncol], ~ {
    frm <- as.formula(paste0(".x ~", paste(covariates, collapse = "+")))
    rm <- lm(frm, data = df)
    ffm <- as.formula(paste0(".x ~ x +", paste(covariates, collapse = "+")))
    fm <- lm(ffm, data = df)
    anova(rm, fm)
  })
  # Correction for Multiple Testing
  # Create a list to store the p-values
  # Extract the p-values of the F-tests from the anova summaries list and store them in p_values
  p_values <- F_tests %>%
    future_map(~ .x[["Pr(>F)"]][[2]])
  # Coerce p_values to dataframe and transpose it
  p_values <- as.data.frame(p_values) %>%
    t() %>%
    data.frame(row.names = colnames(Y))
  # Calculate the FDR-adjusted p-values
  p_values$p_adj <- p.adjust(p_values[, ], method = "fdr", n = 230)
  p_values <- round(p_values, 3)
  # Filter out the non-significant (a=0.05) FDR-adjusted p-values
  sig <- as.data.frame(dplyr::filter(p_values, `.` < 0.05))
  sig <- sig[order(sig$p_adj), ]
  names(sig) <- c("P(>F)", "FDR")
  cat("Significant F-tests: \n")
  print(sig)
  cat("\nDose effects on metabolites:")
  summaries <- furrr::future_map(df[, rownames(sig)], ~ {
    f <- as.formula(paste0(".x ~ x +", paste(covariates, collapse = "+")))
    mdl <- lm(f, data = df)
    summary(mdl)$coefficients[1:length(table(x)), ]
  })
  summaries <- lapply(summaries, function(x) {
    x <- x[, -c(2, 3)]
    x <- t(x)
  })
  summaries <- plyr::ldply(summaries, id = 1:2) %>%
    mutate(across(where(is.numeric), round, digits = 3))
  coeffs <- summaries[c(TRUE, FALSE), ]
  pvalues <- summaries[c(FALSE, TRUE), ]
  t_tests <- cbind.data.frame(coeffs[, 1], coeffs[, 2], pvalues[, 2], coeffs[, 3], pvalues[, 3], coeffs[, 4], pvalues[, 4])
  names(t_tests) <- c("Metabolite", "No ApoE4", "P(>t)", "ApoE4x1", "P(>t)", "ApoE4x2", "P(>t)")
  return(t_tests)
}
```
```{r}
Y <- df[, 1:230]
df$E4dose <- as.numeric(df$E4dose) - 1
df$E4dose[df$E4dose == 0] <- "No ApoE4"
df$E4dose[df$E4dose == 1] <- "1 ApoE4"
df$E4dose[df$E4dose == 2] <- "2 ApoE4"
E4dose <- df$E4dose <- as.factor(df$E4dose)
```
### AD group
Among AD patients, ApoE4 dose seems to have positive effect on several triglycerides, diglycerides, putrescine, 2-ketoglutraric acid, lysophosphatidylcholin, HpH.PAF.-C16.0 and -C18.0
```{r}
load("data/clinical.Rdata")
AD <- subset(df, Diagnosis == "Probable AD")
ADY <- AD[, 1:230]
ADE4dose <- AD$E4dose
ADclinical <- clinical[df$Diagnosis == "Probable AD", ]

## Testing for effects of ApoE4 dose
nested(ADY, relevel(ADE4dose, ref = "No ApoE4"), ADclinical)
```

### SCD group
Among individuals with SCD, lipid metabolites were not affected as much as in the AD group, with only two sphingomyelin species showing a difference. Aminoacids L-serine, tryptophan, glycine,trytptophan, L-homoserine, putrescine were affected in this group.
```{r}
SCD <- df[df$Diagnosis == "Subjectieve klachten", ]
SCDY <- SCD[, 1:230]
SCDE4dose <- SCD$E4dose
SCDclinical <- clinical[df$Diagnosis == "Subjectieve klachten", ]

## Testing for effects of ApoE4 dose
nested(SCDY, relevel(SCDE4dose, ref = "No ApoE4"), SCDclinical)
```
## Multi-class classification of ApoE4 presence and AD
The discriminative potential of the metabolites, correcting for clinical background features, on AD and ApoE4 or not (4-class response ADE4, AD, SCDE4, SCD) is assessed by first fitting a Multi-nomial Logistic Regression (benchmark model). The benchmark model fits only the clinical features. The model's AUC is compared with the same model (with altered hyperparameters) fitting clinical variables + 230 metabolites. The metabolites are then projected to 6 latent ML-estimated factors and 3 more models are fitted: Multi-nomial Logistic Regression, Decision Tree and XGBoost fitting clinical variables + 6 factors. The models' multi-class classification performance is compared using confusion matrices and AUCs.
```{r multifit}
multifit <- function(X,
                     y,
                     model,
                     ctrl = NULL,
                     grid = NULL,
                     seed = 87654, ...) {
  set.seed(seed)
  # Merge X and y into df
  df <- cbind.data.frame(X, y)
  # Train the model
  model <- caret::train(df[, 1:ncol(X)], df$y,
    method = model,
    tuneGrid = grid,
    trControl = ctrl,
    metric = "logLoss",
    # importance=TRUE,
    ...
  )
  obs <- model$pred$obs
  preds <- model$pred$pred

  # Get the multi-class ROC curves
  ys <- as.numeric(obs) - 1
  yhats <- as.numeric(preds) - 1
  roc <- multiclass.roc(
    response = ys,
    predictor = yhats, quiet = T, legacy.axes = T
  )
  print(roc$auc)
  names <- c(
    "ADE4 vs. AD", "ADE4 vs. SCDE4",
    "ADE4 vs. SCD", "AD vs SCDE4", "AD vs. SCD", "SCDE4 vs. SCD"
  )
  for (i in 1:length(names)) {
    names[[i]] <- paste0(
      names[[i]], ", AUC=",
      paste(round(auc(roc$rocs[[i]]), 2))
    )
  }
  names(roc$rocs) <- names
  cat("\n")
  # Create a confusion matrix and get performance metrics from caret
  cm <- confusionMatrix(reference = obs, data = preds, mode = "everything")
  print(cm)
  out <- list("model"= model, "roc"= roc)
  return(out)
}
```

```{r multi_preparation}
load("data/data.Rdata")
X <- scale(df[, 1:230])
y <- df$target

Xclin <- cbind.data.frame(X, clin_dummy)


ctrl <- trainControl(
  method = "repeatedcv",
  number = 10,
  repeats = 100,
  savePredictions = "final",
  classProbs = T,
  summaryFunction = multiClassSummary,
  selectionFunction = best,
  search = "random",
  sampling = "smote"
)
```

### Fitting a benchmark model: clinical characteristics only

```{r benchmark}
benchmark <- multifit(
  X = clin_dummy,
  y = y,
  model = "multinom",
  ctrl = ctrl,
  grid = expand.grid(decay = 0),
  trace = F
)

g <- ggroc(benchmark$roc$rocs) + scale_color_tableau() + theme_tufte() + guides(color = guide_legend(title = "ROC curves"))
plotly::ggplotly(g)
 ```

### Regularized Multinomial Regression adding 230 metabolites

```{r mlr230}
mlr <- multifit(
  X = Xclin,
  y = y,
  model = "multinom",
  ctrl = ctrl,
  grid = expand.grid(decay = 9.187724),
  trace = F,
  importance=T
)
g <- ggroc(mlr$roc$rocs) + scale_color_tableau() + theme_tufte() + guides(color = guide_legend(title = "ROC curves"))
plotly::ggplotly(g)
# Get the feature importance scores
fimp <- varImp(mlr$model)

fimp$importance %>%
  arrange(desc(Overall)) %>%
  top_n(20)
```

### Projection to Latent Factors

```{r get_thomson, results = FALSE}
project <- function(X, m, seed) {
  set.seed(seed)
  X <- scale(as.matrix(df[, 1:230]))
  cov <- cor(X)

  # Find redundant features
  filter <- RF(cov)

  # Filter out redundant features
  filtered <- subSet(X, filter)

  # Regularized correlation matrix estimation
  M <- regcor(filtered)

  # Get the regularized correlation matrix of the filtered dataset
  R <- M$optCor

  mlfa <- mlFA(R, m = 6)
  thomson <- facScore(filtered, mlfa$Loadings, mlfa$Uniqueness)
  return(thomson)
}
```

```{r, results = FALSE}
thomson <- project(X, 6, seed = 1234)
X <- cbind(clin_dummy, thomson)
```

### Multinomial Logistic Regression adding 6 factors
```{r ml6}
mlrf <- multifit(
  X = X,
  y = y,
  model = "multinom",
  ctrl = ctrl,
  trace = FALSE,
  grid = expand.grid(decay = 0)
)
g <- ggroc(mlrf$roc$rocs) + scale_color_tableau() + theme_tufte() + guides(color = guide_legend(title = "ROC curves"))
plotly::ggplotly(g)
```

### Decision Tree

```{r tree}
tree <- multifit(
  X = X,
  y = y,
  model = "rpart2",
  ctrl = ctrl,
  grid = expand.grid(maxdepth = 3)
)
g <- ggroc(tree$roc$rocs) + scale_color_tableau() + theme_tufte() + guides(color = guide_legend(title = "ROC curves"))
plotly::ggplotly(g)
```

### XGBoost Forest

```{r xgb}
xgb <- multifit(
  X = X,
  y = y,
  model = "xgbTree",
  ctrl = ctrl,
  grid = xgb.grid
)
g <- ggroc(xgb$roc$rocs) + scale_color_tableau() + theme_tufte() + guides(color = guide_legend(title = "ROC curves"))
plotly::ggplotly(g)
# cvms::plot_confusion_matrix(cmxgb, palette = "Oranges" , theme_fn = ggthemes::theme_tufte)
```
### Model Comparison
```{r}
aucs <- data.frame(AUC = c("Clinical features only" = benchmark$roc$auc,
"Clinical features + 230 metabolites" = mlr$roc$auc,
"Clinical features + 6 latent factors" = mlrf$roc$auc,
"Decision Tree" = tree$roc$auc, "XGBoost" = xgb$roc$auc))

aucs
```
Observations:

1. Adding serum metabolite information (either the full 230-metabolite matrix or its 6-factor projection) seems to slightly increase the discriminatory power of the models. 

2. Fitting 6 ML-estimated factors obtained by the FMradio package (cummulatively explaining 30% of variace) yields slightly increased classification performance.

3. Looking at the confusion matrices and individual ROC curves, all models were able to discriminate better among certain classes (AD+E4/SCD+E4, AD+E4/SCD, AD-E4/SCD+E4 and AD-E4/SCD-E4) compared to others (AD+E4/AD-E4 and SCD+E4/SCD-E4).


## Metabolite Covariance Network Analysis in AD
The covariance network analysis shows distinct metabolic "images" among ApoE4 carriers and non-carriers in AD. The network vertices and edges are distinct. 
```{r preparation}
# Store all observations of ApoE4 non-carriers in C1
C1 <- scale(AD[AD$E4 == 0, 1:230])
# Store all observations of ApoE4 carriers in C2
C2 <- scale(AD[AD$E4 == 1, 1:230])

# Get the covariance matrices of C1 and C2
S1 <- covML(C1)
S2 <- covML(C2)

# Store them in a list
S <- list(S1 = S1, S2 = S2)

# Get the total number of samples
n <- c(nrow(S1), nrow(S2))

# Create a list of fused covariance matrices T
Ts <- default.target.fused(Slist = S, ns = n, type = "DUPV")
```

```{r, eval=FALSE}
# Get the optimal lambdas per class and fused with 10-fold CV
# set.seed(8910)
# optf <- optPenalty.fused(
#   Ylist = Ys,
#   Tlist = Ts,
#   lambda = default.penalty(Ys),
#   cv.method = "kCV",
#   k = 10,
#   verbose = FALSE
# )
# save(optf, file= "data/optf.Rdata")
```

```{r sparsify, echo = FALSE}
# Create a list of the two-class data Y
load("data/optf.Rdata")
Ys <- list(C1 = C1, C2 = C2)
Ps <- optf$Plist

Ps <- ridgeP.fused(S, n, Ts, lambda = optf$lambda, maxit = 1000)

# Get the sparsified high precision matrices, correcting FDR at .001
P0s <- sparsify.fused(Ps,
  threshold = "localFDR",
  FDRcut = 0.9999
)
```
The Gaussian graphical modelling between ApoE4 carriers and non-
carriers with AD unveiled distinct correlations between the metabolites. The ApoE4 carriers’
network was more sparse, having less edges compared to the non-carriers’. Moreover, the metabolites
could be clearly separated in communities based on their covariances in ApoE4 carriers, as shown in
The communities found contain similar compounds (phosphatidylcholines, triglycerides, sphin-
gomyelines, prostaglandines) or molecules that share metabolic pathways (amines, amines-organic acids,
aspartic-glutamic acid, oxidative stress compounds-lipids). The ApoE4 non-carrier’s network appeared less
cohesive, with only a few small aminic and organic acid communities.
```{r GGMs per class and diff, message=FALSE, results = FALSE}
# Merge the sparse high precision matrices
TST <- Union(P0s$S1$sparseParCor, P0s$S2$sparseParCor)
PCE4NO <- TST$M1subset
PCE4YES <- TST$M2subset

# Create a color map per metabolite class
Colors <- rownames(PCE4YES)
Colors[grep("Am", rownames(PCE4YES))] <- "pink"
Colors[grep("OA", rownames(PCE4YES))] <- "lightblue"
Colors[grep("Lip", rownames(PCE4YES))] <- "yellow"
Colors[grep("OS", rownames(PCE4YES))] <- "grey"
# httpgd::hgd( width = getOption("httpgd.width", 700),
#   height = getOption("httpgd.height", 700),)
set.seed(111213)
# Plot he sparsified ridge matrix of ApoE4 carriers with AD
Coords <- Ugraph(PCE4YES,
  type = "fancy", lay = "layout_with_fr",
  Vcolor = Colors, prune = FALSE, Vsize = 5, Vcex = 0.3, cut = 0.5,
  main = "ApoE4 carriers with AD"
)
# Plot the sparsified ridge matrix of ApoE4 non-carriers with AD
Ugraph(PCE4NO,
  type = "fancy", lay = NULL, coords = Coords,
  Vcolor = Colors, prune = FALSE, Vsize = 5, Vcex = 0.3, cut = 0.5,
  main = "ApoE4 non-carriers with AD"
)
# Plot the differential network
diffgraph <- DiffGraph(PCE4NO, PCE4YES,
  lay = NULL, coords = Coords,
  Vcolor = Colors, Vsize = 5, Vcex = 0.3, main = "Differential Network"
)
```
Distinct metabolite communities among ApoE4 carriers/non carriers
```{r communities}
# Get the communities per class
set.seed(141516)
CommC1 <- Communities(PCE4NO,
  Vcolor = Colors, Vsize = 5, Vcex = 0.3,
  main = "ApoE4 non-carriers with AD"
)

CommC2 <- Communities(PCE4YES,
  Vcolor = Colors, Vsize = 5, Vcex = 0.3,
  main = "ApoE4 carriers with AD"
)
```
```{r centrality}
PC0list <- list(PCE4NO = PCE4NO, PCE4YES = PCE4YES)
# Get the network statistics
NetStats <- GGMnetworkStats.fused(PC0list)

NetStatsE4yes <- NetStats[, 10:18]
NetStatsE4no <- NetStats[, 1:9]

NetStatsE4yes[, c(3, 9)] <- NetStatsE4no[, c(3, 9)] <- NULL
DegreesE4yes <- as.data.frame(NetStatsE4yes[, 1], rownames(NetStats))
DegreesE4no <- as.data.frame(NetStatsE4no[, 1], row.names = row.names(NetStatsE4no))
head(DegreesE4yes[order(DegreesE4yes[, 1], decreasing = TRUE), , drop = FALSE], 20)
head(DegreesE4no[order(DegreesE4no[, 1], decreasing = TRUE), , drop = FALSE], 20)
```
### Network centrality measures
The Wilcoxon Signed Rank tests revealed most network statistics to be larger in ApoE4 carriers compared
to non-carriers (p-values <0.05) with AD (Table 7). In the ApoE4 positive group, top central metabolites
were amines (dopamine and citrulline) (Table 8). ApoE4 non-carriers had 3-methoxytyrosine as top central
metabolite, followed by the oxidative stress compound nitro fatty acid (C18:3)
```{r wilcoxon rank sum }
# Perform a Wilcoxon signed rank test
w <- furrr::future_map2(NetStatsE4yes[, 1:7], NetStatsE4no[, 1:7], ~ {
  wilcox.test(.x, .y, paired = TRUE, alternative = "greater")
})
p.values <- names <- list(1:length(w))
for (i in 1:length(w)) {
  names[[i]] <- names(w[[i]])[7:length(names(w[[i]]))]
  p.values[[i]] <- w[[i]][["p.value"]]
}
names(p.values) <- names(w)
d <- as.data.frame(p.values)
d <- t(d)
d <- cbind.data.frame(d, p.adjust(d, method = "fdr"))
rownames(d) <- sub("PCE4YES.", "", row.names(d))
names(d) <- c("p-value", "FDR")
d
```