s5_3_signature_based_deconvolution_brain.Rmd

---
title: "Signature_based_deconvolution_of_brain_data"
output: html_notebook
---
# Repeat of partial figure Brain


## Imports



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

dir.create(file.path(".", dir_to_save_fig), showWarnings = F, recursive = T)


source('../R/setup.R') # import main package
source("../R/figure_utils.R") # for same colors
```



```{r}
library(CellMix) # to compare with
library(Seurat)
library(digest)
library(biomaRt)
library(tibble)
library(preprocessCore) # needed for CellMix/DSA
library(parallel) # needed for CellMix/DSA
library(e1071) # needed for CellMix/DSA
library(reshape2)
library(ggplot2)
library(ggrastr)
library(svglite)
source("../R/CIBERSORT.R") # to compare with
library(digest)
library(matrixStats)
library(Metrics)
library(dplyr)
```


## Load the full data and proportions
```{r}

mart <- useDataset("hsapiens_gene_ensembl", useMart("ensembl"))
ex <- readRDS(paste0(data_dir, "brain_merge2_grad.rds"))
props <- t(readRDS(paste0(data_dir, "brain_merge2_grad_true_props.rds")))
rownames(props)[rownames(props) == "OPC"] <- "OPCs"

```

## Load external markers
```{r}
## VL, NG, CA, ## possibly (LK, TS, F5, IP, DM, MM) but check column names carefully 
DS_NAME <- "VL"
n_markers <- 20
load(paste0(data_dir, "sigsBrain.rda"))
bas <- sigsBrain[[DS_NAME]]
```


```{r}
colnames(bas)
```
```{r}
print(paste("Dim for the whole data:", toString(dim(ex))))
print(paste("Dim for proportions (H):", toString(dim(props))))
print(paste("Dim for basis (W):", toString(dim(bas))))
```

## Change gene names, transform expression, get markers list

util to process markers, set gene names

```{r}
extract_markers_remove_not_present_genes <- function (ex, bas, props) {
  data_full <- ex
  genes <- rownames(bas)
  gene_names <- getBM(
  filters = "ensembl_gene_id",
  attributes = c("ensembl_gene_id", "hgnc_symbol"),
  values = genes,
  mart = mart)
  rownames(gene_names) <- gene_names$ensembl_gene_id
  new_rnames <- gene_names[rownames(bas), "hgnc_symbol"]
  
  # remove duplicated and not present in ensembl
  bas <- bas[!(is.na(new_rnames) | (new_rnames == "") | duplicated(new_rnames)),]
  rownames(bas) <- gene_names[rownames(bas), "hgnc_symbol"]
  # ensure basis have same columns as rows in props ( Cell types the same)
  bas <- bas[, rownames(props)]
  #select  highly variable genes  in total data (by mad)
  data_full <- data_full[log2(apply(data_full, 1, mad)) > 2, ]
  data_full <- data_full[rownames(data_full) %in% rownames(bas), ]
  # ensure basis have same rows as rows in expression data (Genes the same)
  bas <- bas[rownames(data_full),]
  cat(dim(data_full), "\n")
  colnames(bas) <- gsub(" |-", "_", colnames(bas))
  # Get markers by log fold change
  marker_list <- get_signature_markers(bas, n_marker_genes = 60)
  for (name in names(marker_list)) {
    marker_list[[name]] <- marker_list[[name]][marker_list[[name]] %in% rownames(data_full)]
  }
  # Remove duplicated markers
  duplicated_markers <- unlist(marker_list)[duplicated(unlist(marker_list))]
  # Ensure markers are in the data
  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
  }
  return (list(data_full=data_full, props=props, bas=bas, marker_list=marker_list, rev_mlist=rev_mlist))
}
```

```{r}

r <- extract_markers_remove_not_present_genes(ex, bas, props)
data_full <- r$data_full
props <- r$props
bas <- r$bas
marker_list <-  r$marker_list
rev_mlist <- r$rev_mlist
ntypes <- length(marker_list)
```
## Print dimensions of data 

```{r}
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))))
```


## Visualize mean expression of markers colored by 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)
```

## Make make reduced data (in original and projected space)
```{r}


data_reduced <- data_full[intersect(unlist(marker_list), rownames(data_full)), ]

temp_dso <-  DualSimplexSolver$new()
temp_dso$set_data(data_reduced)
data_reduced_scaled <- temp_dso$get_V_row()
proj_reduced <-  temp_dso$st$proj

temp_dso$set_data(data_full)
data_full_scaled <- temp_dso$get_V_row()
proj_full <-  temp_dso$st$proj



```


