3_Quality control

Using scrublet

This code shows how to set up conda and download scrublet on the cluster. This can be adapted to run locally.

# run in terminal
ml fhPython/3.8.6-foss-2020b-Python-3.8.6 
ml Anaconda3/2022.05
conda create -n py3 python=3
source activate py3                                                                                                                  

Sys.setenv(RETICULATE_PYTHON = "~/.conda/envs/py3/bin/python3")
library("SeuratDisk")
library("Seurat")
library("reticulate")
use_python("~/.conda/envs/py3/bin/python3")
use_condaenv("~/.conda/envs/py3")
py_config()
py_install("numpy")
py_install("scipy")
py_install("scikit-learn") # Scrublet needs this
py_install("scikit-image") # Scrublet needs this
py_install("matplotlib") # Scrublet needs this
py_install("scrublet")

#Make Object read cell ranger output. i like to use the h5 file

# Load the dataset
seul <- Read10X_h5(file.path(DATA_DIR, "AY_cBIR_1", "outs/per_sample_outs/AY_cBIR_1/count/sample_filtered_feature_bc_matrix.h5"))
rna_mat<-seul$`Gene Expression`
seu <- CreateSeuratObject(counts = seul$`Gene Expression`, project = "cBIR", min.cells = 0, min.features = 0)
exprData <- Matrix(as.matrix(rna_mat), sparse = TRUE)
writeMM(exprData, file.path(DATA_DIR, "AY_cBIR_1", "matrix.mtx"))

import scrublet as scr
import scipy.io
import numpy as np
import os

#Load raw counts matrix and gene list
input_dir = '../data'
counts_matrix = scipy.io.mmread(input_dir + '/matrix.mtx').T.tocsc()

#Initialize Scrublet object
scrub = scr.Scrublet(counts_matrix,
                     expected_doublet_rate=0.1,
                     sim_doublet_ratio=2,
                     n_neighbors = 8)


#Run the default pipeline
doublet_scores, predicted_doublets = scrub.scrub_doublets(min_counts=1, 
                                                          min_cells=3, 
                                                          min_gene_variability_pctl=85, 
                                                          n_prin_comps=25)

# Import scrublet's doublet score
seu$Doubletscore <- py$doublet_scores

 

# Plot doublet score
ggplot(seu@meta.data, aes(x = Doubletscore, stat(ndensity))) +
  geom_histogram(bins = 200, colour ="lightgrey")+
  geom_vline(xintercept = 0.23, colour = "red", linetype = 2)+
  geom_vline(xintercept = 0.15, colour = "green", linetype = 2) # Manually set threshold

 


# Manually set threshold at doublet score to 0.15
seu@meta.data$Predicted_doublets <- ifelse(py$doublet_scores > 0.2, "Doublet","Singlet" )
table(seu@meta.data$Predicted_doublets)