TSS shift calculation sets shift_score_FDR=1

[krbe9202]

Hi,

I've been trying to run the example code from your FEATURE_SHIFT.pdf in the documentation, and this error appears when running the section with the plot_shift_rank function:

Error in `combine_vars()`:
! Faceting variables must have at least one value
Run `rlang::last_error()` to see where the error occurred.

traceback():

13: stop(fallback)
12: signal_abort(cnd, .file)
11: rlang::abort(message, ..., call = call, use_cli_format = TRUE, 
        .frame = .frame)
10: cli::cli_abort("Faceting variables must have at least one value")
9: combine_vars(data, params$plot_env, vars, drop = params$drop)
8: compute_layout(..., self = self)
7: self$facet$compute_layout(data, self$facet_params)
6: setup(..., self = self)
5: layout$setup(data, plot$data, plot$plot_env)
4: ggplot_build.ggplot(x)
3: ggplot_build(x)
2: print.ggplot(x)
1: (function (x, ...) 
   UseMethod("print"))(x)

From my understanding, this happens because the shift_score_FDR values are always set to 1, meaning everything gets annotated as 'n.s' when using fdr_cutoff= 0.05, which in turn is filtered out by the plotting function leaving nothing to plot.

My question is, why does the tss_shift() function only report shift_score_FDR =1, and not actual FDR values? This is not consistent with the example plot you have the in documentation, even though I'm apparently running the exact same code on the same data.

I do also get the same issue when running the same code on our own data.

This is the example code I'm trying to run using your toy data:

library("TSRexploreR")
library("ggplot2")
# Load example TSSs
data(TSSs)

# Load genome assembly and annotation
annotation <- system.file("extdata", "S288C_Annotation.gtf", package = "TSRexploreR")
assembly <- system.file("extdata", "S288C_Assembly.fasta", package = "TSRexploreR")

# Generate sample sheet
samples <- data.frame(
  sample_name=c(sprintf("S288C_D_%s", seq_len(3)), sprintf("S288C_WT_%s", seq_len(3))),
  file_1=rep(NA, 6), file_2=rep(NA, 6),
  condition=c(rep("Diamide", 3), rep("Untreated", 3))
)

# Create the TSRexploreR object
exp_test <- TSSs %>%
  tsr_explorer(
    genome_annotation=annotation, genome_assembly=assembly,
    sample_sheet=samples
  ) %>%
  format_counts(data_type = "tss")

# Cluster the TSSs into TSR
exp_test <- tss_clustering(exp_test, threshold=3, n_samples=1) %>%
  merge_samples(data_type = "tss", merge_group="condition") %>%
  merge_samples(data_type = "tsr", merge_group="condition")

# Calculate shift for TSR
exp_test <- tss_shift(
  exp_test,
  sample_1=c(TSS="Untreated", TSR="Untreated"),
  sample_2=c(TSS="Diamide", TSR="Diamide"),
  comparison_name="Untreated_vs_Diamide"
)

plot_shift_rank(exp_test) +
  scale_fill_viridis_c() +
  theme_bw() +
  geom_hline(yintercept = 0)

And my sessionInfo():

R version 4.2.1 (2022-06-23)
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

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

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

other attached packages:
[1] ggplot2_3.4.0   TSRexploreR_0.1

loaded via a namespace (and not attached):
 [1] viridis_0.6.2               MatrixGenerics_1.10.0       Biobase_2.58.0              httr_1.4.4                  tidyr_1.2.1                
 [6] viridisLite_0.4.1           bit64_4.0.5                 assertthat_0.2.1            stats4_4.2.0                BiocFileCache_2.6.0        
[11] blob_1.2.3                  GenomeInfoDbData_1.2.9      Rsamtools_2.14.0            yaml_2.3.6                  progress_1.2.2             
[16] pillar_1.8.1                RSQLite_2.2.19              lattice_0.20-45             glue_1.6.2                  digest_0.6.30              
[21] GenomicRanges_1.50.1        XVector_0.38.0              colorspace_2.0-3            Matrix_1.5-3                XML_3.99-0.13              
[26] pkgconfig_2.0.3             biomaRt_2.54.0              zlibbioc_1.44.0             purrr_0.3.5                 scales_1.2.1               
[31] BiocParallel_1.32.4         tibble_3.1.8                KEGGREST_1.38.0             generics_0.1.3              IRanges_2.32.0             
[36] ellipsis_0.3.2              withr_2.5.0                 cachem_1.0.6                SummarizedExperiment_1.28.0 GenomicFeatures_1.50.2     
[41] BiocGenerics_0.44.0         cli_3.4.1                   magrittr_2.0.3              crayon_1.5.2                memoise_2.0.1              
[46] fansi_1.0.3                 forcats_0.5.2               xml2_1.3.3                  tools_4.2.0                 data.table_1.14.6          
[51] prettyunits_1.1.1           hms_1.1.2                   GlobalOptions_0.1.2         BiocIO_1.8.0                lifecycle_1.0.3            
[56] matrixStats_0.63.0          stringr_1.5.0               plyranges_1.16.0            S4Vectors_0.36.1            munsell_0.5.0              
[61] DelayedArray_0.24.0         AnnotationDbi_1.60.0        Biostrings_2.66.0           compiler_4.2.0              GenomeInfoDb_1.34.4        
[66] rlang_1.0.6                 grid_4.2.0                  RCurl_1.98-1.9              rstudioapi_0.14             rjson_0.2.21               
[71] rappdirs_0.3.3              circlize_0.4.15             bitops_1.0-7                restfulr_0.0.15             gtable_0.3.1               
[76] codetools_0.2-18            DBI_1.1.3                   curl_4.3.3                  R6_2.5.1                    gridExtra_2.3              
[81] GenomicAlignments_1.34.0    dplyr_1.0.10                rtracklayer_1.58.0          fastmap_1.1.0               bit_4.0.5                  
[86] utf8_1.2.2                  filelock_1.0.2              shape_1.4.6                 stringi_1.7.8               parallel_4.2.0             
[91] Rcpp_1.0.9                  vctrs_0.5.1                 png_0.1-8                   dbplyr_2.2.1                tidyselect_1.2.0