Rsamtools::indexTabix: Capture console output

[bschilder]

I recently posted an error I've been encountering with seqminer (here) but realized it seems to trace all the way back to tabix itself.

Regardless, in cases where I get the following kind of error...

[E::hts_idx_push] Unsorted positions on sequence #2: 71999802 followed by 7

...I'd like to wrap the function in atryCatch and capture the messages so that I rerun it with the offending lines in the skip argument. Something like:

out <- tryCatch({  
               Rsamtools::indexTabix(file = bgz_file,
                                     seq = 2,
                                     start = 3,
                                     end = 3, 
                                     comment ="SNP"
               ) 
           })

However, I can't seem to get the output into R. I've tried:

• tryCatch (above)
• capture.output
• sink()

Do you know a way to extract this info?

Thanks,
Brian

Session info

R version 4.1.0 (2021-05-18)
Platform: x86_64-apple-darwin17.0 (64-bit)
Running under: macOS Big Sur 11.4

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

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

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

other attached packages:
[1] arrow_7.0.0      ggimage_0.3.0    ggplot2_3.3.5    dplyr_1.0.8      hexSticker_0.4.9
[6] echotabix_0.99.3

loaded via a namespace (and not attached):
  [1] AnnotationHub_3.2.2           BiocFileCache_2.2.1           systemfonts_1.0.4            
  [4] igraph_1.2.11                 BiocParallel_1.28.3           GenomeInfoDb_1.30.1          
  [7] digest_0.6.29                 yulab.utils_0.0.4             htmltools_0.5.2              
 [10] magick_2.7.3                  fansi_1.0.2                   magrittr_2.0.2               
 [13] memoise_2.0.1                 BSgenome_1.62.0               echoverseTemplate_0.99.0     
 [16] ontologyPlot_1.6              openxlsx_4.2.5                Biostrings_2.62.0            
 [19] matrixStats_0.61.0            R.utils_2.11.0                sysfonts_0.8.5               
 [22] prettyunits_1.1.1             colorspace_2.0-3              blob_1.2.2                   
 [25] rappdirs_0.3.3                textshaping_0.3.6             xfun_0.30                    
 [28] crayon_1.5.0                  RCurl_1.98-1.6                echodata_0.99.6              
 [31] jsonlite_1.8.0                hexbin_1.28.2                 graph_1.72.0                 
 [34] VariantAnnotation_1.40.0      glue_1.6.2                    gtable_0.3.0                 
 [37] zlibbioc_1.40.0               XVector_0.34.0                DelayedArray_0.20.0          
 [40] Rgraphviz_2.38.0              BiocGenerics_0.40.0           scales_1.1.1                 
 [43] DBI_1.1.2                     Rcpp_1.0.8.2                  showtextdb_3.0               
 [46] xtable_1.8-4                  progress_1.2.2                gridGraphics_0.5-1           
 [49] bit_4.0.4                     clisymbols_1.2.0              stats4_4.1.0                 
 [52] DT_0.21                       htmlwidgets_1.5.4             httr_1.4.2                   
 [55] ontologyIndex_2.7             ellipsis_0.3.2                pkgconfig_2.0.3              
 [58] XML_3.99-0.9                  R.methodsS3_1.8.1             farver_2.1.0                 
 [61] seqminer_8.4                  dbplyr_2.1.1                  utf8_1.2.2                   
 [64] here_1.0.1                    ggplotify_0.1.0               tidyselect_1.1.2             
 [67] labeling_0.4.2                rlang_1.0.2                   later_1.3.0                  
 [70] AnnotationDbi_1.56.2          BiocVersion_3.14.0            munsell_0.5.0                
 [73] tools_4.1.0                   cachem_1.0.6                  cli_3.2.0                    
 [76] generics_0.1.2                RSQLite_2.2.10                evaluate_0.15                
 [79] stringr_1.4.0                 fastmap_1.1.0                 yaml_2.3.5                   
 [82] ragg_1.2.2                    knitr_1.37                    bit64_4.0.5                  
 [85] fs_1.5.2                      zip_2.2.0                     purrr_0.3.4                  
 [88] KEGGREST_1.34.0               gh_1.3.0                      showtext_0.9-5               
 [91] mime_0.12                     R.oo_1.24.0                   xml2_1.3.3                   
 [94] biomaRt_2.50.3                brio_1.1.3                    compiler_4.1.0               
 [97] rstudioapi_0.13               interactiveDisplayBase_1.32.0 filelock_1.0.2               
