7a_gse19830.Rmd

---
title: "Figure 6B: Deconvolution of GSE19830"
output: html_document
---

```{r}
source("../R/figure_utils.R")
source("../R/setup.R")
dir_to_save_fig <- "../out/GSE19830/"
dir.create(file.path(".", dir_to_save_fig), showWarnings = F, recursive = T)
```


## Dependencies
```{r}
if (!require("rat2302.db", quietly = TRUE))
  BiocManager::install("rat2302.db")
if (!require("hgu133plus2.db", quietly = TRUE))
  BiocManager::install("hgu133plus2.db")
if (!require("ComplexHeatmap", quietly = TRUE))
  BiocManager::install("ComplexHeatmap")
if (!require("nnls", quietly = TRUE))
  install.packages("nnls")
if (!require("irlba", quietly = TRUE))
  install.packages("irlba")
if (!require("reshape", quietly = TRUE))
  install.packages("reshape")
if (!require("linseed", quietly = TRUE))
  install.packages("linseed")
if (!require("linseed", quietly = TRUE))
  devtools::install_github("ctlab/linseed")

library('GEOquery')
library(biomaRt)
library(hgu133plus2.db)
library(ComplexHeatmap)
library(progress)
library(reshape)
library(plotly)
library(ggrastr)
library(ggpubr)
library(linseed)
library(matrixStats)
library(ComplexHeatmap)
```


# Prepare Data
## Get it from GEO
```{r}
n_ct <- 3
dataset <- "GSE19830"
gse <- getGEO(dataset, AnnotGPL = T)
gse <- gse[[1]]
data_raw <- Biobase::exprs(gse)
#linearize the data!
data_raw <- linearize_dataset(data_raw)
```

## Map probe IDs to genes
probe id will be converted and collapsed to gene names
```{r}
# probes which are present in data
probe_mapping <- biomaRt::select(rat2302.db::rat2302.db, rownames(data_raw), c("SYMBOL"))
probe_mapping <-
  probe_mapping[probe_mapping$PROBEID %in% rownames(data_raw),]
# not empty result gene id (now we  many-to-many in in genes to probes)
probe_mapping <- probe_mapping[!is.na(probe_mapping$SYMBOL),]
# This dataset is strange, 1 probe could be many genes, which is not correct (but you can try without this)
probe_mapping <- probe_mapping[!duplicated(probe_mapping$PROBEID),]
## want to collapse gene for each gene and take mean value of data
data_raw <- data_raw[probe_mapping$PROBEID, ]
rownames(data_raw) <- probe_mapping$SYMBOL
# should be unique gene name as row name after this
data_raw <- tapply(data_raw, list(row.names(data_raw)[row(data_raw)], colnames(data_raw)[col(data_raw)]), FUN = median)
```


## Extract true proportions (H) and basis (W) from the data
Info about true samples is included for this dataset. We use this info for validation.

```{r}
# H
component_names <- c("Liver", "Brain", "Lung")
pdata <- pData(gse)
parsed_proportions <- strsplit(pdata$source_name_ch1,split = '/')
parsed_proportions <- lapply(parsed_proportions, function(sample_props){
  return(strtoi(gsub("\\D*(\\d+)\\D*","\\1",sample_props)))
})
true_proportions <- matrix(unlist(parsed_proportions), ncol=3, byrow=TRUE)/100
rownames(true_proportions) <- rownames(pdata)
colnames(true_proportions) <- component_names
true_proportions <- t(true_proportions)

# W
components <- lapply(c(1:n_ct), function(comp_num){
  component_subset_columns <- colnames(true_proportions[,true_proportions[comp_num, ] == 1])
  component_subset <- data_raw[, component_subset_columns]
  component_vector <- as.matrix(rowMedians(component_subset))
  rownames(component_vector) <- rownames(component_subset)
  return(component_vector)
  
})
true_basis <- do.call(cbind, components)
colnames(true_basis) <- component_names
```


## Remove pure samples
(Optional, to make the task harder)