## Genes UMAP
Sometimes UMAP tend to place markers separetely, play siwh set.seed() with different numbers to get better picture
```{r}
# This option makes all dots as a single layer reducing result image size
options("dualsimplex-rasterize"=T)
```


```{r fig.width=10, fig.height=7}

data_normalized <- t(apply(data_full, 1, function(x) x / sum(x)))
pca <- prcomp(data_normalized)
plot(
  1:length(pca$sdev),
  cumsum(pca$sdev^2) / sum(pca$sdev^2),
  type = "b",
  xlab = "Principal Component",
  ylab = "Cumulative Variance Explained"
)

n_components <- 7
subset_pca <- pca$x[, 1:n_components]
```


```{r}
# Perform UMAP
set.seed(5)
um <- uwot::umap(subset_pca, ret_model = TRUE)
genes_umap <- um$embedding


```


```{r}
ct_labels <- rep(NA, length(rownames(genes_umap)))
names(ct_labels) <- rownames(genes_umap)
ct_labels[names(rev_mlist)] <- rev_mlist
ct_labels <- unlist(ct_labels[rownames(genes_umap)])

to_plot <- as.data.frame(genes_umap)
to_plot[["cell_type"]] <- brain_ct_map_long[ct_labels]
res_plot <-  ggplot(to_plot[is.na(ct_labels),], aes_string(x="V2", y="V1", col="cell_type")) +
    rasterise(geom_point(size=1.5), dpi=600) +
    scale_color_manual(values = colors_v, na.value = "grey70") + 
      rasterise(geom_point(data=to_plot[!is.na(ct_labels),], 
                           size=3, 
                           aes_string(col="cell_type")), dpi=600) + 
    labs(col = "Marker Cell Type", x = "UMAP1", y = "UMAP2") +
    theme_bw(base_family = "sans", base_size = 18) +
    theme(axis.ticks = element_blank(),
          axis.text = element_blank()) +  ggtitle("Genes PCA 1-7 UMAP, markers highlighted")

filename <- paste0(dir_to_save_fig,"7_2_1_partial_umap.svg")
ggsave(file=filename, plot=res_plot, width=8, height=6, device=svglite)

res_plot
```


## Create Dual Simplex object

```{r}
dso <- DualSimplexSolver$new()
dso$set_save_dir(objects_save_dir)
dso$set_data(data_full)
dso$project(ntypes)
dso$plot_svd_history()


```
## Check umap of projected points

```{r fig.width=10, fig.height=7}
set.seed(13)
dso$run_umap(neighbors_X = 15, neighbors_Omega = 20)
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())
```

## Single cell enrichment of markers

```{r}
load(paste0(data_dir, "SeuratObjects.rda"))
so <- obj$VL
rm(obj)
so <- ScaleData(so)
so <- RunPCA(so)
so <- RunUMAP(so, dims = 1:5)
```

```{r}
so$orig.celltype[so$orig.celltype == "AST-FB" | so$orig.celltype == "AST-PP"] <- "Astrocytes"
so$orig.celltype[startsWith(so$orig.celltype, "IN")] <- "Inhibitory"
so$orig.celltype[so$orig.celltype %in% c("L2/3", "L4", "L5/6", "L5/6-CC", "Neu-mat")] <- "Excitatory"
so$orig.celltype[so$orig.celltype %in% c("Neu-NRGN-I", "Neu-NRGN-II")] <- NA
```

### Write annotation to picture with rasterization
```{r fig.width=5, fig.height=5}
res_plot <- rasterise(DimPlot(so, group.by = "orig.celltype"), dpi=300) + 
  NoAxes() + 
  NoLegend() + 
  ggtitle(NULL) + 
  scale_colour_manual(values = alpha(colors_m, 1), na.value = alpha("grey70", 1))

filename <- paste0(dir_to_save_fig, "5_3_2_annotation_sc.svg")
ggsave(file=filename, plot=res_plot, width=2.5, height=2.5, device=svglite)
res_plot


```

### Write each individual cell type enrichment

```{r fig.height=5, fig.width=5}
ml_sorted <- marker_list[sort(names(marker_list))]
so <- add_list_markers(so, ml_sorted)

list_of_plots <- plot_marker_enrichment(so, names(ml_sorted), ncol=7, limits = c(-0.2, 1), function(plt,i) {
  rasterise(plt, dpi=300) + NoLegend() + ggtitle(NULL) + NoAxes()
}, wrap = F)


list_of_plots2 <- plot_marker_enrichment(so, names(ml_sorted), ncol=7, limits = c(-0.2, 1), function(plt,i) {
  rasterise(plt, dpi=300) + ggtitle(NULL) + NoAxes()
}, wrap = F)
```



