FractionCountsInRegion() returns values > 1

[Almr95]

Hello!

I’m currently learning to use Signac for analysing single-cell chromatin data from histone marks, and while exploring quality metrics like fraction of reads in enhancer and promoter regions I noticed that FractionCountsInRegion() sometimes returns values greater than 1.

I initially thought this might reflect some normalization step, but when I looked into the source code I saw that the helper function CountsInRegion() uses queryHits(findOverlaps(...)) to subset the count matrix. I wonder if this might accidentally count the same bin multiple times when a single bin overlaps several regions (e.g., enhancers). Of course, it’s possible that I’m misunderstanding the intended behaviour — if this design is correct for a reason I’ve missed, I’d really appreciate clarification.

In the toy example below, all five enhancer intervals overlap the same bin. FractionCountsInRegion() seems to sum that bin five times (once per overlap), producing fractions > 1. When I collapse to unique bin indices (unique(queryHits)), the result becomes ≤ 1, which in my mind matches the idea of “fraction of reads in regions.”

Thank you very much for your time and for maintaining this great package

suppressPackageStartupMessages({
  library(Signac)
  library(Seurat)
  library(GenomicRanges)
  library(IRanges)
  library(Matrix)
})

# Bins
bins <- GRanges(
  seqnames = "chr1",
  ranges = IRanges(start = c(1, 101, 201), width = 100),
  strand = "*"
)
names(bins) <- paste0("bin", 1:3)

# Counts
counts <- Matrix(c(
  10,  5,   # bin1
  50, 60,   # bin2  (dominant bin)
   5,  5    # bin3
), nrow = 3, byrow = TRUE, sparse = TRUE)
rownames(counts) <- names(bins)
colnames(counts) <- c("cellA", "cellB")

# Build assay & object
chrom_assay <- CreateChromatinAssay(counts = counts, ranges = bins, genome = "hg38")
obj <- CreateSeuratObject(counts = chrom_assay, assay = "peaks")

# Enhancer regions – all five fall inside bin2 (chr1:101–200)
enhancers <- GRanges(
  seqnames = rep("chr1", 5),
  ranges = IRanges(start = c(120, 130, 140, 150, 160), width = 10),
  strand = "*"
)

# Fraction of Reads in Enhancers
FractionCountsInRegion(obj, enhancers)
# cellA 3.846154   cellB 4.285714  → >1

# Internal check
obj.gr <- GetAssayData(obj, assay = "peaks", layer = "ranges")
ov <- findOverlaps(query = obj.gr, subject = enhancers, ignore.strand = TRUE)
hit.regions.dup <- queryHits(ov)          # c(2, 2, 2, 2, 2)
hit.regions.uni <- unique(hit.regions.dup)  # 2

num_dup <- colSums(GetAssayData(obj, assay = "peaks", layer = "counts")[hit.regions.dup, , drop = FALSE])
num_uni <- colSums(GetAssayData(obj, assay = "peaks", layer = "counts")[hit.regions.uni, , drop = FALSE])
den <- colSums(GetAssayData(obj, assay = "peaks", layer = "counts"))

FR_dup <- num_dup / den
FR_uni <- num_uni / den

FR_dup
# cellA 3.846154   cellB 4.285714
FR_uni
# cellA 0.7692308  cellB 0.8571429

sessionInfo()

R version 4.5.1 (2025-06-13)
Platform: aarch64-apple-darwin20
Running under: macOS Sequoia 15.6

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.5-arm64/Resources/lib/libRlapack.dylib;  LAPACK version 3.12.1

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: Europe/Paris
tzcode source: internal

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

other attached packages:
 [1] Matrix_1.7-4                      future.apply_1.20.0               future_1.67.0                     writexl_1.5.4                    
 [5] readxl_1.4.5                      mclust_6.1.1                      slingshot_2.16.0                  TrajectoryUtils_1.16.1           
 [9] SingleCellExperiment_1.30.1       SummarizedExperiment_1.38.1       MatrixGenerics_1.20.0             matrixStats_1.5.0                
