debug.Rmd

---
title: "ScrubletR"
output: html_document
date: "2024-02-02"
always_allow_html: true
editor_options: 
  
  chunk_output_type: console
---

## Installing scrubletR

You will need to have the devtools package installed...

```{r, eval=F}
devtools::install_github("furlan-lab/scrubletR")
```


## Load data

```{r, dpi=300, fig.height=4, fig.width = 6, warning=F, message=F}

suppressPackageStartupMessages({
  library(viewmastRust)
  library(Seurat)
  library(scCustomize)
})

if(grepl("^gizmo", Sys.info()["nodename"])){
  ROOT_DIR2<-"/fh/fast/furlan_s/grp/data/ddata/BM_data"
} else {
  ROOT_DIR2<-"/Users/sfurlan/Library/CloudStorage/OneDrive-SharedLibraries-FredHutchinsonCancerCenter/Furlan_Lab - General/experiments/patient_marrows/aggr/cds/indy"
}

#query dataset
seuP<-readRDS(file.path(ROOT_DIR2, "220831_WC1.RDS"))
DimPlot_scCustom(seuP)

counts_matrix<-t(seuP@assays$RNA@counts[,1:2000])
counts_matrix<-t(seuP@assays$RNA@counts)
```

## Check Python Scrublet doublet scores

```{python}
import numpy as np
import scrublet as scr
scrub = scr.Scrublet(r.counts_matrix)
scrub._E_obs.data.shape

##set params
synthetic_doublet_umi_subsampling=1.0
use_approx_neighbors=True
distance_metric='euclidean'
get_doublet_neighbor_parents=False
min_counts=3
min_cells=3
min_gene_variability_pctl=85
log_transform=False
mean_center=True,
normalize_variance=True
n_prin_comps=30
svd_solver='arpack'
verbose=True

# doublet_scores, predicted_doublets = scrub.scrub_doublets()
scrub._E_sim = None
scrub._E_obs_norm = None
scrub._E_sim_norm = None
scrub._gene_filter = np.arange(scrub._E_obs.shape[1])
scr.pipeline_normalize(scrub)
r.py_E_obs_norm = scrub._E_obs_norm.data

scr.pipeline_get_gene_filter(scrub, min_counts=min_counts, min_cells=min_cells, min_gene_variability_pctl=min_gene_variability_pctl)
scrub._gene_filter.shape
pipeline_apply_gene_filter(scrub)

        print_optional('Simulating doublets...', verbose)
        scrub.simulate_doublets(sim_doublet_ratio=scrub.sim_doublet_ratio, synthetic_doublet_umi_subsampling=synthetic_doublet_umi_subsampling)
        pipeline_normalize(scrub, postnorm_total=1e6)
        if log_transform:
            pipeline_log_transform(scrub)
        if mean_center and normalize_variance:
            pipeline_zscore(scrub)
        elif mean_center:
            pipeline_mean_center(scrub)
        elif normalize_variance: 
            pipeline_normalize_variance(scrub)

        if mean_center:
            print_optional('Embedding transcriptomes using PCA...', verbose)
            pipeline_pca(scrub, n_prin_comps=n_prin_comps, random_state=scrub.random_state, svd_solver=svd_solver)
        else:
            print_optional('Embedding transcriptomes using Truncated SVD...', verbose)
            pipeline_truncated_svd(scrub, n_prin_comps=n_prin_comps, random_state=scrub.random_state, algorithm=svd_solver)            

        print_optional('Calculating doublet scores...', verbose)
        scrub.calculate_doublet_scores(
            use_approx_neighbors=use_approx_neighbors,
            distance_metric=distance_metric,
            get_doublet_neighbor_parents=get_doublet_neighbor_parents
            )
        scrub.call_doublets(verbose=verbose)

        t1=time.time()
        print_optional('Elapsed time: {:.1f} seconds'.format(t1 - t0), verbose)

```

## Check scrubletR

```{r}
scr<-Scrublet$new(counts_matrix = counts_matrix)

scr$pipeline_normalize()

ix<-sample(1:length(scr$E_obs_norm@x), 20000)
ggplot(data.frame(x=scr$E_obs_norm@x[ix], y = py_E_obs_norm[ix]), aes(x=x, y=y))+geom_point()
#debug(scr$pipeline_pca)
#debug(scr$scrub_doublets)
ds<-scr$scrub_doublets()

length(scr$gene_filter)


seuP$ds<-ds$doublet_scores

seuP$dsP <- dsP
seuP$dsPr <- dsPr

FeaturePlot_scCustom(seuP, features = "ds")
FeaturePlot_scCustom(seuP, features = "dsP")
DimPlot(seuP, group.by = "dsPr")
ggplot(data.frame(ds = seuP$dsP), aes(x=ds))+geom_density()

```

```{r}

library("mixtools")

#' Plot a Mixture Component
#' 
#' @param x Input data
#' @param mu Mean of component
#' @param sigma Standard deviation of component
#' @param lam Mixture weight of component
plot_mix_comps <- function(x, mu, sigma, lam) {
  lam * dnorm(x, mu, sigma)
}

set.seed(1)
wait <- seuP$dsP
mixmdl <- normalmixEM(wait, k = 2)
data.frame(x = mixmdl$x) %>%
  ggplot() +
  geom_histogram(aes(x, ..density..), binwidth = 0.01, colour = "black", 
                 fill = "white") +
  stat_function(geom = "line", fun = plot_mix_comps,
                args = list(mixmdl$mu[1], mixmdl$sigma[1], lam = mixmdl$lambda[1]),
                colour = "red", lwd = 1.5) +
  stat_function(geom = "line", fun = plot_mix_comps,
                args = list(mixmdl$mu[2], mixmdl$sigma[2], lam = mixmdl$lambda[2]),
                colour = "blue", lwd = 1.5) +
  ylab("Density")


seuP$dp_sf<-mixmdl$posterior[,2]>0.85 & mixmdl$posterior[,1]<0.15
DimPlot(seuP, group.by = "dp_sf")
```


## Appendix
```{r Appendix}
sessionInfo()
getwd()
```