Update tutorials

Author: jaydu1

README.md

@@ -10,13 +10,22 @@ Advances in single-cell sequencing and CRISPR technologies have enabled detailed
 
 # Usage
 
+We recommend using `causarray` in a conda environment:
+```cmd
+# create a new conda environment and install the necessary packages
+conda create -n causarray python=3.12 -y
+
+# activate the environment
+conda activate causarray
+```
 
 The module can be installed via PyPI:
 ```cmd
 pip install causarray
 ```
 
-The documentation is available at [causarray.readthedocs.io](https://causarray.readthedocs.io/en/latest/).
+For `R` users, `reticulate` can be used to call `causarray` from `R`.
+The documentation and tutorials using both `Python` and `R` are available at [causarray.readthedocs.io](https://causarray.readthedocs.io/en/latest/).
 
 
 
@@ -60,6 +69,9 @@ mkdir docs
 sphinx-quickstart
 cd docs
 make html # sphinx-build source build
+
+
+rmarkdown::render("perturbseq.Rmd", rmarkdown::md_document(variant = "markdown_github"))
 ```
 -->
 

docs/requirements.txt

@@ -4,7 +4,7 @@ numpydoc
 nbsphinx
 sphinxcontrib-bibtex
 sphinxcontrib-svg2pdfconverter
-
+myst-parser
 
 h5py
 joblib

docs/source/conf.py

@@ -10,13 +10,13 @@
 
 
 
-
 # -- General configuration ---------------------------------------------------
 # https://www.sphinx-doc.org/en/master/usage/configuration.html#general-configuration
 
 extensions = [
     'numpydoc',
     'nbsphinx',
+    'myst_parser',
     'sphinx.ext.mathjax',  # for math equations
     'sphinxcontrib.bibtex',  # for bibliographic references
     'sphinxcontrib.rsvgconverter',  # for SVG->PDF conversion in LaTeX output

docs/source/index.rst

@@ -11,10 +11,24 @@ Advances in single-cell sequencing and CRISPR technologies have enabled detailed
 
 
 .. toctree::
-   :maxdepth: 4
+   :maxdepth: 1
    :glob:
-   :caption: Contents:
+   :caption: Main functions
 
-   gcate
+   main_function/gcate
+   main_function/lfc
 
-   lfc
+
+.. toctree::
+   :maxdepth: 1
+   :glob:
+   :caption: Tutorials (Python)
+
+   tutorial/perturbseq-py.ipynb
+
+.. toctree::
+   :maxdepth: 1
+   :glob:
+   :caption: Tutorials (R)
+
+   tutorial/perturbseq-r.md   

docs/source/gcate.rst docs/source/main_function/gcate.rstdocs/source/gcate.rst renamed to docs/source/main_function/gcate.rst

@@ -0,0 +1,107 @@
+---
+title: "Joint Pertub-seq perturbations analysis"
+output:
+  md_document:
+    variant: markdown_github
+date: "`r format(Sys.time(), '%d %B, %Y')`"
+---
+
+```{r setup, include=FALSE}
+knitr::opts_chunk$set(echo = TRUE)
+setwd("/Users/dujinhong/Documents/study/Kathryn Roeder/Project/causarray/causarray/docs/source/tutorial")
+```
+
+
+The original data of Jin et al 2020 can be downloaded from the Broad single cell portal (https://singlecell.broadinstitute.org/single_cell/study/SCP1184).
+Here, we just use a subset of the data to demonstrate the workflow of the analysis.
+
+
+
+```{r}
+library(Seurat)
+library(caret)
+
+sc.seurat <- readRDS("perturbseq-exneu.rds")
+
+Y <- data.frame(t(sc.seurat[['RNA']]$counts)) # cell-by-gene matrix
+metadata <- sc.seurat@meta.data
+
+perturb <- metadata
+colnames(perturb) <- gsub("Perturbation", "trt_", colnames(perturb))
+perturb$trt_ <- relevel(as.factor(perturb$trt_), ref = "GFP")
+dmy <- dummyVars(" ~ trt_", data = perturb)
+A <- data.frame(predict(dmy, newdata = perturb))[,-1] # cell-by-trt matrix
+```
+
+For running causarray, we require the following inputs:
+
+- `Y`: the cell-by-gene gene expression matrix.
+- `A`: the cell-by-condition binary matrix of the perturbation/treatment conditions.
+- `X, X_A`: (optional) the cell-by-covariate matrix of the covariates of interest for outcome and propensity models.
+
+Here, `Y` and `A` can be dataframes.
+
+
+
+Use R package `reticulate` to load the Python package `causarray`.
+
+```{r}
+require(reticulate)
+Sys.setenv(PYTHONUNBUFFERED = TRUE)
+use_condaenv('causarray')
+causarray <- import("causarray")
+cat(causarray$`__version__`)
+
+# (Y, A) should be either data.frame or matrix
+# optional covariates can be provided as matrices
+dat <- causarray$prep_causarray_data(Y, A)
+names(dat) <- c("Y", "A", "X", "X_A")
+list2env(dat, .GlobalEnv)
+```
+
+
+We first apply gcate to estimate unmeasured confounders.
+
+```{r}
+r <- 10
+res_gate <- causarray$fit_gcate(Y, X, A, r, verbose=TRUE) # a list of results from 2 stages optimization
+U <- res_gate[[2]]$U
+```
+
+Next, we apply causarray to estimate the causal effects of perturbations on gene expression.
+
+```{r}
+offsets <- log(res_gate[[2]][['kwargs_glm']][['size_factor']]) # use the precomputed size factors
+res <- causarray$LFC(Y, cbind(X, U), A, cbind(X_A, U), offset=offsets, verbose=TRUE)
+```
+
+```{r}
+names(res) <- c("df_res", "estimation")
+list2env(res, .GlobalEnv)
+```
+
+
+```{r, fig.width=12, fig.height=6}
+library(dplyr)
+library(ggplot2)
+
+# Filter the results for significant discoveries
+significant_discoveries <- df_res[df_res$padj < 0.1, ]
+
+# Count the number of discoveries for each perturbation condition
+discovery_counts <- as.data.frame(table(significant_discoveries$trt))
+colnames(discovery_counts) <- c('Perturbation', 'Count')
+
+# Order the discovery_counts by Count in descending order
+discovery_counts <- discovery_counts %>% arrange(desc(Count))
+
+# Set the factor levels of Perturbation to ensure ggplot respects the order
+discovery_counts$Perturbation <- factor(discovery_counts$Perturbation, levels = discovery_counts$Perturbation)
+
+# Plot the number of discoveries for each perturbation condition
+ggplot(discovery_counts, aes(x = Perturbation, y = Count)) +
+  geom_bar(stat = "identity", fill = "royalblue") +  theme(axis.text.x = element_text(angle = 90, hjust = 1)) +
+  ggtitle('Number of Discoveries (padj < 0.1) for Each Perturbation Condition') +
+  xlab('Perturbation Condition') +
+  ylab('Number of Discoveries')
+```