[13] princurve_2.1.6                   tidygraph_1.3.1                   igraph_2.1.4                      fgsea_1.34.2                     
[17] decoupleR_2.14.0                  msigdbr_25.1.1                    sctransform_0.4.2                 Signac_1.16.0                    
[21] SeuratDisk_0.0.0.9021             glmGamPoi_1.20.0                  presto_1.0.0                      data.table_1.17.8                
[25] BPCells_0.3.1                     Seurat_5.3.0                      SeuratObject_5.2.0                sp_2.2-0                         
[29] biovizBase_1.56.0                 ensembldb_2.32.0                  AnnotationFilter_1.32.0           GenomicFeatures_1.60.0           
[33] AnnotationDbi_1.70.0              Biobase_2.68.0                    AnnotationHub_3.16.1              BiocFileCache_2.16.2             
[37] dbplyr_2.5.1                      BSgenome.Hsapiens.UCSC.hg38_1.4.5 BSgenome_1.76.0                   rtracklayer_1.68.0               
[41] BiocIO_1.18.0                     Rsamtools_2.24.1                  Biostrings_2.76.0                 XVector_0.48.0                   
[45] GenomicRanges_1.60.0              GenomeInfoDb_1.44.3               IRanges_2.42.0                    S4Vectors_0.46.0                 
[49] BiocGenerics_0.54.0               generics_0.1.4                    clustree_0.5.1                    ggraph_2.2.2                     
[53] gridExtra_2.3                     ggtree_3.16.3                     paletteer_1.6.0                   viridis_0.6.5                    
[57] viridisLite_0.4.2                 pheatmap_1.0.13                   circlize_0.4.16                   ComplexHeatmap_2.24.1            
[61] ComplexUpset_1.3.3                ggupset_0.4.1                     UpSetR_1.4.0                      ggVennDiagram_1.5.4              
[65] patchwork_1.3.2                   magrittr_2.0.4                    lubridate_1.9.4                   forcats_1.0.1                    
[69] stringr_1.5.2                     purrr_1.1.0                       readr_2.1.5                       tidyr_1.3.1                      
[73] tibble_3.3.0                      ggplot2_4.0.0                     tidyverse_2.0.0                   dplyr_1.1.4                      
[77] proxy_0.4-27                      Rfast_2.1.5.2                     RcppParallel_5.1.11-1             zigg_0.0.2                       
[81] Rcpp_1.1.0                        qs2_0.1.5                         fs_1.6.6                         

loaded via a namespace (and not attached):
  [1] dichromat_2.0-0.1        nnet_7.3-20              goftest_1.2-3            vctrs_0.6.5              spatstat.random_3.4-2   
  [6] digest_0.6.37            png_0.1-8                shape_1.4.6.1            ggrepel_0.9.6            deldir_2.0-4            
 [11] parallelly_1.45.1        MASS_7.3-65              reshape2_1.4.4           foreach_1.5.2            httpuv_1.6.16           
 [16] withr_3.0.2              xfun_0.53                ggfun_0.2.0              survival_3.8-3           memoise_2.0.1           
 [21] tidytree_0.4.6           zoo_1.8-14               GlobalOptions_0.1.2      pbapply_1.7-4            Formula_1.2-5           
 [26] rematch2_2.1.2           KEGGREST_1.48.1          promises_1.3.3           httr_1.4.7               restfulr_0.0.16         
 [31] globals_0.18.0           fitdistrplus_1.2-4       stringfish_0.17.0        rstudioapi_0.17.1        UCSC.utils_1.4.0        
 [36] miniUI_0.1.2             base64enc_0.1-3          babelgene_22.9           curl_7.0.0               polyclip_1.10-7         
 [41] GenomeInfoDbData_1.2.14  SparseArray_1.8.1        xtable_1.8-4             doParallel_1.0.17        evaluate_1.0.5          
 [46] S4Arrays_1.8.1           hms_1.1.3                irlba_2.3.5.1            colorspace_2.1-2         filelock_1.0.3          
 [51] hdf5r_1.3.12             ROCR_1.0-11              reticulate_1.43.0        spatstat.data_3.1-8      lmtest_0.9-40           
 [56] later_1.4.4              lattice_0.22-7           spatstat.geom_3.6-0      scattermore_1.2          XML_3.99-0.19           
 [61] cowplot_1.2.0            RcppAnnoy_0.0.22         Hmisc_5.2-4              pillar_1.11.1            nlme_3.1-168            
 [66] iterators_1.0.14         compiler_4.5.1           beachmat_2.24.0          RSpectra_0.16-2          stringi_1.8.7           
 [71] tensor_1.5.1             GenomicAlignments_1.44.0 plyr_1.8.9               crayon_1.5.3             abind_1.4-8             
 [76] gridGraphics_0.5-1       graphlayouts_1.2.2       bit_4.6.0                fastmatch_1.1-6          codetools_0.2-20        
 [81] GetoptLong_1.0.5         plotly_4.11.0            mime_0.13                splines_4.5.1            fastDummies_1.7.5       
 [86] prismatic_1.1.2          cellranger_1.1.0         utf8_1.2.6               knitr_1.50               blob_1.2.4              
 [91] clue_0.3-66              BiocVersion_3.21.1       listenv_0.9.1            checkmate_2.3.3          ggplotify_0.1.3         
 [96] tzdb_0.5.0               tweenr_2.0.3             pkgconfig_2.0.3          tools_4.5.1              cachem_1.1.0            
