s6_2_signature_based_deconvolution_GSE11058.Rmd

---
title: "R Notebook"
output: html_notebook
---

# Repeat of partial figure GSE11058
Here we assume that data is already generated (see another file)
You need to install CellMix to use dsa

# Tricky installation
Note: maximum R version supported by CellMix is 3.6. 
(docker for this https://hub.docker.com/r/aladyevae/cellmix)
Original package is not available for later version of R
But we found a fork of this repo that is working
```{r}
# Before: copy CellMix folder we provided in our repo to .libPaths() folder
if (!require("CellMix", quietly = TRUE))
  devtools::install_github('denklewer/cellmix')
library("CellMix")
library(digest)
library(matrixStats)
library(dplyr)
```



# Imports 

```{r}
# path to results and to save everything
dir_to_save_fig <-  "../out/signature_based_S6a/"
objects_save_dir <- "../out/signature_based_S6a/dso_objects"
data_dir <- "../data/large/signature_based/GSE11058_benchmark/"

dir.create(file.path(".", dir_to_save_fig), showWarnings = F, recursive = T)
dir.create(file.path(".", objects_save_dir), showWarnings = F, recursive = T)
source('../R/setup.R') # import main package
source("../R/figure_utils.R") # for same colors
```


```{r}
library(tibble)
library(preprocessCore)
library(parallel)
library(e1071)
library(reshape2)
library(ggplot2)
library(svglite)
source("../R/CIBERSORT.R")
```

## Download the data

```{r}

data_full <- readRDS(paste0(data_dir, "single_run_data_full.rds"))
bas <- readRDS(paste0(data_dir,"single_run_bas.rds"))
props <- readRDS(paste0(data_dir,"single_run_props.rds"))

print(paste("Dim for the whole data:", toString(dim(data_full))))
print(paste("Dim for proportions (H):", toString(dim(props))))
print(paste("Dim for basis (W):", toString(dim(bas))))
```

## Prepare markers
```{r}
DS_NAME <- "GSE11058"
n_markers <- 20
# extract markers from basis based on log fold change
marker_list <- get_signature_markers(bas, n_marker_genes = 60)
fc <- log(get_fold_change(bas))
for (name in names(marker_list)) {
  marker_list[[name]] <- marker_list[[name]][marker_list[[name]] %in% rownames(data_full)]
}
# remove intersected markers
duplicated_markers <- unlist(marker_list)[duplicated(unlist(marker_list))]
for (name in names(marker_list)) {
  marker_list[[name]] <- marker_list[[name]][
    (marker_list[[name]] %in% rownames(data_full)) &
      !(marker_list[[name]] %in% duplicated_markers)
  ]
  marker_list[[name]] <- marker_list[[name]][1:min(n_markers, length(marker_list[[name]]))]
}

for (name in names(marker_list)) {
  cat(name, ": ", length(marker_list[[name]]), "\n")
}

rev_mlist <- list()
for (name in names(marker_list)) {
  for (gene in marker_list[[name]])
    rev_mlist[[gene]] <- name
}

ntypes <- length(marker_list)
```

## Observe that markers are enriched in respective cell type 

```{r}
res <- as.data.frame(lapply(marker_list, function(x) apply(as.matrix(log2(bas[x,] + 1)), 2, mean)))
res$markers_for <- rownames(res)
res <- melt(res, id = "markers_for")
colnames(res) <- c("markers_for", "enrichment_in", "mean_log_expression")
res$markers_for <- as.factor(res$markers_for)
res$enrichment_in <- as.factor(res$enrichment_in)
ggplot(
  res,
  aes(x = enrichment_in, y = mean_log_expression, group = markers_for, col = markers_for)
) + geom_line() + geom_point() + scale_color_manual(values = colors_v)
```
# Prepare linseed object 

```{r}
data_reduced <- data_full[intersect(unlist(marker_list), rownames(data_full)), ]
print(paste("Dim for the full data:", toString(dim(data_full))))
print(paste("Dim for the reduced data:", toString(dim(data_reduced))))
```


```{r}
dso <- DualSimplexSolver$new()
dso$set_save_dir(objects_save_dir)
dso$set_data(data_full)
dso$project(4)
dso$plot_svd_history()
```
## This is how we usually do the optimization
Start in random solution points and optimize

```{r}
dso$init_solution("random")
dso$optim_solution(5000, 
                   optim_config(
                                coef_der_X = 0.001, 
                                coef_der_Omega = 0.001
                                ))
## Extract solution for markers in original space
solution <- dso$finalize_solution()
dso$plot_error_history()
```


# Work with predefined markers

## Projected picture
Here we just take prepared markers and plot them in projected space. 
Expectation is that these previously selected markers will be in the corners.

```{r fig.width=9, fig.height=7}
this_projection <- dso$st$proj
## Alternative methods to get full/reduced projections
#data_reduced_scaled <- sinkhorn_scale(data_reduced, 20)
#data_full_scaled <-  sinkhorn_scale(data_full, 20)
#proj_reduced <- svd_project(data_reduced_scaled, 1:ntypes)
#proj_full <- svd_project(data_full_scaled, 1:ntypes)


dso$st$marker_genes <- marker_list
colors <-  which_marker(rownames(fData(dso$st$data)), dso$st$marker_genes)

plot_projection_points(this_projection, 
                       spaces = c("X"), 
                       use_dims=c(2:3),
                       pt_size = 1.5,
                       pt_opacity = 0.7, 
                       color=colors) +
  scale_color_manual(values = colors_v, na.value = "grey60") + 
  labs(col = "Marker Cell Type", x = "R2", y = "R3") +
  theme_bw(base_family = "sans", base_size = 18) +
  theme(legend.position = "right", axis.ticks = element_blank(),
        axis.text = element_blank())
```
## 3D picture

```{r}
library(plotly)

dso$st$marker_genes <- marker_list
colors <-  which_marker(rownames(fData(dso$st$data)), dso$st$marker_genes)

font <- list(
  family = "Helvetica",
  size = 30)
ax_settings <- list(title = '', 
                    autorange = T, 
                    showgrid = T, 
                    zeroline = F, 
                    showline = T, 
                    autotick = T, 
                    ticks = '', 
                    showticklabels = F, 
                    titlefont = font)
to_plot <- as.data.frame(dso$st$proj$X)
to_plot$color <- as.character(colors)
to_plot$color[is.na(to_plot$color)] <- "grey70"
plotly::plot_ly(
  to_plot,
  x = ~dim_2,
  y = ~dim_3,
  z = ~dim_4,
  color = ~color,
  colors = c("grey70", colors_v),
  opacity = 0.9,
  size = 1,
  type = "scatter3d",
  mode = "markers"
) %>%
  plotly::layout(showlegend = FALSE,
    scene = list(
      xaxis = c(title = "R2", ax_settings),
      yaxis = c(title = "R3", ax_settings),
      zaxis = c(title = "R4", ax_settings)
    )
  )  %>%  config(toImageButtonOptions = list(format = "svg", width = 600,
                                                              height = 600))

```
## Projected picture using dso object method. (color by zero distance)
We can initialize solution in the coordinates of markers

```{r}
dso$init_solution("marker_means", marker_list=marker_list)
dso$plot_projected("zero_distance", "zero_distance", use_dims = list(3:4), with_solution = T)
```
## Make umap plot (only for visualization purposes)

```{r}
dso$run_umap(neighbors_X = 20, neighbors_Omega = 20)
```

```{r}
dso$plot_projected("zero_distance", "zero_distance", use_dims = list(NULL,3:4), with_solution = T)
```


```{r fig.width=9, fig.height=7}
dso$st$marker_genes <- marker_list
colors <-  which_marker(rownames(fData(dso$st$data)), dso$st$marker_genes)
plot_projection_points(dso$st$proj,  spaces = c("X"), pt_size = 2, color=colors) +
  scale_color_manual(values = colors_v, na.value = adjustcolor("grey70", alpha.f = 0.7)) + 
  labs(col = "Marker Cell Type", x = "UMAP1", y = "UMAP2") +
  theme_bw(base_family = "sans", base_size = 18) +
  theme(legend.position = "right", axis.ticks = element_blank(),
        axis.text = element_blank())
```



## Multiple runs

```{r}
save_pt <- function(
  this_H, this_W, model_n, rep_n, dir = objects_save_dir, other_H = props, other_W = bas
) {
  other_H <- other_H[names(marker_list), ]
  other_W <- other_W[, names(marker_list)]
  ptp <- coerce_pred_true_props(this_H, other_H)
  dir.create(file.path(dir, model_n))
  saveRDS(ptp, paste0(dir, "/", model_n , "/", rep_n, "_ptp_", DS_NAME, ".rds"))
  
  if (!is.null(this_W)) {
    ptb <- coerce_pred_true_basis(this_W, other_W)
    # ptb[[1]] <- t(ptb[[1]])
    # ptb[[2]] <- t(ptb[[2]])
    # ptb[[1]][ptb[[1]] < 0] <- 0
    # ptb[[1]] <- log2(ptb[[1]] + 1)
    # ptb[[2]] <- log2(ptb[[2]] + 1)
    
    saveRDS(ptb, paste0(dir, "/", model_n, "/", rep_n, "_ptb_", DS_NAME, ".rds"))
  }
}

read_pt <- function(model_n, rep_n, sig = DS_NAME, dir = objects_save_dir) {
  ptp_name <- paste0(dir, "/", model_n , "/", rep_n, "_ptp_", sig, ".rds")
  ptp <- readRDS(ptp_name)
  ptb_name <- paste0(dir, "/", model_n, "/", rep_n, "_ptb_", sig, ".rds")
  ptb <- if (file.exists(ptb_name)&& !dir.exists(ptb_name)) readRDS(ptb_name) else NULL
  list(ptp = ptp, ptb = ptb)
}
```




```{r fig.width=20, fig.height=5}
for (rep_n in seq(1, 5)) {
  # Read data
  data_full <- readRDS(paste0(data_dir,rep_n,"_run_data_full.rds")) 
  bas <- readRDS(paste0(data_dir, rep_n, "_run_bas.rds"))
  props <- readRDS(paste0(data_dir, rep_n, "_run_props.rds"))
  
  # Get Signature markers from 'bas'
  marker_list <- get_signature_markers(bas, n_marker_genes = 60)
    fc <- log(get_fold_change(bas))
  for (name in names(marker_list)) {
    marker_list[[name]] <- marker_list[[name]][marker_list[[name]] %in% rownames(data_full)]
  }
  duplicated_markers <- unlist(marker_list)[duplicated(unlist(marker_list))]
  for (name in names(marker_list)) {
    marker_list[[name]] <- marker_list[[name]][
      (marker_list[[name]] %in% rownames(data_full)) &
        !(marker_list[[name]] %in% duplicated_markers)
    ]
    marker_list[[name]] <- marker_list[[name]][1:min(n_markers, length(marker_list[[name]]))]
  }
  for (name in names(marker_list)) {
    cat(name, ": ", length(marker_list[[name]]), "\n")
  }
  rev_mlist <- list()
  for (name in names(marker_list)) {
    for (gene in marker_list[[name]])
      rev_mlist[[gene]] <- name
  }
  
  ntypes <- length(marker_list)
  # Get reduced data (only marker rows present)
  data <- data_full[intersect(unlist(marker_list), rownames(data_full)), ]
  
  # Run linseed2
  dso <- DualSimplexSolver$new()
  dso$set_save_dir(objects_save_dir)
  dso$set_data(data)
  dso$project(4)
  dso$init_solution("marker_means", marker_list=marker_list)
  dso$optim_solution(3, 
                     optim_config(coef_hinge_W = 0.0001,
                                  coef_der_X = 0, 
                                  coef_der_Omega = 0
                                  ))
  ## Extract solution for markers in original space
  solution <- dso$finalize_solution()
  H_res <- solution$H
  W_mark <- solution$W
  X_mark <- dso$st$solution_proj$X  
  Omega_mark <- t(dso$st$solution_proj$Omega)
  
  ## Get W for the whole data based on proportions
  res <- nnls_C__(t(H_res), t(data_full))
  W_full <- t(res)
  colnames(W_full) <- paste0("Cell_type_", 1:ntypes)
  
  ## Save result
  save_pt(H_res, W_full, "dualsimplex", rep_n, dir=objects_save_dir)
  
  # Run CIBERSORT
  cib <- CIBERSORT(
    as.data.frame(bas[rownames(data_full), ]),
    as.data.frame(data_full),
    perm = 0,
    QN = T,
    absolute = F
    )
  H_res <- t(cib)
  H_res <- H_res[names(marker_list), ]
  save_pt(H_res, NULL, "cib", rep_n,  dir=objects_save_dir)
  
  # Run DSA (comment if no CellMix installed)
  dsa <- CellMix::ged(as.matrix(data_full), 
                      x = CellMix::MarkerList(marker_list), 
                      "DSA", 
                      verbose = F) 
  H_res <- coef(dsa)
  W_full <- basis(dsa)
  save_pt(H_res, W_full, "dsa", rep_n,  dir=objects_save_dir)
  
}
```


##  Read results and produce plots 

```{r fig.height=7, fig.width=5}
models <- c("dualsimplex", "cib", "dsa")
signatures <- list("GSE11058")

total_rmses <- data.frame()
for (rep_n in seq(1, 5)) {
  this_dfs <- lapply(models, function(model) {
    results <- lapply(signatures, function(sig) {
      read_pt(model, rep_n, sig, dir=objects_save_dir)$ptp
    })
    rmses <- sapply(results, function(res) {
      sapply(rownames(res[[1]]), function(ct) Metrics::rmse(res[[1]][ct, ], res[[2]][ct, ]))
    })
    colnames(rmses) <- signatures
    rmses <- melt(rmses)
    colnames(rmses) <- c("cell_type", "signature", "rmse")
    if (model == "cib") {
      rmses$model <- "CIBERSORT"
    } else if (model == "dsa") {
      rmses$model <- "DSA"
    } else if (model == "dualsimplex") {
      rmses$model <- "DualSimplex"
    }
    rmses$rep <- rep_n
    rmses
  })
  total_rmses <- rbind(total_rmses, do.call(rbind, this_dfs))
}
rmses <- total_rmses
rmses$cell_type <- factor(rmses$cell_type, levels = levels(rmses$cell_type)[order(levels(rmses$cell_type))])

agg_data <- rmses %>%
  group_by(cell_type, signature, model) %>%
  summarize(mean_rmse = mean(rmse), 
            sd_rmse = sd(rmse), 
            n = n())

res_plot <- ggplot(agg_data, aes(x = cell_type, y = mean_rmse, fill = cell_type)) +
  geom_bar(stat = "identity", position = position_dodge()) +
  geom_errorbar(aes(ymin = mean_rmse - sd_rmse/sqrt(n), 
                    ymax = mean_rmse + sd_rmse/sqrt(n)), 
                position = position_dodge(width = 0.9), 
                width = 0.25) +
  labs(y = label_wrap_gen(width = 25)("Mean RMSE, true vs. predicted proportions")) +
  xlab("Cell type") +
  scale_fill_manual(values = colors_m) +
  theme_bw(base_family = "sans", base_size = 18) +
  theme(legend.position = "none", strip.text.y = element_text(angle = 0), strip.background = element_blank()) +
  facet_grid(rows = vars(model))


filename <- paste0(dir_to_save_fig,"7_2_2_rmse_barplot.svg")

ggsave(file=filename, plot=res_plot, width=5, height=7, device=svg)
res_plot
```
```{r}
res <- read_pt("dualsimplex", 2,dir=objects_save_dir,)
ptp <- res$ptp
ptb <- res$ptb
```

```{r}
library(grid)
library(gridExtra)
library(svglite)
```

```{r, fig.width=20, fig.height=5}
options("dualsimplex-rasterize" = FALSE)

proportions_plot <- plot_ptp_scatter(ptp)
filename <- paste0(dir_to_save_fig, "7_2_3_proportions_plot.svg")

ggsave(file=filename, plot=proportions_plot, width=20, height=5, device=svglite)


```

```{r, fig.width=20, fig.height=5}
options("dualsimplex-rasterize" = TRUE)
expressions_plot <- plot_ptb_scatter(ptb, max_expr = 18) 

filename <- paste0(dir_to_save_fig, "7_2_4_expressions_plot.svg") 

ggsave(file=filename, plot=expressions_plot, width=16, height=4, device=svglite)

```