```{r}
#sample_is_mixed <- apply(true_proportions, 2, mad) != 0
#true_proportions <- true_proportions[, sample_is_mixed]
#data_raw <- data_raw[, sample_is_mixed]
```


## Print data dimensions
```{r}
print(paste("Dim for the whole data:", toString(dim(data_raw))))
print(paste("Dim for hidden proportions (H):", toString(dim(true_proportions))))
print(paste("Dim for hidden basis (W):", toString(dim(true_basis))))
```


# Dual simplex method
## Create DualSimplexSolver object and prefilter

keep coding genes only

```{r}
dso <- DualSimplexSolver$new()
dso$set_data(data_raw)
dso$basic_filter(remove_true_cols_default = c())
dso$project(n_ct)

print(paste("Dim for the prefiltered data data:", toString(dim(dso$get_data()))))
```

```{r fig.height = 5, fig.width = 11}
dso$plot_svd(dims=1:15)
dso$plot_projected("zero_distance", "zero_distance", use_dims = list(2:3))
```


## Additional filters

### MAD filtering

check distribution of median absolute deviation to see which genes will be deleted

```{r}
dummy_threshold <- 1
data <- dso$get_data()
anno <- get_anno(data)
anno$PASS_FILTER <- FALSE
anno[anno$log_mad > dummy_threshold,]$PASS_FILTER <- TRUE
data <- set_anno(anno, data)
plot_feature(data, feature = "log_mad", col_by ='PASS_FILTER')
```

Filter the data with selected threshold
```{r fig.height = 4, fig.width=14}
dso$basic_filter(log_mad_gt = dummy_threshold)
dso$project(3)

svd_plot_2 <- dso$plot_svd(1:10) + theme_minimal(base_size = 8)
points_2 <- dso$plot_projected("black", "black",use_dims = 2:3, show_plots=F) +  theme_minimal(base_size = 8)
plotlist = list(points_2, svd_plot_2)
cowplot::plot_grid(plotlist=plotlist, rel_widths = c(0.66, 0.33)) 
```

Repeat previous cell with increasing thresholds to observe how SVD singular values will change. (in a perfect case should see only n_ct non-zero singular values)



### Denoising the data

Plot distance of genes from the projection hyperplane
```{r fig.height = 4, fig.width=14}
distance_dist_plot_noizy <-
  plot_feature(dso$get_data(), feature = "plane_distance") + theme_minimal(base_size = 15)
color = fData(dso$get_data())$plane_distance
X_points_noizy <-
  plot_projection_points(
    dso$st$proj,
    use_dims = 2:3,
    color = color,
    color_name = "plane_distance",
    spaces = c("X"),
    pt_size = 2
  ) + theme_minimal(base_size = 15)

distances_distribution <-
  plot_feature_pair(dso$get_data(), "plane_distance", "zero_distance", T, size = 0.1)


plotlist = list(X_points_noizy,
                distance_dist_plot_noizy,
                distances_distribution)

cowplot::plot_grid(
  plotlist = plotlist,
  nrow = 1,
  rel_widths = c(0.4, 0.3, 0.3)
) 
```

Remove genes which are to far from projection hyperplane
```{r}
dso$distance_filter(plane_d_lt = 0.015, zero_d_lt = NULL, genes = T)
dso$project(n_ct)
```


Color points by distance
```{r fig.height = 4, fig.width=14}
current_svd_plot <- dso$plot_svd(1:10) + theme_minimal(base_size = 8)
current_points_plot <-
  dso$plot_projected("plane_distance",
                     "plane_distance",
                     use_dims = 2:3,
                     show_plots = F) +  theme_minimal(base_size = 8)
plotlist = list(current_points_plot, current_svd_plot)
cowplot::plot_grid(plotlist = plotlist, rel_widths = c(0.66, 0.33)) 
```

Repeat previous cell with decreasing thresholds to observe how SVD singular values will change. (in a perfect case should see only n_ct non-zero singular values)

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