[101] RSQLite_2.4.3            DBI_1.2.3                fastmap_1.2.0            rmarkdown_2.30           scales_1.4.0            
[106] ica_1.0-3                BiocManager_1.30.26      dotCall64_1.2            VariantAnnotation_1.54.1 RANN_2.6.2              
[111] rpart_4.1.24             farver_2.1.2             yaml_2.3.10              foreign_0.8-90           cli_3.6.5               
[116] lifecycle_1.0.4          uwot_0.2.3               backports_1.5.0          BiocParallel_1.42.2      timechange_0.3.0        
[121] gtable_0.3.6             rjson_0.2.23             ggridges_0.5.7           progressr_0.16.0         parallel_4.5.1          
[126] ape_5.8-1                jsonlite_2.0.0           RcppHNSW_0.6.0           bitops_1.0-9             bit64_4.6.0-1           
[131] assertthat_0.2.1         Rtsne_0.17               yulab.utils_0.2.1        spatstat.utils_3.2-0     spatstat.univar_3.1-4   
[136] lazyeval_0.2.2           shiny_1.11.1             htmltools_0.5.8.1        rappdirs_0.3.3           glue_1.8.0              
[141] spam_2.11-1              RCurl_1.98-1.17          treeio_1.32.0            R6_2.6.1                 labeling_0.4.3          
[146] RcppRoll_0.3.1           cluster_2.1.8.1          pkgload_1.4.1            aplot_0.2.9              DelayedArray_0.34.1     
[151] tidyselect_1.2.1         ProtGenerics_1.40.0      htmlTable_2.4.3          ggforce_0.5.0            KernSmooth_2.23-26      
[156] S7_0.2.0                 htmlwidgets_1.6.4        RColorBrewer_1.1-3       rlang_1.1.6              spatstat.sparse_3.1-0   
[161] spatstat.explore_3.5-3  

[shayanjl]

Yes, I think your interpretation is correct. CountsInRegion() is documented as counting reads per cell overlapping a set of regions, and the current Signac source uses findOverlaps(...), takes queryHits(overlaps), subsets the count matrix with those hits, and then sums them. Because those hits are not made unique, the same bin/feature can be counted multiple times if it overlaps multiple supplied regions. That means FractionCountsInRegion() can exceed 1 when the region set itself contains overlapping intervals.

So in practice, if you want a true “fraction of reads in regions,” it is safer to first collapse overlapping regions, for example with reduce(regions), or otherwise deduplicate the overlapping feature indices before summing. For disjoint region sets, the current function should behave as expected.