NARES v. recount2 LV comparsion

J. Taroni 2018

Our goal is to map LVs between different PLIER models based on the gene loadings (Z matrix). The rationale is that LVs with the most similar gene loadings map to each other. In reality, a clear one-to-one mapping may not exist. Here, we’re looking the NARES PLIER model and the recount2 PLIER model.

We’d like to examine whether or not this transfer learning approach (e.g., using the recount2 PLIER model with the NARES data) is likely to give us very different LV (B matrix) values. If we get very different values, that should give us pause about using the transfer learning approach (e.g., multi-PLIER). We’re particularly concerned with LVs we consider to be “interpretable”, i.e., they have some non-random association with a pathway. For instance, do the trends we see in the “Neutrophil LV” from the NARES model hold true in the “Neutrophil” LV from the recount2 model?

Functions and directory setup

`%>%` <- dplyr::`%>%`
source(file.path("util", "plier_util.R"))

# plot and result directory setup for this notebook
plot.dir <- file.path("plots", "13")
dir.create(plot.dir, recursive = TRUE, showWarnings = FALSE)
results.dir <- file.path("results", "13")
dir.create(results.dir, recursive = TRUE, showWarnings = FALSE)

Load PLIER models, Z matrices

NARES PLIER

# Load NARES PLIER model
nares.plier.file <- file.path("results", "12", "NARES_PLIER_model.RDS")
nares.plier <- readRDS(nares.plier.file)
nares.summary <- nares.plier$summary

recount2 PLIER

# recount2 model
# load recount2 PLIER model, rename
recount.plier <- readRDS(file.path("data", "recount2_PLIER_data", 
                                   "recount_PLIER_model.RDS"))

Z matrices

# get Z matrices from both models
recount.z <- as.data.frame(recount.plier$Z)
nares.z <- as.data.frame(nares.plier$Z)
# we'll need to add the gene identifiers (symbols in this case) to a column
# rather than as rownames -- this will facilitate joining
recount.z <- tibble::rownames_to_column(recount.z, var = "Gene")
colnames(recount.z)[2:ncol(recount.z)] <- paste0("recountLV", 
                                                 1:(ncol(recount.z) - 1))
nares.z <- tibble::rownames_to_column(nares.z, var = "Gene")
colnames(nares.z)[2:ncol(nares.z)] <- paste0("naresLV", 
                                                 1:(ncol(nares.z) - 1))
# join -- only genes present in both models
z.df <- dplyr::inner_join(recount.z, nares.z, by = "Gene")
# need matrix to calculate correlation
z.matrix <- as.matrix(z.df[, 2:ncol(z.df)])
rownames(z.matrix) <- z.df$Gene

head(z.df)

Calculate correlation between the Z matrices

# calculate pearson correlation between LVs -- can we map between models using
# this distance metric?
cor.z.mat <- cor(z.matrix, method = "pearson")
# set diagonal to 0 to help find max correlation between LVs
diag(cor.z.mat) <- 0
# indices for each model
nares.indx <- grep("nares", rownames(cor.z.mat))
recount.indx <- grep("recount", rownames(cor.z.mat))
# pertinent indices
impt.cor.z.mat <- cor.z.mat[recount.indx, nares.indx] 
# for each NARES model LV, want the highest correlated LV from recount
mapping.df <- reshape2::melt(impt.cor.z.mat, 
                              varnames = c("recount_LV", "NARES_LV"),
                              value.name = "pearson_Z") %>%
                dplyr::group_by(NARES_LV) %>%
                dplyr::top_n(1, pearson_Z)
mapping.df

NARES data in recount2 latent space

We’ll need to calculate this to compare B matrices.

# get NARES expression matrix
exprs.file <- file.path("data", "expression_data", 
                        "NARES_SCANfast_ComBat_with_GeneSymbol.pcl")
exprs.df <- data.table::fread(exprs.file, data.table = FALSE)
exprs.mat <- dplyr::select(exprs.df, -(EntrezID:GeneSymbol))
rownames(exprs.mat) <- exprs.df$GeneSymbol
rm(exprs.df)

recount.nares.b <- GetNewDataB(exprs.mat = as.matrix(exprs.mat),
                               plier.model = recount.plier)

Loading required package: PLIER
Loading required package: RColorBrewer
Loading required package: gplots

Attaching package: ‘gplots’

The following object is masked from ‘package:stats’:

    lowess

Loading required package: pheatmap
Loading required package: glmnet
Loading required package: Matrix
Loading required package: foreach
Loaded glmnet 2.0-13

Loading required package: knitr
Loading required package: rsvd
Loading required package: qvalue

# save to file
recount.b.file <- file.path(results.dir, "NARES_recount2_B.RDS")
saveRDS(recount.nares.b, file = recount.b.file)

Compare B matrices

# get the NARES B matrices from both models, relabel the latent variables to
# indicate which model they are from
nares.b <- nares.plier$B
rownames(nares.b) <- paste0("naresLV", 1:nrow(nares.b))
rownames(recount.nares.b) <- paste0("recountLV", 1:nrow(recount.nares.b))
# correlation
cor.b.mat <- cor(t(recount.nares.b), t(nares.b))
# melt
cor.b.df <- reshape2::melt(cor.b.mat, varnames = c("recount_LV", "NARES_LV"),
                           value.name = "pearson_B")
mapping.df <- dplyr::inner_join(mapping.df, cor.b.df, 
                                by = c("recount_LV" = "recount_LV", 
                                       "NARES_LV" = "NARES_LV"))
