search

js2264 / HiContacts

Analyzing Hi-C data in R with HiCExperiment objects
https://js2264.github.io/OHCA/
Other
10 stars 1 forks source link

v4C - rowSum ERROR #4

Closed matteozoia4 closed 1 year ago

matteozoia4 commented 1 year ago

Dear all,

I am trying to extract the v4C from a .cool file. It works with another view point, but with this one "chr8:57483000-57493000", it gives an error, why?

v4C_Sap30_Prom_10kb <- virtual4C(hic, viewpoint = GRanges("chr8:57483000-57493000")) _Error in rowSums(cm[, regions_inviewpoint], na.rm = TRUE) : 'x' must be an array of at least two dimensions

How can I fix this?

Best regards and many thanks!

Link to the .cool file : https://we.tl/t-oi2vbQR8SW

js2264 commented 1 year ago

I believe this error might be due to the fact that your query GRanges is rather narrow. If this is the issue, this bug could be fixed by a recent update in the package (https://github.com/js2264/HiContacts/commit/9c76595ca4daab4e22eee6f578f5f37b76fd405e). You can update HiContacts to the development version from GitHub using the following command:

remotes::install_github(js2264/HiContacts)

I do not have access to a computer right now, but I will check this as soon as I have a machine at my disposal.

matteozoia4 commented 1 year ago

Many thanks for the prompt response!

The error persists unfortunately.

js2264 commented 1 year ago

I am able to compute a v4C profile from the file and command you provided me: 

> hic <- import('GSE161259_ES_chr8_55400000_59200000_MAPQ30.40kb.cool')
> v4C <- virtual4C(hic, viewpoint = GenomicRanges::GRanges("chr8:57483000-57493000"))
Warning message:
In asMethod(object) :
  sparse->dense coercion: allocating vector of size 34.8 GiB

> v4C
GRanges object with 68301 ranges and 4 metadata columns:
                      seqnames        ranges strand |     score              viewpoint    center in_viewpoint
                         <Rle>     <IRanges>  <Rle> | <numeric>            <character> <numeric>    <logical>
      [1]                 chr1       1-40000      * |         0 chr8:57483000-57493000   20000.5        FALSE
      [2]                 chr1   40001-80000      * |         0 chr8:57483000-57493000   60000.5        FALSE
      [3]                 chr1  80001-120000      * |         0 chr8:57483000-57493000  100000.5        FALSE
      [4]                 chr1 120001-160000      * |         0 chr8:57483000-57493000  140000.5        FALSE
      [5]                 chr1 160001-200000      * |         0 chr8:57483000-57493000  180000.5        FALSE
      ...                  ...           ...    ... .       ...                    ...       ...          ...
  [68297]       chrUn_GL456396       1-21240      * |         0 chr8:57483000-57493000   10620.5        FALSE
  [68298]       chrUn_GL456368       1-20208      * |         0 chr8:57483000-57493000   10104.5        FALSE
  [68299]                 chrM       1-16299      * |         0 chr8:57483000-57493000    8150.0        FALSE
  [68300] chr4_JH584292_random       1-14945      * |         0 chr8:57483000-57493000    7473.0        FALSE
  [68301] chr4_JH584295_random        1-1976      * |         0 chr8:57483000-57493000     988.5        FALSE
  -------
  seqinfo: 66 sequences from an unspecified genome; no seqlengths

My sessionInfo(): 

R version 4.3.0 (2023-04-21)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Ubuntu 18.04.6 LTS

Matrix products: default
BLAS/LAPACK: /usr/lib/x86_64-linux-gnu/libopenblasp-r0.2.20.so;  LAPACK version 3.7.1

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

time zone: Europe/Paris
tzcode source: system (glibc)

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

other attached packages:
[1] HiContacts_1.3.0    HiCExperiment_1.1.0

loaded via a namespace (and not attached):
 [1] SummarizedExperiment_1.30.2 gtable_0.3.3                beeswarm_0.4.0
 [4] ggplot2_3.4.2               rhdf5_2.44.0                Biobase_2.60.0
 [7] lattice_0.21-8              tzdb_0.4.0                  rhdf5filters_1.12.1
[10] vctrs_0.6.2                 tools_4.3.0                 bitops_1.0-7
[13] generics_0.1.3              stats4_4.3.0                parallel_4.3.0
[16] tibble_3.2.1                fansi_1.0.4                 pkgconfig_2.0.3
[19] Matrix_1.5-4.1              S4Vectors_0.38.1            lifecycle_1.0.3
[22] GenomeInfoDbData_1.2.10     compiler_4.3.0              stringr_1.5.0
[25] munsell_0.5.0               codetools_0.2-19            InteractionSet_1.29.1
[28] vipor_0.4.5                 GenomeInfoDb_1.36.0         RCurl_1.98-1.12
[31] pillar_1.9.0                crayon_1.5.2                tidyr_1.3.0
[34] BiocParallel_1.34.2         DelayedArray_0.26.3         RSpectra_0.16-1
[37] tidyselect_1.2.0            stringi_1.7.12              dplyr_1.1.2
[40] purrr_1.0.1                 grid_4.3.0                  colorspace_2.1-0
[43] cli_3.6.1                   magrittr_2.0.3              S4Arrays_1.0.4
[46] utf8_1.2.3                  withr_2.5.0                 readr_2.1.4
[49] scales_1.2.1                bit64_4.0.5                 ggbeeswarm_0.7.2
[52] XVector_0.40.0              matrixStats_1.0.0           bit_4.0.5
[55] hms_1.1.3                   ggrastr_1.0.2               GenomicRanges_1.52.0
[58] IRanges_2.34.0              BiocIO_1.10.0               rlang_1.1.1
[61] Rcpp_1.0.10                 glue_1.6.2                  BiocGenerics_0.46.0
[64] vroom_1.6.3                 strawr_0.0.91               R6_2.5.1
[67] Rhdf5lib_1.22.0             MatrixGenerics_1.12.2       zlibbioc_1.46.0

Could you please share your sessionInfo()?

Also, I filled an issue to InteractionSet, as a hot-fix had not been included in the latest release version (1.28.0). You may have to install 1.29.1 from github: remotes::install_github("LTLA/InteractionSet") to sort all the issues (though the one you reported should be fixed in

matteozoia4 commented 1 year ago

Many thank,

I installed the InteractionSet tool many thanks!

Here my session:

R version 4.2.1 (2022-06-23)
Platform: x86_64-apple-darwin17.0 (64-bit)
Running under: macOS Monterey 12.0.1

Matrix products: default
LAPACK: /Library/Frameworks/R.framework/Versions/4.2/Resources/lib/libRlapack.dylib

locale:
[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8

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

other attached packages:
 [1] SummarizedExperiment_1.28.0 Biobase_2.58.0              MatrixGenerics_1.10.0      
 [4] matrixStats_1.0.0           remotes_2.4.2               devtools_2.4.5             
 [7] usethis_2.2.0               HiCExperiment_1.1.0         rtracklayer_1.58.0         
[10] ggplot2_3.4.2               HiContactsData_1.0.0        ExperimentHub_2.6.0        
[13] AnnotationHub_3.6.0         BiocFileCache_2.6.1         dbplyr_2.3.2               
[16] GenomicRanges_1.50.2        GenomeInfoDb_1.34.9         IRanges_2.32.0             
[19] S4Vectors_0.36.2            BiocGenerics_0.44.0         HiContacts_1.3.0           

loaded via a namespace (and not attached):
  [1] backports_1.4.1               Hmisc_5.1-0                  
  [3] igraph_1.4.3                  lazyeval_0.2.2               
  [5] BiocParallel_1.32.6           digest_0.6.31                
  [7] ensembldb_2.22.0              htmltools_0.5.5              
  [9] fansi_1.0.4                   magrittr_2.0.3               
 [11] checkmate_2.2.0               memoise_2.0.1                
 [13] BSgenome_1.66.3               cluster_2.1.4                
 [15] tzdb_0.4.0                    InteractionSet_1.29.1        
 [17] Biostrings_2.66.0             vroom_1.6.3                  
 [19] prettyunits_1.1.1             jpeg_0.1-10                  
 [21] colorspace_2.1-0              blob_1.2.4                   
 [23] rappdirs_0.3.3                xfun_0.39                    
 [25] dplyr_1.1.2                   jsonlite_1.8.5               
 [27] callr_3.7.3                   crayon_1.5.2                 
 [29] RCurl_1.98-1.12               VariantAnnotation_1.44.1     
 [31] glue_1.6.2                    gtable_0.3.3                 
 [33] zlibbioc_1.44.0               XVector_0.38.0               
 [35] strawr_0.0.91                 DelayedArray_0.24.0          
 [37] pkgbuild_1.4.0                Rhdf5lib_1.20.0              
 [39] scales_1.2.1                  DBI_1.1.3                    
 [41] miniUI_0.1.1.1                Rcpp_1.0.10                  
 [43] xtable_1.8-4                  progress_1.2.2               
 [45] htmlTable_2.4.1               reticulate_1.30              
 [47] foreign_0.8-84                bit_4.0.5                    
 [49] Formula_1.2-5                 profvis_0.3.8                
 [51] htmlwidgets_1.6.2             httr_1.4.6                   
 [53] RColorBrewer_1.1-3            ellipsis_0.3.2               
 [55] urlchecker_1.0.1              pkgconfig_2.0.3              
 [57] XML_3.99-0.14                 Gviz_1.42.1                  
 [59] nnet_7.3-19                   deldir_1.0-9                 
 [61] utf8_1.2.3                    tidyselect_1.2.0             
 [63] rlang_1.1.1                   later_1.3.1                  
 [65] AnnotationDbi_1.60.2          BiocVersion_3.16.0           
 [67] munsell_0.5.0                 tools_4.2.1                  
 [69] cachem_1.0.8                  cli_3.6.1                    
 [71] generics_0.1.3                RSQLite_2.3.1                
 [73] evaluate_0.21                 stringr_1.5.0                
 [75] fastmap_1.1.1                 yaml_2.3.7                   
 [77] processx_3.8.1                knitr_1.43                   
 [79] bit64_4.0.5                   fs_1.6.2                     
 [81] purrr_1.0.1                   KEGGREST_1.38.0              
 [83] AnnotationFilter_1.22.0       mime_0.12                    
 [85] ggrastr_1.0.2                 xml2_1.3.4                   
 [87] biomaRt_2.54.1                compiler_4.2.1               
 [89] rstudioapi_0.14               beeswarm_0.4.0               
 [91] interactiveDisplayBase_1.36.0 filelock_1.0.2               
 [93] curl_5.0.1                    png_0.1-8                    
 [95] tibble_3.2.1                  stringi_1.7.12               
 [97] ps_1.7.5                      desc_1.4.2                   
 [99] GenomicFeatures_1.50.4        lattice_0.21-8               
[101] ProtGenerics_1.30.0           Matrix_1.5-4.1               
[103] vctrs_0.6.2                   rhdf5filters_1.10.1          
[105] pillar_1.9.0                  GenomicInteractions_1.32.0   
[107] lifecycle_1.0.3               BiocManager_1.30.21          
[109] data.table_1.14.8             bitops_1.0-7                 
[111] httpuv_1.6.11                 R6_2.5.1                     
[113] BiocIO_1.8.0                  latticeExtra_0.6-30          
[115] promises_1.2.0.1              gridExtra_2.3                
[117] vipor_0.4.5                   sessioninfo_1.2.2            
[119] codetools_0.2-19              dichromat_2.0-0.1            
[121] pkgload_1.3.2                 rhdf5_2.42.1                 
[123] rprojroot_2.0.3               rjson_0.2.21                 
[125] withr_2.5.0                   GenomicAlignments_1.34.1     
[127] Rsamtools_2.14.0              GenomeInfoDbData_1.2.9       
[129] parallel_4.2.1                hms_1.1.3                    
[131] grid_4.2.1                    rpart_4.1.19                 
[133] tidyr_1.3.0                   rmarkdown_2.22               
[135] biovizBase_1.46.0             shiny_1.7.4                  
[137] base64enc_0.1-3               ggbeeswarm_0.7.2             
[139] interp_1.1-4                  restfulr_0.0.15
js2264 commented 1 year ago

This is odd. For troubleshooting, could you start a new R instance, paste the following commands and send me the output? 

library(HiCExperiment)
library(GenomicRanges)
hic <- import('GSE161259_ES_chr8_55400000_59200000_MAPQ30.40kb.cool')

##

print(HiContacts::virtual4C)
HiContacts::virtual4C(hic, GRanges("chr8:57483000-57493000"))

##

gis <- interactions(hic)
cm <- Matrix::as.matrix(gi2cm(gis, 'balanced'))
cm <- base::as.matrix(cm)
regions <- regions(gis)
regions_in_viewpoint <- seq_along(regions) %in% queryHits(
    findOverlaps(regions, GRanges("chr8:57483000-57493000"))
)
score <- cm[, regions_in_viewpoint]

##

table(regions_in_viewpoint)
head(score)

Thanks