FeatureMatrix silently generates empty matrix when used with merged Seurat object

[aabarug]

Hi there,

I was attempting to use a combination of merging and integration in a workflow with Signac. I have two datasets D1 (samples 1 to 4) and D2 (samples 5 to 6). There multiple different un-replicated samples per dataset. Additionally, the two datasets have different experimental setups. I was originally working on the two datasets separately, having merged the samples together to create two Seurat objects, until my collaborator expressed interest in seeing them together.

I wasn't quite sure how to attempt the integration workflow (take the intersection of all individual samples or take the intersection of the merged datasets). I attempted to use the workflow on the merged objects but that failed, so I ended up merging all samples together for both datasets for the sake of moving on with the analysis. When I was investigating why integration was failing, I noticed that the FeatureMatrix step returns an empty matrix when used on a merged object:

library(Signac)
library(Seurat)

# unifiedpeaks is the assay I generated when taking the union of peaks as per the merge workflow
combined_d1 <- readRDS("data/combined_d1.RDS")
DefaultAssay(combined_d1) <- "unifiedpeaks"
combined_d2 <- readRDS("data/combined_d2.RDS")
DefaultAssay(combined_nem) <- "unifiedpeaks"

# find peaks that intersect in both datasets
intersecting_peaks <- findOverlaps(query = combined_d1, subject = combined_d2)

# find the coordinates of peaks in the intersecting peaks
intersections_loc <- unique(queryHits(intersecting_peaks)) %>% sort()

# choose a subset of intersecting peaks
subset_peaks <- sort(granges(combined_d1)[sample(intersections_loc, size = 10000, replace = FALSE)])

## count fragments per cell overlapping the intersecting peaks
peaks_mat_d2 <- FeatureMatrix(
  fragments = Fragments(combined_d2),
  features = subset_peaks,
  cells = colnames(combined_d2)
)

FeatureMatrix runs without an error message but peaks_mat_d2 is completely empty.

> peaks_mat_d2 %>% sum()
[1] 0

I've checked the same step with one of the unmerged Seurat objects with the subset_peaks and it ran fine and produced a matrix with counts in it, so I know that there should be counts for subset_peaks found for the merged object.

# where filtered_samp5 is one of my individual Seurat objects
frag_counts <- FeatureMatrix(
  fragments = Fragments(filtered_samp5),
  features = subset_peaks,
  cells = colnames(filtered_samp5)
)

> frag_counts %>% sum()
[1] 2507048

Both the combined_d1 and combined_d2 objects work fine - I have been able to use them for label transfers and to generate coverage plots.

I tried to take a look at what the FeatureMatrix function itself was doing and each matrix inside mat.list was empty. So something is going wrong with the SingleFeatureMatrix step and I haven't investigated further. I am using the dev version of Signac as of last week, when I was having a similar problem with the GeneActivity step for the merged object and saw someone had filed #205, which fixed that issue for me.

My sessionInfo:

R version 4.0.2 (2020-06-22)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Ubuntu 16.04.7 LTS

Matrix products: default
BLAS:   /mnt/software/apps/R/4.0.2/lib/R/lib/libRblas.so
LAPACK: /mnt/software/apps/R/4.0.2/lib/R/lib/libRlapack.so

locale:
 [1] LC_CTYPE=en_AU.UTF-8       LC_NUMERIC=C               LC_TIME=en_AU.UTF-8        LC_COLLATE=en_AU.UTF-8    
 [5] LC_MONETARY=en_AU.UTF-8    LC_MESSAGES=en_AU.UTF-8    LC_PAPER=en_AU.UTF-8       LC_NAME=C                 
 [9] LC_ADDRESS=C               LC_TELEPHONE=C             LC_MEASUREMENT=en_AU.UTF-8 LC_IDENTIFICATION=C       

attached base packages:
[1] stats4    parallel  stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
 [1] Matrix_1.2-18              pbapply_1.4-2              future_1.18.0              Rsamtools_2.4.0           
 [5] Biostrings_2.56.0          XVector_0.28.0             magrittr_1.5               ggplot2_3.3.2             
 [9] EnsDb.Mmusculus.v79_2.99.0 ensembldb_2.12.1           AnnotationFilter_1.12.0    GenomicFeatures_1.40.1    
[13] AnnotationDbi_1.50.3       Biobase_2.48.0             GenomicRanges_1.40.0       GenomeInfoDb_1.24.2       
[17] IRanges_2.22.2             S4Vectors_0.26.1           BiocGenerics_0.34.0        Seurat_3.2.0              
[21] Signac_1.0.0              

