ClosestFeature function returning fewer ranges than queried

[sabrinacamp2]

Hi!

Thank you for making this analytical toolset.

I am attempting to use the ClosestFeature function to annotate my FindMarkers result with gene names. However, I keep receiving this error

Error in $<-.data.frame(*tmp*, "query_region", value = c("chr1-208176523-208178576", : replacement has 87552 rows, data has 87532

Digging into the function, I found that the Genomic Ranges distanceToNearest function is not returning the same amount of rows as the number of regions I am querying. This function is returning 87532 rows, while the number of regions I queried is 87552. This seems to cause issues later in the function when a column is trying to be added to the dataframe adding my queried regions, which now have more rows than the dataframe. I would assume distanceToNearest should be filling in regions where there is no result with "NA" to keep the rows the same length, but it's not.

## example commands
ranges_da_peaks <- rownames(da_peaks[da_peaks$p_val_adj < 0.05, ])
ClosestFeature(object, regions = ranges_da_peaks )

Session info:
R version 4.1.0 (2021-05-18)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Ubuntu 18.04.5 LTS

Matrix products: default
BLAS: /usr/lib/x86_64-linux-gnu/openblas/libblas.so.3
LAPACK: /usr/lib/x86_64-linux-gnu/libopenblasp-r0.2.20.so

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

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

other attached packages:
[1] RColorBrewer_1.1-2 future_1.21.0
[3] hdf5r_1.3.3 patchwork_1.1.1
[5] ggplot2_3.3.5 EnsDb.Hsapiens.v86_2.99.0
[7] ensembldb_2.16.2 AnnotationFilter_1.16.0
[9] GenomicFeatures_1.44.0 AnnotationDbi_1.54.1
[11] Biobase_2.52.0 GenomicRanges_1.44.0
[13] GenomeInfoDb_1.28.1 IRanges_2.26.0
[15] S4Vectors_0.30.0 BiocGenerics_0.38.0
[17] Signac_1.3.0 SeuratObject_4.0.2
[19] Seurat_4.0.3

loaded via a namespace (and not attached):
[1] utf8_1.2.1 reticulate_1.20
[3] tidyselect_1.1.1 RSQLite_2.2.7
[5] htmlwidgets_1.5.3 grid_4.1.0
[7] docopt_0.7.1 BiocParallel_1.26.1
[9] Rtsne_0.15 munsell_0.5.0
[11] codetools_0.2-18 ica_1.0-2
[13] pbdZMQ_0.3-5 miniUI_0.1.1.1
[15] withr_2.4.2 colorspace_2.0-2
[17] filelock_1.0.2 uuid_0.1-4
[19] ROCR_1.0-11 tensor_1.5
[21] listenv_0.8.0 MatrixGenerics_1.4.0
[23] labeling_0.4.2 slam_0.1-48
[25] repr_1.1.3 GenomeInfoDbData_1.2.6
[27] polyclip_1.10-0 bit64_4.0.5
[29] farver_2.1.0 parallelly_1.26.1
[31] vctrs_0.3.8 generics_0.1.0
[33] BiocFileCache_2.0.0 lsa_0.73.2
[35] ggseqlogo_0.1 R6_2.5.0
[37] bitops_1.0-7 spatstat.utils_2.2-0
[39] cachem_1.0.5 DelayedArray_0.18.0
[41] assertthat_0.2.1 promises_1.2.0.1
[43] BiocIO_1.2.0 scales_1.1.1
[45] gtable_0.3.0 globals_0.14.0
[47] goftest_1.2-2 rlang_0.4.11
[49] RcppRoll_0.3.0 splines_4.1.0
[51] rtracklayer_1.52.0 lazyeval_0.2.2
[53] spatstat.geom_2.2-0 yaml_2.2.1
[55] reshape2_1.4.4 abind_1.4-5
[57] httpuv_1.6.1 tools_4.1.0
[59] ellipsis_0.3.2 spatstat.core_2.2-0
[61] ggridges_0.5.3 Rcpp_1.0.7
[63] plyr_1.8.6 base64enc_0.1-3
[65] progress_1.2.2 zlibbioc_1.38.0
[67] purrr_0.3.4 RCurl_1.98-1.3
[69] prettyunits_1.1.1 rpart_4.1-15
[71] deldir_0.2-10 pbapply_1.4-3
[73] cowplot_1.1.1 zoo_1.8-9
[75] SummarizedExperiment_1.22.0 ggrepel_0.9.1
[77] cluster_2.1.2 magrittr_2.0.1
[79] data.table_1.14.0 scattermore_0.7
[81] lmtest_0.9-38 RANN_2.6.1
[83] SnowballC_0.7.0 ProtGenerics_1.24.0
[85] fitdistrplus_1.1-5 matrixStats_0.59.0
[87] hms_1.1.0 mime_0.11
[89] evaluate_0.14 xtable_1.8-4
[91] XML_3.99-0.6 sparsesvd_0.2
[93] gridExtra_2.3 compiler_4.1.0
[95] biomaRt_2.48.2 tibble_3.1.2
[97] KernSmooth_2.23-20 crayon_1.4.1
[99] htmltools_0.5.1.1 mgcv_1.8-36
[101] later_1.2.0 tidyr_1.1.3
[103] DBI_1.1.1 tweenr_1.0.2
[105] dbplyr_2.1.1 MASS_7.3-54
[107] rappdirs_0.3.3 Matrix_1.3-4
[109] igraph_1.2.6 pkgconfig_2.0.3
[111] GenomicAlignments_1.28.0 IRdisplay_1.0
[113] plotly_4.9.4.1 spatstat.sparse_2.0-0
[115] xml2_1.3.2 XVector_0.32.0
[117] stringr_1.4.0 digest_0.6.27
[119] sctransform_0.3.2 RcppAnnoy_0.0.18
[121] spatstat.data_2.1-0 Biostrings_2.60.1
[123] leiden_0.3.8 fastmatch_1.1-3
[125] uwot_0.1.10 restfulr_0.0.13
[127] curl_4.3.2 shiny_1.6.0
[129] Rsamtools_2.8.0 rjson_0.2.20
[131] lifecycle_1.0.0 nlme_3.1-152
[133] jsonlite_1.7.2 viridisLite_0.4.0
[135] fansi_0.5.0 pillar_1.6.1
[137] lattice_0.20-44 KEGGREST_1.32.0
[139] fastmap_1.1.0 httr_1.4.2
[141] survival_3.2-11 glue_1.4.2
[143] qlcMatrix_0.9.7 png_0.1-7
[145] bit_4.0.4 ggforce_0.3.3
[147] stringi_1.6.2 blob_1.2.1
[149] memoise_2.0.0 IRkernel_1.2
[151] dplyr_1.0.7 irlba_2.3.3
[153] future.apply_1.7.0

[strawberry789]

I'm having the same issue

[timoast]

Are you running FindMarkers() or FindAllMarkers(). If you're using FindAllMarkers(), you may have some duplicate peaks that were differential in multiple comparisons. Since you use the row names of the dataframe, these will have been made unique by adding a suffix (make.unique() function), and would not be able to be turned back into genomic ranges. You can try using the peak name extracted from the gene column instead to avoid this.

[strawberry789]

I used FindAllMarkers(), but aren't the peak names the same as the rownames?

Anyhow, I tried using the peak names from the gene column, like open_c0 <- c0_peaks$gene, and then closest_genes_c0 <- ClosestFeature(integrated_res0.5, regions = open_c0) and head(closest_genes_c0), but I still got the same error.

[timoast]

Peak names will not be the same as row names if there are duplicate peaks in the dataframe (see my comment above).

Can you share the open_c0 dataframe that gives this issue, and the gene annotations stored in the integrated_res0.5 object? You can send a link to tstuart@nygenome.org

[timoast]

Thanks for sending the information, the reason for this error was that there were peaks supplied on contigs not present in the annotation (eg, GL456216.1). I have now updated ClosestFeature() to check if any regions are on contigs not present in the annotation, and remove them and show a warning message. This is available now on the develop branch.

[menage95]

Hello,

I keep getting this error even though I use 'FindMarkers' and have installed the develop version of Signac.

da_peaks <- FindMarkers(data, 
                        ident.1 = rownames(data[[]][data$dataset == "condition_1",]),
                        ident.2 = rownamesdata[[]][data$dataset == "condition_2",]),
                        min.pct = 0.05,
                        test.use = 'LR',
                        latent.vars = 'atac_peak_region_fragments'
)
da_peaks$closest_gene <- ClosestFeature(data, regions = rownames(da_peaks))$gene_name

Error in $<-.data.frame(*tmp*, closest_gene, value = c("Pogz", "AI314180", :
replacement has 2593 rows, data has 2596
In addition: Warning message:
In ClosestFeature(data, regions = rownames(da_peaks)) :
The following seqlevels present in query regions are not present
in the supplied gene annotations and will be removed: GL456216.1, JH584304.1

sessionInfo()
R version 4.3.1 (2023-06-16 ucrt)
Platform: x86_64-w64-mingw32/x64 (64-bit)
Running under: Windows 11 x64 (build 22621)

Matrix products: default

locale:
[1] LC_COLLATE=English_United Kingdom.utf8 LC_CTYPE=English_United Kingdom.utf8
[3] LC_MONETARY=English_United Kingdom.utf8 LC_NUMERIC=C
[5] LC_TIME=English_United Kingdom.utf8

time zone: Europe/London
tzcode source: internal

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

other attached packages:
[1] biovizBase_1.48.0 hdf5r_1.3.8 patchwork_1.1.3
[4] EnsDb.Mmusculus.v79_2.99.0 ensembldb_2.24.0 AnnotationFilter_1.24.0
[7] GenomicFeatures_1.52.1 AnnotationDbi_1.62.2 Biobase_2.60.0
[10] GenomicRanges_1.52.0 GenomeInfoDb_1.36.1 IRanges_2.34.1
[13] S4Vectors_0.38.1 BiocGenerics_0.46.0 ggplot2_3.4.3
[16] readr_2.1.4 Seurat_4.3.0.1 Signac_1.10.9001
[19] SeuratObject_4.1.3 sp_2.0-0

loaded via a namespace (and not attached):
[1] RcppAnnoy_0.0.20 splines_4.3.1 later_1.3.1
[4] BiocIO_1.10.0 bitops_1.0-7 filelock_1.0.2
[7] tibble_3.2.1 polyclip_1.10-4 rpart_4.1.19
[10] XML_3.99-0.14 lifecycle_1.0.3 globals_0.16.2
[13] lattice_0.21-8 MASS_7.3-60 backports_1.4.1
[16] magrittr_2.0.3 rmarkdown_2.24 Hmisc_5.1-0
[19] plotly_4.10.2 yaml_2.3.7 httpuv_1.6.11
[22] sctransform_0.3.5 spatstat.sparse_3.0-2 reticulate_1.30
[25] cowplot_1.1.1 pbapply_1.7-2 DBI_1.1.3
[28] RColorBrewer_1.1-3 abind_1.4-5 zlibbioc_1.46.0
[31] Rtsne_0.16 purrr_1.0.1 RCurl_1.98-1.12
[34] nnet_7.3-19 VariantAnnotation_1.46.0 rappdirs_0.3.3
[37] GenomeInfoDbData_1.2.10 ggrepel_0.9.3 irlba_2.3.5.1
[40] listenv_0.9.0 spatstat.utils_3.0-3 goftest_1.2-3
[43] spatstat.random_3.1-5 fitdistrplus_1.1-11 parallelly_1.36.0
[46] DelayedArray_0.26.7 leiden_0.4.3 codetools_0.2-19
[49] RcppRoll_0.3.0 xml2_1.3.5 tidyselect_1.2.0
[52] base64enc_0.1-3 matrixStats_1.0.0 BiocFileCache_2.8.0
[55] spatstat.explore_3.2-1 GenomicAlignments_1.36.0 jsonlite_1.8.5
[58] Formula_1.2-5 ellipsis_0.3.2 progressr_0.14.0
[61] ggridges_0.5.4 survival_3.5-5 tools_4.3.1
[64] progress_1.2.2 ica_1.0-3 Rcpp_1.0.10
[67] glue_1.6.2 gridExtra_2.3 xfun_0.40
[70] MatrixGenerics_1.12.3 dplyr_1.1.2 withr_2.5.0
[73] fastmap_1.1.1 fansi_1.0.4 digest_0.6.32
[76] R6_2.5.1 mime_0.12 colorspace_2.1-0
[79] scattermore_1.2 tensor_1.5 dichromat_2.0-0.1
[82] spatstat.data_3.0-1 biomaRt_2.56.1 RSQLite_2.3.1
[85] utf8_1.2.3 tidyr_1.3.0 generics_0.1.3
[88] data.table_1.14.8 rtracklayer_1.60.1 S4Arrays_1.0.5
[91] prettyunits_1.1.1 httr_1.4.7 htmlwidgets_1.6.2
[94] uwot_0.1.15 pkgconfig_2.0.3 gtable_0.3.4
[97] blob_1.2.4 lmtest_0.9-40 XVector_0.40.0
[100] htmltools_0.5.5 ProtGenerics_1.32.0 scales_1.2.1
[103] png_0.1-8 knitr_1.43 rstudioapi_0.15.0
[106] tzdb_0.4.0 reshape2_1.4.4 rjson_0.2.21
[109] checkmate_2.2.0 nlme_3.1-162 curl_5.0.1
[112] zoo_1.8-12 cachem_1.0.8 stringr_1.5.0
[115] KernSmooth_2.23-21 parallel_4.3.1 miniUI_0.1.1.1
[118] foreign_0.8-84 restfulr_0.0.15 pillar_1.9.0
[121] grid_4.3.1 vctrs_0.6.3 RANN_2.6.1
[124] promises_1.2.0.1 dbplyr_2.3.3 xtable_1.8-4
[127] cluster_2.1.4 htmlTable_2.4.1 evaluate_0.21
[130] cli_3.6.1 compiler_4.3.1 Rsamtools_2.16.0
[133] rlang_1.1.1 crayon_1.5.2 future.apply_1.11.0
[136] plyr_1.8.8 stringi_1.7.12 viridisLite_0.4.2
[139] deldir_1.0-9 BiocParallel_1.34.2 munsell_0.5.0
[142] Biostrings_2.68.1 lazyeval_0.2.2 spatstat.geom_3.2-1
[145] Matrix_1.6-0 BSgenome_1.68.0 hms_1.1.3
[148] bit64_4.0.5 future_1.33.0 KEGGREST_1.40.0
[151] shiny_1.7.5 SummarizedExperiment_1.30.2 ROCR_1.0-11
[154] igraph_1.5.0 memoise_2.0.1 fastmatch_1.1-4
[157] bit_4.0.5

[Jyng321]

Hi guys,
I am not a bioinformatic, but I tried to figure out the problem by looking for the raw code of ClosestFeature in signac/R/utilities.R. The problem is distanceToNearest function doesn't add NA in the list. So I just exported what are actually in the output of distanceToNearest. My manual code is:

Load data

object <-obj[["ATAC"]]
regions <- top.da.peak.all #row.name from FindAllMarkers
annotations <- GetGRangesFromEnsDb(ensdb = EnsDb.Hsapiens.v86)
seqlevelsStyle(annotations) <- 'UCSC'
genome(annotations) <- "hg38"

regions <- StringToGRanges(regions = regions,sep = c("-", "-"))
if (!is(object = regions, class2 = 'GRanges')) {
regions <- StringToGRanges(regions = regions) }

Remove the NA

annotation <- SetIfNull(x = annotations, y = Annotation(object = object))
missing_seqlevels <- setdiff(
x = seqlevels(x = regions), y = seqlevels(x = annotation))
regions <- dropSeqlevels(
x = regions,
value = missing_seqlevels,
pruning.mode = "coarse"
)

After nearest_feature is calculated

nearest_feature <- distanceToNearest(x = regions, subject = annotation)

Ensure that only matched regions are retained

feature_hits <- annotation[subjectHits(x = nearest_feature)]
filtered_regions <- regions[queryHits(x = nearest_feature)]

Convert to data frames with aligned sizes

df <- as.data.frame(x = mcols(x = feature_hits))
df$closest_region <- GRangesToString(grange = feature_hits)
df$query_region <- GRangesToString(grange = filtered_regions)
df$distance <- mcols(x = nearest_feature)$distance

Save

write.csv(df, file = here("Results_Analysis/Closest_Features_Top_Da_peak.csv"), quote = F, row.names = F)
print(df)
I hope it can help.
Best.