search

stuart-lab / signac

R toolkit for the analysis of single-cell chromatin data
https://stuartlab.org/signac/
Other
327 stars 88 forks source link

TSSEnrichment fails only on fast = F #1395

Closed mossconfuse closed 1 year ago

mossconfuse commented 1 year ago

I am running into an issue with TSSEnrichment. When I run TSSEnrichment on my multimodal seurat object:

> TSSEnrichment(seuratObject, fast = T)
Extracting TSS positions
Extracting fragments at TSSs
  |++++++++++++++++++++++++++++++++++++++++++++++++++| 100% elapsed=03m 22s

Computing TSS enrichment score
An object of class Seurat 
50661 features across 2648 samples within 2 assays 
Active assay: ATAC (17333 features, 0 variable features)
 2 layers present: counts, data
 1 other assay present: RNA
> TSSEnrichment(seuratObject, fast = F)
Extracting TSS positions
Finding + strand cut sites
Finding - strand cut sites
Error: non-conformable matrix dimensions in dimCheck(e1, e2)
> table(seuratObject@assays$ATAC@positionEnrichment[1:50])
< table of extent 0 x 0 x 0 x 0 x 0 x 0 x 0 x 0 x 0 x 0 x 0 x 0 x 0 x 0 x 0 x 0 x 0 x 0 x 0 x 0 x 0 x 0 x 0 x 0 x 0 x 0 x 0 x 0 x 0 x 0 x 0 x 0 x 0 x 0 x 0 x 0 x 0 x 0 x 0 x 0 x 0 x 0 x 0 x 0 x 0 x 0 x 0 x 0 x 0 x 0 >
> seuratObject@meta.data$TSS.enrichment
   [1]  2.501846  2.236894  3.187248  2.960618  3.141062  2.987967  2.031775  2.869239  3.314771
  [10]  3.849963  3.264457  2.866023  2.936598  3.570835  3.050826  3.617816  3.675272  3.512585
  [19]  3.005801  3.606243  2.338743  3.626644  3.096224  3.521044  1.760992  3.295038  3.092574

So the fast = T version works fine, but not the slow version, which means I can't plot the TSSPlot. I am not sure what Error: non-conformable matrix dimensions in dimCheck(e1, e2) points to, as I don't get any search hits for it.

Thanks again for all the help.

timoast commented 1 year ago

Can you include the output of sessionInfo()?

mossconfuse commented 1 year ago

Sorry for the delayed response. Here is more details on how I got exactly my output:

> load("Saved RObjects/PreprocessedSeurat.RData")
> ppsn
An object of class Seurat 
172071 features across 11667 samples within 4 assays 
Active assay: peaks (88220 features, 88220 variable features)
 3 other assays present: RNA, ATAC, SCT
 7 dimensional reductions calculated: pca, lsi, sct_harmony, sct_harmony_UMAP, lsi_harmony, umap.atac, wnn.umap