```{r fig.width=5, fig.height=5}
names(list_of_plots) <- names(ml_sorted)
list_of_plots <- lapply(names(list_of_plots), function(current_name){
  res_plot <-  list_of_plots[[current_name]]
  filename <- paste0(dir_to_save_fig, "5_3_2_", current_name,"_sc.svg")

  ggsave(file=filename, plot=res_plot, width=1.5, height=1.5, device=svglite)
  res_plot
  

  return (res_plot)
})
names(list_of_plots2) <- names(ml_sorted)
list_of_plots2 <- lapply(names(list_of_plots2), function(current_name){
  res_plot <-  list_of_plots2[[current_name]]
  filename <- paste0(dir_to_save_fig, "5_3_2_", current_name,"_sc_with_legend.svg")

  ggsave(file=filename, plot=res_plot, width=1.5, height=1.5, device=svglite)
  res_plot
  return (res_plot)
})

list_of_plots
```
```{r}
names(list_of_plots2)
```

## Perform single run of signature-based deconvolution

```{r}
dso <- DualSimplexSolver$new()
dso$set_save_dir(objects_save_dir)
dso$set_data(data_reduced)
dso$project(ntypes)
dso$init_solution("marker_means", marker_list=marker_list)
dso$default_optimization()


## 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)
  
```

## Make projected plot

```{r}
dso$plot_projected()
```

## Compare result proportions


```{r fig.width=20, fig.height=5}
ptp <- coerce_pred_true_props(H_res, props)
plot_ptp_scatter(ptp)
```

```{r}
ptp_lines <-  plotProportions(ptp[[1]], 
                                       ptp[[2]], 
                                       pnames = c("predicted", "true"),
                                       point_size = 1,
                                       line_size = 0.7
                                       ) + theme_bw(base_size=12) + theme(legend.title = element_blank(), legend.position = "bottom", axis.title.x = element_blank(), 
    axis.text.x = element_text(angle = 45, hjust = 1))

filename <- paste0(dir_to_save_fig,"7_3_5_ptp_lines.svg")
ggsave(filename, ptp_lines, width=3.5, height=3, device=svglite)
ptp_lines
```


## Multiple runs of the same. Save plots

 Util methods to save proportion comparison
```{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)
}
```

## Do 5 independent runs on the same data 

```{r}

for (rep_n in seq(1, 5)) {
  print(rep_n)
  # Read the data
  ex <- readRDS(paste0(data_dir, "brain_merge2_grad.rds"))
  props <- t(readRDS(paste0(data_dir, "brain_merge2_grad_true_props.rds")))
  rownames(props)[rownames(props) == "OPC"] <- "OPCs"
  ## VL, NG, CA, LK, TS, F5, IP, DM, MM
  DS_NAME <- "CA"
  n_markers <- 20
  bas <- sigsBrain[[DS_NAME]]
  
  r <- extract_markers_remove_not_present_genes(ex, bas, props)
  data_full <- r$data_full
  props <- r$props
  bas <- r$bas
  marker_list <-  r$marker_list
  rev_mlist <- r$rev_mlist
  ntypes <- length(marker_list)
    
  
  # Work with reduced data below
  data_reduced <- data_full[intersect(unlist(marker_list), rownames(data_full)), ]
  
  ## Solve with Linseed2
  
  # Run linseed2
  dso <- DualSimplexSolver$new()
  dso$set_save_dir(objects_save_dir)
  dso$set_data(data_reduced)
  dso$project(ntypes)
  dso$init_solution("marker_means", marker_list=marker_list)
  dso$default_optimization()
  ## 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)
  print("CIBERSORT")
  # 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)
  
  
  print("DSA")
  # 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)

  }
```

## Make comparison plot for multiple methods

```{r fig.height=10, fig.width=11}
library(svglite)
#save_dir <- "./7_partial_figure_save/brain"
filename <- paste0(dir_to_save_fig, "7_3_3_barplot.svg")

models <- c("dualsimplex", "dsa"
            , "cib"
            )
signatures <- list(
  #"VL", "NG",
  "CA")
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) 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))])
rmses$cell_type <- factor(brain_ct_map_short[rmses$cell_type], levels = brain_ct_map_short)

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) +  # Adjust width as needed
  ylab("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") +
  facet_grid(rows = vars(signature), cols = vars(model))

ggsave(file=filename, plot=res_plot, width=11, height=10, device=svglite)
res_plot
```