This report reads the outputs produced by
scripts/test_decoupled_pipeline_star_fullsize.sh (or the
BWA variant) and compares the joint -W 1 pipeline’s BSJ/FSJ
matrices against the decoupled
SCAN1 → BUILD_UNIVERSE → SCAN2 → FINALIZE matrices.
To rerender after a new test run:
Rscript -e "rmarkdown::render('scripts/validation_report.Rmd')"
or, to point at a different output directory:
Rscript -e "rmarkdown::render('scripts/validation_report.Rmd', params=list(out_root='/path/to/out_root'))"
Test configuration
out_root <- params$out_root
orig_dir <- file.path(out_root, "original")
decoupled_dir <- file.path(out_root, "decoupled")
finalize_dir <- file.path(decoupled_dir, "finalize")
bench_dir <- file.path(out_root, "bench")
required <- c(
"original/result.BSJ_Matrix",
"original/result.FSJ_Matrix",
"decoupled/finalize/result.BSJ_Matrix",
"decoupled/finalize/result.FSJ_Matrix"
)
missing <- required[!file.exists(file.path(out_root, required))]
if (length(missing) > 0) {
stop("Missing expected input files:\n ",
paste(file.path(out_root, missing), collapse = "\n "))
}
tibble(key = c("out_root", "orig_dir", "decoupled_dir"),
value = c(out_root, orig_dir, decoupled_dir)) |>
kable(caption = "Paths used by this report")
Paths used by this report
| out_root |
/pastel/tools/circRNA_tools/test_data/decoupled_comparison |
| orig_dir |
/pastel/tools/circRNA_tools/test_data/decoupled_comparison/original |
| decoupled_dir |
/pastel/tools/circRNA_tools/test_data/decoupled_comparison/decoupled |
Load matrices
joint_BSJ <- read_tsv(file.path(orig_dir, "result.BSJ_Matrix"))
joint_FSJ <- read_tsv(file.path(orig_dir, "result.FSJ_Matrix"))
decoupled_BSJ <- read_tsv(file.path(finalize_dir, "result.BSJ_Matrix"))
decoupled_FSJ <- read_tsv(file.path(finalize_dir, "result.FSJ_Matrix"))
# Joint writes sample columns by input-file basename (e.g. "bwa.sam");
# decoupled uses the sampleName column from finalize_samples.tsv. They mean
# the same thing in the same order, so align positionally using the decoupled
# names — those are more informative.
align_cols <- function(df, ref_names) {
stopifnot(ncol(df) == length(ref_names) + 1)
colnames(df) <- c("circRNA_ID", ref_names)
df
}
ref_names <- colnames(decoupled_BSJ)[-1]
joint_BSJ <- align_cols(joint_BSJ, ref_names)
joint_FSJ <- align_cols(joint_FSJ, ref_names)
tibble(
matrix = c("BSJ joint", "BSJ decoupled", "FSJ joint", "FSJ decoupled"),
nrow = c(nrow(joint_BSJ), nrow(decoupled_BSJ),
nrow(joint_FSJ), nrow(decoupled_FSJ)),
n_samples = length(ref_names)
) |> kable(caption = "Matrix dimensions (sample columns aligned positionally)")
Matrix dimensions (sample columns aligned
positionally)
| BSJ joint |
35869 |
4 |
| BSJ decoupled |
35869 |
4 |
| FSJ joint |
35869 |
4 |
| FSJ decoupled |
35869 |
4 |
cat("Sample column names used:\n")
## Sample column names used:
cat(paste0(" ", ref_names, collapse = "\n"), "\n")
## Div_100_S91
## Div_101_S92
## PARDOS_1_S1
## PARDOS_2_S2
circRNA set agreement
both <- intersect(joint_BSJ$circRNA_ID, decoupled_BSJ$circRNA_ID)
only_j <- setdiff(joint_BSJ$circRNA_ID, decoupled_BSJ$circRNA_ID)
only_d <- setdiff(decoupled_BSJ$circRNA_ID, joint_BSJ$circRNA_ID)
tibble(
set = c("circRNAs in both pipelines",
"joint only",
"decoupled only",
"Jaccard (intersection / union)"),
n = c(length(both),
length(only_j),
length(only_d),
round(length(both) /
length(union(joint_BSJ$circRNA_ID,
decoupled_BSJ$circRNA_ID)), 4))
) |> kable(caption = "circRNA ID set agreement")
circRNA ID set agreement
| circRNAs in both pipelines |
35869 |
| joint only |
0 |
| decoupled only |
0 |
| Jaccard (intersection / union) |
1 |
Sample joint-only IDs (first 20)
head(only_j, 20) |> as_tibble() |> kable(col.names = "joint-only circRNA_ID")
Sample decoupled-only IDs (first 20)
head(only_d, 20) |> as_tibble() |> kable(col.names = "decoupled-only circRNA_ID")
BSJ matrix agreement
bsj_merged <- full_join(joint_BSJ, decoupled_BSJ,
by = "circRNA_ID",
suffix = c(".joint", ".decoupled"))
bsj_long <- bsj_merged |>
pivot_longer(cols = -circRNA_ID,
names_to = c("Sample", "Run"),
names_pattern = "(.*)\\.(joint|decoupled)$") |>
pivot_wider(names_from = Run, values_from = value) |>
replace_na(list(joint = 0, decoupled = 0))
bsj_long |>
group_by(Sample) |>
summarise(
n_rows = n(),
exactly_equal = sum(joint == decoupled),
joint_greater = sum(joint > decoupled),
decoupled_greater = sum(decoupled > joint),
max_abs_diff = max(abs(joint - decoupled)),
pct_exact = round(100 * mean(joint == decoupled), 2)
) |> kable(caption = "BSJ matrix: per-sample cell-level agreement")
BSJ matrix: per-sample cell-level agreement
| Div_100_S91 |
35869 |
35869 |
0 |
0 |
0 |
100 |
| Div_101_S92 |
35869 |
35869 |
0 |
0 |
0 |
100 |
| PARDOS_1_S1 |
35869 |
35869 |
0 |
0 |
0 |
100 |
| PARDOS_2_S2 |
35869 |
35869 |
0 |
0 |
0 |
100 |
BSJ scatter: joint vs. decoupled
ggplot(bsj_long, aes(x = joint, y = decoupled)) +
geom_point(alpha = 0.4, size = 0.6) +
geom_abline(slope = 1, intercept = 0, linetype = "dashed", color = "red") +
facet_wrap(~ Sample, scales = "free") +
labs(x = "Joint pipeline BSJ count",
y = "Decoupled pipeline BSJ count",
title = "BSJ counts per circRNA per sample") +
theme_bw()
BSJ cell-level discordance (top 20)
bsj_long |>
mutate(abs_diff = abs(joint - decoupled)) |>
filter(abs_diff > 0) |>
arrange(desc(abs_diff)) |>
head(20) |> kable(caption = "Largest per-cell BSJ disagreements")
Largest per-cell BSJ disagreements
FSJ matrix agreement
fsj_merged <- full_join(joint_FSJ, decoupled_FSJ,
by = "circRNA_ID",
suffix = c(".joint", ".decoupled"))
fsj_long <- fsj_merged |>
pivot_longer(cols = -circRNA_ID,
names_to = c("Sample", "Run"),
names_pattern = "(.*)\\.(joint|decoupled)$") |>
pivot_wider(names_from = Run, values_from = value) |>
replace_na(list(joint = 0, decoupled = 0))
fsj_long |>
group_by(Sample) |>
summarise(
n_rows = n(),
exactly_equal = sum(joint == decoupled),
joint_greater = sum(joint > decoupled),
decoupled_greater = sum(decoupled > joint),
max_abs_diff = max(abs(joint - decoupled)),
pct_exact = round(100 * mean(joint == decoupled), 2)
) |> kable(caption = "FSJ matrix: per-sample cell-level agreement")
FSJ matrix: per-sample cell-level agreement
| Div_100_S91 |
35869 |
35869 |
0 |
0 |
0 |
100 |
| Div_101_S92 |
35869 |
35869 |
0 |
0 |
0 |
100 |
| PARDOS_1_S1 |
35869 |
35869 |
0 |
0 |
0 |
100 |
| PARDOS_2_S2 |
35869 |
35869 |
0 |
0 |
0 |
100 |
FSJ scatter: joint vs. decoupled
ggplot(fsj_long, aes(x = joint, y = decoupled)) +
geom_point(alpha = 0.4, size = 0.6) +
geom_abline(slope = 1, intercept = 0, linetype = "dashed", color = "red") +
facet_wrap(~ Sample, scales = "free") +
labs(x = "Joint pipeline FSJ count",
y = "Decoupled pipeline FSJ count",
title = "FSJ counts per circRNA per sample") +
theme_bw()
FSJ cell-level discordance (top 20)
fsj_long |>
mutate(abs_diff = abs(joint - decoupled)) |>
filter(abs_diff > 0) |>
arrange(desc(abs_diff)) |>
head(20) |> kable(caption = "Largest per-cell FSJ disagreements")
Largest per-cell FSJ disagreements
Detection-correlation per sample
Per-sample Pearson and Spearman correlation between joint and
decoupled counts — a value of 1.00 means per-cell agreement is perfect;
values under ~0.99 flag a real divergence.
sample_corr <- function(df, what) {
df |>
group_by(Sample) |>
summarise(
pearson = round(cor(joint, decoupled, method = "pearson"), 4),
spearman = round(cor(joint, decoupled, method = "spearman"), 4)
) |>
mutate(matrix = what, .before = 1)
}
bind_rows(sample_corr(bsj_long, "BSJ"),
sample_corr(fsj_long, "FSJ")) |>
kable(caption = "Joint vs. decoupled correlation (per sample)")
Joint vs. decoupled correlation (per sample)
| BSJ |
Div_100_S91 |
1 |
1 |
| BSJ |
Div_101_S92 |
1 |
1 |
| BSJ |
PARDOS_1_S1 |
1 |
1 |
| BSJ |
PARDOS_2_S2 |
1 |
1 |
| FSJ |
Div_100_S91 |
1 |
1 |
| FSJ |
Div_101_S92 |
1 |
1 |
| FSJ |
PARDOS_1_S1 |
1 |
1 |
| FSJ |
PARDOS_2_S2 |
1 |
1 |
Benchmark (if present)
bench_files <- list.files(bench_dir, pattern = "\\.bench$", full.names = TRUE)
if (length(bench_files) == 0) {
cat("No benchmark files found in", bench_dir, "\n")
} else {
parse_bench <- function(f) {
txt <- readLines(f, warn = FALSE)
if (length(txt) == 0) return(NULL)
kv <- strsplit(txt[1], " ")[[1]]
as_tibble(setNames(
lapply(kv, function(x) strsplit(x, "=")[[1]][2]),
sapply(kv, function(x) strsplit(x, "=")[[1]][1])
)) |> mutate(stage = sub("\\.bench$", "", basename(f)), .before = 1)
}
bench_df <- bind_rows(lapply(bench_files, parse_bench))
# Best-effort numeric coercion
for (col in c("wall_sec", "user_sec", "sys_sec", "max_rss_kb")) {
if (col %in% names(bench_df)) {
bench_df[[col]] <- suppressWarnings(as.numeric(bench_df[[col]]))
}
}
bench_df <- bench_df |> arrange(stage)
if ("max_rss_kb" %in% names(bench_df)) {
bench_df$max_rss_gb <- round(bench_df$max_rss_kb / 1024 / 1024, 2)
}
bench_df |> kable(caption = "Per-stage benchmark")
}
Per-stage benchmark
| 00_original_joint |
250.16 |
12447.89 |
561.28 |
13782240 |
5200% |
0 |
13.14 |
| 10_scan1_01_Div_100_S91 |
119.94 |
1886.27 |
104.28 |
11848004 |
1659% |
0 |
11.30 |
| 10_scan1_02_Div_101_S92 |
123.21 |
3735.27 |
80.45 |
11722928 |
3096% |
0 |
11.18 |
| 10_scan1_03_PARDOS_1_S1 |
32.97 |
462.75 |
54.18 |
12309132 |
1567% |
0 |
11.74 |
| 10_scan1_04_PARDOS_2_S2 |
34.56 |
499.20 |
78.16 |
11574668 |
1670% |
0 |
11.04 |
| 20_build_universe |
20.15 |
30.43 |
20.60 |
6727544 |
253% |
0 |
6.42 |
| 30_scan2_01_Div_100_S91 |
158.10 |
2813.38 |
241.65 |
20013064 |
1932% |
0 |
19.09 |
| 30_scan2_02_Div_101_S92 |
131.88 |
3113.51 |
264.95 |
22656412 |
2561% |
0 |
21.61 |
| 30_scan2_03_PARDOS_1_S1 |
90.29 |
2580.47 |
185.58 |
12933104 |
3063% |
0 |
12.33 |
| 30_scan2_04_PARDOS_2_S2 |
95.84 |
2065.58 |
200.04 |
21966236 |
2363% |
0 |
20.95 |
| 40_finalize |
28.15 |
48.67 |
33.87 |
9593792 |
293% |
0 |
9.15 |
Session info
sessionInfo()
## R version 4.1.2 (2021-11-01)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: CentOS Stream 8
##
## Matrix products: default
## BLAS/LAPACK: /usr/lib64/libopenblasp-r0.3.15.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] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] knitr_1.51 lubridate_1.9.5 forcats_1.0.1 stringr_1.6.0
## [5] dplyr_1.2.0 purrr_1.2.1 readr_2.2.0 tidyr_1.3.2
## [9] tibble_3.3.1 ggplot2_4.0.2 tidyverse_2.0.0
##
## loaded via a namespace (and not attached):
## [1] bslib_0.10.0 compiler_4.1.2 pillar_1.11.1 RColorBrewer_1.1-3
## [5] jquerylib_0.1.4 tools_4.1.2 bit_4.6.0 digest_0.6.39
## [9] timechange_0.4.0 jsonlite_2.0.0 evaluate_1.0.5 lifecycle_1.0.5
## [13] gtable_0.3.6 pkgconfig_2.0.3 rlang_1.1.7 cli_3.6.5
## [17] rstudioapi_0.18.0 parallel_4.1.2 yaml_2.3.12 xfun_0.57
## [21] fastmap_1.2.0 withr_3.0.2 hms_1.1.4 generics_0.1.4
## [25] sass_0.4.10 vctrs_0.7.2 bit64_4.6.0-1 grid_4.1.2
## [29] tidyselect_1.2.1 glue_1.8.0 R6_2.6.1 otel_0.2.0
## [33] vroom_1.7.1 rmarkdown_2.31 tzdb_0.5.0 farver_2.1.2
## [37] magrittr_2.0.4 scales_1.4.0 htmltools_0.5.9 dichromat_2.0-0.1
## [41] labeling_0.4.3 S7_0.2.1 stringi_1.8.7 cachem_1.1.0
## [45] crayon_1.5.3