> DefaultAssay(ppsn) <- "ATAC"
> rawDataPath <- "~/Library/xxx/Raw Data/"
> frags <- Fragments(ppsn)  # get list of fragment objects
> frags[[1]] <- UpdatePath(frags[[1]], new.path = paste0(rawDataPath,"Wt_T00","/atac_fragments.tsv.gz")) 
Computing hash
> frags[[2]] <- UpdatePath(frags[[2]], new.path = paste0(rawDataPath,"Wt_T06","/atac_fragments.tsv.gz")) 
Computing hash
> frags[[3]] <- UpdatePath(frags[[3]], new.path = paste0(rawDataPath,"Wt_T12","/atac_fragments.tsv.gz")) 
Computing hash
> frags[[4]] <- UpdatePath(frags[[4]], new.path = paste0(rawDataPath,"ste_T00","/atac_fragments.tsv.gz")) 
Computing hash
> frags[[5]] <- UpdatePath(frags[[5]], new.path = paste0(rawDataPath,"ste_T06","/atac_fragments.tsv.gz")) 
Computing hash
> frags[[6]] <- UpdatePath(frags[[6]], new.path = paste0(rawDataPath,"ste_T12","/atac_fragments.tsv.gz"))  # update path. Do this for any/all fragment objects in the list
Computing hash
> Fragments(ppsn) <- NULL  # remove fragment information from assay
> Fragments(ppsn) <- frags  # assign update list of fragment objects back to the assay
> ppsn <- NucleosomeSignal(ppsn, n = ncol(ppsn)*5000,verbose = TRUE)
Found 1838 cell barcodes
Done Processing 58 million linesFound 2596 cell barcodes
Done Processing 58 million linesFound 2579 cell barcodes
Done Processing 58 million linesFound 2717 cell barcodes
Done Processing 58 million linesFound 927 cell barcodes
Done Processing 58 million linesFound 1010 cell barcodes
Done Processing 44 million lines                  
> ppsn <- TSSEnrichment(ppsn, fast = TRUE) 
Extracting TSS positions
Extracting fragments at TSSs
  |++++++++++++++++++++++++++++++++++++++++++++++++++| 100% elapsed=02m 34srning: longer object length is not a multiple of shorter object lengthWarning: number of items to replace is not a multiple of replacement lengthWarning: longer object length is not a multiple of shorter object lengthWarning: number of items to replace is not a multiple of replacement length
  |++++++++++++++++++++++++++++++++++++++++++++++++++| 100% elapsed=02m 42srning: longer object length is not a multiple of shorter object lengthWarning: number of items to replace is not a multiple of replacement lengthWarning: longer object length is not a multiple of shorter object lengthWarning: number of items to replace is not a multiple of replacement length
  |++++++++++++++++++++++++++++++++++++++++++++++++++| 100% elapsed=03m 43srning: longer object length is not a multiple of shorter object lengthWarning: number of items to replace is not a multiple of replacement lengthWarning: longer object length is not a multiple of shorter object lengthWarning: number of items to replace is not a multiple of replacement length
  |++++++++++++++++++++++++++++++++++++++++++++++++++| 100% elapsed=02m 37srning: longer object length is not a multiple of shorter object lengthWarning: number of items to replace is not a multiple of replacement lengthWarning: longer object length is not a multiple of shorter object lengthWarning: number of items to replace is not a multiple of replacement length
  |++++++++++++++++++++++++++++++++++++++++++++++++++| 100% elapsed=02m 40srning: longer object length is not a multiple of shorter object lengthWarning: number of items to replace is not a multiple of replacement lengthWarning: longer object length is not a multiple of shorter object lengthWarning: number of items to replace is not a multiple of replacement length
  |++++++++++++++++++++++++++++++++++++++++++++++++++| 100% elapsed=01m 42srning: longer object length is not a multiple of shorter object lengthWarning: number of items to replace is not a multiple of replacement lengthWarning: longer object length is not a multiple of shorter object lengthWarning: number of items to replace is not a multiple of replacement length

Computing TSS enrichment score
> ppsn@meta.data$TSS.enrichment
   [1] 2.236894 3.187248 2.960618 3.141062 2.987967 2.031775 3.314771 3.849963 3.264457 3.570835 3.050826 3.617816 3.512585

> ppsn <- TSSEnrichment(ppsn, fast = FALSE)
Extracting TSS positions
Finding + strand cut sites
Finding - strand cut sites
Error: non-conformable matrix dimensions in dimCheck(e1, e2)

> session_info()
─ Session info ─────────────────────────────────────────────────────────────────────────────────────────────────────────────────
 setting  value
 version  R version 4.3.0 (2023-04-21)
 os       macOS Ventura 13.2.1
 system   aarch64, darwin20
 ui       RStudio
 language (EN)
 collate  en_US.UTF-8
 ctype    en_US.UTF-8
 tz       Asia/Tokyo
 date     2023-05-17
 rstudio  2023.03.1+446 Cherry Blossom (desktop)
 pandoc   NA