loaded via a namespace (and not attached):
  [1] reticulate_1.16             tidyselect_1.1.0            RSQLite_2.2.0               htmlwidgets_1.5.1          
  [5] grid_4.0.2                  BiocParallel_1.22.0         Rtsne_0.15                  munsell_0.5.0              
  [9] codetools_0.2-16            ica_1.0-2                   miniUI_0.1.1.1              withr_2.2.0                
 [13] colorspace_1.4-1            OrganismDbi_1.30.0          knitr_1.29                  rstudioapi_0.11            
 [17] ROCR_1.0-11                 tensor_1.5                  listenv_0.8.0               GenomeInfoDbData_1.2.3     
 [21] polyclip_1.10-0             bit64_0.9-7.1               vctrs_0.3.2                 generics_0.0.2             
 [25] xfun_0.16                   biovizBase_1.36.0           BiocFileCache_1.12.1        lsa_0.73.2                 
 [29] ggseqlogo_0.1               R6_2.4.1                    rsvd_1.0.3                  bitops_1.0-6               
 [33] spatstat.utils_1.17-0       reshape_0.8.8               DelayedArray_0.14.1         assertthat_0.2.1           
 [37] promises_1.1.1              scales_1.1.1                nnet_7.3-14                 gtable_0.3.0               
 [41] globals_0.12.5              goftest_1.2-2               ggbio_1.36.0                rlang_0.4.7                
 [45] RcppRoll_0.3.0              splines_4.0.2               rtracklayer_1.48.0          lazyeval_0.2.2             
 [49] acepack_1.4.1               dichromat_2.0-0             checkmate_2.0.0             BiocManager_1.30.10        
 [53] reshape2_1.4.4              abind_1.4-5                 backports_1.1.8             httpuv_1.5.4               
 [57] qvalue_2.20.0               Hmisc_4.4-0                 RBGL_1.64.0                 tools_4.0.2                
 [61] ellipsis_0.3.1              RColorBrewer_1.1-2          ggridges_0.5.2              Rcpp_1.0.5                 
 [65] plyr_1.8.6                  base64enc_0.1-3             progress_1.2.2              zlibbioc_1.34.0            
 [69] purrr_0.3.4                 RCurl_1.98-1.2              prettyunits_1.1.1           rpart_4.1-15               
 [73] openssl_1.4.2               deldir_0.1-28               cowplot_1.0.0               zoo_1.8-8                  
 [77] SummarizedExperiment_1.18.2 ggrepel_0.8.2               cluster_2.1.0               fs_1.4.2                   
 [81] data.table_1.12.8           lmtest_0.9-37               RANN_2.6.1                  SnowballC_0.7.0            
 [85] packrat_0.5.0               ProtGenerics_1.20.0         fitdistrplus_1.1-1          matrixStats_0.56.0         
 [89] hms_0.5.3                   patchwork_1.0.1             mime_0.9                    evaluate_0.14              
 [93] xtable_1.8-4                XML_3.99-0.5                jpeg_0.1-8.1                gridExtra_2.3              
 [97] compiler_4.0.2              biomaRt_2.44.1              tibble_3.0.3                KernSmooth_2.23-17         
[101] crayon_1.3.4                htmltools_0.5.0             mgcv_1.8-31                 later_1.1.0.1              
[105] Formula_1.2-3               tidyr_1.1.0                 DBI_1.1.0                   dbplyr_1.4.4               
[109] MASS_7.3-51.6               rappdirs_0.3.1              igraph_1.2.5                pkgconfig_2.0.3            
[113] GenomicAlignments_1.24.0    foreign_0.8-80              plotly_4.9.2.1              stringr_1.4.0              
[117] VariantAnnotation_1.34.0    digest_0.6.25               sctransform_0.2.1           RcppAnnoy_0.0.16           
[121] graph_1.66.0                spatstat.data_1.4-3         rmarkdown_2.3               leiden_0.3.3               
[125] fastmatch_1.1-0             htmlTable_2.0.1             uwot_0.1.8                  curl_4.3                   
[129] shiny_1.5.0                 lifecycle_0.2.0             nlme_3.1-148                jsonlite_1.7.0             
[133] viridisLite_0.3.0           askpass_1.1                 BSgenome_1.56.0             pillar_1.4.6               
[137] lattice_0.20-41             GGally_2.0.0                fastmap_1.0.1               httr_1.4.2                 
[141] survival_3.1-12             glue_1.4.1                  spatstat_1.64-1             png_0.1-7                  
[145] bit_1.1-15.2                stringi_1.4.6               blob_1.2.1                  latticeExtra_0.6-29        
[149] memoise_1.1.0               dplyr_1.0.0                 irlba_2.3.3                 future.apply_1.6.0         
[153] ape_5.4 

[ninal88]

Commenting because I have the same issue trying to follow the merge vignette. FeatureMatrix works fine when following the analyzing adult mouse vignette. From there I saved the filtered files as a .RDS and fed them into the merge objects steps and FeatureMatrix returns NULL.

[timoast]

This is fixed on the develop branch (see #205)

[ninal88]

Hi @timoast, I updated to the develop version of Signac successfully and the problem still persists.

[timoast]

@ninal88 try restarting R after updating Signac and verify that running packageVersion("Signac") gives 1.0.0.9000

If you can verify that you still have this issue after restarting R please open a new issue and include your full code and the output of sessionInfo()

[aabarug]

@timoast - bit puzzled why I wasn't on the dev version because I had restarted R a couple of times since updating and the GeneActivity function was working ¯_(ツ)_/¯ so I assumed I was running the dev version. Anyway FeatureMatrix is also working on merged objects. Thank you for the help!