mapping.df

Plotting results

First, let’s identify which of the NARES LVs are significantly associated with pathways. We’ll want to highlight these in our plots, as we expect that these are most likely to be “preserved” between models because they are related to the input gene sets (which were constant between models).

sig.summary <- nares.summary %>%
                dplyr::filter(FDR < 0.05)
sig.summary

sig.lvs <- as.integer(unique(sig.summary$`LV index`))

Heatmap of correlation values

# annotation bar indicating significant association with a pathway
row.color.bar <- rep("#FFFFFF", nrow(mapping.df))
row.color.bar[sig.lvs] <- "#0000FF"

cor.map.mat <- as.matrix(mapping.df[, c("pearson_Z", "pearson_B")])
rownames(cor.map.mat) <- mapping.df$NARES_LV
gplots::heatmap.2(cor.map.mat,
                  Rowv = NA, Colv = NA, trace = "none", dendrogram = "none",
                  col = colorRampPalette(colors = c("#FFFFFF", "#FF0000")),
                  cellnote = round(x = cor.map.mat, 2), notecol = "#000000",
                  notecex = 0.5, margins = c(5, 9),
                  cexCol = 0.75, cexRow = 0.7,
                  RowSideColors = row.color.bar, key = FALSE,
                  labRow = sub("nares", "NARES ", rownames(cor.map.mat)),
                  labCol = c("Z", "B"),
                  srtCol = 45)
legend("topleft", 
       legend = c("significant"),
       col = c("#0000FF"), cex = 0.75,
       lty = 1, lwd = 6)

# save plot to file
pdf(file.path(plot.dir, "NARES_Z_B_correlation_heatmap.pdf"), width = 5,
    height = 7)
gplots::heatmap.2(cor.map.mat,
                  Rowv = NA, Colv = NA, trace = "none", dendrogram = "none",
                  col = colorRampPalette(colors = c("#FFFFFF", "#FF0000")),
                  cellnote = round(x = cor.map.mat, 2), notecol = "#000000",
                  notecex = 0.5, margins = c(5, 9),
                  cexCol = 0.75, cexRow = 0.7,
                  RowSideColors = row.color.bar, key = FALSE,
                  labRow = sub("nares", "NARES ", rownames(cor.map.mat)),
                  labCol = c("Z", "B"))
legend("topleft", 
       legend = c("significant"),
       col = c("#0000FF"), cex = 0.75,
       lty = 1, lwd = 6)
dev.off()

null device 
          1 

Density plots

Let’s look at the entire distribution of cor.b.mat. It’s possible that the correlation of B between the “best match” (as determined by Z matrix correlation) LVs is not particularly notable.

# add LV type information
cor.b.df$LV_type <- rep("all", nrow(cor.b.df))

# add best match information
best.match.df <- dplyr::select(mapping.df, -pearson_Z) %>%
                  dplyr::mutate(LV_type = "best match")
# we'll add points for the significant lvs in blue
best.match.df$Significance <- NA
best.match.df$Significance[sig.lvs] <- "Significant"
b.cor.df <- dplyr::bind_rows(cor.b.df, best.match.df) 

b.cor.df %>%
  ggplot2::ggplot(ggplot2::aes(x = pearson_B,
                               group = LV_type,
                               fill = LV_type)) +
  ggplot2::geom_density(alpha = 0.3) +
  ggplot2::scale_fill_manual(values = c("#FFFFFF", "#5E5E5E")) +
  ggplot2::theme_bw() +
  ggplot2::geom_point(ggplot2::aes(y = 0, colour = Significance),
                      shape = 15, size = 3, alpha = 0.7) +
  ggplot2::labs(title = "Latent Space Agreement",
                subtitle = "NARES-MultiPLIER correlation",
                x = "Correlation") +
  ggplot2::theme(text = ggplot2::element_text(size = 15),
                 plot.title = ggplot2::element_text(hjust = 0.5, 
                                                    face = "bold"),
                 plot.subtitle = ggplot2::element_text(hjust = 0.5)) +
  ggplot2::scale_colour_manual(breaks = "Significant", values = c("#0000FF")) +
  ggplot2::guides(fill = ggplot2::guide_legend(override.aes = list(shape = NA),
                                               title = "LV Type"))

plot.file <- file.path(plot.dir, "NARES_B_correlation_density.pdf")
ggplot2::ggsave(plot.file, plot = ggplot2::last_plot())

Saving 7 x 7 in image

Summary

Some ambiguity in mapping latent variables across models suggests that the transfer learning using a model trained on gene expression data from diverse biological conditions (e.g., recount2 data) could be valuable for meta-analyses where curating a large enough compendium specific to one’s biological domain of interest is infeasible (e.g., in ANCA-associated vasculitis).

This analysis suggests that LVs that can be mapped across models (e.g., “best match”) have reasonably correlated sample values, much more so than random pairs of LVs. Pathway-associated LVs are particularly well-correlated.

The neutrophil (NARES LV3) and ECM (NARES LV14) LVs results from the NARES model suggests that these are among the strongest signals in this expression data (and this makes a lot of biological sense given the study at hand). These LVs had correlation values > 0.9 with their best matches in the recount2 PLIER model.

NARES LV2 values have no association with their best match in the recount2 model. The gene sets significantly associated with NARES LV2 are macrophage, dendritic cell, and mast cell gene sets. This is a broad swath of the myeloid lineage, and may suggest that the lack of specificity contributes to the poor correlation.