Benchmarking Isosceles on Pacbio, Nanopore and Illumina GM12878 data
Preparing the data
Load the required packages:
library(tidyverse)
library(glue)
library(Matrix)
library(SingleCellExperiment)Global parameters:
result_dir <- file.path("..", "pacbio_analysis", "report_data")
illumina_sample_id <- "ENCFF485OUK"
pacbio_sample_id <- "ENCFF450VAU"
ont_sample_id <- "NA12878_ONT"
min_tpm <- 1Prepare Isosceles transcript ID map data:
se <- readRDS(file.path(
result_dir, glue("isosceles_{pacbio_sample_id}_se_transcript.rds")
))
tx_df <- rowData(se)[, c("transcript_id", "compatible_tx")]
tx_df <- as.data.frame(tx_df)
rownames(tx_df) <- NULL
tx_df$compatible_tx <- strsplit(tx_df$compatible_tx, ",")
tx_df <- unchop(tx_df, compatible_tx)
tx2isosceles <- setNames(tx_df$transcript_id, tx_df$compatible_tx)
all_isosceles_ids <- sort(unique(tx2isosceles))
all_tx_ids <- sort(names(tx2isosceles))Prepare Illumina TPM data:
illumina_df <- read_delim(file.path(
result_dir, glue("illumina_{illumina_sample_id}.tsv")
))
illumina_df <- filter(illumina_df, grepl("^ENST", transcript_id), TPM > min_tpm)
illumina_tpm <- illumina_df$TPM
names(illumina_tpm) <- sapply(strsplit(illumina_df$transcript_id, "\\."), "[", 1)
illumina_tpm <- illumina_tpm[names(illumina_tpm) %in% all_tx_ids]
illumina_tpm <- illumina_tpm / sum(illumina_tpm) * 1e6
illumina_tx_ids <- names(illumina_tpm)
illumina_isosceles_tpm <- tapply(
illumina_tpm, tx2isosceles[names(illumina_tpm)], sum
)
illumina_isosceles_tpm <- illumina_isosceles_tpm /
sum(illumina_isosceles_tpm) * 1e6
illumina_isosceles_tx_ids <- names(illumina_isosceles_tpm)Prepare Isosceles TPM data:
se_pacbio <- readRDS(file.path(
result_dir, glue("isosceles_{pacbio_sample_id}_se_transcript.rds")
))
isosceles_pacbio_tpm <- assay(se_pacbio, "tpm")[, 1]
isosceles_pacbio_tpm <- isosceles_pacbio_tpm[illumina_isosceles_tx_ids]
isosceles_pacbio_tpm <- isosceles_pacbio_tpm / sum(isosceles_pacbio_tpm) * 1e6
se_ont <- readRDS(file.path(
result_dir, glue("isosceles_{ont_sample_id}_se_transcript.rds")
))
isosceles_ont_tpm <- assay(se_ont, "tpm")[, 1]
isosceles_ont_tpm <- isosceles_ont_tpm[illumina_isosceles_tx_ids]
isosceles_ont_tpm <- isosceles_ont_tpm / sum(isosceles_ont_tpm) * 1e6Prepare IsoQuant TPM data:
isoquant_pacbio_df <- read_delim(file.path(
result_dir, glue("isoquant_{pacbio_sample_id}_transcript_tpm.tsv")
))
isoquant_pacbio_tpm <- isoquant_pacbio_df$TPM
names(isoquant_pacbio_tpm) <- isoquant_pacbio_df[, 1, drop = TRUE]
isoquant_pacbio_tpm <- isoquant_pacbio_tpm[illumina_tx_ids]
isoquant_pacbio_tpm[is.na(isoquant_pacbio_tpm)] <- 0
names(isoquant_pacbio_tpm) <- illumina_tx_ids
isoquant_pacbio_tpm <- isoquant_pacbio_tpm / sum(isoquant_pacbio_tpm) * 1e6
isoquant_ont_df <- read_delim(file.path(
result_dir, glue("isoquant_{ont_sample_id}_transcript_tpm.tsv")
))
isoquant_ont_tpm <- isoquant_ont_df$TPM
names(isoquant_ont_tpm) <- isoquant_ont_df[, 1, drop = TRUE]
isoquant_ont_tpm <- isoquant_ont_tpm[illumina_tx_ids]
isoquant_ont_tpm[is.na(isoquant_ont_tpm)] <- 0
names(isoquant_ont_tpm) <- illumina_tx_ids
isoquant_ont_tpm <- isoquant_ont_tpm / sum(isoquant_ont_tpm) * 1e6Prepare bambu TPM data:
bambu_pacbio_df <- read_delim(file.path(
result_dir, glue("bambu_{pacbio_sample_id}_counts_transcript.txt")
))
bambu_pacbio_tpm <- bambu_pacbio_df[, 3, drop = TRUE]
names(bambu_pacbio_tpm) <- bambu_pacbio_df$TXNAME
bambu_pacbio_tpm <- bambu_pacbio_tpm / sum(bambu_pacbio_tpm) * 1e6
bambu_pacbio_tpm <- bambu_pacbio_tpm[illumina_tx_ids]
bambu_pacbio_tpm <- bambu_pacbio_tpm / sum(bambu_pacbio_tpm) * 1e6
bambu_ont_df <- read_delim(file.path(
result_dir, glue("bambu_{ont_sample_id}_counts_transcript.txt")
))
bambu_ont_tpm <- bambu_ont_df[, 3, drop = TRUE]
names(bambu_ont_tpm) <- bambu_ont_df$TXNAME
bambu_ont_tpm <- bambu_ont_tpm / sum(bambu_ont_tpm) * 1e6
bambu_ont_tpm <- bambu_ont_tpm[illumina_tx_ids]
bambu_ont_tpm <- bambu_ont_tpm / sum(bambu_ont_tpm) * 1e6Results
Prepare Spearman correlation data:
spearman_cor_df <- data.frame(
tool = c("Isosceles", "IsoQuant", "Bambu"),
pacbio_illumina = c(
cor(isosceles_pacbio_tpm, illumina_isosceles_tpm, method = "spearman"),
cor(isoquant_pacbio_tpm, illumina_tpm, method = "spearman"),
cor(bambu_pacbio_tpm, illumina_tpm, method = "spearman")
),
ont_illumina = c(
cor(isosceles_ont_tpm, illumina_isosceles_tpm, method = "spearman"),
cor(isoquant_ont_tpm, illumina_tpm, method = "spearman"),
cor(bambu_ont_tpm, illumina_tpm, method = "spearman")
),
pacbio_ont = c(
cor(isosceles_pacbio_tpm, isosceles_ont_tpm, method = "spearman"),
cor(isoquant_pacbio_tpm, isoquant_ont_tpm, method = "spearman"),
cor(bambu_pacbio_tpm, bambu_ont_tpm, method = "spearman")
)
)Spearman correlation data:
spearman_cor_dfPlot Spearman correlation (Pacbio vs Illumina):
spearman_cor_df$value <- spearman_cor_df$pacbio_illumina
spearman_cor_df$tool <- fct_reorder(spearman_cor_df$tool,
spearman_cor_df$value)
ggplot(spearman_cor_df, mapping = aes(x = value,
y = tool,
fill = tool)) +
geom_col() +
geom_text(
aes(label = format(round(value, 2), nsmall = 2)),
size = 3,
hjust = -0.4
) +
scale_fill_manual(values = c(Bambu = "mediumseagreen",
Isosceles = "black",
IsoQuant = "gold")) +
lims(x = c(0, 1)) +
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())
Plot Spearman correlation (Nanopore vs Illumina):
spearman_cor_df$value <- spearman_cor_df$ont_illumina
spearman_cor_df$tool <- fct_reorder(spearman_cor_df$tool,
spearman_cor_df$value)
ggplot(spearman_cor_df, mapping = aes(x = value,
y = tool,
fill = tool)) +
geom_col() +
geom_text(
aes(label = format(round(value, 2), nsmall = 2)),
size = 3,
hjust = -0.4
) +
scale_fill_manual(values = c(Bambu = "mediumseagreen",
Isosceles = "black",
IsoQuant = "gold")) +
lims(x = c(0, 1)) +
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())
Plot Spearman correlation (Pacbio vs Nanopore):
spearman_cor_df$value <- spearman_cor_df$pacbio_ont
spearman_cor_df$tool <- fct_reorder(spearman_cor_df$tool,
spearman_cor_df$value)
ggplot(spearman_cor_df, mapping = aes(x = value,
y = tool,
fill = tool)) +
geom_col() +
geom_text(
aes(label = format(round(value, 2), nsmall = 2)),
size = 3,
hjust = -0.4
) +
scale_fill_manual(values = c(Bambu = "mediumseagreen",
Isosceles = "black",
IsoQuant = "gold")) +
lims(x = c(0, 1)) +
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())
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] Matrix_1.5-1 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] pkgconfig_2.0.3 broom_1.0.1 haven_2.5.1
## [25] zlibbioc_1.42.0 scales_1.2.1 tzdb_0.3.0
## [28] googledrive_2.0.0 farver_2.1.1 generics_0.1.3
## [31] ellipsis_0.3.2 cachem_1.0.6 withr_2.5.0
## [34] cli_3.4.1 magrittr_2.0.3 crayon_1.5.2
## [37] readxl_1.4.1 evaluate_0.17 fs_1.5.2
## [40] fansi_1.0.3 xml2_1.3.3 tools_4.2.1
## [43] hms_1.1.2 gargle_1.2.1 lifecycle_1.0.3
## [46] munsell_0.5.0 reprex_2.0.2 DelayedArray_0.22.0
## [49] compiler_4.2.1 jquerylib_0.1.4 rlang_1.0.6
## [52] grid_4.2.1 RCurl_1.98-1.9 labeling_0.4.2
## [55] bitops_1.0-7 rmarkdown_2.17 gtable_0.3.1
## [58] DBI_1.1.3 R6_2.5.1 lubridate_1.8.0
## [61] knitr_1.40 bit_4.0.4 fastmap_1.1.0
## [64] utf8_1.2.2 stringi_1.7.8 parallel_4.2.1
## [67] vctrs_0.5.0 dbplyr_2.2.1 tidyselect_1.2.0
## [70] xfun_0.34