Percent_Expressing doesn't work with chromatin assay
Closed this issue · 4 comments
Hi, thank you for this incredibly useful package! I hesitate to post this as a bug/issue since I don't know if scCustomize was designed to handle chromatin assays, and not sure if this will be possible or feasible. I have a single-cell multiome dataset (RNA-seq and ATAC-seq). Running Percent_Expressing on the RNA assay works as expected, but running it on the ATAC assay produces the following error:
Error in if (!(layer %in% slotNames(x = object))) { :
argument is of length zero
I am still using Seurat v.4.4.0, but I'm wondering if this is an issue that was introduced with the update for Seurat 5 since I think some things changed with GetAssayData() (seems to be where it's getting stuck; see traceback below).
I can reproduce this behavior using the pbmcMultiome data package in SeuratData (datasets pbmc.rna and pbmc.atac).
InstallData("pbmcMultiome")
pbmc.rna <- LoadData("pbmcMultiome", "pbmc.rna")
pbmc.atac <- LoadData("pbmcMultiome", "pbmc.atac")
# To find names of features that definitely exist in dataset
head(rownames(pbmc.rna@assays$RNA))
## "MIR1302-2HG" "FAM138A" "OR4F5" "AL627309.1" "AL627309.3" "AL627309.2"
head(rownames(pbmc.atac@assays$ATAC))
## "chr1-10109-10357" "chr1-180730-181630" "chr1-191491-191736" "chr1-267816-268196" "chr1-586028-586373" "chr1-629721-630172"
# Set idents
Idents(pbmc.rna) <- pbmc.rna$seurat_annotations
Idents(pbmc.atac) <- pbmc.atac$seurat_annotations
# Try RNA assay
Percent_Expressing(pbmc.rna, "AL627309.1", assay = "RNA")
## CD4.Naive CD4.TCM filtered CD8.Naive CD16.Mono NK Treg
## AL627309.1 0.07047216 0.435161 0.6012024 0.2836879 0.7782101 0.2136752 0
## CD14.Mono cDC CD8.TEM_1 Intermediate.B Naive.B Plasma CD4.TEM MAIT
## AL627309.1 1.88478 0 0 0 0 0 0.3355705 0.729927
## Memory.B gdT pDC CD8.TEM_2 HSPC
## AL627309.1 0 0.6849315 1.886792 0.2793296 0
# Try ATAC assay
Percent_Expressing(pbmc.atac, "chr1-191491-191736", assay = "ATAC")
## Error in if (!(layer %in% slotNames(x = object))) { :
## argument is of length zero
traceback()
## 7: GetAssayData.ChromatinAssay(object = object[[assay]], layer = layer)
## 6: GetAssayData(object = object[[assay]], layer = layer)
## 5: GetAssayData.Seurat(object = data, assay = assay)
## 4: GetAssayData(object = data, assay = assay)
## 3: rownames(x = GetAssayData(object = data, assay = assay))
## 2: Gene_Present(data = seurat_object, gene_list = features, print_msg = FALSE,
## case_check = TRUE, seurat_assay = assay)
## 1: Percent_Expressing(pbmc.atac, "chr1-191491-191736")sessionInfo()
sessionInfo() output
R version 4.3.1 (2023-06-16)
Platform: aarch64-apple-darwin20 (64-bit)
Running under: macOS Ventura 13.4
Matrix products: default
BLAS: /System/Library/Frameworks/Accelerate.framework/Versions/A/Frameworks/vecLib.framework/Versions/A/libBLAS.dylib
LAPACK: /Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/lib/libRlapack.dylib; LAPACK version 3.11.0
locale:
[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
time zone: America/Chicago
tzcode source: internal
attached base packages:
[1] stats graphics grDevices utils datasets methods base
other attached packages:
[1] pbmcMultiome.SeuratData_0.1.4 pbmc3k.SeuratData_3.1.4 ifnb.SeuratData_3.1.0 SeuratData_0.2.2
[5] scCustomize_1.1.3 Signac_1.12.9003 SeuratObject_5.0.1 Seurat_4.4.0
loaded via a namespace (and not attached):
[1] RColorBrewer_1.1-3 shape_1.4.6 rstudioapi_0.15.0 jsonlite_1.8.7 magrittr_2.0.3
[6] ggbeeswarm_0.7.2 spatstat.utils_3.0-3 GlobalOptions_0.1.2 zlibbioc_1.46.0 vctrs_0.6.4
[11] ROCR_1.0-11 spatstat.explore_3.2-3 Rsamtools_2.16.0 paletteer_1.5.0 RCurl_1.98-1.12
[16] janitor_2.2.0 RcppRoll_0.3.0 forcats_1.0.0 htmltools_0.5.6.1 sctransform_0.4.0
[21] parallelly_1.36.0 KernSmooth_2.23-22 htmlwidgets_1.6.2 ica_1.0-3 plyr_1.8.8
[26] lubridate_1.9.2 plotly_4.10.2 zoo_1.8-12 igraph_1.5.1 mime_0.12
[31] lifecycle_1.0.3 pkgconfig_2.0.3 Matrix_1.6-4 R6_2.5.1 fastmap_1.1.1
[36] snakecase_0.11.1 GenomeInfoDbData_1.2.10 fitdistrplus_1.1-11 future_1.33.0 shiny_1.7.5
[41] digest_0.6.33 colorspace_2.1-0 rematch2_2.1.2 patchwork_1.1.3 S4Vectors_0.38.1
[46] tensor_1.5 grr_0.9.5 irlba_2.3.5.1 GenomicRanges_1.52.0 progressr_0.14.0
[51] timechange_0.2.0 fansi_1.0.4 spatstat.sparse_3.0-2 httr_1.4.7 polyclip_1.10-4
[56] abind_1.4-5 compiler_4.3.1 BiocParallel_1.34.2 MASS_7.3-60 rappdirs_0.3.3
[61] tools_4.3.1 vipor_0.4.5 lmtest_0.9-40 beeswarm_0.4.0 httpuv_1.6.11
[66] future.apply_1.11.0 goftest_1.2-3 glue_1.6.2 nlme_3.1-163 promises_1.2.1
[71] grid_4.3.1 Rtsne_0.16 cluster_2.1.4 reshape2_1.4.4 generics_0.1.3
[76] gtable_0.3.4 spatstat.data_3.0-1 tidyr_1.3.0 data.table_1.14.8 sp_2.0-0
[81] utf8_1.2.3 XVector_0.40.0 BiocGenerics_0.46.0 spatstat.geom_3.2-5 RcppAnnoy_0.0.21
[86] ggrepel_0.9.3 RANN_2.6.1 pillar_1.9.0 stringr_1.5.0 ggprism_1.0.4
[91] spam_2.9-1 later_1.3.1 circlize_0.4.15 splines_4.3.1 dplyr_1.1.3
[96] lattice_0.21-8 survival_3.5-7 deldir_1.0-9 tidyselect_1.2.0 Biostrings_2.68.1
[101] Matrix.utils_0.9.8 miniUI_0.1.1.1 pbapply_1.7-2 gridExtra_2.3 IRanges_2.34.1
[106] scattermore_1.2 stats4_4.3.1 matrixStats_1.0.0 stringi_1.7.12 lazyeval_0.2.2
[111] codetools_0.2-19 tibble_3.2.1 cli_3.6.1 uwot_0.1.16 xtable_1.8-4
[116] reticulate_1.32.0 munsell_0.5.0 Rcpp_1.0.11 GenomeInfoDb_1.36.3 globals_0.16.2
[121] spatstat.random_3.1-6 png_0.1-8 ggrastr_1.0.2 parallel_4.3.1 ellipsis_0.3.2
[126] ggplot2_3.4.3 dotCall64_1.0-2 bitops_1.0-7 listenv_0.9.0 viridisLite_0.4.2
[131] scales_1.2.1 ggridges_0.5.4 crayon_1.5.2 leiden_0.4.3 purrr_1.0.2
[136] rlang_1.1.1 fastmatch_1.1-4 cowplot_1.1.1Hi @cathynewman,
Thanks so much for kind words!! I've tried to make it so it will work with any Seurat object but I don't explicitly test on ATAC assays (but maybe something in future). Thanks so much for reproducible example too!!
So I think the error is actually due to the way some of the feature checking functions within scCustomize were originally written. I believe the issue should be solved if you update to current CRAN release of scCustomize v2.1.2 or to current "develop" branch v2.1.2.9003. When I run your code provided (using either Seurat 4.4.0 or Seurat 5+) it does return value for me.
The updated scCustomize will not require updating to Seurat v5 but does require SeuratObject V5+ because that ensures that a lot of the getter/setter functions are cross compatible between v5 and v3/4 style objects. If you are unable to update to current CRAN release or dev branch (there shouldn't be any breaking or reproducibility issues outside of some function names; only improvements) let me know and I can work on alternative fix for you.
Also if this doesn't fix your issue please let me know and I'll reopen the issue here.
Best,
Sam
Updating to scCustomize v.2.1.2 fixed the issue, thanks! Could've sworn I checked first to see if there was a newer version...
Great! No worries, glad it fixed the issue! :)
Also If there is anything Chromatin assay related that doesn't work in future in scCustomize or you think would be helpful to have please let me know!