## Finaly train the model
### Initialze random points
```{r}
set.seed(1)
dso$init_solution("random")
dso$plot_projected(
  "zero_distance",
  "zero_distance",
  with_solution = TRUE,
  use_dims = list(2:3)
)
```

### Train multiple times with changing hyperparameters

We make learning smaller for later iterations We increase contribution of term related to positivity constraint for later iterations

```{r}
start.time <- Sys.time()

LR_DECAY_STEPS = 2
PARAMETERS_INCREASE_STEPS = 3
lr_decay <- 0.5
params_increase <- 2
original_lambda_term <- 0.1  #coef_hinge_H
original_beta_term <- 1 #coef_hinge_W
lr_x <- 0.01
lr_omega <- 0.01

RUNS_EACH_STEP <- 2000
pb <-
  progress_bar$new(total = LR_DECAY_STEPS * PARAMETERS_INCREASE_STEPS)

for (lr_step in 1:LR_DECAY_STEPS) {
  lambda_term <- original_lambda_term
  beta_term <- original_beta_term
  for (x in 1:PARAMETERS_INCREASE_STEPS) {
    # Main training method, you can just run this
    dso$optim_solution(
      RUNS_EACH_STEP,
      optim_config(
        coef_hinge_H = lambda_term,
        coef_hinge_W = beta_term,
        coef_der_X = lr_x,
        coef_der_Omega = lr_omega
      )
    )
    lambda_term <- lambda_term * params_increase
    beta_term <- beta_term * params_increase
    pb$tick()
  }
  lr_x <- lr_x * lr_decay
  lr_omega <- lr_omega * lr_decay
  
}

end.time <- Sys.time()
time.taken <- round(end.time - start.time, 2)
print(paste("Total time", time.taken))
```


### Plot errors
```{r}
plot_errors(dso$st$solution_proj, variables = c(
    "deconv_error",
    "lamdba_error",
    "beta_error",
    "total_error"
  ))
```

## Prepare plots as they are in the paper
```{r}
# This option makes all dots as a single layer reducing result image size
options("dualsimplex-rasterize"=T)
```

### History of solution
```{r fig.width=5, fig.height=5}
colnames(dso$st$proj$X) <- c("R1", "R2", "R3")
colnames(dso$st$proj$Omega) <- c("S1", "S2", "S3")
proj_solution_history <-
  dso$plot_projected(use_dims = (2:3),
                     wrap = F,
                     with_legend = F)
proj_solution_history[[1]] # genes in a space of samples
proj_solution_history[[2]] # samples in a space of genes
```

Decorate it how it was done in the paper

```{r fig.width=10, fig.height=5}
pltX <- proj_solution_history[[1]] +
  scale_fill_manual(values = colors_v[3:5]) +
  theme_minimal(base_size = 12) +
  theme(legend.position = "none")

filename <-
  paste0(dir_to_save_fig, "proj_solution_history_X", ".svg")

ggsave(
  file = filename,
  plot = pltX,
  width = 3,
  height = 3,
  device = svg
)
pltOmega <-
  proj_solution_history[[2]] + scale_fill_manual(values = colors_v[3:5])
filename <-
  paste0(dir_to_save_fig, "proj_solution_history_Omega", ".svg")
ggsave(
  file = filename,
  plot = pltOmega,
  width = 3,
  height = 3,
  device = svg
)
ggarrange(pltX, pltOmega)
```

## Solution check plots

### Get solution (W and H)
```{r}
solution <- dso$finalize_solution()
```

### Compare with true proportions/basis

Check ratio of true proportion or component expression to predicted
```{r fig.width=20, fig.height=5}
cur_H <- solution$H

norm_res_list <- lapply(1:dim(cur_H)[[2]], function(col_ind) {
  cur_col <- cur_H[, col_ind]
  normalized_column <-
    (cur_col - min(cur_col)) / (max(cur_col) - min(cur_col))
  normalized_column <- cur_col / sum(cur_col)
  return(normalized_column)
})
normalized_H <- do.call(cbind, norm_res_list)

ptp <- coerce_pred_true_props(normalized_H, true_proportions)
plot_ptp_scatter(ptp)
ptb <-
  coerce_pred_true_basis(solution$W, true_basis[rownames(solution$W),])
plot_ptb_scatter(ptb)
```