─ Packages ─────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
 package          * version   date (UTC) lib source
 abind              1.4-5     2016-07-21 [1] CRAN (R 4.3.0)
 AnnotationDbi    * 1.62.1    2023-05-08 [1] Bioconductor
 Biobase          * 2.60.0    2023-05-08 [1] Bioconductor
 BiocGenerics     * 0.46.0    2023-04-25 [1] Bioconductor
 BiocParallel       1.34.1    2023-05-08 [1] Bioconductor
 Biostrings         2.68.0    2023-05-08 [1] Bioconductor
 bit                4.0.5     2022-11-15 [1] CRAN (R 4.3.0)
 bit64              4.0.5     2020-08-30 [1] CRAN (R 4.3.0)
 bitops             1.0-7     2021-04-24 [1] CRAN (R 4.3.0)
 blob               1.2.4     2023-03-17 [1] CRAN (R 4.3.0)
 cachem             1.0.8     2023-05-01 [1] CRAN (R 4.3.0)
 callr              3.7.3     2022-11-02 [1] CRAN (R 4.3.0)
 cli                3.6.1     2023-03-23 [1] CRAN (R 4.3.0)
 cluster            2.1.4     2022-08-22 [1] CRAN (R 4.3.0)
 codetools          0.2-19    2023-02-01 [1] CRAN (R 4.3.0)
 colorspace         2.1-0     2023-01-23 [1] CRAN (R 4.3.0)
 cowplot            1.1.1     2020-12-30 [1] CRAN (R 4.3.0)
 crayon             1.5.2     2022-09-29 [1] CRAN (R 4.3.0)
 data.table         1.14.8    2023-02-17 [1] CRAN (R 4.3.0)
 DBI                1.1.3     2022-06-18 [1] CRAN (R 4.3.0)
 deldir             1.0-6     2021-10-23 [1] CRAN (R 4.3.0)
 devtools         * 2.4.5     2022-10-11 [1] CRAN (R 4.3.0)
 digest             0.6.31    2022-12-11 [1] CRAN (R 4.3.0)
 dplyr            * 1.1.2     2023-04-20 [1] CRAN (R 4.3.0)
 ellipsis           0.3.2     2021-04-29 [1] CRAN (R 4.3.0)
 fansi              1.0.4     2023-01-22 [1] CRAN (R 4.3.0)
 fastmap            1.1.1     2023-02-24 [1] CRAN (R 4.3.0)
 fastmatch          1.1-3     2021-07-23 [1] CRAN (R 4.3.0)
 fitdistrplus       1.1-11    2023-04-25 [1] CRAN (R 4.3.0)
 forcats          * 1.0.0     2023-01-29 [1] CRAN (R 4.3.0)
 fs                 1.6.2     2023-04-25 [1] CRAN (R 4.3.0)
 future             1.32.0    2023-03-07 [1] CRAN (R 4.3.0)
 future.apply       1.10.0    2022-11-05 [1] CRAN (R 4.3.0)
 generics           0.1.3     2022-07-05 [1] CRAN (R 4.3.0)
 GenomeInfoDb     * 1.36.0    2023-05-08 [1] Bioconductor
 GenomeInfoDbData   1.2.10    2023-05-16 [1] Bioconductor
 GenomicRanges    * 1.52.0    2023-05-08 [1] Bioconductor
 ggplot2          * 3.4.2     2023-04-03 [1] CRAN (R 4.3.0)
 ggrepel            0.9.3     2023-02-03 [1] CRAN (R 4.3.0)
 ggridges           0.5.4     2022-09-26 [1] CRAN (R 4.3.0)
 globals            0.16.2    2022-11-21 [1] CRAN (R 4.3.0)
 glue               1.6.2     2022-02-24 [1] CRAN (R 4.3.0)
 GO.db            * 3.17.0    2023-05-16 [1] Bioconductor
 goftest            1.2-3     2021-10-07 [1] CRAN (R 4.3.0)
 graph            * 1.78.0    2023-05-08 [1] Bioconductor
 gridExtra          2.3       2017-09-09 [1] CRAN (R 4.3.0)
 gtable             0.3.3     2023-03-21 [1] CRAN (R 4.3.0)
 hms                1.1.3     2023-03-21 [1] CRAN (R 4.3.0)
 htmltools          0.5.5     2023-03-23 [1] CRAN (R 4.3.0)
 htmlwidgets        1.6.2     2023-03-17 [1] CRAN (R 4.3.0)
 httpuv             1.6.9     2023-02-14 [1] CRAN (R 4.3.0)
 httr               1.4.5     2023-02-24 [1] CRAN (R 4.3.0)
 ica                1.0-3     2022-07-08 [1] CRAN (R 4.3.0)
 igraph             1.4.2     2023-04-07 [1] CRAN (R 4.3.0)
 IRanges          * 2.34.0    2023-05-08 [1] Bioconductor
 irlba              2.3.5.1   2022-10-03 [1] CRAN (R 4.3.0)
 jsonlite           1.8.4     2022-12-06 [1] CRAN (R 4.3.0)
 KEGGREST           1.40.0    2023-05-08 [1] Bioconductor
 KernSmooth         2.23-20   2021-05-03 [1] CRAN (R 4.3.0)
 knitr              1.42      2023-01-25 [1] CRAN (R 4.3.0)
 later              1.3.1     2023-05-02 [1] CRAN (R 4.3.0)
 lattice            0.21-8    2023-04-05 [1] CRAN (R 4.3.0)
 lazyeval           0.2.2     2019-03-15 [1] CRAN (R 4.3.0)
 leiden             0.4.3     2022-09-10 [1] CRAN (R 4.3.0)
 lifecycle          1.0.3     2022-10-07 [1] CRAN (R 4.3.0)
 listenv            0.9.0     2022-12-16 [1] CRAN (R 4.3.0)
 lmtest             0.9-40    2022-03-21 [1] CRAN (R 4.3.0)
 lubridate        * 1.9.2     2023-02-10 [1] CRAN (R 4.3.0)
 magrittr           2.0.3     2022-03-30 [1] CRAN (R 4.3.0)
 MASS               7.3-58.4  2023-03-07 [1] CRAN (R 4.3.0)
 Matrix             1.5-4     2023-04-04 [1] CRAN (R 4.3.0)
 matrixStats        0.63.0    2022-11-18 [1] CRAN (R 4.3.0)
 memoise            2.0.1     2021-11-26 [1] CRAN (R 4.3.0)
 mime               0.12      2021-09-28 [1] CRAN (R 4.3.0)
 miniUI             0.1.1.1   2018-05-18 [1] CRAN (R 4.3.0)
 munsell            0.5.0     2018-06-12 [1] CRAN (R 4.3.0)
 nlme               3.1-162   2023-01-31 [1] CRAN (R 4.3.0)
 pacman             0.5.1     2019-03-11 [1] CRAN (R 4.3.0)
 parallelly         1.35.0    2023-03-23 [1] CRAN (R 4.3.0)
 patchwork          1.1.2     2022-08-19 [1] CRAN (R 4.3.0)
 pbapply            1.7-0     2023-01-13 [1] CRAN (R 4.3.0)
 pillar             1.9.0     2023-03-22 [1] CRAN (R 4.3.0)
 pkgbuild           1.4.0     2022-11-27 [1] CRAN (R 4.3.0)
 pkgconfig          2.0.3     2019-09-22 [1] CRAN (R 4.3.0)
 pkgload            1.3.2     2022-11-16 [1] CRAN (R 4.3.0)
 plotly             4.10.1    2022-11-07 [1] CRAN (R 4.3.0)
 plyr               1.8.8     2022-11-11 [1] CRAN (R 4.3.0)
 png                0.1-8     2022-11-29 [1] CRAN (R 4.3.0)
 polyclip           1.10-4    2022-10-20 [1] CRAN (R 4.3.0)
 prettyunits        1.1.1     2020-01-24 [1] CRAN (R 4.3.0)
 processx           3.8.1     2023-04-18 [1] CRAN (R 4.3.0)
 profvis            0.3.8     2023-05-02 [1] CRAN (R 4.3.0)
 progressr          0.13.0    2023-01-10 [1] CRAN (R 4.3.0)
 promises           1.2.0.1   2021-02-11 [1] CRAN (R 4.3.0)
 ps                 1.7.5     2023-04-18 [1] CRAN (R 4.3.0)
 purrr            * 1.0.1     2023-01-10 [1] CRAN (R 4.3.0)
 R6                 2.5.1     2021-08-19 [1] CRAN (R 4.3.0)
 RANN               2.6.1     2019-01-08 [1] CRAN (R 4.3.0)
 RColorBrewer       1.1-3     2022-04-03 [1] CRAN (R 4.3.0)
 Rcpp               1.0.10    2023-01-22 [1] CRAN (R 4.3.0)
 RcppAnnoy          0.0.20    2022-10-27 [1] CRAN (R 4.3.0)
 RcppRoll           0.3.0     2018-06-05 [1] CRAN (R 4.3.0)
 RCurl              1.98-1.12 2023-03-27 [1] CRAN (R 4.3.0)
 readr            * 2.1.4     2023-02-10 [1] CRAN (R 4.3.0)
 remotes          * 2.4.2     2021-11-30 [1] CRAN (R 4.3.0)
 reshape2           1.4.4     2020-04-09 [1] CRAN (R 4.3.0)
 reticulate         1.28      2023-01-27 [1] CRAN (R 4.3.0)
 rlang              1.1.1     2023-04-28 [1] CRAN (R 4.3.0)
 ROCR               1.0-11    2020-05-02 [1] CRAN (R 4.3.0)
 Rsamtools          2.16.0    2023-04-25 [1] Bioconductor
 RSQLite            2.3.1     2023-04-03 [1] CRAN (R 4.3.0)
 rstudioapi         0.14      2022-08-22 [1] CRAN (R 4.3.0)
 Rtsne              0.16      2022-04-17 [1] CRAN (R 4.3.0)
 S4Vectors        * 0.38.1    2023-05-08 [1] Bioconductor
 scales             1.2.1     2022-08-20 [1] CRAN (R 4.3.0)
 scattermore        1.0       2023-05-03 [1] CRAN (R 4.3.0)
 sctransform        0.3.5     2022-09-21 [1] CRAN (R 4.3.0)
 sessioninfo        1.2.2     2021-12-06 [1] CRAN (R 4.3.0)
 Seurat           * 4.3.0     2022-11-18 [1] CRAN (R 4.3.0)
 SeuratObject     * 4.1.3     2022-11-07 [1] CRAN (R 4.3.0)
 shiny              1.7.4     2022-12-15 [1] CRAN (R 4.3.0)
 Signac           * 1.9.0     2022-12-08 [1] CRAN (R 4.3.0)
 sp                 1.6-0     2023-01-19 [1] CRAN (R 4.3.0)
 SparseM          * 1.81      2021-02-18 [1] CRAN (R 4.3.0)
 spatstat.data      3.0-1     2023-03-12 [1] CRAN (R 4.3.0)
 spatstat.explore   3.1-0     2023-03-14 [1] CRAN (R 4.3.0)
 spatstat.geom      3.1-0     2023-03-12 [1] CRAN (R 4.3.0)
 spatstat.random    3.1-4     2023-03-13 [1] CRAN (R 4.3.0)
 spatstat.sparse    3.0-1     2023-03-12 [1] CRAN (R 4.3.0)
 spatstat.utils     3.0-2     2023-03-11 [1] CRAN (R 4.3.0)
 stringi            1.7.12    2023-01-11 [1] CRAN (R 4.3.0)
 stringr          * 1.5.0     2022-12-02 [1] CRAN (R 4.3.0)
 survival           3.5-5     2023-03-12 [1] CRAN (R 4.3.0)
 tensor             1.5       2012-05-05 [1] CRAN (R 4.3.0)
 tibble           * 3.2.1     2023-03-20 [1] CRAN (R 4.3.0)
 tidyr            * 1.3.0     2023-01-24 [1] CRAN (R 4.3.0)
 tidyselect         1.2.0     2022-10-10 [1] CRAN (R 4.3.0)
 tidyverse        * 2.0.0     2023-02-22 [1] CRAN (R 4.3.0)
 timechange         0.2.0     2023-01-11 [1] CRAN (R 4.3.0)
 topGO            * 2.52.0    2023-05-08 [1] Bioconductor
 tzdb               0.3.0     2022-03-28 [1] CRAN (R 4.3.0)
 urlchecker         1.0.1     2021-11-30 [1] CRAN (R 4.3.0)
 usethis          * 2.1.6     2022-05-25 [1] CRAN (R 4.3.0)
 utf8               1.2.3     2023-01-31 [1] CRAN (R 4.3.0)
 uwot               0.1.14    2022-08-22 [1] CRAN (R 4.3.0)
 vctrs              0.6.2     2023-04-19 [1] CRAN (R 4.3.0)
 viridisLite        0.4.2     2023-05-02 [1] CRAN (R 4.3.0)
 withr              2.5.0     2022-03-03 [1] CRAN (R 4.3.0)
 xfun               0.39      2023-04-20 [1] CRAN (R 4.3.0)
 xtable             1.8-4     2019-04-21 [1] CRAN (R 4.3.0)
 XVector            0.40.0    2023-05-08 [1] Bioconductor
 zlibbioc           1.46.0    2023-05-08 [1] Bioconductor
 zoo                1.8-12    2023-04-13 [1] CRAN (R 4.3.0)

 [1] /Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library
