scMAGeCK is a computational model to identify genes associated with multiple expression phenotypes from CRISPR screening coupled with single-cell RNA sequencing data (e.g., CROP-seq).
scMAGeCK is based on our previous MAGeCK and MAGeCK-VISPR models for pooled CRISPR screens, but further extends to scRNA-seq as the readout of the screening experiment. scMAGeCK consists of two modules: scMAGeCK-Robust Rank Aggregation (RRA), a sensitive and precise algorithm to detect genes whose perturbation links to one single marker expression; and scMAGeCK-LR, a linear-regression based approach that unravels the perturbation effects on thousands of gene expressions, especially from cells undergo multiple perturbations.
library(scMAGeCK)
### BARCODE file contains cell identity information, generated from the cell identity collection step
BARCODE <- system.file("extdata","barcode_rec.txt",package = "scMAGeCK")
### RDS can be a Seurat object or local RDS file path that contains the scRNA-seq dataset
RDS <- system.file("extdata","singles_dox_mki67_v3.RDS",package = "scMAGeCK")
### Set RRA executable file path.
### You can generate RRA executable file by following commands:
### wget https://bitbucket.org/weililab/scmageck/downloads/RRA_0.5.9.zip
### unzip RRA_0.5.9.zip
### cd RRA_0.5.9
### make
RRAPATH <- "/Library/RRA_0.5.9/bin/RRA"
target_gene <- "MKI67"
rra_result <- scmageck_rra(BARCODE=BARCODE, RDS=RDS, GENE=target_gene, RRAPATH=RRAPATH,
LABEL='dox_mki67', NEGCTRL=NULL, KEEPTMP=FALSE, PATHWAY=FALSE, SAVEPATH=NULL)## Testing Rcpp:## Provided RRAPATH not found; falling back to the built-in RRA.## keep_tmp: FALSE## ispathway: FALSE## Total barcode records: 8425## Unique barcode records: 6704## Reading RDS file: /tmp/RtmpDQ3ihF/Rinst15f334fcf062/scMAGeCK/extdata/singles_dox_mki67_v3.RDS## Validating object structure## Updating object slots## Ensuring keys are in the proper structure## Updating matrix keys for DimReduc 'pca'## Updating matrix keys for DimReduc 'tsne'## Updating matrix keys for DimReduc 'umap'## Ensuring keys are in the proper structure## Ensuring feature names don't have underscores or pipes## Updating slots in RNA## Updating slots in pca## Updating slots in tsne## Setting tsne DimReduc to global## Updating slots in umap## Setting umap DimReduc to global## Validating object structure for Assay 'RNA'## Validating object structure for DimReduc 'pca'## Validating object structure for DimReduc 'tsne'## Validating object structure for DimReduc 'umap'## Object representation is consistent with the most current Seurat version## Cell names in expression matrix and barcode file do not match. Try to remove possible trailing "-1"s...## Target gene: MKI67## Testing gene MKI67 ...## Calling RRA parameters:
## RRA -i sample_9342.7506705441_rra_low.txt -o sample_9342.7506705441_rra_low.out -p 0.3 --max-sgrnapergene-permutation 10000
## Welcome to RRA v 0.5.9.
## Reading input file...
## Summary: 5698 sgRNAs, 30 genes, 1 lists; skipped sgRNAs:0
## Increase the number of permutations to 3334 to get precise p values. To avoid this, use the --permutation option.
## Permuting genes with 118 sgRNAs...
## Permuting genes with 119 sgRNAs...
## Permuting genes with 131 sgRNAs...
## Permuting genes with 143 sgRNAs...
## Permuting genes with 145 sgRNAs...
## Permuting genes with 146 sgRNAs...
## Permuting genes with 149 sgRNAs...
## Permuting genes with 150 sgRNAs...
## Permuting genes with 155 sgRNAs...
## Permuting genes with 157 sgRNAs...
## Permuting genes with 160 sgRNAs...
## Permuting genes with 161 sgRNAs...
## Permuting genes with 163 sgRNAs...
## Permuting genes with 164 sgRNAs...
## Permuting genes with 165 sgRNAs...
## Permuting genes with 172 sgRNAs...
## Permuting genes with 176 sgRNAs...
## Permuting genes with 179 sgRNAs...
## Permuting genes with 182 sgRNAs...
## Permuting genes with 184 sgRNAs...
## Permuting genes with 212 sgRNAs...
## Permuting genes with 252 sgRNAs...
## Permuting genes with 277 sgRNAs...
## Permuting genes with 302 sgRNAs...
## Permuting genes with 368 sgRNAs...
## Permuting genes with 496 sgRNAs...
## Number of genes under FDR adjustment: 30
## Calling RRA parameters:
## RRA -i sample_9342.7506705441_rra_high.txt -o sample_9342.7506705441_rra_high.out -p 0.3 --max-sgrnapergene-permutation 10000
## Welcome to RRA v 0.5.9.
## Reading input file...
## Summary: 5698 sgRNAs, 30 genes, 1 lists; skipped sgRNAs:0
## Increase the number of permutations to 3334 to get precise p values. To avoid this, use the --permutation option.
## Permuting genes with 118 sgRNAs...
## Permuting genes with 119 sgRNAs...
## Permuting genes with 131 sgRNAs...
## Permuting genes with 143 sgRNAs...
## Permuting genes with 145 sgRNAs...
## Permuting genes with 146 sgRNAs...
## Permuting genes with 149 sgRNAs...
## Permuting genes with 150 sgRNAs...
## Permuting genes with 155 sgRNAs...
## Permuting genes with 157 sgRNAs...
## Permuting genes with 160 sgRNAs...
## Permuting genes with 161 sgRNAs...
## Permuting genes with 163 sgRNAs...
## Permuting genes with 164 sgRNAs...
## Permuting genes with 165 sgRNAs...
## Permuting genes with 172 sgRNAs...
## Permuting genes with 176 sgRNAs...
## Permuting genes with 179 sgRNAs...
## Permuting genes with 182 sgRNAs...
## Permuting genes with 184 sgRNAs...
## Permuting genes with 212 sgRNAs...
## Permuting genes with 252 sgRNAs...
## Permuting genes with 277 sgRNAs...
## Permuting genes with 302 sgRNAs...
## Permuting genes with 368 sgRNAs...
## Permuting genes with 496 sgRNAs...
## Number of genes under FDR adjustment: 30 library(scMAGeCK)
### BARCODE file contains cell identity information, generated from the cell identity collection step
BARCODE <- system.file("extdata","barcode_rec.txt",package = "scMAGeCK")
### RDS can be a Seurat object or local RDS file path that contains the scRNA-seq dataset
RDS <- system.file("extdata","singles_dox_mki67_v3.RDS",package = "scMAGeCK")
lr_result <- scmageck_lr(BARCODE=BARCODE, RDS=RDS, LABEL='dox_scmageck_lr',
NEGCTRL = 'NonTargetingControlGuideForHuman', PERMUTATION = 1000, SAVEPATH=NULL, LAMBDA=0.01)## run_signature: FALSE## Total barcode records: 8425## Neg Ctrl guide: NonTargetingControlGuideForHuman## Reading RDS file: /tmp/RtmpDQ3ihF/Rinst15f334fcf062/scMAGeCK/extdata/singles_dox_mki67_v3.RDS## Validating object structure## Updating object slots## Ensuring keys are in the proper structure## Updating matrix keys for DimReduc 'pca'## Updating matrix keys for DimReduc 'tsne'## Updating matrix keys for DimReduc 'umap'## Ensuring keys are in the proper structure## Ensuring feature names don't have underscores or pipes## Updating slots in RNA## Updating slots in pca## Updating slots in tsne## Setting tsne DimReduc to global## Updating slots in umap## Setting umap DimReduc to global## Validating object structure for Assay 'RNA'## Validating object structure for DimReduc 'pca'## Validating object structure for DimReduc 'tsne'## Validating object structure for DimReduc 'umap'## Object representation is consistent with the most current Seurat version## Cell names in expression matrix and barcode file do not match. Try to remove possible trailing "-1"s...## 6704 ...## 6229 ...## Index matrix dimension: 5698 , 30## Filter genes whose expression is greater than 0 in raw read count in less than 0.01 single-cell populations.## Selected genes: 25## Selected cells: 5698## Permutation: 100 / 1000 ...## Permutation: 200 / 1000 ...## Permutation: 300 / 1000 ...## Permutation: 400 / 1000 ...## Permutation: 500 / 1000 ...## Permutation: 600 / 1000 ...## Permutation: 700 / 1000 ...## Permutation: 800 / 1000 ...## Permutation: 900 / 1000 ...## Permutation: 1000 / 1000 ...The scmageck_rra function will output the ranking and p values of each perturbed genes, using the RRA program in MAGeCK. Users familiar with the MAGeCK program may find it similar with the gene_summary output in MAGeCK.
Here is the example output of scMAGeCK-RRA:
Row.names items_in_group.low lo_value.low p.low FDR.low goodsgrna.low items_in_group.high lo_value.high p.high FDR.high goodsgrna.high
TP53 271 0.11832 0.95619 1 48 271 1.014e-83 4.9975e-06 0.00015 184Explanations of each column are below:
| Column | Content |
|---|---|
| Row.names | Perturbed gene name |
| items_in_group.low | The number of single-cells with each gene perturbed |
| lo_value.low | The RRA score in negative selection (reducing the marker expression if this gene is perturbed). The RRA score uses a p value from rank order statistics to measure the degree of selection; the smaller score, the stronger the selection is. More information on the calculation of RRA score can be found in our original MAGeCK paper. |
| p.low | The raw p-value (using permutation) of this gene in negative selection |
| FDR.low | The false discovery rate of this gene in negative selection |
| goodsgrna.low | The number of single-cells that passes the threshold and is considered in the RRA score calculation in negative selection |
| items_in_group.high | The same as items_in_group.low: the number of single-cells with each gene perturbed) |
| lo_value.high | The RRA score in positive selection (increasing the marker expression if this gene is perturbed) |
| p.high | The raw p-value (using permutation) of this gene in positive selection |
| FDR.high | The false discovery rate of this gene in positive selection |
| goodsgrna.high | The number of single-cells that passes the threshold and is considered in the RRA score calculation in positive selection |
The scmageck_lr function will generate several files below:
| File | Description |
|---|---|
| lr_score | The score (similar with log fold change) of each perturbed gene (rows) on each marker gene (columns) |
| lr_score.pval | The associated p values of each score |
| LR.RData | An R object to store scores and p values |
The format of score.txt and score.pval.txt is a simple table file with rows corresponding to perturbed genes and columns corresponding to marker genes. For example in the score.txt,
Perturbedgene APC ARID1A TP53 MKI67
APC 0.138075836476524 -0.0343441660045313 0.214449590551132 -0.150287676553705This row records the effects of perturbing APC gene on the expressions of APC, ARID1A, TP53 and MKI67.
Questions? Comments? Join the MAGeCK Google group or email us ([email protected]) directly.
Any advice and suggestions will be greatly appreciated.
## R version 4.6.1 (2026-06-24)
## Platform: x86_64-pc-linux-gnu
## Running under: Ubuntu 26.04 LTS
##
## Matrix products: default
## BLAS: /usr/lib/x86_64-linux-gnu/openblas-pthread/libblas.so.3
## LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/libopenblasp-r0.3.32.so; LAPACK version 3.12.0
##
## locale:
## [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
## [3] LC_TIME=en_US.UTF-8 LC_COLLATE=en_US.UTF-8
## [5] LC_MONETARY=en_US.UTF-8 LC_MESSAGES=en_US.UTF-8
## [7] LC_PAPER=en_US.UTF-8 LC_NAME=C
## [9] LC_ADDRESS=C LC_TELEPHONE=C
## [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C
##
## time zone: Etc/UTC
## tzcode source: system (glibc)
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] scMAGeCK_0.99.1 rmarkdown_2.31
##
## loaded via a namespace (and not attached):
## [1] deldir_2.0-4 pbapply_1.7-4 gridExtra_2.3.1
## [4] rlang_1.3.0 magrittr_2.0.5 RcppAnnoy_0.0.23
## [7] otel_0.2.0 spatstat.geom_3.8-1 matrixStats_1.5.0
## [10] ggridges_0.5.7 compiler_4.6.1 png_0.1-9
## [13] vctrs_0.7.3 reshape2_1.4.5 stringr_1.6.0
## [16] pkgconfig_2.0.3 fastmap_1.2.0 promises_1.5.0
## [19] purrr_1.2.2 xfun_0.60 cachem_1.1.0
## [22] jsonlite_2.0.0 goftest_1.2-3 later_1.4.8
## [25] spatstat.utils_3.2-3 irlba_2.3.7 parallel_4.6.1
## [28] cluster_2.1.8.2 R6_2.6.1 ica_1.0-3
## [31] spatstat.data_3.1-9 bslib_0.11.0 stringi_1.8.7
## [34] RColorBrewer_1.1-3 reticulate_1.46.0 spatstat.univar_3.2-0
## [37] parallelly_1.48.0 lmtest_0.9-40 jquerylib_0.1.4
## [40] scattermore_1.2 Rcpp_1.1.2 knitr_1.51
## [43] tensor_1.5.1 future.apply_1.20.2 zoo_1.8-15
## [46] sctransform_0.4.3 httpuv_1.6.17 Matrix_1.7-5
## [49] splines_4.6.1 igraph_2.3.3 tidyselect_1.2.1
## [52] abind_1.4-8 yaml_2.3.12 spatstat.random_3.5-0
## [55] codetools_0.2-20 miniUI_0.1.2 spatstat.explore_3.8-1
## [58] listenv_1.0.0 lattice_0.22-9 tibble_3.3.1
## [61] plyr_1.8.9 shiny_1.14.0 S7_0.2.2
## [64] ROCR_1.0-12 evaluate_1.0.5 Rtsne_0.17
## [67] future_1.70.0 fastDummies_1.7.6 survival_3.8-9
## [70] polyclip_1.10-7 fitdistrplus_1.2-6 pillar_1.11.1
## [73] Seurat_5.5.1 KernSmooth_2.23-26 plotly_4.12.0
## [76] generics_0.1.4 RcppHNSW_0.7.0 sp_2.2-1
## [79] ggplot2_4.0.3 scales_1.4.0 globals_0.19.1
## [82] xtable_1.8-8 glue_1.8.1 lazyeval_0.2.3
## [85] maketools_1.3.2 tools_4.6.1 sys_3.4.3
## [88] data.table_1.18.4 RSpectra_0.16-2 RANN_2.6.2
## [91] buildtools_1.0.0 dotCall64_1.2 cowplot_1.2.0
## [94] grid_4.6.1 tidyr_1.3.2 nlme_3.1-169
## [97] patchwork_1.3.2 cli_3.6.6 spatstat.sparse_3.2-0
## [100] spam_2.11-4 viridisLite_0.4.3 dplyr_1.2.1
## [103] uwot_0.2.4 gtable_0.3.6 sass_0.4.10
## [106] digest_0.6.39 progressr_1.0.0 ggrepel_0.9.8
## [109] htmlwidgets_1.6.4 SeuratObject_5.4.0 farver_2.1.2
## [112] htmltools_0.5.9 lifecycle_1.0.5 httr_1.4.8
## [115] mime_0.13 MASS_7.3-65