[100] curl_4.3.2                    png_0.1-7                     testthat_3.1.2               
[103] paintmap_1.0                  tibble_3.1.6                  stringi_1.7.6                
[106] GenomicFeatures_1.46.5        desc_1.4.1                    lattice_0.20-45              
[109] Matrix_1.4-0                  vctrs_0.3.8                   pillar_1.7.0                 
[112] lifecycle_1.0.1               BiocManager_1.30.16           data.table_1.14.2            
[115] bitops_1.0-7                  httpuv_1.6.5                  rtracklayer_1.54.0           
[118] GenomicRanges_1.46.1          R6_2.5.1                      BiocIO_1.4.0                 
[121] promises_1.2.0.1              IRanges_2.28.0                ontoProc_1.16.0              
[124] assertthat_0.2.1              pkgload_1.2.4                 SummarizedExperiment_1.24.0  
[127] rprojroot_2.0.2               rjson_0.2.21                  withr_2.5.0                  
[130] GenomicAlignments_1.30.0      Rsamtools_2.10.0              S4Vectors_0.32.3             
[133] GenomeInfoDbData_1.2.7        parallel_4.1.0                hms_1.1.1                    
[136] grid_4.1.0                    ggfun_0.0.5                   tidyr_1.2.0                  
[139] rmarkdown_2.13                MatrixGenerics_1.6.0          piggyback_0.1.1              
[142] Biobase_2.54.0                shiny_1.7.1                   lubridate_1.8.0              
[145] restfulr_0.0.13          

[mtmorgan]

These errors are written to C language 'stderr' but unfortunately this is not directly available within R. From the command line, on Linux, one could run a script script.R and redirect stderr to a file

Rscript -f script.R 2> stderr.log

but probably your intention is something different.

Previously, Rsamtools used a 'hack' to remap fprintf(stderr, ...) to R's error stream, with something like

DFLAGS = \
  -Dfprintf=_samtools_fprintf \
  -Dexit=_samtools_exit \
  -Dabort=_samtools_abort

This would have captured samtools errors as R warnings, but unfortunately this functionality was lost when we shifted to using Rhtslib quite some time ago; perhaps @hpages has additional comment.

I don't think there's an easy way to capture what you'd like to do.

[hpages]

Hi Martin,

HTSlib uses a more sophisticated mechanism than Samtools for logging events: now everything goes thru the hts_log() function defined in header file hts_log.h. The function handles various levels of event severity: error, warning, debug, and trace. There's a macro defined for each level (hts_log_error(), hts_log_warning(), etc...). Each macro is a simple wrapper around hts_log().

Note that hts_log() itself is not a macro, it's a real function, and it's not an inline function either, so I'm not sure why it's defined in a header file.

Anyways, this new mechanism and its implementation differ substantially from the plain fprintf(stderr, ...) calls used all over the place in Samtools. I didn't investigate so I don't know how easy/hard it would be to come up with a hack similar to yours where all the events logged thru hts_log_error() get captured as R warnings.

Ideas, suggestions, patches are welcome.

Thanks,
H.

[hpages]

Personally I'm not planning to tweak Rhtslib to make the low-level C code in htslib spit output that is capturable from R but maybe someone else wants to give it a shot? If not, we can close this.

Thanks,
H.