timoast commented 1 year ago

Unfortunately I wasn't able to reproduce this error, could you show the output of traceback() after the error occurs?

mossconfuse commented 1 year ago 
> ppsn <- TSSEnrichment(ppsn, fast = FALSE)
Extracting TSS positions
Finding + strand cut sites
Finding - strand cut sites
Error: non-conformable matrix dimensions in dimCheck(e1, e2)
> traceback()
8: stop(gettextf("non-conformable matrix dimensions in %s", deparse(sys.call(sys.parent()))), 
       call. = FALSE, domain = NA)
7: checkDim(dim(a), dim(b))
6: dimCheck(e1, e2)
5: .Arith.Csparse(e1, e2, .Generic, class. = "dgCMatrix")
4: cut.matrix.plus + cut.matrix.minus[, rev(x = colnames(x = cut.matrix.minus))]
3: cut.matrix.plus + cut.matrix.minus[, rev(x = colnames(x = cut.matrix.minus))]
2: CreateRegionPileupMatrix(object = object, regions = tss.positions, 
       assay = assay, cells = cells, verbose = verbose)
1: TSSEnrichment(ppsn, fast = FALSE)
timoast commented 1 year ago

Would you mind sharing the seurat object and associated fragment file that reproduces this error? Then I can look into it more. You can email a download link to stuartt@gis.a-star.edu.sg

mossconfuse commented 1 year ago

Thank you for responding. I solved the problem, partially. I reran cellranger-arc mkref on a gtf annotation with only chromosomal annotations (no mitochondrial, chloroplast or scaffold annotations) and reran cellranger-arc count. Now Signac::TSSEnrichment(object, fast = F) runs smoothly. I had a look at the code for TSSEnrichment and there is some removal mechanism for mitochondrial features, but my scaffold and chloroplast features may have been what caused the hiccups.

timoast commented 1 year ago

Ok that's good to know, thanks