Visualize the same as lines
```{r fig.height = 5, fig.width = 7}
ptp_lines <- linseed::plotProportions(
  as.data.frame(ptp[[1]]),
  as.data.frame(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, "ptp_lines.svg")
ggsave(
  filename,
  ptp_lines,
  width = 3.5,
  height = 3,
  device = svg
)
ptp_lines
```


### Select and visualize marker genes
```{r fig.height = 5, fig.width = 11}
markers  <- dso$get_marker_genes()

to_plot <- data_raw[
  unlist(markers),
  apply(
    true_proportions,
    2,
    function(x) any(x == 1)
  )
]

dso$plot_projected(rownames(to_plot), use_dims = (2:3))
```


### More decorated markers plot (as it is in the paper)
```{r fig.height = 5, fig.width = 5}
colors <-
  which_marker(rownames(fData(dso$st$data)), dso$st$marker_genes)
plot_markers <-
  plot_projection_points(
    dso$st$proj,
    use_dims = (2:3),
    spaces = c("X"),
    pt_size = 1,
    color = colors
  ) +
  scale_color_manual(values = colors_v[3:5],
                     na.value = adjustcolor("grey70", alpha.f = 0.7)) +
  labs(col = "Marker Cell Type", x = "R2", y = "R3") +
  theme_bw(base_family = "sans", base_size = 12) +
  theme(
    legend.position = "right",
    axis.ticks = element_blank(),
    axis.text = element_blank()
  )

filename <- paste0(dir_to_save_fig, "marker_highlight_with_legend.svg")
ggsave(
  filename,
  plot_markers,
  width = 4.5,
  height = 3,
  device = svg
)

colors <- which_marker(rownames(fData(dso$st$data)), dso$st$marker_genes)
plot_markers <-
  plot_projection_points(
    dso$st$proj,
    use_dims = (2:3),
    spaces = c("X"),
    pt_size = 1,
    color = colors
  ) +
  scale_color_manual(values = colors_v[3:5],
                     na.value = adjustcolor("grey70", alpha.f = 0.7)) +
  labs(col = "Marker Cell Type", x = "R2", y = "R3") +
  theme_bw(base_family = "sans", base_size = 12) +
  theme(
    legend.position = "none",
    axis.ticks = element_blank(),
    axis.text = element_blank()
  )

filename <- paste0(dir_to_save_fig, "marker_highlight_without_legend.svg")
ggsave(
  filename,
  plot_markers,
  width = 3,
  height = 3,
  device = svg
)

plot_markers
```


## Plot heatmap of marker genes expressions per component
```{r fig.height = 5, fig.width = 2}
# Minmax scaling for visualization only
to_plot <-
  as.data.frame(apply(to_plot, 1, function(x)
    (x - min(x)) / (max(x) - min(x))))


nrows <- dim(to_plot)[[1]]
ncols <- dim(to_plot)[[2]]
factor_rows <-
  factor(rep(c("Liver", "Brain", "Lung"), each = nrows / 3), levels = c("Liver", "Brain", "Lung"))
factor_columns <-
  factor(rep(c("1", "2", "3"), each = ncols / 3), levels = c("1", "2", "3"))
marker_map <- Heatmap(
  to_plot,
  name = "W",
  cluster_rows = FALSE,
  cluster_columns = FALSE,
  show_column_names = FALSE,
  show_row_names = FALSE,
  show_heatmap_legend = F,
  column_split = factor_columns,
  row_split = factor_rows,
  use_raster = T,
  column_title_gp = gpar(
    fill = colors_v[3:5],
    border = "black",
    fontsize = 20
  ),
  row_title_gp = gpar(fontsize = 20),
  column_title_side  = "bottom"
)
filename <- paste0(dir_to_save_fig, "heatmap.svg")
svg(filename, width = 2, height = 5)
draw(marker_map)
dev.off()
marker_map
```