Nanopore bulk RNA-Seq and scRNA-Seq data benchmarks (top highly variable transcripts results summary)
Preparing the data
Load the required packages:
library(tidyverse)
library(glue)Global parameters:
csv_dir <- "nanopore_bulk_sc_benchmarks_hvts_data"
top_n_hvts_values <- c(500, 1000, 2000, 4000, 6000, 10000)Prepare plot data:
plot_df_list <- lapply(top_n_hvts_values, function(top_n_hvts) {
plot_df <- read.csv(file.path(
csv_dir, glue("data_{top_n_hvts}_hvts.csv")
), row.names = 1)
plot_df <- as.data.frame(t(plot_df))
plot_df$tool <- rownames(plot_df)
rownames(plot_df) <- NULL
plot_df$top_n_hvts <- top_n_hvts
plot_df <- plot_df %>%
select("tool", "top_n_hvts", "corr_spearman",
"corr_spearman_upper_ci_95", "corr_spearman_lower_ci_95",
"corr_spearman_upper_ci_90", "corr_spearman_lower_ci_90",
"corr_spearman_upper_ci_68", "corr_spearman_lower_ci_68",
"mean_rel_diff",
"mean_rel_diff_upper_ci_95", "mean_rel_diff_lower_ci_95",
"mean_rel_diff_upper_ci_90", "mean_rel_diff_lower_ci_90",
"mean_rel_diff_upper_ci_68", "mean_rel_diff_lower_ci_68")
return(plot_df)
})
plot_df <- do.call(rbind, plot_df_list)
cor_spearman_hvts_fc_df <- plot_df %>%
select("tool", "top_n_hvts", "corr_spearman") %>%
pivot_wider(names_from = tool, values_from = corr_spearman) %>%
mutate(
IsoQuant = Isosceles / IsoQuant,
FLAMES = Isosceles / FLAMES,
Sicelore = Isosceles / Sicelore
) %>%
select(-"Isosceles")
cor_spearman_fc_df <- plot_df %>%
select("tool", "top_n_hvts", "corr_spearman") %>%
pivot_wider(names_from = tool, values_from = corr_spearman) %>%
transmute(
IsoQuant = Isosceles / IsoQuant,
FLAMES = Isosceles / FLAMES,
Sicelore = Isosceles / Sicelore
)
cor_spearman_fc_df <- data.frame(
median_fc = apply(cor_spearman_fc_df, 2, median),
max_fc = apply(cor_spearman_fc_df, 2, max)
)
mean_rel_diff_hvts_fc_df <- plot_df %>%
select("tool", "top_n_hvts", "mean_rel_diff") %>%
pivot_wider(names_from = tool, values_from = mean_rel_diff) %>%
mutate(
IsoQuant = Isosceles / IsoQuant,
FLAMES = Isosceles / FLAMES,
Sicelore = Isosceles / Sicelore
) %>%
select(-"Isosceles")
mean_rel_diff_fc_df <- plot_df %>%
select("tool", "top_n_hvts", "mean_rel_diff") %>%
pivot_wider(names_from = tool, values_from = mean_rel_diff) %>%
transmute(
IsoQuant = Isosceles / IsoQuant,
FLAMES = Isosceles / FLAMES,
Sicelore = Isosceles / Sicelore
)
mean_rel_diff_fc_df <- data.frame(
median_fc = apply(mean_rel_diff_fc_df, 2, median),
max_fc = apply(mean_rel_diff_fc_df, 2, max)
)
test_df_list <- lapply(top_n_hvts_values, function(top_n_hvts) {
test_df <- read.csv(file.path(
csv_dir, glue("test_values_{top_n_hvts}_hvts.csv")
))
test_df$top_n_hvts <- top_n_hvts
return(test_df)
})
test_df <- do.call(rbind, test_df_list)
corr_spearman_isosceles <- filter(test_df, tool == "Isosceles") %>%
pull(corr_spearman)
mean_rel_diff_isosceles <- filter(test_df, tool == "Isosceles") %>%
pull(mean_rel_diff)
wilcox_pvalue <- function(x, y) {
wilcox.test(x, y, paired = TRUE)$p.value
}
test_df <- test_df %>%
filter(tool != "Isosceles") %>%
group_by(tool) %>%
summarise(
corr_spearman_pvalue = wilcox_pvalue(corr_spearman,
corr_spearman_isosceles),
mean_rel_diff_pvalue = wilcox_pvalue(mean_rel_diff,
mean_rel_diff_isosceles)
) %>%
as.data.frame()
rownames(test_df) <- test_df$tool
test_df$tool <- NULL
test_df <- test_df[c("IsoQuant", "FLAMES", "Sicelore"),]
cor_spearman_fc_df$pvalue <- test_df$corr_spearman_pvalue
mean_rel_diff_fc_df$pvalue <- test_df$mean_rel_diff_pvalueResults
Matched - decoy fold-change values (Spearman correlation):
cor_spearman_hvts_fc_dfMatched - decoy fold-change statistics (Spearman correlation):
cor_spearman_fc_dfMatched - decoy fold-change values (mean relative difference):
mean_rel_diff_hvts_fc_dfMatched - decoy fold-change statistics (mean relative difference):
mean_rel_diff_fc_dfMatched - decoy plot (Spearman correlation, standard error interval):
ggplot(plot_df, mapping = aes(x = top_n_hvts)) +
geom_line(mapping = aes(y = corr_spearman,
col = tool),
size = 1) +
geom_ribbon(mapping = aes(ymin = corr_spearman_lower_ci_68,
ymax = corr_spearman_upper_ci_68,
fill = tool),
alpha = 0.1) +
scale_color_manual(values = c(Isosceles = "black",
IsoQuant = "gold",
FLAMES = "tomato",
Sicelore = "royalblue")) +
scale_fill_manual(values = c(Isosceles = "black",
IsoQuant = "gold",
FLAMES = "tomato",
Sicelore = "royalblue")) +
scale_x_continuous(
breaks = c(500, 1000, 2000, 4000, 6000, 10000),
labels = c("500", "1k", "2k", "4k", "6k", "10k")
) +
coord_cartesian(ylim = c(0, NA)) +
labs(
x = "Top N variable transcripts",
y = "Matched - decoy (correlation)"
) +
theme_bw() +
theme(aspect.ratio = 1,
legend.title = element_blank(),
axis.text.x = element_text(size = 7, hjust = 0.85))
Matched - decoy plot (mean relative difference, standard error interval):
ggplot(plot_df, mapping = aes(x = top_n_hvts)) +
geom_line(mapping = aes(y = mean_rel_diff,
col = tool),
size = 1) +
geom_ribbon(mapping = aes(ymin = mean_rel_diff_lower_ci_68,
ymax = mean_rel_diff_upper_ci_68,
fill = tool),
alpha = 0.1) +
scale_color_manual(values = c(Isosceles = "black",
IsoQuant = "gold",
FLAMES = "tomato",
Sicelore = "royalblue")) +
scale_fill_manual(values = c(Isosceles = "black",
IsoQuant = "gold",
FLAMES = "tomato",
Sicelore = "royalblue")) +
scale_x_continuous(
breaks = c(500, 1000, 2000, 4000, 6000, 10000),
labels = c("500", "1k", "2k", "4k", "6k", "10k")
) +
coord_cartesian(ylim = c(0, NA)) +
labs(
x = "Top N variable transcripts",
y = "Matched - decoy (mean rel. diff.)"
) +
theme_bw() +
theme(aspect.ratio = 1,
legend.title = element_blank(),
axis.text.x = element_text(size = 7, hjust = 0.85))
Session Info
sessionInfo()## R version 4.2.1 (2022-06-23)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 20.04.4 LTS
##
## Matrix products: default
## BLAS: /usr/lib/x86_64-linux-gnu/openblas-pthread/libblas.so.3
## LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/liblapack.so.3
##
## locale:
## [1] C
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] glue_1.6.2 forcats_0.5.2 stringr_1.4.1 dplyr_1.0.10
## [5] purrr_0.3.5 readr_2.1.3 tidyr_1.2.1 tibble_3.1.8
## [9] ggplot2_3.3.6 tidyverse_1.3.2
##
## loaded via a namespace (and not attached):
## [1] tidyselect_1.2.0 xfun_0.34 bslib_0.4.0
## [4] haven_2.5.1 gargle_1.2.1 colorspace_2.0-3
## [7] vctrs_0.5.0 generics_0.1.3 htmltools_0.5.3
## [10] yaml_2.3.6 utf8_1.2.2 rlang_1.0.6
## [13] jquerylib_0.1.4 pillar_1.8.1 withr_2.5.0
## [16] DBI_1.1.3 dbplyr_2.2.1 readxl_1.4.1
## [19] modelr_0.1.9 lifecycle_1.0.3 munsell_0.5.0
## [22] gtable_0.3.1 cellranger_1.1.0 rvest_1.0.3
## [25] evaluate_0.17 labeling_0.4.2 knitr_1.40
## [28] tzdb_0.3.0 fastmap_1.1.0 fansi_1.0.3
## [31] highr_0.9 broom_1.0.1 backports_1.4.1
## [34] scales_1.2.1 googlesheets4_1.0.1 cachem_1.0.6
## [37] jsonlite_1.8.3 farver_2.1.1 fs_1.5.2
## [40] hms_1.1.2 digest_0.6.30 stringi_1.7.8
## [43] grid_4.2.1 cli_3.4.1 tools_4.2.1
## [46] magrittr_2.0.3 sass_0.4.2 crayon_1.5.2
## [49] pkgconfig_2.0.3 ellipsis_0.3.2 xml2_1.3.3
## [52] reprex_2.0.2 googledrive_2.0.0 lubridate_1.8.0
## [55] assertthat_0.2.1 rmarkdown_2.17 httr_1.4.4
## [58] R6_2.5.1 compiler_4.2.1