Benchmarking Isosceles on SIRV data
Preparing the data
Load the required packages:
library(tidyverse)
library(glue)
library(pheatmap)
library(viridis)
library(SingleCellExperiment)Global parameters:
result_dir <- file.path("..", "sirv_analysis", "report_data")
sample_ids <- "cdna"Over-annotations analysis
Prepare reference data:
annotations_correct_df <- readRDS(file.path(
result_dir, "annotations_correct.rds"
))$tx_df
correct_ref_tx_ids <- annotations_correct_df$compatible_tx %>%
strsplit(",") %>%
unlist() %>%
sort()Prepare Isosceles TPM data:
isosceles_tpm_list <- lapply(sample_ids, function(sample_id) {
se_transcript <- readRDS(file.path(
result_dir, glue("isosceles_{sample_id}_over_se_transcript.rds")
))
tx_df <- as.data.frame(rowData(se_transcript))
tpm_values <- assay(se_transcript, "tpm")[, 1]
names(tpm_values) <- sapply(strsplit(tx_df$compatible_tx, ","), "[", 1)
return(tpm_values)
})
names(isosceles_tpm_list) <- paste0("Isosceles.", sample_ids)
isosceles_correct_tpm_list <- lapply(isosceles_tpm_list, function(tpms) {
tpms[names(tpms) %in% correct_ref_tx_ids]
})
isosceles_fractions <- sapply(isosceles_correct_tpm_list, sum) / sapply(isosceles_tpm_list, sum)Prepare IsoQuant TPM data:
isoquant_tpm_list <- lapply(sample_ids, function(sample_id) {
tx_df <- read_delim(file.path(
result_dir, glue("isoquant_{sample_id}_over_transcript_tpm.tsv")
))
tpm_values <- tx_df$TPM
names(tpm_values) <- tx_df$`#feature_id`
return(tpm_values)
})
names(isoquant_tpm_list) <- paste0("IsoQuant.", sample_ids)
isoquant_correct_tpm_list <- lapply(isoquant_tpm_list, function(tpms) {
tpms[names(tpms) %in% correct_ref_tx_ids]
})
isoquant_fractions <- sapply(isoquant_correct_tpm_list, sum) / sapply(isoquant_tpm_list, sum)Prepare bambu TPM data:
bambu_tpm_list <- lapply(sample_ids, function(sample_id) {
tx_df <- read_delim(file.path(
result_dir, glue("bambu_{sample_id}_over_counts_transcript.txt")
))
tpm_values <- tx_df[, 3, drop = TRUE]
tpm_values <- tpm_values / sum(tpm_values) * 1e6
names(tpm_values) <- tx_df$TXNAME
return(tpm_values)
})
names(bambu_tpm_list) <- paste0("bambu.", sample_ids)
bambu_correct_tpm_list <- lapply(bambu_tpm_list, function(tpms) {
tpms[names(tpms) %in% correct_ref_tx_ids]
})
bambu_fractions <- sapply(bambu_correct_tpm_list, sum) / sapply(bambu_tpm_list, sum)Plot fractions of TPMs from correct transcripts:
plot_df <- data.frame(
tool = c("Isosceles", "IsoQuant", "Bambu"),
value = c(
isosceles_fractions["Isosceles.cdna"],
isoquant_fractions["IsoQuant.cdna"],
bambu_fractions["bambu.cdna"]
)
)
plot_df$tool <- fct_reorder(plot_df$tool, plot_df$value)
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",
Isosceles = "black",
IsoQuant = "gold")) +
lims(x = c(0, 1.1)) +
labs(
x = "Read precision (expression from correct transcripts)",
y = ""
) +
theme_classic() +
theme(legend.position = "none",
axis.text = element_text(size = 7),
axis.title.y = element_blank(),
axis.title.x = element_text(size = 8),
axis.text.x = element_blank(),
axis.ticks.x = element_blank())
Insufficient annotations analysis
Prepare reference data:
annotations_correct_df <- readRDS(file.path(
result_dir, "annotations_correct.rds"
))$tx_df
annotations_insufficient_df <- readRDS(file.path(
result_dir, "annotations_insufficient.rds"
))$tx_df
correct_intron_positions <- annotations_correct_df$intron_positions
correct_intron_positions <- unique(correct_intron_positions[
!is.na(correct_intron_positions)
])
insufficient_intron_positions <- annotations_insufficient_df$intron_positions
insufficient_intron_positions <- unique(insufficient_intron_positions[
!is.na(insufficient_intron_positions)
])Helper functions:
precision <- function(intron_positions) {
TP <- sum(intron_positions %in% correct_intron_positions)
FP <- sum(!intron_positions %in% correct_intron_positions)
return(TP / (FP + TP))
}
recall <- function(intron_positions) {
TP <- sum(intron_positions %in% correct_intron_positions)
P <- length(correct_intron_positions)
return(TP / P)
}
f1_score <- function(intron_positions) {
TP <- sum(intron_positions %in% correct_intron_positions)
FP <- sum(!intron_positions %in% correct_intron_positions)
P <- length(correct_intron_positions)
FN <- P - TP
return(2 * TP / (2*TP + FP + FN))
}
classify_tx <- function(intron_positions) {
return(c(
TP = sum(intron_positions %in% correct_intron_positions),
TP_KNOWN = sum(intron_positions %in% insufficient_intron_positions),
TP_NOVEL = sum(
(intron_positions %in% correct_intron_positions) &
!(intron_positions %in% insufficient_intron_positions)
),
FN = length(setdiff(correct_intron_positions, intron_positions)),
FN_KNOWN = length(setdiff(
insufficient_intron_positions, intron_positions
)),
FN_NOVEL = length(setdiff(
setdiff(correct_intron_positions, insufficient_intron_positions),
intron_positions
)),
FP = sum(!intron_positions %in% correct_intron_positions)
))
}Prepare Isosceles intron position data:
isosceles_intron_list <- lapply(sample_ids, function(sample_id) {
se_transcript <- readRDS(file.path(
result_dir, glue("isosceles_{sample_id}_insufficient_se_transcript.rds")
))
tpm_values <- assay(se_transcript, "tpm")[, 1]
se_transcript <- se_transcript[tpm_values > 0,]
intron_positions <- rowData(se_transcript)$intron_positions
intron_positions <- unique(intron_positions[!is.na(intron_positions)])
return(intron_positions)
})
names(isosceles_intron_list) <- paste0("Isosceles.", sample_ids)Prepare IsoQuant intron position data:
isoquant_intron_list <- lapply(sample_ids, function(sample_id) {
splicing_df <- Isosceles:::prepare_reference_annotations(file.path(
result_dir, glue("isoquant_{sample_id}_insufficient_extended_annotations.gtf")
))$splicing_df
tpm_df <- read_delim(file.path(
result_dir, glue("isoquant_{sample_id}_insufficient_transcript_tpm.tsv")
))
tpm_df <- tpm_df[tpm_df$TPM > 0,]
expressed_tx_ids <- tpm_df[, 1, drop = TRUE]
tpm_selector <- sapply(strsplit(splicing_df$compatible_tx, ","),
function(x) {any(x %in% expressed_tx_ids)})
splicing_df <- splicing_df[tpm_selector,]
intron_positions <- splicing_df$intron_positions
return(intron_positions)
})
names(isoquant_intron_list) <- paste0("IsoQuant.", sample_ids)Prepare bambu intron position data:
bambu_intron_list <- lapply(sample_ids, function(sample_id) {
splicing_df <- Isosceles:::prepare_reference_annotations(file.path(
result_dir, glue("bambu_{sample_id}_insufficient_extended_annotations.gtf")
))$splicing_df
count_df <- read_delim(file.path(
result_dir, glue("bambu_{sample_id}_insufficient_counts_transcript.txt")
))
count_df <- count_df[count_df[, 3, drop = TRUE] > 0,]
expressed_tx_ids <- count_df$TXNAME
count_selector <- sapply(strsplit(splicing_df$compatible_tx, ","),
function(x) {any(x %in% expressed_tx_ids)})
splicing_df <- splicing_df[count_selector,]
intron_positions <- splicing_df$intron_positions
return(intron_positions)
})
names(bambu_intron_list) <- paste0("bambu.", sample_ids)Plot transcript detection metrics:
tools <- c("Isosceles", "IsoQuant", "Bambu")
plot_df <- data.frame(
tool = rep(tools, each = 3),
metric = rep(c("Precision", "Recall", "F1"), 3),
value = c(
precision(isosceles_intron_list[["Isosceles.cdna"]]),
recall(isosceles_intron_list[["Isosceles.cdna"]]),
f1_score(isosceles_intron_list[["Isosceles.cdna"]]),
precision(isoquant_intron_list[["IsoQuant.cdna"]]),
recall(isoquant_intron_list[["IsoQuant.cdna"]]),
f1_score(isoquant_intron_list[["IsoQuant.cdna"]]),
precision(bambu_intron_list[["bambu.cdna"]]),
recall(bambu_intron_list[["bambu.cdna"]]),
f1_score(bambu_intron_list[["bambu.cdna"]])
)
)
plot_df$tool <- fct_reorder(plot_df$tool, plot_df$value)
plot_df$metric <- fct_rev(factor(
plot_df$metric, levels = c("Precision", "Recall", "F1")
))
ggplot(plot_df, mapping = aes(x = value,
y = metric,
fill = tool)) +
geom_col(position = position_dodge2()) +
geom_text(
aes(label = format(round(value, 2), nsmall = 2)),
position = position_dodge2(width = .9),
size = 3,
hjust = -0.3
) +
lims(x = c(0, 1.1)) +
scale_fill_manual(values = c(Isosceles = "black",
Bambu = "mediumseagreen",
IsoQuant = "gold")) +
theme_classic() +
theme(legend.title = element_blank(),
legend.position = "top",
legend.key.size = unit(0.3, "cm"),
legend.text = element_text(size = 7),
axis.text = element_text(size = 7),
axis.title.x = element_blank(),
axis.title.y = element_blank(),
axis.text.x = element_blank(),
axis.ticks.x = element_blank())
Transcript detection results:
classification_df <- data.frame(
Isosceles = classify_tx(isosceles_intron_list[["Isosceles.cdna"]]),
IsoQuant = classify_tx(isoquant_intron_list[["IsoQuant.cdna"]]),
Bambu = classify_tx(bambu_intron_list[["bambu.cdna"]])
)
classification_dfPlot transcript detection results:
plot_df <- classification_df
plot_df$metric <- rownames(plot_df)
plot_df <- plot_df %>%
pivot_longer(-metric, names_to = "tool", values_to = "value") %>%
filter(metric %in% c("TP_KNOWN", "TP_NOVEL", "FP", "FN_KNOWN", "FN_NOVEL"))
plot_df$metric <- fct_rev(factor(
plot_df$metric,
levels = c(c("TP_KNOWN", "TP_NOVEL", "FN_KNOWN", "FN_NOVEL", "FP"))
))
plot_df$tool <- paste0(
plot_df$tool, " ", ifelse(grepl("^TP", plot_df$metric), "TP", "FP/FN")
)
ggplot(plot_df, mapping = aes(x = value,
y = tool,
fill = metric)) +
geom_col() +
scale_fill_manual(values = c(TP_KNOWN = "darkgreen",
TP_NOVEL = "mediumseagreen",
FN_KNOWN = "royalblue4",
FN_NOVEL = "cornflowerblue",
FP = "darkred")) +
theme_classic() +
theme(legend.title = element_blank(),
legend.position = "top",
legend.key.size = unit(0.3, "cm"),
legend.text = element_text(size = 5),
axis.text = element_text(size = 7),
axis.title.x = element_blank(),
axis.title.y = element_blank())
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] SingleCellExperiment_1.18.1 SummarizedExperiment_1.26.1
## [3] Biobase_2.56.0 GenomicRanges_1.48.0
## [5] GenomeInfoDb_1.32.4 IRanges_2.30.1
## [7] S4Vectors_0.34.0 BiocGenerics_0.42.0
## [9] MatrixGenerics_1.8.1 matrixStats_0.62.0
## [11] viridis_0.6.2 viridisLite_0.4.1
## [13] pheatmap_1.0.12 glue_1.6.2
## [15] forcats_0.5.2 stringr_1.4.1
## [17] dplyr_1.0.10 purrr_0.3.5
## [19] readr_2.1.3 tidyr_1.2.1
## [21] tibble_3.1.8 ggplot2_3.3.6
## [23] tidyverse_1.3.2
##
## loaded via a namespace (and not attached):
## [1] googledrive_2.0.0 colorspace_2.0-3 rjson_0.2.21
## [4] ellipsis_0.3.2 XVector_0.36.0 fs_1.5.2
## [7] farver_2.1.1 bit64_4.0.5 AnnotationDbi_1.58.0
## [10] fansi_1.0.3 lubridate_1.8.0 xml2_1.3.3
## [13] codetools_0.2-18 cachem_1.0.6 knitr_1.40
## [16] jsonlite_1.8.3 Rsamtools_2.12.0 broom_1.0.1
## [19] dbplyr_2.2.1 png_0.1-7 compiler_4.2.1
## [22] httr_1.4.4 backports_1.4.1 assertthat_0.2.1
## [25] Matrix_1.5-1 fastmap_1.1.0 gargle_1.2.1
## [28] cli_3.4.1 htmltools_0.5.3 prettyunits_1.1.1
## [31] tools_4.2.1 igraph_1.3.5 gtable_0.3.1
## [34] GenomeInfoDbData_1.2.8 rappdirs_0.3.3 Rcpp_1.0.9
## [37] cellranger_1.1.0 jquerylib_0.1.4 vctrs_0.5.0
## [40] Biostrings_2.64.1 rtracklayer_1.56.1 xfun_0.34
## [43] rvest_1.0.3 lifecycle_1.0.3 restfulr_0.0.15
## [46] XML_3.99-0.12 googlesheets4_1.0.1 zlibbioc_1.42.0
## [49] scales_1.2.1 vroom_1.6.0 hms_1.1.2
## [52] parallel_4.2.1 RColorBrewer_1.1-3 yaml_2.3.6
## [55] curl_4.3.3 memoise_2.0.1 gridExtra_2.3
## [58] sass_0.4.2 biomaRt_2.52.0 stringi_1.7.8
## [61] RSQLite_2.2.18 highr_0.9 BiocIO_1.6.0
## [64] Isosceles_0.2.0 GenomicFeatures_1.48.4 filelock_1.0.2
## [67] BiocParallel_1.30.4 rlang_1.0.6 pkgconfig_2.0.3
## [70] bitops_1.0-7 evaluate_0.17 lattice_0.20-45
## [73] GenomicAlignments_1.32.1 labeling_0.4.2 bit_4.0.4
## [76] tidyselect_1.2.0 magrittr_2.0.3 R6_2.5.1
## [79] generics_0.1.3 DelayedArray_0.22.0 DBI_1.1.3
## [82] pillar_1.8.1 haven_2.5.1 withr_2.5.0
## [85] KEGGREST_1.36.3 RCurl_1.98-1.9 modelr_0.1.9
## [88] crayon_1.5.2 utf8_1.2.2 BiocFileCache_2.4.0
## [91] tzdb_0.3.0 rmarkdown_2.17 progress_1.2.2
## [94] grid_4.2.1 readxl_1.4.1 blob_1.2.3
## [97] reprex_2.0.2 digest_0.6.30 munsell_0.5.0
## [100] bslib_0.4.0