Benchmarking Isosceles on Sequin data
Preparing the data
Load the required packages:
library(tidyverse)
library(glue)
library(SingleCellExperiment)
library(patchwork)Global parameters:
result_dir <- file.path("..", "sequin_analysis", "report_data")
sample_ids <- c("mixA_ont", "mixB_ont")Helper functions:
mean_rel_diff <- function(x, y) {
mean(abs(x - y) / ((x + y) / 2))
}Prepare reference TPM data:
ref_tx_df <- read_delim(file.path(
result_dir, "sequin_transcripts.tsv"
))
ref_tx_df$MIX_A <- ref_tx_df$MIX_A / sum(ref_tx_df$MIX_A) * 1e6
ref_tx_df$MIX_B <- ref_tx_df$MIX_B / sum(ref_tx_df$MIX_B) * 1e6Prepare Isosceles TPM data:
isosceles_tpm_df <- as.data.frame(sapply(sample_ids, function(sample_id) {
se_transcript <- readRDS(file.path(
result_dir, glue("isosceles_{sample_id}_se_transcript.rds")
))
tx_df <- as.data.frame(rowData(se_transcript))
tpm_values <- assay(se_transcript, "tpm")[, 1]
names(tpm_values) <- tx_df$compatible_tx
tpm_values <- tpm_values[ref_tx_df$NAME]
return(tpm_values)
}))Prepare IsoQuant TPM data:
isoquant_tpm_df <- as.data.frame(sapply(sample_ids, function(sample_id) {
tx_df <- read_delim(file.path(
result_dir, glue("isoquant_{sample_id}_transcript_tpm.tsv")
))
tpm_values <- tx_df$TPM
names(tpm_values) <- tx_df$`#feature_id`
tpm_values <- tpm_values[ref_tx_df$NAME]
tpm_values[is.na(tpm_values)] <- 0
names(tpm_values) <- ref_tx_df$NAME
return(tpm_values)
}))Prepare bambu TPM data:
bambu_tpm_df <- as.data.frame(sapply(sample_ids, function(sample_id) {
tx_df <- read_delim(file.path(
result_dir, glue("bambu_{sample_id}_counts_transcript.txt")
))
tpm_values <- tx_df[, 3, drop = TRUE]
names(tpm_values) <- tx_df$TXNAME
tpm_values <- tpm_values[ref_tx_df$NAME]
tpm_values <- tpm_values / sum(tpm_values) * 1e6
return(tpm_values)
}))Prepare Flair TPM data:
flair_tpm_df <- as.data.frame(sapply(sample_ids, function(sample_id) {
tx_df <- read_delim(file.path(
result_dir, glue("flair_{sample_id}_counts_matrix.tsv")
))
tx_df <- tx_df[grepl("^R\\d+", tx_df$ids),]
tx_df$ids <- sapply(
strsplit(tx_df$ids, "_"),
function(x) {paste0(x[1:3], collapse = "_")}
)
tx_df <- tx_df[!grepl("-", tx_df$ids),]
tpm_values <- tx_df$sample_condition_batch
names(tpm_values) <- tx_df$ids
tpm_values <- tpm_values[ref_tx_df$NAME]
tpm_values[is.na(tpm_values)] <- 0
names(tpm_values) <- ref_tx_df$NAME
tpm_values <- tpm_values / sum(tpm_values) * 1e6
return(tpm_values)
}))Prepare NanoCount TPM data:
nanocount_tpm_df <- as.data.frame(sapply(sample_ids, function(sample_id) {
tx_df <- read_delim(file.path(
result_dir, glue("nanocount_{sample_id}.tsv")
))
tpm_values <- tx_df$tpm
names(tpm_values) <- tx_df$transcript_name
tpm_values <- tpm_values[ref_tx_df$NAME]
tpm_values[is.na(tpm_values)] <- 0
names(tpm_values) <- ref_tx_df$NAME
return(tpm_values)
}))Prepare LIQA TPM data:
liqa_tpm_df <- as.data.frame(sapply(sample_ids, function(sample_id) {
tx_df <- read_delim(file.path(
result_dir, glue("liqa_{sample_id}.tsv")
))
tpm_values <- tx_df$ReadPerGene_corrected
names(tpm_values) <- tx_df$IsoformName
tpm_values <- tpm_values[ref_tx_df$NAME]
tpm_values[is.na(tpm_values)] <- 0
names(tpm_values) <- ref_tx_df$NAME
tpm_values <- tpm_values / sum(tpm_values) * 1e6
return(tpm_values)
}))Prepare ESPRESSO TPM data:
espresso_tpm_df <- as.data.frame(sapply(sample_ids, function(sample_id) {
tx_df <- read_delim(file.path(
result_dir, glue("espresso_{sample_id}_abundance.esp")
))
tpm_values <- tx_df$sample
names(tpm_values) <- tx_df$transcript_ID
tpm_values <- tpm_values[ref_tx_df$NAME]
tpm_values[is.na(tpm_values)] <- 0
names(tpm_values) <- ref_tx_df$NAME
tpm_values <- tpm_values / sum(tpm_values) * 1e6
return(tpm_values)
}))Results
Prepare Spearman correlation data:
tools <- c("Isosceles", "IsoQuant", "Bambu", "Flair", "NanoCount", "LIQA",
"ESPRESSO")
spearman_cor_df <- data.frame(
tool = rep(tools, each = 3),
mix = rep(c("mix_A", "mix_B", "both"), 7),
value = c(
cor(isosceles_tpm_df$mixA_ont, ref_tx_df$MIX_A, method = "spearman"),
cor(isosceles_tpm_df$mixB_ont, ref_tx_df$MIX_B, method = "spearman"),
cor(c(isosceles_tpm_df$mixA_ont, isosceles_tpm_df$mixB_ont),
c(ref_tx_df$MIX_A, ref_tx_df$MIX_B), method = "spearman"),
cor(isoquant_tpm_df$mixA_ont, ref_tx_df$MIX_A, method = "spearman"),
cor(isoquant_tpm_df$mixB_ont, ref_tx_df$MIX_B, method = "spearman"),
cor(c(isoquant_tpm_df$mixA_ont, isoquant_tpm_df$mixB_ont),
c(ref_tx_df$MIX_A, ref_tx_df$MIX_B), method = "spearman"),
cor(bambu_tpm_df$mixA_ont, ref_tx_df$MIX_A, method = "spearman"),
cor(bambu_tpm_df$mixB_ont, ref_tx_df$MIX_B, method = "spearman"),
cor(c(bambu_tpm_df$mixA_ont, bambu_tpm_df$mixB_ont),
c(ref_tx_df$MIX_A, ref_tx_df$MIX_B), method = "spearman"),
cor(flair_tpm_df$mixA_ont, ref_tx_df$MIX_A, method = "spearman"),
cor(flair_tpm_df$mixB_ont, ref_tx_df$MIX_B, method = "spearman"),
cor(c(flair_tpm_df$mixA_ont, flair_tpm_df$mixB_ont),
c(ref_tx_df$MIX_A, ref_tx_df$MIX_B), method = "spearman"),
cor(nanocount_tpm_df$mixA_ont, ref_tx_df$MIX_A, method = "spearman"),
cor(nanocount_tpm_df$mixB_ont, ref_tx_df$MIX_B, method = "spearman"),
cor(c(nanocount_tpm_df$mixA_ont, nanocount_tpm_df$mixB_ont),
c(ref_tx_df$MIX_A, ref_tx_df$MIX_B), method = "spearman"),
cor(liqa_tpm_df$mixA_ont, ref_tx_df$MIX_A, method = "spearman"),
cor(liqa_tpm_df$mixB_ont, ref_tx_df$MIX_B, method = "spearman"),
cor(c(liqa_tpm_df$mixA_ont, liqa_tpm_df$mixB_ont),
c(ref_tx_df$MIX_A, ref_tx_df$MIX_B), method = "spearman"),
cor(espresso_tpm_df$mixA_ont, ref_tx_df$MIX_A, method = "spearman"),
cor(espresso_tpm_df$mixB_ont, ref_tx_df$MIX_B, method = "spearman"),
cor(c(espresso_tpm_df$mixA_ont, espresso_tpm_df$mixB_ont),
c(ref_tx_df$MIX_A, ref_tx_df$MIX_B), method = "spearman")
)
)Spearman correlation data:
spearman_cor_dfPrepare mean relative difference data:
tools <- c("Isosceles", "IsoQuant", "Bambu", "Flair", "NanoCount", "LIQA",
"ESPRESSO")
mean_rel_diff_df <- data.frame(
tool = rep(tools, each = 3),
mix = rep(c("mix_A", "mix_B", "both"), 7),
value = c(
mean_rel_diff(isosceles_tpm_df$mixA_ont, ref_tx_df$MIX_A),
mean_rel_diff(isosceles_tpm_df$mixB_ont, ref_tx_df$MIX_B),
mean_rel_diff(c(isosceles_tpm_df$mixA_ont, isosceles_tpm_df$mixB_ont),
c(ref_tx_df$MIX_A, ref_tx_df$MIX_B)),
mean_rel_diff(isoquant_tpm_df$mixA_ont, ref_tx_df$MIX_A),
mean_rel_diff(isoquant_tpm_df$mixB_ont, ref_tx_df$MIX_B),
mean_rel_diff(c(isoquant_tpm_df$mixA_ont, isoquant_tpm_df$mixB_ont),
c(ref_tx_df$MIX_A, ref_tx_df$MIX_B)),
mean_rel_diff(bambu_tpm_df$mixA_ont, ref_tx_df$MIX_A),
mean_rel_diff(bambu_tpm_df$mixB_ont, ref_tx_df$MIX_B),
mean_rel_diff(c(bambu_tpm_df$mixA_ont, bambu_tpm_df$mixB_ont),
c(ref_tx_df$MIX_A, ref_tx_df$MIX_B)),
mean_rel_diff(flair_tpm_df$mixA_ont, ref_tx_df$MIX_A),
mean_rel_diff(flair_tpm_df$mixB_ont, ref_tx_df$MIX_B),
mean_rel_diff(c(flair_tpm_df$mixA_ont, flair_tpm_df$mixB_ont),
c(ref_tx_df$MIX_A, ref_tx_df$MIX_B)),
mean_rel_diff(nanocount_tpm_df$mixA_ont, ref_tx_df$MIX_A),
mean_rel_diff(nanocount_tpm_df$mixB_ont, ref_tx_df$MIX_B),
mean_rel_diff(c(nanocount_tpm_df$mixA_ont, nanocount_tpm_df$mixB_ont),
c(ref_tx_df$MIX_A, ref_tx_df$MIX_B)),
mean_rel_diff(liqa_tpm_df$mixA_ont, ref_tx_df$MIX_A),
mean_rel_diff(liqa_tpm_df$mixB_ont, ref_tx_df$MIX_B),
mean_rel_diff(c(liqa_tpm_df$mixA_ont, liqa_tpm_df$mixB_ont),
c(ref_tx_df$MIX_A, ref_tx_df$MIX_B)),
mean_rel_diff(espresso_tpm_df$mixA_ont, ref_tx_df$MIX_A),
mean_rel_diff(espresso_tpm_df$mixB_ont, ref_tx_df$MIX_B),
mean_rel_diff(c(espresso_tpm_df$mixA_ont, espresso_tpm_df$mixB_ont),
c(ref_tx_df$MIX_A, ref_tx_df$MIX_B))
)
)Mean relative difference data:
mean_rel_diff_dfHelper function for plotting Spearman correlation data:
plot_spearman_cor <- function(mix_value, title = FALSE, max_x = 1) {
plot_df <- spearman_cor_df %>%
filter(mix == mix_value)
plot_df$tool <- fct_reorder(plot_df$tool, plot_df$value)
bar_plot <- ggplot(plot_df, mapping = aes(x = value,
y = tool,
fill = tool)) +
geom_col() +
geom_text(
aes(label = format(round(value, 2), nsmall = 2)),
size = 3,
nudge_x = 0.03,
hjust = "outward"
) +
scale_fill_manual(values = c(Bambu = "mediumseagreen",
ESPRESSO = "lightsalmon",
Flair = "mediumpurple",
Isosceles = "black",
IsoQuant = "gold",
LIQA = "deeppink",
NanoCount = "cornflowerblue")) +
lims(x = c(0, max_x)) +
labs(
x = "Correlation value",
y = ""
) +
theme_classic() +
theme(legend.position = "none",
axis.text = element_text(size = 7),
axis.title.y = element_blank(),
axis.title.x = element_text(size = 9),
axis.text.x = element_blank(),
axis.ticks.x = element_blank())
if (title) {
bar_plot <- bar_plot +
labs(title = mix_value)
}
return(bar_plot)
}Plot Spearman correlation (mix_A):
plot_spearman_cor("mix_A", max_x = 1.1)
Plot Spearman correlation (mix_B):
plot_spearman_cor("mix_B", max_x = 1.1)
Plot Spearman correlation (both mixes):
plot_spearman_cor("both", max_x = 1.1)
Plot Spearman correlation (mixA & mix_B):
plot_spearman_cor("mix_A", title = TRUE, max_x = 1.25) +
plot_spearman_cor("mix_B", title = TRUE, max_x = 1.25)
Helper function for plotting mean relative difference data:
plot_mean_rel_diff <- function(mix_value, title = FALSE, max_x = 2) {
plot_df <- mean_rel_diff_df %>%
filter(mix == mix_value)
plot_df$tool <- fct_reorder(plot_df$tool, plot_df$value)
plot_df$tool <- fct_rev(plot_df$tool)
bar_plot <- ggplot(plot_df, mapping = aes(x = value,
y = tool,
fill = tool)) +
geom_col() +
geom_text(
aes(label = format(round(value, 2), nsmall = 2)),
size = 3,
nudge_x = 0.03,
hjust = "outward"
) +
scale_fill_manual(values = c(Bambu = "mediumseagreen",
ESPRESSO = "lightsalmon",
Flair = "mediumpurple",
Isosceles = "black",
IsoQuant = "gold",
LIQA = "deeppink",
NanoCount = "cornflowerblue")) +
lims(x = c(0, max_x)) +
labs(
x = "Mean rel. diff.",
y = ""
) +
theme_classic() +
theme(legend.position = "none",
axis.text = element_text(size = 7),
axis.title.y = element_blank(),
axis.title.x = element_text(size = 9),
axis.text.x = element_blank(),
axis.ticks.x = element_blank())
if (title) {
bar_plot <- bar_plot +
labs(title = mix_value)
}
return(bar_plot)
}Plot mean relative difference (mix_A):
plot_mean_rel_diff("mix_A", max_x = 1.3)
Plot mean relative difference (mix_B):
plot_mean_rel_diff("mix_B", max_x = 1.3)
Plot mean relative difference (both mixes):
plot_mean_rel_diff("both", max_x = 1.3)
Plot mean relative difference (mixA & mix_B):
plot_mean_rel_diff("mix_A", title = TRUE, max_x = 1.47) +
plot_mean_rel_diff("mix_B", title = TRUE, max_x = 1.47)
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] stats4 stats graphics grDevices utils datasets methods
## [8] base
##
## other attached packages:
## [1] patchwork_1.1.2 SingleCellExperiment_1.18.1
## [3] SummarizedExperiment_1.26.1 Biobase_2.56.0
## [5] GenomicRanges_1.48.0 GenomeInfoDb_1.32.4
## [7] IRanges_2.30.1 S4Vectors_0.34.0
## [9] BiocGenerics_0.42.0 MatrixGenerics_1.8.1
## [11] matrixStats_0.62.0 glue_1.6.2
## [13] forcats_0.5.2 stringr_1.4.1
## [15] dplyr_1.0.10 purrr_0.3.5
## [17] readr_2.1.3 tidyr_1.2.1
## [19] tibble_3.1.8 ggplot2_3.3.6
## [21] tidyverse_1.3.2
##
## loaded via a namespace (and not attached):
## [1] httr_1.4.4 sass_0.4.2 bit64_4.0.5
## [4] vroom_1.6.0 jsonlite_1.8.3 modelr_0.1.9
## [7] bslib_0.4.0 assertthat_0.2.1 highr_0.9
## [10] GenomeInfoDbData_1.2.8 googlesheets4_1.0.1 cellranger_1.1.0
## [13] yaml_2.3.6 pillar_1.8.1 backports_1.4.1
## [16] lattice_0.20-45 digest_0.6.30 XVector_0.36.0
## [19] rvest_1.0.3 colorspace_2.0-3 htmltools_0.5.3
## [22] Matrix_1.5-1 pkgconfig_2.0.3 broom_1.0.1
## [25] haven_2.5.1 zlibbioc_1.42.0 scales_1.2.1
## [28] tzdb_0.3.0 googledrive_2.0.0 farver_2.1.1
## [31] generics_0.1.3 ellipsis_0.3.2 cachem_1.0.6
## [34] withr_2.5.0 cli_3.4.1 magrittr_2.0.3
## [37] crayon_1.5.2 readxl_1.4.1 evaluate_0.17
## [40] fs_1.5.2 fansi_1.0.3 xml2_1.3.3
## [43] tools_4.2.1 hms_1.1.2 gargle_1.2.1
## [46] lifecycle_1.0.3 munsell_0.5.0 reprex_2.0.2
## [49] DelayedArray_0.22.0 compiler_4.2.1 jquerylib_0.1.4
## [52] rlang_1.0.6 grid_4.2.1 RCurl_1.98-1.9
## [55] labeling_0.4.2 bitops_1.0-7 rmarkdown_2.17
## [58] gtable_0.3.1 DBI_1.1.3 R6_2.5.1
## [61] lubridate_1.8.0 knitr_1.40 bit_4.0.4
## [64] fastmap_1.1.0 utf8_1.2.2 stringi_1.7.8
## [67] parallel_4.2.1 vctrs_0.5.0 dbplyr_2.2.1
## [70] tidyselect_1.2.0 xfun_0.34