Merge branch 'r-test-infra-only' into 'dev'
[CW-7161] Add workflow-local R test infrastructure See merge request epi2melabs/workflows/wf-transcriptomes!236
This commit is contained in:
Chris Wright
commit
0633a6304d
@ -13,6 +13,8 @@ variables:
|
|||||||
CI_FLAVOUR: "new"
|
CI_FLAVOUR: "new"
|
||||||
PYTEST_CONTAINER_NAME: "wf-common"
|
PYTEST_CONTAINER_NAME: "wf-common"
|
||||||
PYTEST_CONTAINER_CONFIG_KEY: "common_sha"
|
PYTEST_CONTAINER_CONFIG_KEY: "common_sha"
|
||||||
|
RTEST_CONTAINER_NAME: "wf-transcriptomes-core"
|
||||||
|
RTEST_CONTAINER_CONFIG_KEY: "container_sha"
|
||||||
|
|
||||||
macos-run:
|
macos-run:
|
||||||
# Let's avoid those ARM64 runners for now
|
# Let's avoid those ARM64 runners for now
|
||||||
|
|||||||
449
bin/run_bambu.R
449
bin/run_bambu.R
@ -1,449 +0,0 @@
|
|||||||
#!/usr/bin/env Rscript
|
|
||||||
|
|
||||||
# Set seed for reproducibility
|
|
||||||
set.seed(42)
|
|
||||||
|
|
||||||
suppressPackageStartupMessages({
|
|
||||||
library(argparser)
|
|
||||||
library(bambu)
|
|
||||||
library(Rsamtools)
|
|
||||||
library(SummarizedExperiment)
|
|
||||||
library(jsonlite)
|
|
||||||
})
|
|
||||||
|
|
||||||
parser <- arg_parser("Run bambu transcript discovery and quantification.")
|
|
||||||
parser <- add_argument(parser, "--bam_dir", help = "Directory containing BAM files.")
|
|
||||||
parser <- add_argument(parser, "--bam_path", help = "Path to a single BAM file.")
|
|
||||||
parser <- add_argument(parser, "--sample_alias", help = "Alias to use for a single BAM file.")
|
|
||||||
parser <- add_argument(parser, "--sample_sheet", help = "Optional sample sheet CSV.")
|
|
||||||
parser <- add_argument(parser, "--annotation", help = "Reference annotation GTF/GFF.")
|
|
||||||
parser <- add_argument(parser, "--genome", help = "Reference genome FASTA.")
|
|
||||||
parser <- add_argument(parser, "--transcriptome_mode", help = "discover or fixed_annotation.", default = "discover")
|
|
||||||
parser <- add_argument(parser, "--threads", help = "Number of worker threads.", type = "numeric", default = 1)
|
|
||||||
parser <- add_argument(parser, "--ndr", help = "Optional novel discovery rate.", type = "numeric")
|
|
||||||
parser <- add_argument(parser, "--out_dir", help = "Output directory.")
|
|
||||||
argv <- parse_args(parser)
|
|
||||||
|
|
||||||
arg_missing <- function(value) {
|
|
||||||
if (is.null(value) || length(value) == 0 || all(is.na(value))) {
|
|
||||||
return(TRUE)
|
|
||||||
}
|
|
||||||
if (is.character(value)) {
|
|
||||||
return(all(!nzchar(value)))
|
|
||||||
}
|
|
||||||
FALSE
|
|
||||||
}
|
|
||||||
|
|
||||||
required_args <- c("annotation", "genome", "out_dir")
|
|
||||||
missing_args <- required_args[vapply(required_args, function(arg_name) {
|
|
||||||
value <- argv[[arg_name]]
|
|
||||||
arg_missing(value)
|
|
||||||
}, logical(1))]
|
|
||||||
if (length(missing_args) > 0) {
|
|
||||||
stop(sprintf(
|
|
||||||
"Missing required arguments: %s",
|
|
||||||
paste(sprintf("--%s", missing_args), collapse = ", ")
|
|
||||||
))
|
|
||||||
}
|
|
||||||
|
|
||||||
if (arg_missing(argv$bam_dir) == arg_missing(argv$bam_path)) {
|
|
||||||
stop("Provide exactly one of --bam_dir or --bam_path.")
|
|
||||||
}
|
|
||||||
|
|
||||||
dir.create(argv$out_dir, showWarnings = FALSE, recursive = TRUE)
|
|
||||||
|
|
||||||
sample_df <- NULL
|
|
||||||
if (!arg_missing(argv$sample_sheet)) {
|
|
||||||
sample_df <- read.csv(argv$sample_sheet, check.names = FALSE, stringsAsFactors = FALSE)
|
|
||||||
if (!"alias" %in% names(sample_df)) {
|
|
||||||
stop("Sample sheet must contain an 'alias' column.")
|
|
||||||
}
|
|
||||||
}
|
|
||||||
|
|
||||||
strip_alias <- function(path) {
|
|
||||||
name <- basename(path)
|
|
||||||
name <- sub("\\.aligned\\.sorted\\.bam$", "", name)
|
|
||||||
name <- tools::file_path_sans_ext(name)
|
|
||||||
name
|
|
||||||
}
|
|
||||||
|
|
||||||
if (!arg_missing(argv$bam_dir)) {
|
|
||||||
bam_paths <- sort(list.files(argv$bam_dir, pattern = "\\.bam$", full.names = TRUE))
|
|
||||||
if (length(bam_paths) < 1) {
|
|
||||||
stop("No BAM files were found in bam_dir.")
|
|
||||||
}
|
|
||||||
aliases <- vapply(bam_paths, strip_alias, character(1))
|
|
||||||
} else {
|
|
||||||
bam_paths <- argv$bam_path
|
|
||||||
aliases <- if (!arg_missing(argv$sample_alias)) argv$sample_alias else strip_alias(argv$bam_path)
|
|
||||||
}
|
|
||||||
|
|
||||||
if (!is.null(sample_df)) {
|
|
||||||
missing_aliases <- setdiff(aliases, sample_df$alias)
|
|
||||||
if (length(missing_aliases) > 0) {
|
|
||||||
stop(sprintf(
|
|
||||||
"Sample sheet is missing alias rows for BAM files: %s",
|
|
||||||
paste(missing_aliases, collapse = ", ")
|
|
||||||
))
|
|
||||||
}
|
|
||||||
sample_df <- sample_df[match(aliases, sample_df$alias), , drop = FALSE]
|
|
||||||
} else {
|
|
||||||
sample_df <- data.frame(alias = aliases, stringsAsFactors = FALSE)
|
|
||||||
}
|
|
||||||
|
|
||||||
annotation_obj <- prepareAnnotations(argv$annotation)
|
|
||||||
reads <- if (length(bam_paths) == 1) bam_paths else BamFileList(bam_paths, yieldSize = 1000000)
|
|
||||||
|
|
||||||
# Handle NDR parameter with validation and documentation
|
|
||||||
default_ndr <- 0.1
|
|
||||||
ndr_value <- default_ndr
|
|
||||||
|
|
||||||
if (!arg_missing(argv$ndr)) {
|
|
||||||
if (argv$ndr < 0 || argv$ndr > 1) {
|
|
||||||
stop("NDR (Novel Discovery Rate) must be between 0 and 1")
|
|
||||||
}
|
|
||||||
ndr_value <- argv$ndr
|
|
||||||
message(sprintf("Using user-specified NDR = %.3f", ndr_value))
|
|
||||||
} else {
|
|
||||||
message(sprintf("Using default NDR = %.3f", default_ndr))
|
|
||||||
}
|
|
||||||
|
|
||||||
if (identical(argv$transcriptome_mode, "discover")) {
|
|
||||||
message("Novel Discovery Rate (NDR) controls transcript discovery stringency:")
|
|
||||||
message(" Lower NDR (e.g., 0.05) = fewer false positive transcripts, may miss real ones")
|
|
||||||
message(" Higher NDR (e.g., 0.2) = more sensitive discovery, more false positives")
|
|
||||||
message(sprintf(" Current NDR = %.3f balances precision and recall", ndr_value))
|
|
||||||
}
|
|
||||||
|
|
||||||
bambu_args <- list(
|
|
||||||
reads = reads,
|
|
||||||
annotations = annotation_obj,
|
|
||||||
genome = argv$genome,
|
|
||||||
ncore = as.integer(argv$threads),
|
|
||||||
discovery = identical(argv$transcriptome_mode, "discover")
|
|
||||||
)
|
|
||||||
|
|
||||||
if (identical(argv$transcriptome_mode, "discover")) {
|
|
||||||
bambu_args$NDR <- ndr_value
|
|
||||||
}
|
|
||||||
|
|
||||||
message("Running bambu...")
|
|
||||||
se <- do.call(bambu, bambu_args)
|
|
||||||
message("Bambu completed successfully")
|
|
||||||
colnames(se) <- aliases
|
|
||||||
|
|
||||||
counts_mat <- assays(se)$counts
|
|
||||||
full_length_mat <- assays(se)$fullLengthCounts
|
|
||||||
|
|
||||||
# Collect QC statistics before filtering
|
|
||||||
qc_stats <- list()
|
|
||||||
qc_stats$total_transcripts_before_filter <- nrow(se)
|
|
||||||
qc_stats$total_genes_before_filter <- length(unique(rowData(se)$GENEID))
|
|
||||||
qc_stats$samples <- ncol(se)
|
|
||||||
|
|
||||||
# Filter low-count transcripts
|
|
||||||
if (is.null(full_length_mat)) {
|
|
||||||
keep_idx <- rowSums(counts_mat) > 0
|
|
||||||
} else {
|
|
||||||
keep_idx <- rowSums(full_length_mat) > 0
|
|
||||||
}
|
|
||||||
if (!any(keep_idx)) {
|
|
||||||
keep_idx <- rowSums(counts_mat) >= 0
|
|
||||||
}
|
|
||||||
|
|
||||||
qc_stats$transcripts_filtered <- sum(!keep_idx)
|
|
||||||
message(sprintf("Filtering: keeping %d / %d transcripts", sum(keep_idx), length(keep_idx)))
|
|
||||||
|
|
||||||
se <- se[keep_idx, ]
|
|
||||||
|
|
||||||
# Library size statistics and warnings
|
|
||||||
lib_sizes <- colSums(assays(se)$counts)
|
|
||||||
qc_stats$library_sizes <- as.list(lib_sizes)
|
|
||||||
qc_stats$min_library_size <- min(lib_sizes)
|
|
||||||
qc_stats$max_library_size <- max(lib_sizes)
|
|
||||||
qc_stats$median_library_size <- median(lib_sizes)
|
|
||||||
|
|
||||||
if (length(lib_sizes) > 1) {
|
|
||||||
lib_size_ratio <- max(lib_sizes) / min(lib_sizes)
|
|
||||||
qc_stats$library_size_ratio <- lib_size_ratio
|
|
||||||
|
|
||||||
if (lib_size_ratio > 3) {
|
|
||||||
warning(sprintf(
|
|
||||||
"Large library size variation detected (%.1fx difference).\n Min: %d, Max: %d reads.\n CPM normalization may not be appropriate for such variation.",
|
|
||||||
lib_size_ratio, min(lib_sizes), max(lib_sizes)
|
|
||||||
))
|
|
||||||
qc_stats$library_size_warning <- sprintf("%.1fx variation (>3x threshold)", lib_size_ratio)
|
|
||||||
}
|
|
||||||
}
|
|
||||||
|
|
||||||
# Per-sample detection statistics
|
|
||||||
qc_stats$transcripts_detected_per_sample <- as.list(colSums(assays(se)$counts > 0))
|
|
||||||
qc_stats$median_transcripts_detected <- median(colSums(assays(se)$counts > 0))
|
|
||||||
|
|
||||||
qc_stats$total_transcripts_after_filter <- nrow(se)
|
|
||||||
qc_stats$total_genes_after_filter <- length(unique(rowData(se)$GENEID))
|
|
||||||
|
|
||||||
row_ranges <- rowRanges(se)
|
|
||||||
writeToGTF(row_ranges, file = file.path(argv$out_dir, "transcripts.gtf"))
|
|
||||||
|
|
||||||
gene_se <- transcriptToGeneExpression(se)
|
|
||||||
colnames(gene_se) <- aliases
|
|
||||||
|
|
||||||
saveRDS(se, file.path(argv$out_dir, "bambu_transcripts.rds"))
|
|
||||||
saveRDS(gene_se, file.path(argv$out_dir, "bambu_genes.rds"))
|
|
||||||
write.csv(sample_df, file.path(argv$out_dir, "samples.csv"), row.names = FALSE, quote = FALSE)
|
|
||||||
|
|
||||||
tx_meta <- as.data.frame(rowData(se))
|
|
||||||
if (!"TXNAME" %in% names(tx_meta)) {
|
|
||||||
tx_meta$TXNAME <- rownames(se)
|
|
||||||
}
|
|
||||||
if (!"GENEID" %in% names(tx_meta)) {
|
|
||||||
tx_meta$GENEID <- NA_character_
|
|
||||||
}
|
|
||||||
gene_meta <- as.data.frame(rowData(gene_se))
|
|
||||||
if (!"GENEID" %in% names(gene_meta)) {
|
|
||||||
gene_meta$GENEID <- rownames(gene_se)
|
|
||||||
}
|
|
||||||
|
|
||||||
matrix_to_df <- function(se_obj, assay_name, id_col, meta_df) {
|
|
||||||
assay_df <- as.data.frame(assays(se_obj)[[assay_name]])
|
|
||||||
assay_df[[id_col]] <- rownames(se_obj)
|
|
||||||
assay_df <- assay_df[, c(id_col, setdiff(names(assay_df), id_col)), drop = FALSE]
|
|
||||||
merge(meta_df, assay_df, by.x = id_col, by.y = id_col, all.y = TRUE, sort = FALSE)
|
|
||||||
}
|
|
||||||
|
|
||||||
normalise_tsv_value <- function(value) {
|
|
||||||
if (length(value) == 0 || all(is.na(value))) {
|
|
||||||
return(NA_character_)
|
|
||||||
}
|
|
||||||
if (is.list(value)) {
|
|
||||||
value <- unlist(value, recursive = TRUE, use.names = FALSE)
|
|
||||||
}
|
|
||||||
if (length(value) == 0 || all(is.na(value))) {
|
|
||||||
return(NA_character_)
|
|
||||||
}
|
|
||||||
paste(as.character(value), collapse = ";")
|
|
||||||
}
|
|
||||||
|
|
||||||
normalise_tsv_df <- function(df) {
|
|
||||||
as.data.frame(
|
|
||||||
lapply(df, function(column) {
|
|
||||||
if (is.list(column)) {
|
|
||||||
vapply(column, normalise_tsv_value, character(1))
|
|
||||||
} else {
|
|
||||||
column
|
|
||||||
}
|
|
||||||
}),
|
|
||||||
stringsAsFactors = FALSE,
|
|
||||||
check.names = FALSE
|
|
||||||
)
|
|
||||||
}
|
|
||||||
|
|
||||||
extract_gtf_attribute <- function(attr_field, key) {
|
|
||||||
match <- regexec(sprintf('%s "([^"]*)";', key), attr_field, perl = TRUE)
|
|
||||||
captures <- regmatches(attr_field, match)[[1]]
|
|
||||||
if (length(captures) < 2) {
|
|
||||||
return(NULL)
|
|
||||||
}
|
|
||||||
captures[2]
|
|
||||||
}
|
|
||||||
|
|
||||||
normalise_gtf_attribute_value <- function(value) {
|
|
||||||
if (is.null(value)) {
|
|
||||||
return(NULL)
|
|
||||||
}
|
|
||||||
value <- gsub('[";]', "", value)
|
|
||||||
value <- trimws(gsub("\\s+", " ", value))
|
|
||||||
if (!nzchar(value)) {
|
|
||||||
return(NULL)
|
|
||||||
}
|
|
||||||
value
|
|
||||||
}
|
|
||||||
|
|
||||||
sanitise_gtf_file <- function(path) {
|
|
||||||
lines <- readLines(path, warn = FALSE)
|
|
||||||
cleaned_lines <- vapply(lines, function(line) {
|
|
||||||
if (!nzchar(line) || startsWith(line, "#")) {
|
|
||||||
return(line)
|
|
||||||
}
|
|
||||||
|
|
||||||
fields <- strsplit(line, "\t", fixed = TRUE)[[1]]
|
|
||||||
if (length(fields) < 9) {
|
|
||||||
return(line)
|
|
||||||
}
|
|
||||||
|
|
||||||
attr_field <- fields[9]
|
|
||||||
transcript_id <- normalise_gtf_attribute_value(
|
|
||||||
extract_gtf_attribute(attr_field, "transcript_id")
|
|
||||||
)
|
|
||||||
gene_id <- normalise_gtf_attribute_value(
|
|
||||||
extract_gtf_attribute(attr_field, "gene_id")
|
|
||||||
)
|
|
||||||
|
|
||||||
if (!is.null(gene_id) && grepl("\\btranscript_id\\b", gene_id)) {
|
|
||||||
gene_id <- transcript_id
|
|
||||||
}
|
|
||||||
if (is.null(gene_id)) {
|
|
||||||
gene_id <- transcript_id
|
|
||||||
}
|
|
||||||
|
|
||||||
if (!is.null(gene_id)) {
|
|
||||||
attr_field <- sub(
|
|
||||||
'gene_id "([^"]*)";',
|
|
||||||
sprintf('gene_id "%s";', gene_id),
|
|
||||||
attr_field,
|
|
||||||
perl = TRUE
|
|
||||||
)
|
|
||||||
}
|
|
||||||
if (!is.null(transcript_id)) {
|
|
||||||
attr_field <- sub(
|
|
||||||
'transcript_id "([^"]*)";',
|
|
||||||
sprintf('transcript_id "%s";', transcript_id),
|
|
||||||
attr_field,
|
|
||||||
perl = TRUE
|
|
||||||
)
|
|
||||||
}
|
|
||||||
|
|
||||||
fields[9] <- attr_field
|
|
||||||
paste(fields, collapse = "\t")
|
|
||||||
}, character(1))
|
|
||||||
|
|
||||||
writeLines(cleaned_lines, path)
|
|
||||||
}
|
|
||||||
|
|
||||||
tx_meta <- normalise_tsv_df(tx_meta)
|
|
||||||
gene_meta <- normalise_tsv_df(gene_meta)
|
|
||||||
|
|
||||||
write.table(
|
|
||||||
tx_meta,
|
|
||||||
file = file.path(argv$out_dir, "transcript_metadata.tsv"),
|
|
||||||
sep = "\t",
|
|
||||||
quote = FALSE,
|
|
||||||
row.names = FALSE
|
|
||||||
)
|
|
||||||
write.table(
|
|
||||||
gene_meta,
|
|
||||||
file = file.path(argv$out_dir, "gene_metadata.tsv"),
|
|
||||||
sep = "\t",
|
|
||||||
quote = FALSE,
|
|
||||||
row.names = FALSE
|
|
||||||
)
|
|
||||||
|
|
||||||
tx_counts <- matrix_to_df(se, "counts", "TXNAME", tx_meta)
|
|
||||||
tx_cpm <- matrix_to_df(se, "CPM", "TXNAME", tx_meta)
|
|
||||||
gene_counts <- matrix_to_df(gene_se, "counts", "GENEID", gene_meta)
|
|
||||||
gene_cpm <- matrix_to_df(gene_se, "CPM", "GENEID", gene_meta)
|
|
||||||
|
|
||||||
write.table(
|
|
||||||
tx_counts,
|
|
||||||
file = file.path(argv$out_dir, "transcript_counts.tsv"),
|
|
||||||
sep = "\t",
|
|
||||||
quote = FALSE,
|
|
||||||
row.names = FALSE
|
|
||||||
)
|
|
||||||
write.table(
|
|
||||||
tx_cpm,
|
|
||||||
file = file.path(argv$out_dir, "transcript_cpm.tsv"),
|
|
||||||
sep = "\t",
|
|
||||||
quote = FALSE,
|
|
||||||
row.names = FALSE
|
|
||||||
)
|
|
||||||
write.table(
|
|
||||||
gene_counts,
|
|
||||||
file = file.path(argv$out_dir, "gene_counts.tsv"),
|
|
||||||
sep = "\t",
|
|
||||||
quote = FALSE,
|
|
||||||
row.names = FALSE
|
|
||||||
)
|
|
||||||
write.table(
|
|
||||||
gene_cpm,
|
|
||||||
file = file.path(argv$out_dir, "gene_cpm.tsv"),
|
|
||||||
sep = "\t",
|
|
||||||
quote = FALSE,
|
|
||||||
row.names = FALSE
|
|
||||||
)
|
|
||||||
|
|
||||||
sanitise_gtf_file(file.path(argv$out_dir, "transcripts.gtf"))
|
|
||||||
|
|
||||||
# Write QC statistics as JSON for HTML report
|
|
||||||
qc_stats$transcriptome_mode <- argv$transcriptome_mode
|
|
||||||
qc_stats$ndr_used <- if (identical(argv$transcriptome_mode, "discover")) ndr_value else "N/A"
|
|
||||||
qc_stats$timestamp <- format(Sys.time(), "%Y-%m-%d %H:%M:%S")
|
|
||||||
|
|
||||||
write_json(
|
|
||||||
qc_stats,
|
|
||||||
file.path(argv$out_dir, "bambu_qc_stats.json"),
|
|
||||||
pretty = TRUE,
|
|
||||||
auto_unbox = TRUE
|
|
||||||
)
|
|
||||||
|
|
||||||
format_count <- function(value) {
|
|
||||||
if (length(value) == 0 || all(is.na(value))) {
|
|
||||||
return("NA")
|
|
||||||
}
|
|
||||||
format(
|
|
||||||
round(as.numeric(value), 0),
|
|
||||||
scientific = FALSE,
|
|
||||||
trim = TRUE,
|
|
||||||
big.mark = ","
|
|
||||||
)
|
|
||||||
}
|
|
||||||
|
|
||||||
# Write human-readable QC summary
|
|
||||||
qc_summary <- c(
|
|
||||||
"Bambu Quantification QC Summary",
|
|
||||||
"================================",
|
|
||||||
"",
|
|
||||||
sprintf("Timestamp: %s", qc_stats$timestamp),
|
|
||||||
sprintf("Mode: %s", argv$transcriptome_mode),
|
|
||||||
if (identical(argv$transcriptome_mode, "discover")) sprintf("NDR: %.3f", ndr_value) else NULL,
|
|
||||||
"",
|
|
||||||
"Sample Statistics:",
|
|
||||||
sprintf(" Samples analyzed: %s", format_count(qc_stats$samples)),
|
|
||||||
sprintf(
|
|
||||||
" Median library size: %s reads",
|
|
||||||
format_count(qc_stats$median_library_size)
|
|
||||||
),
|
|
||||||
sprintf(
|
|
||||||
" Library size range: %s - %s reads",
|
|
||||||
format_count(qc_stats$min_library_size),
|
|
||||||
format_count(qc_stats$max_library_size)
|
|
||||||
),
|
|
||||||
if (!is.null(qc_stats$library_size_warning)) sprintf(" WARNING: %s", qc_stats$library_size_warning) else NULL,
|
|
||||||
"",
|
|
||||||
"Transcript Discovery:",
|
|
||||||
sprintf(
|
|
||||||
" Transcripts before filtering: %s",
|
|
||||||
format_count(qc_stats$total_transcripts_before_filter)
|
|
||||||
),
|
|
||||||
sprintf(
|
|
||||||
" Transcripts after filtering: %s",
|
|
||||||
format_count(qc_stats$total_transcripts_after_filter)
|
|
||||||
),
|
|
||||||
sprintf(
|
|
||||||
" Transcripts removed: %s",
|
|
||||||
format_count(qc_stats$transcripts_filtered)
|
|
||||||
),
|
|
||||||
sprintf(
|
|
||||||
" Median transcripts detected per sample: %s",
|
|
||||||
format_count(qc_stats$median_transcripts_detected)
|
|
||||||
),
|
|
||||||
"",
|
|
||||||
"Gene-Level Summary:",
|
|
||||||
sprintf(
|
|
||||||
" Unique genes (before filter): %s",
|
|
||||||
format_count(qc_stats$total_genes_before_filter)
|
|
||||||
),
|
|
||||||
sprintf(
|
|
||||||
" Unique genes (after filter): %s",
|
|
||||||
format_count(qc_stats$total_genes_after_filter)
|
|
||||||
),
|
|
||||||
""
|
|
||||||
)
|
|
||||||
|
|
||||||
writeLines(qc_summary, file.path(argv$out_dir, "bambu_qc_summary.txt"))
|
|
||||||
message("QC statistics written to bambu_qc_stats.json and bambu_qc_summary.txt")
|
|
||||||
|
|
||||||
# Save session info for reproducibility
|
|
||||||
writeLines(capture.output(sessionInfo()), file.path(argv$out_dir, "session_info.txt"))
|
|
||||||
message("Session info saved for reproducibility")
|
|
||||||
@ -1,717 +0,0 @@
|
|||||||
#!/usr/bin/env Rscript
|
|
||||||
|
|
||||||
# Set seed for reproducibility
|
|
||||||
set.seed(42)
|
|
||||||
|
|
||||||
suppressPackageStartupMessages({
|
|
||||||
library(argparser)
|
|
||||||
library(DESeq2)
|
|
||||||
library(DEXSeq)
|
|
||||||
library(SummarizedExperiment)
|
|
||||||
library(jsonlite)
|
|
||||||
})
|
|
||||||
|
|
||||||
parser <- arg_parser("Run DESeq2 and DEXSeq on bambu output.")
|
|
||||||
parser <- add_argument(parser, "--transcript_rds", help = "bambu transcript RDS.")
|
|
||||||
parser <- add_argument(parser, "--gene_rds", help = "bambu gene RDS.")
|
|
||||||
parser <- add_argument(parser, "--sample_sheet", help = "Sample sheet CSV.")
|
|
||||||
parser <- add_argument(parser, "--condition_column", help = "Primary condition column.", default = "condition")
|
|
||||||
parser <- add_argument(parser, "--covariates", help = "Comma-separated nuisance covariates.")
|
|
||||||
parser <- add_argument(parser, "--reference_level", help = "Reference level for the condition column.")
|
|
||||||
parser <- add_argument(parser, "--out_dir", help = "Output directory.", default = "de_analysis")
|
|
||||||
argv <- parse_args(parser)
|
|
||||||
|
|
||||||
arg_missing <- function(value) {
|
|
||||||
if (is.null(value) || length(value) == 0 || all(is.na(value))) {
|
|
||||||
return(TRUE)
|
|
||||||
}
|
|
||||||
if (is.character(value)) {
|
|
||||||
return(all(!nzchar(value)))
|
|
||||||
}
|
|
||||||
FALSE
|
|
||||||
}
|
|
||||||
|
|
||||||
required_args <- c("transcript_rds", "gene_rds", "sample_sheet")
|
|
||||||
missing_args <- required_args[vapply(required_args, function(arg_name) {
|
|
||||||
value <- argv[[arg_name]]
|
|
||||||
arg_missing(value)
|
|
||||||
}, logical(1))]
|
|
||||||
if (length(missing_args) > 0) {
|
|
||||||
stop(sprintf(
|
|
||||||
"Missing required arguments: %s",
|
|
||||||
paste(sprintf("--%s", missing_args), collapse = ", ")
|
|
||||||
))
|
|
||||||
}
|
|
||||||
|
|
||||||
dir.create(argv$out_dir, showWarnings = FALSE, recursive = TRUE)
|
|
||||||
|
|
||||||
tx_se <- readRDS(argv$transcript_rds)
|
|
||||||
gene_se <- readRDS(argv$gene_rds)
|
|
||||||
sample_df <- read.csv(argv$sample_sheet, check.names = FALSE, stringsAsFactors = FALSE)
|
|
||||||
|
|
||||||
if (!"alias" %in% names(sample_df)) {
|
|
||||||
stop("Sample sheet must contain an 'alias' column.")
|
|
||||||
}
|
|
||||||
if (!(argv$condition_column %in% names(sample_df))) {
|
|
||||||
stop(sprintf("Sample sheet must contain the '%s' column.", argv$condition_column))
|
|
||||||
}
|
|
||||||
|
|
||||||
covariates <- character(0)
|
|
||||||
if (!arg_missing(argv$covariates)) {
|
|
||||||
covariates <- trimws(strsplit(argv$covariates, ",", fixed = TRUE)[[1]])
|
|
||||||
covariates <- covariates[nzchar(covariates)]
|
|
||||||
}
|
|
||||||
missing_covariates <- setdiff(covariates, names(sample_df))
|
|
||||||
if (length(missing_covariates) > 0) {
|
|
||||||
stop(sprintf("Missing covariate columns: %s", paste(missing_covariates, collapse = ", ")))
|
|
||||||
}
|
|
||||||
|
|
||||||
sample_df <- sample_df[match(colnames(tx_se), sample_df$alias), , drop = FALSE]
|
|
||||||
if (any(is.na(sample_df$alias))) {
|
|
||||||
stop("Sample sheet aliases do not match the bambu output sample names.")
|
|
||||||
}
|
|
||||||
|
|
||||||
condition_values <- unique(sample_df[[argv$condition_column]])
|
|
||||||
if (length(condition_values) < 2) {
|
|
||||||
stop("Differential analysis requires at least two condition levels.")
|
|
||||||
}
|
|
||||||
|
|
||||||
reference_level <- argv$reference_level
|
|
||||||
if (arg_missing(reference_level)) {
|
|
||||||
if ("control" %in% condition_values) {
|
|
||||||
reference_level <- "control"
|
|
||||||
} else {
|
|
||||||
stop("Provide --reference_level when the condition column does not contain 'control'.")
|
|
||||||
}
|
|
||||||
}
|
|
||||||
if (!(reference_level %in% condition_values)) {
|
|
||||||
stop("The requested reference level is not present in the condition column.")
|
|
||||||
}
|
|
||||||
|
|
||||||
sample_df[[argv$condition_column]] <- factor(sample_df[[argv$condition_column]])
|
|
||||||
for (covariate in covariates) {
|
|
||||||
sample_df[[covariate]] <- factor(sample_df[[covariate]])
|
|
||||||
}
|
|
||||||
|
|
||||||
run_deseq_with_fallback <- function(dds, contrast_name = "unknown") {
|
|
||||||
tryCatch(
|
|
||||||
DESeq(dds, quiet = TRUE),
|
|
||||||
error = function(err) {
|
|
||||||
if (!grepl(
|
|
||||||
"all gene-wise dispersion estimates are within 2 orders of magnitude",
|
|
||||||
conditionMessage(err),
|
|
||||||
fixed = TRUE
|
|
||||||
)) {
|
|
||||||
stop(err)
|
|
||||||
}
|
|
||||||
|
|
||||||
warning(
|
|
||||||
"STATISTICAL POWER REDUCED: DESeq2 dispersion estimation failed for ", contrast_name, ".\n",
|
|
||||||
"This usually indicates:\n",
|
|
||||||
" 1. Too few replicates (recommend n>=3 per group)\n",
|
|
||||||
" 2. High biological variability\n",
|
|
||||||
" 3. Poor data quality\n",
|
|
||||||
"Falling back to gene-wise dispersion (no information sharing).\n",
|
|
||||||
"Results will have reduced power and wider confidence intervals."
|
|
||||||
)
|
|
||||||
|
|
||||||
dds <- estimateSizeFactors(dds)
|
|
||||||
dds <- estimateDispersionsGeneEst(dds)
|
|
||||||
dispersions(dds) <- mcols(dds)$dispGeneEst
|
|
||||||
|
|
||||||
# Write diagnostic file
|
|
||||||
diag_content <- c(
|
|
||||||
"DESeq2 Dispersion Estimation Fallback Applied",
|
|
||||||
"==============================================",
|
|
||||||
"",
|
|
||||||
sprintf("Timestamp: %s", format(Sys.time(), "%Y-%m-%d %H:%M:%S")),
|
|
||||||
sprintf("Contrast: %s", contrast_name),
|
|
||||||
sprintf("Samples: %d", ncol(dds)),
|
|
||||||
sprintf("Genes tested: %d", nrow(dds)),
|
|
||||||
sprintf("Dispersion range: %.3f to %.3f", min(dispersions(dds)), max(dispersions(dds))),
|
|
||||||
"",
|
|
||||||
"WHAT HAPPENED:",
|
|
||||||
" Curve fitting failed. Using gene-wise dispersion estimates.",
|
|
||||||
"",
|
|
||||||
"IMPLICATIONS:",
|
|
||||||
" - No information sharing across genes",
|
|
||||||
" - Reduced statistical power",
|
|
||||||
" - Wider confidence intervals",
|
|
||||||
" - More conservative results (fewer discoveries)",
|
|
||||||
"",
|
|
||||||
"LIKELY CAUSES:",
|
|
||||||
" 1. Too few replicates (recommend n>=3 per group)",
|
|
||||||
" 2. High biological variability",
|
|
||||||
" 3. Poor data quality or outlier samples",
|
|
||||||
"",
|
|
||||||
"RECOMMENDATIONS:",
|
|
||||||
" - Add more biological replicates if possible",
|
|
||||||
" - Check sample quality metrics",
|
|
||||||
" - Consider filtering low-count genes more stringently"
|
|
||||||
)
|
|
||||||
|
|
||||||
diag_file <- file.path(argv$out_dir, sprintf("DESeq2_dispersion_fallback_%s.txt", gsub("[^A-Za-z0-9_-]", "_", contrast_name)))
|
|
||||||
writeLines(diag_content, diag_file)
|
|
||||||
|
|
||||||
nbinomWaldTest(dds)
|
|
||||||
}
|
|
||||||
)
|
|
||||||
}
|
|
||||||
|
|
||||||
estimate_dispersions_with_fallback <- function(object, context_label, allow_gene_est = TRUE) {
|
|
||||||
tryCatch(
|
|
||||||
estimateDispersions(object),
|
|
||||||
error = function(err) {
|
|
||||||
if (!grepl(
|
|
||||||
"all gene-wise dispersion estimates are within 2 orders of magnitude",
|
|
||||||
conditionMessage(err),
|
|
||||||
fixed = TRUE
|
|
||||||
)) {
|
|
||||||
stop(err)
|
|
||||||
}
|
|
||||||
|
|
||||||
message(
|
|
||||||
context_label,
|
|
||||||
" dispersion fitting failed; ",
|
|
||||||
"retrying with fitType='local'."
|
|
||||||
)
|
|
||||||
tryCatch(
|
|
||||||
estimateDispersions(object, fitType = "local"),
|
|
||||||
error = function(local_err) {
|
|
||||||
if (!grepl(
|
|
||||||
"all gene-wise dispersion estimates are within 2 orders of magnitude",
|
|
||||||
conditionMessage(local_err),
|
|
||||||
fixed = TRUE
|
|
||||||
)) {
|
|
||||||
stop(local_err)
|
|
||||||
}
|
|
||||||
|
|
||||||
message(
|
|
||||||
context_label,
|
|
||||||
" local-fit dispersion retry failed; ",
|
|
||||||
"retrying with fitType='mean'."
|
|
||||||
)
|
|
||||||
tryCatch(
|
|
||||||
estimateDispersions(object, fitType = "mean"),
|
|
||||||
error = function(mean_err) {
|
|
||||||
if (!grepl(
|
|
||||||
"all gene-wise dispersion estimates are within 2 orders of magnitude",
|
|
||||||
conditionMessage(mean_err),
|
|
||||||
fixed = TRUE
|
|
||||||
)) {
|
|
||||||
stop(mean_err)
|
|
||||||
}
|
|
||||||
if (!allow_gene_est) {
|
|
||||||
stop(mean_err)
|
|
||||||
}
|
|
||||||
|
|
||||||
message(
|
|
||||||
context_label,
|
|
||||||
" mean-fit dispersion retry failed; ",
|
|
||||||
"falling back to gene-wise dispersion estimates."
|
|
||||||
)
|
|
||||||
object <- estimateDispersionsGeneEst(object)
|
|
||||||
dispersions(object) <- mcols(object)$dispGeneEst
|
|
||||||
object
|
|
||||||
}
|
|
||||||
)
|
|
||||||
}
|
|
||||||
)
|
|
||||||
}
|
|
||||||
)
|
|
||||||
}
|
|
||||||
|
|
||||||
normalise_tsv_value <- function(value) {
|
|
||||||
if (length(value) == 0 || all(is.na(value))) {
|
|
||||||
return(NA_character_)
|
|
||||||
}
|
|
||||||
if (is.list(value)) {
|
|
||||||
value <- unlist(value, recursive = TRUE, use.names = FALSE)
|
|
||||||
}
|
|
||||||
if (length(value) == 0 || all(is.na(value))) {
|
|
||||||
return(NA_character_)
|
|
||||||
}
|
|
||||||
paste(as.character(value), collapse = ";")
|
|
||||||
}
|
|
||||||
|
|
||||||
normalise_tsv_df <- function(df) {
|
|
||||||
as.data.frame(
|
|
||||||
lapply(df, function(column) {
|
|
||||||
if (is.list(column)) {
|
|
||||||
vapply(column, normalise_tsv_value, character(1))
|
|
||||||
} else {
|
|
||||||
column
|
|
||||||
}
|
|
||||||
}),
|
|
||||||
stringsAsFactors = FALSE,
|
|
||||||
check.names = FALSE
|
|
||||||
)
|
|
||||||
}
|
|
||||||
|
|
||||||
is_recoverable_dexseq_error <- function(message_text) {
|
|
||||||
grepl(
|
|
||||||
"all gene-wise dispersion estimates are within 2 orders of magnitude",
|
|
||||||
message_text,
|
|
||||||
fixed = TRUE
|
|
||||||
) || grepl(
|
|
||||||
"model matrix is not full rank",
|
|
||||||
message_text,
|
|
||||||
fixed = TRUE
|
|
||||||
) || grepl(
|
|
||||||
"replacement has 1 row, data has 0",
|
|
||||||
message_text,
|
|
||||||
fixed = TRUE
|
|
||||||
)
|
|
||||||
}
|
|
||||||
|
|
||||||
empty_tsv <- function(columns) {
|
|
||||||
out <- as.data.frame(matrix(nrow = 0, ncol = length(columns)))
|
|
||||||
names(out) <- columns
|
|
||||||
out
|
|
||||||
}
|
|
||||||
|
|
||||||
write_placeholder_pdf <- function(path, label) {
|
|
||||||
pdf(path)
|
|
||||||
plot.new()
|
|
||||||
text(0.5, 0.5, label, cex = 0.9)
|
|
||||||
dev.off()
|
|
||||||
}
|
|
||||||
|
|
||||||
run_deseq2 <- function(count_mat, coldata, target_level, contrast_name) {
|
|
||||||
design_terms <- c(covariates, argv$condition_column)
|
|
||||||
design_formula <- as.formula(paste("~", paste(design_terms, collapse = " + ")))
|
|
||||||
dds <- DESeqDataSetFromMatrix(
|
|
||||||
countData = round(count_mat),
|
|
||||||
colData = coldata,
|
|
||||||
design = design_formula
|
|
||||||
)
|
|
||||||
dds <- run_deseq_with_fallback(dds, contrast_name)
|
|
||||||
results(dds, contrast = c(argv$condition_column, target_level, reference_level), independentFiltering = TRUE)
|
|
||||||
}
|
|
||||||
|
|
||||||
run_dexseq <- function(tx_counts, tx_meta, coldata, active_covariates = covariates) {
|
|
||||||
coldata$sample <- factor(coldata$alias)
|
|
||||||
coldata[[argv$condition_column]] <- factor(coldata[[argv$condition_column]])
|
|
||||||
for (covariate in active_covariates) {
|
|
||||||
coldata[[covariate]] <- factor(coldata[[covariate]])
|
|
||||||
}
|
|
||||||
|
|
||||||
covariate_exon_terms <- if (length(active_covariates) > 0) {
|
|
||||||
paste0(active_covariates, ":exon")
|
|
||||||
} else {
|
|
||||||
character(0)
|
|
||||||
}
|
|
||||||
design_terms <- c("sample", "exon", covariate_exon_terms, paste0(argv$condition_column, ":exon"))
|
|
||||||
reduced_terms <- c("sample", "exon", covariate_exon_terms)
|
|
||||||
full_formula <- as.formula(paste("~", paste(design_terms, collapse = " + ")))
|
|
||||||
reduced_formula <- as.formula(paste("~", paste(reduced_terms, collapse = " + ")))
|
|
||||||
|
|
||||||
tryCatch({
|
|
||||||
dxd <- DEXSeqDataSet(
|
|
||||||
countData = round(tx_counts),
|
|
||||||
sampleData = as.data.frame(coldata),
|
|
||||||
design = full_formula,
|
|
||||||
featureID = tx_meta$TXNAME,
|
|
||||||
groupID = tx_meta$GENEID
|
|
||||||
)
|
|
||||||
dxd <- estimateSizeFactors(dxd)
|
|
||||||
dxd <- estimate_dispersions_with_fallback(dxd, "DEXSeq", allow_gene_est = TRUE)
|
|
||||||
dxd <- testForDEU(dxd, reducedModel = reduced_formula)
|
|
||||||
dxd <- estimateExonFoldChanges(dxd, fitExpToVar = argv$condition_column)
|
|
||||||
dxr <- DEXSeqResults(dxd, independentFiltering = FALSE)
|
|
||||||
list(dxd = dxd, dxr = dxr)
|
|
||||||
}, error = function(err) {
|
|
||||||
if (length(active_covariates) == 0 || !grepl(
|
|
||||||
"model matrix is not full rank",
|
|
||||||
conditionMessage(err),
|
|
||||||
fixed = TRUE
|
|
||||||
)) {
|
|
||||||
stop(err)
|
|
||||||
}
|
|
||||||
|
|
||||||
dropped_covariate <- tail(active_covariates, 1)
|
|
||||||
kept_covariates <- head(active_covariates, -1)
|
|
||||||
message(
|
|
||||||
"DEXSeq design was not full rank with covariate '",
|
|
||||||
dropped_covariate,
|
|
||||||
"'; retrying without it."
|
|
||||||
)
|
|
||||||
run_dexseq(tx_counts, tx_meta, coldata, kept_covariates)
|
|
||||||
})
|
|
||||||
}
|
|
||||||
|
|
||||||
tx_meta <- as.data.frame(rowData(tx_se))
|
|
||||||
if (!"TXNAME" %in% names(tx_meta)) {
|
|
||||||
tx_meta$TXNAME <- rownames(tx_se)
|
|
||||||
}
|
|
||||||
if (!"GENEID" %in% names(tx_meta)) {
|
|
||||||
stop("Transcript rowData must contain GENEID for DEXSeq.")
|
|
||||||
}
|
|
||||||
|
|
||||||
gene_meta <- as.data.frame(rowData(gene_se))
|
|
||||||
if (!"GENEID" %in% names(gene_meta)) {
|
|
||||||
gene_meta$GENEID <- rownames(gene_se)
|
|
||||||
}
|
|
||||||
|
|
||||||
targets <- setdiff(as.character(condition_values), reference_level)
|
|
||||||
|
|
||||||
# Initialize QC statistics collector
|
|
||||||
de_qc_stats <- list()
|
|
||||||
de_qc_stats$timestamp <- format(Sys.time(), "%Y-%m-%d %H:%M:%S")
|
|
||||||
de_qc_stats$total_samples <- nrow(sample_df)
|
|
||||||
de_qc_stats$condition_column <- argv$condition_column
|
|
||||||
de_qc_stats$reference_level <- reference_level
|
|
||||||
de_qc_stats$covariates <- if (length(covariates) > 0) covariates else "none"
|
|
||||||
de_qc_stats$num_contrasts <- length(targets)
|
|
||||||
de_qc_stats$contrasts <- list()
|
|
||||||
|
|
||||||
# Check sample sizes and warn if underpowered
|
|
||||||
n_per_group <- table(sample_df[[argv$condition_column]])
|
|
||||||
de_qc_stats$samples_per_group <- as.list(n_per_group)
|
|
||||||
|
|
||||||
sample_size_warnings <- c()
|
|
||||||
if (any(n_per_group < 3)) {
|
|
||||||
warning(
|
|
||||||
"WARNING: Some condition groups have fewer than 3 replicates.\n",
|
|
||||||
"Recommended minimum for DGE: n=3 per group\n",
|
|
||||||
"Current sample sizes: ", paste(names(n_per_group), "=", n_per_group, collapse=", "), "\n",
|
|
||||||
"Results may have reduced statistical power."
|
|
||||||
)
|
|
||||||
sample_size_warnings <- c(sample_size_warnings, "Some groups have n<3 (recommended minimum)")
|
|
||||||
}
|
|
||||||
|
|
||||||
if (any(n_per_group < 2)) {
|
|
||||||
stop("ERROR: Some condition groups have fewer than 2 replicates. Cannot perform statistical testing.")
|
|
||||||
}
|
|
||||||
|
|
||||||
de_qc_stats$sample_size_warnings <- if (length(sample_size_warnings) > 0) sample_size_warnings else "none"
|
|
||||||
|
|
||||||
# Multiple testing warning
|
|
||||||
if (length(targets) > 1) {
|
|
||||||
fwer <- (1 - (1-0.05)^length(targets)) * 100
|
|
||||||
mt_warning <- sprintf(
|
|
||||||
"Multiple contrasts tested (%d). Per-contrast FDR < 0.05 yields family-wise error rate of ~%.1f%%",
|
|
||||||
length(targets), fwer
|
|
||||||
)
|
|
||||||
message("WARNING: ", mt_warning)
|
|
||||||
de_qc_stats$multiple_testing_note <- mt_warning
|
|
||||||
|
|
||||||
mt_content <- c(
|
|
||||||
"Multiple Testing Across Contrasts",
|
|
||||||
"==================================",
|
|
||||||
"",
|
|
||||||
sprintf("Timestamp: %s", format(Sys.time(), "%Y-%m-%d %H:%M:%S")),
|
|
||||||
sprintf("Number of contrasts tested: %d", length(targets)),
|
|
||||||
sprintf("Contrasts: %s", paste(sprintf("%s vs %s", targets, reference_level), collapse=", ")),
|
|
||||||
"",
|
|
||||||
"PER-CONTRAST FDR THRESHOLD: 0.05",
|
|
||||||
sprintf("FAMILY-WISE ERROR RATE: ~%.1f%%", fwer),
|
|
||||||
"",
|
|
||||||
"WHAT THIS MEANS:",
|
|
||||||
" Each contrast uses FDR < 0.05 independently.",
|
|
||||||
" When testing multiple contrasts, the overall false positive rate increases.",
|
|
||||||
sprintf(" Expected: %.1f%% chance of at least one false positive across all contrasts", fwer),
|
|
||||||
"",
|
|
||||||
"RECOMMENDATIONS:",
|
|
||||||
" 1. Use stricter per-contrast threshold:",
|
|
||||||
sprintf(" Bonferroni correction: 0.05 / %d = %.4f", length(targets), 0.05/length(targets)),
|
|
||||||
" 2. Focus on pre-specified contrasts of interest",
|
|
||||||
" 3. Treat results as exploratory and validate key findings",
|
|
||||||
" 4. Consider using hierarchical testing procedures",
|
|
||||||
"",
|
|
||||||
"INTERPRETATION:",
|
|
||||||
" - Results passing FDR < 0.05 in each contrast are discoveries for that contrast",
|
|
||||||
" - But the overall false discovery burden is higher than 5%",
|
|
||||||
" - Prioritize genes significant across multiple contrasts",
|
|
||||||
" - Validate top findings experimentally"
|
|
||||||
)
|
|
||||||
writeLines(mt_content, file.path(argv$out_dir, "MULTIPLE_TESTING_WARNING.txt"))
|
|
||||||
}
|
|
||||||
|
|
||||||
for (target_level in targets) {
|
|
||||||
contrast_name <- sprintf("%s_%s_vs_%s", argv$condition_column, target_level, reference_level)
|
|
||||||
contrast_dir <- file.path(argv$out_dir, contrast_name)
|
|
||||||
dir.create(contrast_dir, showWarnings = FALSE, recursive = TRUE)
|
|
||||||
|
|
||||||
keep_samples <- sample_df[[argv$condition_column]] %in% c(reference_level, target_level)
|
|
||||||
contrast_samples <- droplevels(sample_df[keep_samples, , drop = FALSE])
|
|
||||||
contrast_samples[[argv$condition_column]] <- relevel(
|
|
||||||
factor(contrast_samples[[argv$condition_column]]),
|
|
||||||
ref = reference_level
|
|
||||||
)
|
|
||||||
|
|
||||||
# Collect per-contrast QC stats
|
|
||||||
contrast_qc <- list()
|
|
||||||
contrast_qc$name <- contrast_name
|
|
||||||
contrast_qc$target_level <- target_level
|
|
||||||
contrast_qc$reference_level <- reference_level
|
|
||||||
contrast_qc$n_samples <- nrow(contrast_samples)
|
|
||||||
contrast_qc$n_target <- sum(contrast_samples[[argv$condition_column]] == target_level)
|
|
||||||
contrast_qc$n_reference <- sum(contrast_samples[[argv$condition_column]] == reference_level)
|
|
||||||
|
|
||||||
# DTU power warning
|
|
||||||
if (nrow(contrast_samples) < 6) {
|
|
||||||
dtu_warning <- sprintf(
|
|
||||||
"DTU analysis may be underpowered (n=%d, recommend n>=6 with >=3 per group)",
|
|
||||||
nrow(contrast_samples)
|
|
||||||
)
|
|
||||||
warning(dtu_warning)
|
|
||||||
contrast_qc$dtu_power_warning <- dtu_warning
|
|
||||||
}
|
|
||||||
|
|
||||||
gene_counts <- assays(gene_se)$counts[, contrast_samples$alias, drop = FALSE]
|
|
||||||
tx_counts <- assays(tx_se)$counts[, contrast_samples$alias, drop = FALSE]
|
|
||||||
|
|
||||||
contrast_qc$genes_tested <- nrow(gene_counts)
|
|
||||||
contrast_qc$transcripts_tested <- nrow(tx_counts)
|
|
||||||
|
|
||||||
dge_res <- as.data.frame(run_deseq2(gene_counts, contrast_samples, target_level, contrast_name))
|
|
||||||
dge_res$GENEID <- rownames(dge_res)
|
|
||||||
dge_res <- merge(gene_meta, dge_res, by = "GENEID", all.y = TRUE, sort = FALSE)
|
|
||||||
dge_res <- normalise_tsv_df(dge_res)
|
|
||||||
|
|
||||||
# Collect DGE statistics
|
|
||||||
contrast_qc$dge_total_genes <- nrow(dge_res)
|
|
||||||
contrast_qc$dge_significant_fdr05 <- sum(dge_res$padj < 0.05, na.rm = TRUE)
|
|
||||||
contrast_qc$dge_significant_fdr01 <- sum(dge_res$padj < 0.01, na.rm = TRUE)
|
|
||||||
contrast_qc$dge_upregulated <- sum(dge_res$padj < 0.05 & dge_res$log2FoldChange > 0, na.rm = TRUE)
|
|
||||||
contrast_qc$dge_downregulated <- sum(dge_res$padj < 0.05 & dge_res$log2FoldChange < 0, na.rm = TRUE)
|
|
||||||
|
|
||||||
write.table(
|
|
||||||
dge_res[order(dge_res$padj), ],
|
|
||||||
file = file.path(contrast_dir, "results_dge.tsv"),
|
|
||||||
sep = "\t",
|
|
||||||
quote = FALSE,
|
|
||||||
row.names = FALSE
|
|
||||||
)
|
|
||||||
|
|
||||||
pdf(file.path(contrast_dir, "results_dge.pdf"))
|
|
||||||
dds_plot <- DESeqDataSetFromMatrix(
|
|
||||||
countData = round(gene_counts),
|
|
||||||
colData = contrast_samples,
|
|
||||||
design = as.formula(paste("~", paste(c(covariates, argv$condition_column), collapse = " + ")))
|
|
||||||
)
|
|
||||||
dds_plot <- run_deseq_with_fallback(dds_plot, contrast_name)
|
|
||||||
plotMA(results(dds_plot, contrast = c(argv$condition_column, target_level, reference_level), independentFiltering = TRUE))
|
|
||||||
dev.off()
|
|
||||||
|
|
||||||
dex_res <- tryCatch(
|
|
||||||
run_dexseq(tx_counts, tx_meta, contrast_samples),
|
|
||||||
error = function(err) {
|
|
||||||
message_text <- conditionMessage(err)
|
|
||||||
if (!is_recoverable_dexseq_error(message_text)) {
|
|
||||||
stop(err)
|
|
||||||
}
|
|
||||||
|
|
||||||
warning(
|
|
||||||
"DEXSeq failed for contrast ", target_level, " vs ", reference_level, "\n",
|
|
||||||
"Error: ", message_text
|
|
||||||
)
|
|
||||||
|
|
||||||
# Write explicit failure report
|
|
||||||
failure_content <- c(
|
|
||||||
"DTU Analysis Failed",
|
|
||||||
"===================",
|
|
||||||
"",
|
|
||||||
sprintf("Timestamp: %s", format(Sys.time(), "%Y-%m-%d %H:%M:%S")),
|
|
||||||
sprintf("Contrast: %s vs %s", target_level, reference_level),
|
|
||||||
sprintf("Samples: %d (%d %s, %d %s)",
|
|
||||||
nrow(contrast_samples),
|
|
||||||
sum(contrast_samples[[argv$condition_column]] == target_level), target_level,
|
|
||||||
sum(contrast_samples[[argv$condition_column]] == reference_level), reference_level),
|
|
||||||
sprintf("Transcripts: %d", nrow(tx_counts)),
|
|
||||||
"",
|
|
||||||
"ERROR MESSAGE:",
|
|
||||||
sprintf(" %s", message_text),
|
|
||||||
"",
|
|
||||||
"DTU RESULTS CANNOT BE INTERPRETED",
|
|
||||||
"",
|
|
||||||
"This failure is likely due to:",
|
|
||||||
" 1. Insufficient samples (need >=3 per group, recommend >=6 total for DTU)",
|
|
||||||
" 2. Too few transcripts with sufficient counts",
|
|
||||||
" 3. Design matrix not full rank (covariate confounding)",
|
|
||||||
" 4. Extreme count distributions",
|
|
||||||
"",
|
|
||||||
"RECOMMENDATIONS:",
|
|
||||||
" - Use gene-level DGE results (less power required)",
|
|
||||||
" - Add more biological replicates",
|
|
||||||
" - Filter transcripts more stringently",
|
|
||||||
" - Simplify experimental design (remove problematic covariates)",
|
|
||||||
"",
|
|
||||||
"NOTE: Empty DTU result files indicate analysis failure, not 'no DTU detected'"
|
|
||||||
)
|
|
||||||
writeLines(failure_content, file.path(contrast_dir, "DTU_ANALYSIS_FAILED.txt"))
|
|
||||||
|
|
||||||
NULL
|
|
||||||
}
|
|
||||||
)
|
|
||||||
|
|
||||||
if (is.null(dex_res)) {
|
|
||||||
dex_df <- empty_tsv(c(
|
|
||||||
"featureID",
|
|
||||||
"groupID",
|
|
||||||
"log2fold",
|
|
||||||
"pvalue",
|
|
||||||
"padj",
|
|
||||||
"exonBaseMean"
|
|
||||||
))
|
|
||||||
tx_dtu <- dex_df
|
|
||||||
gene_dtu <- empty_tsv(c("GENEID", "qval"))
|
|
||||||
write_placeholder_pdf(
|
|
||||||
file.path(contrast_dir, "results_dtu.pdf"),
|
|
||||||
"DEXSeq did not converge for this contrast.\nSee DTU_ANALYSIS_FAILED.txt for details."
|
|
||||||
)
|
|
||||||
contrast_qc$dtu_status <- "FAILED"
|
|
||||||
contrast_qc$dtu_significant_transcripts <- 0
|
|
||||||
contrast_qc$dtu_significant_genes <- 0
|
|
||||||
} else {
|
|
||||||
dxr <- dex_res$dxr
|
|
||||||
dxd <- dex_res$dxd
|
|
||||||
dex_df <- as.data.frame(dxr)
|
|
||||||
dex_df <- normalise_tsv_df(dex_df)
|
|
||||||
tx_dtu <- dex_df[, intersect(
|
|
||||||
c("featureID", "groupID", "log2fold", "pvalue", "padj", "exonBaseMean"),
|
|
||||||
names(dex_df)
|
|
||||||
), drop = FALSE]
|
|
||||||
tx_dtu <- normalise_tsv_df(tx_dtu)
|
|
||||||
|
|
||||||
gene_q <- perGeneQValue(dxr)
|
|
||||||
gene_dtu <- data.frame(
|
|
||||||
GENEID = names(gene_q),
|
|
||||||
qval = unname(gene_q),
|
|
||||||
row.names = NULL
|
|
||||||
)
|
|
||||||
|
|
||||||
# Collect DTU statistics
|
|
||||||
contrast_qc$dtu_status <- "SUCCESS"
|
|
||||||
contrast_qc$dtu_significant_transcripts <- sum(tx_dtu$padj < 0.05, na.rm = TRUE)
|
|
||||||
contrast_qc$dtu_significant_genes <- sum(gene_dtu$qval < 0.05, na.rm = TRUE)
|
|
||||||
|
|
||||||
pdf(file.path(contrast_dir, "results_dtu.pdf"))
|
|
||||||
plotMA(dxr, cex = 0.8, alpha = 0.05)
|
|
||||||
plotDispEsts(dxd)
|
|
||||||
dev.off()
|
|
||||||
}
|
|
||||||
|
|
||||||
write.table(
|
|
||||||
dex_df,
|
|
||||||
file = file.path(contrast_dir, "results_dexseq.tsv"),
|
|
||||||
sep = "\t",
|
|
||||||
quote = FALSE,
|
|
||||||
row.names = FALSE
|
|
||||||
)
|
|
||||||
|
|
||||||
write.table(
|
|
||||||
tx_dtu[order(tx_dtu$padj), ],
|
|
||||||
file = file.path(contrast_dir, "results_dtu_transcript.tsv"),
|
|
||||||
sep = "\t",
|
|
||||||
quote = FALSE,
|
|
||||||
row.names = FALSE
|
|
||||||
)
|
|
||||||
|
|
||||||
write.table(
|
|
||||||
gene_dtu[order(gene_dtu$qval), ],
|
|
||||||
file = file.path(contrast_dir, "results_dtu_gene.tsv"),
|
|
||||||
sep = "\t",
|
|
||||||
quote = FALSE,
|
|
||||||
row.names = FALSE
|
|
||||||
)
|
|
||||||
|
|
||||||
write.table(
|
|
||||||
contrast_samples,
|
|
||||||
file = file.path(contrast_dir, "samples_used.tsv"),
|
|
||||||
sep = "\t",
|
|
||||||
quote = FALSE,
|
|
||||||
row.names = FALSE
|
|
||||||
)
|
|
||||||
|
|
||||||
# Write per-contrast QC summary
|
|
||||||
contrast_qc_summary <- c(
|
|
||||||
sprintf("Contrast QC Summary: %s", contrast_name),
|
|
||||||
paste(rep("=", 50), collapse = ""),
|
|
||||||
"",
|
|
||||||
"Sample Information:",
|
|
||||||
sprintf(" Target level (%s): %d samples", target_level, contrast_qc$n_target),
|
|
||||||
sprintf(" Reference level (%s): %d samples", reference_level, contrast_qc$n_reference),
|
|
||||||
sprintf(" Total samples: %d", contrast_qc$n_samples),
|
|
||||||
"",
|
|
||||||
"DGE Results:",
|
|
||||||
sprintf(" Genes tested: %d", contrast_qc$genes_tested),
|
|
||||||
sprintf(" Significant (FDR < 0.05): %d", contrast_qc$dge_significant_fdr05),
|
|
||||||
sprintf(" Significant (FDR < 0.01): %d", contrast_qc$dge_significant_fdr01),
|
|
||||||
sprintf(" Upregulated: %d", contrast_qc$dge_upregulated),
|
|
||||||
sprintf(" Downregulated: %d", contrast_qc$dge_downregulated),
|
|
||||||
"",
|
|
||||||
"DTU Results:",
|
|
||||||
sprintf(" Status: %s", contrast_qc$dtu_status),
|
|
||||||
sprintf(" Transcripts tested: %d", contrast_qc$transcripts_tested),
|
|
||||||
if (contrast_qc$dtu_status == "SUCCESS") {
|
|
||||||
c(
|
|
||||||
sprintf(" Significant transcripts (FDR < 0.05): %d", contrast_qc$dtu_significant_transcripts),
|
|
||||||
sprintf(" Genes with DTU (q < 0.05): %d", contrast_qc$dtu_significant_genes)
|
|
||||||
)
|
|
||||||
} else {
|
|
||||||
" See DTU_ANALYSIS_FAILED.txt for details"
|
|
||||||
},
|
|
||||||
if (!is.null(contrast_qc$dtu_power_warning)) paste0(" WARNING: ", contrast_qc$dtu_power_warning) else NULL,
|
|
||||||
""
|
|
||||||
)
|
|
||||||
writeLines(contrast_qc_summary, file.path(contrast_dir, "contrast_qc_summary.txt"))
|
|
||||||
|
|
||||||
# Add to overall QC stats
|
|
||||||
de_qc_stats$contrasts[[contrast_name]] <- contrast_qc
|
|
||||||
}
|
|
||||||
|
|
||||||
# Write overall DE/DTU QC statistics as JSON for HTML report
|
|
||||||
write_json(
|
|
||||||
de_qc_stats,
|
|
||||||
file.path(argv$out_dir, "de_qc_stats.json"),
|
|
||||||
pretty = TRUE,
|
|
||||||
auto_unbox = TRUE
|
|
||||||
)
|
|
||||||
|
|
||||||
# Write human-readable overall summary
|
|
||||||
overall_summary <- c(
|
|
||||||
"Differential Expression/Usage Analysis Summary",
|
|
||||||
paste(rep("=", 50), collapse = ""),
|
|
||||||
"",
|
|
||||||
sprintf("Timestamp: %s", de_qc_stats$timestamp),
|
|
||||||
sprintf("Total samples: %d", de_qc_stats$total_samples),
|
|
||||||
sprintf("Condition column: %s", de_qc_stats$condition_column),
|
|
||||||
sprintf("Reference level: %s", de_qc_stats$reference_level),
|
|
||||||
sprintf("Covariates: %s", paste(de_qc_stats$covariates, collapse = ", ")),
|
|
||||||
"",
|
|
||||||
"Sample Sizes:",
|
|
||||||
sapply(names(de_qc_stats$samples_per_group), function(grp) {
|
|
||||||
sprintf(" %s: %d samples", grp, de_qc_stats$samples_per_group[[grp]])
|
|
||||||
}),
|
|
||||||
if (de_qc_stats$sample_size_warnings != "none") paste0(" WARNING: ", de_qc_stats$sample_size_warnings) else NULL,
|
|
||||||
"",
|
|
||||||
sprintf("Number of contrasts tested: %d", de_qc_stats$num_contrasts),
|
|
||||||
if (!is.null(de_qc_stats$multiple_testing_note)) paste0(" NOTE: ", de_qc_stats$multiple_testing_note) else NULL,
|
|
||||||
"",
|
|
||||||
"Per-Contrast Results:",
|
|
||||||
sapply(names(de_qc_stats$contrasts), function(cname) {
|
|
||||||
cqc <- de_qc_stats$contrasts[[cname]]
|
|
||||||
c(
|
|
||||||
"",
|
|
||||||
sprintf(" %s:", cname),
|
|
||||||
sprintf(" Samples: %d (%d vs %d)", cqc$n_samples, cqc$n_target, cqc$n_reference),
|
|
||||||
sprintf(" DGE significant: %d genes (FDR<0.05)", cqc$dge_significant_fdr05),
|
|
||||||
sprintf(" DTU status: %s", cqc$dtu_status),
|
|
||||||
if (cqc$dtu_status == "SUCCESS") sprintf(" DTU significant: %d genes", cqc$dtu_significant_genes) else NULL
|
|
||||||
)
|
|
||||||
}),
|
|
||||||
"",
|
|
||||||
"For detailed per-contrast statistics, see:",
|
|
||||||
" - <contrast>/contrast_qc_summary.txt",
|
|
||||||
" - <contrast>/results_dge.tsv",
|
|
||||||
" - <contrast>/results_dtu_gene.tsv",
|
|
||||||
""
|
|
||||||
)
|
|
||||||
writeLines(unlist(overall_summary), file.path(argv$out_dir, "de_overall_summary.txt"))
|
|
||||||
|
|
||||||
# Save session info for reproducibility
|
|
||||||
writeLines(capture.output(sessionInfo()), file.path(argv$out_dir, "session_info.txt"))
|
|
||||||
message("QC statistics written to de_qc_stats.json and de_overall_summary.txt")
|
|
||||||
message("Session info saved for reproducibility")
|
|
||||||
24
bin/supeRglue
Executable file
24
bin/supeRglue
Executable file
@ -0,0 +1,24 @@
|
|||||||
|
#!/usr/bin/env Rscript
|
||||||
|
|
||||||
|
supeRglue_script_dir <- function() {
|
||||||
|
file_arg <- grep("^--file=", commandArgs(trailingOnly = FALSE), value = TRUE)
|
||||||
|
if (length(file_arg) > 0) {
|
||||||
|
return(dirname(normalizePath(sub("^--file=", "", file_arg[[1]]), mustWork = TRUE)))
|
||||||
|
}
|
||||||
|
|
||||||
|
frame_files <- vapply(sys.frames(), function(frame) {
|
||||||
|
path <- frame$ofile
|
||||||
|
if (is.null(path)) "" else path
|
||||||
|
}, character(1))
|
||||||
|
frame_files <- frame_files[nzchar(frame_files)]
|
||||||
|
if (length(frame_files) > 0) {
|
||||||
|
return(dirname(normalizePath(frame_files[[length(frame_files)]], mustWork = TRUE)))
|
||||||
|
}
|
||||||
|
|
||||||
|
stop("Unable to locate supeRglue on disk.", call. = FALSE)
|
||||||
|
}
|
||||||
|
|
||||||
|
script_dir <- supeRglue_script_dir()
|
||||||
|
source(file.path(script_dir, "workflow_glue_r", "load.R"))
|
||||||
|
workflow_glue_r_load(file.path(script_dir, "workflow_glue_r"))
|
||||||
|
workflow_glue_r_cli()
|
||||||
20
bin/workflow_glue_r/DESCRIPTION
Normal file
20
bin/workflow_glue_r/DESCRIPTION
Normal file
@ -0,0 +1,20 @@
|
|||||||
|
Package: workflowGlueR
|
||||||
|
Title: Local R Helpers for wf-transcriptomes
|
||||||
|
Version: 0.0.1.0000
|
||||||
|
Authors@R:
|
||||||
|
person("Oxford Nanopore Technologies", role = c("aut", "cre"))
|
||||||
|
Description: Local R helper functions used by wf-transcriptomes.
|
||||||
|
License: Proprietary
|
||||||
|
Encoding: UTF-8
|
||||||
|
LazyData: false
|
||||||
|
Imports:
|
||||||
|
argparser,
|
||||||
|
bambu,
|
||||||
|
DESeq2,
|
||||||
|
DEXSeq,
|
||||||
|
jsonlite,
|
||||||
|
Rsamtools,
|
||||||
|
S4Vectors,
|
||||||
|
SummarizedExperiment
|
||||||
|
Suggests:
|
||||||
|
testthat
|
||||||
26
bin/workflow_glue_r/NAMESPACE
Normal file
26
bin/workflow_glue_r/NAMESPACE
Normal file
@ -0,0 +1,26 @@
|
|||||||
|
export(bambu_arg_parser)
|
||||||
|
export(bambu_build_args)
|
||||||
|
export(bambu_discovery_enabled)
|
||||||
|
export(bambu_filter_transcripts)
|
||||||
|
export(bambu_normalise_tsv_df)
|
||||||
|
export(bambu_resolve_inputs)
|
||||||
|
export(bambu_resolve_ndr)
|
||||||
|
export(bambu_strip_alias)
|
||||||
|
export(bambu_validate_args)
|
||||||
|
export(bambu_write_outputs)
|
||||||
|
export(de_analysis_arg_parser)
|
||||||
|
export(de_build_contrast_name)
|
||||||
|
export(de_parse_covariates)
|
||||||
|
export(de_validate_inputs)
|
||||||
|
export(main_run_bambu)
|
||||||
|
export(main_run_de_analysis)
|
||||||
|
export(run_bambu_cli)
|
||||||
|
export(run_de_analysis_cli)
|
||||||
|
export(workflow_glue_r_cli)
|
||||||
|
export(workflow_glue_r_components)
|
||||||
|
export(workflow_glue_r_arg_missing)
|
||||||
|
export(workflow_glue_r_empty_tsv)
|
||||||
|
export(workflow_glue_r_normalise_tsv_df)
|
||||||
|
export(workflow_glue_r_parse_csv_list)
|
||||||
|
export(workflow_glue_r_read_csv)
|
||||||
|
export(workflow_glue_r_require_args)
|
||||||
538
bin/workflow_glue_r/R/bambu.R
Normal file
538
bin/workflow_glue_r/R/bambu.R
Normal file
@ -0,0 +1,538 @@
|
|||||||
|
bambu_arg_parser <- function() {
|
||||||
|
parser <- argparser::arg_parser("Run bambu transcript discovery and quantification.")
|
||||||
|
parser <- argparser::add_argument(parser, "--bam_dir", help = "Directory containing BAM files.")
|
||||||
|
parser <- argparser::add_argument(parser, "--bam_path", help = "Path to a single BAM file.")
|
||||||
|
parser <- argparser::add_argument(parser, "--sample_alias", help = "Alias to use for a single BAM file.")
|
||||||
|
parser <- argparser::add_argument(parser, "--sample_sheet", help = "Optional sample sheet CSV.")
|
||||||
|
parser <- argparser::add_argument(parser, "--annotation", help = "Reference annotation GTF/GFF.")
|
||||||
|
parser <- argparser::add_argument(parser, "--genome", help = "Reference genome FASTA.")
|
||||||
|
parser <- argparser::add_argument(
|
||||||
|
parser,
|
||||||
|
"--transcriptome_mode",
|
||||||
|
help = "discover or fixed_annotation.",
|
||||||
|
default = "discover"
|
||||||
|
)
|
||||||
|
parser <- argparser::add_argument(
|
||||||
|
parser,
|
||||||
|
"--threads",
|
||||||
|
help = "Number of worker threads.",
|
||||||
|
type = "numeric",
|
||||||
|
default = 1
|
||||||
|
)
|
||||||
|
parser <- argparser::add_argument(
|
||||||
|
parser,
|
||||||
|
"--ndr",
|
||||||
|
help = "Optional novel discovery rate.",
|
||||||
|
type = "numeric"
|
||||||
|
)
|
||||||
|
argparser::add_argument(parser, "--out_dir", help = "Output directory.")
|
||||||
|
}
|
||||||
|
|
||||||
|
bambu_validate_args <- function(argv) {
|
||||||
|
workflow_glue_r_require_args(argv, c("annotation", "genome", "out_dir"))
|
||||||
|
|
||||||
|
if (workflow_glue_r_arg_missing(argv$bam_dir) == workflow_glue_r_arg_missing(argv$bam_path)) {
|
||||||
|
stop("Provide exactly one of --bam_dir or --bam_path.", call. = FALSE)
|
||||||
|
}
|
||||||
|
|
||||||
|
if (!argv$transcriptome_mode %in% c("discover", "fixed_annotation")) {
|
||||||
|
stop(
|
||||||
|
sprintf(
|
||||||
|
"transcriptome_mode must be one of: %s",
|
||||||
|
paste(c("discover", "fixed_annotation"), collapse = ", ")
|
||||||
|
),
|
||||||
|
call. = FALSE
|
||||||
|
)
|
||||||
|
}
|
||||||
|
|
||||||
|
if (!workflow_glue_r_arg_missing(argv$ndr) && (argv$ndr < 0 || argv$ndr > 1)) {
|
||||||
|
stop("NDR (Novel Discovery Rate) must be between 0 and 1", call. = FALSE)
|
||||||
|
}
|
||||||
|
|
||||||
|
invisible(argv)
|
||||||
|
}
|
||||||
|
|
||||||
|
bambu_strip_alias <- function(path) {
|
||||||
|
name <- basename(path)
|
||||||
|
name <- sub("\\.aligned\\.sorted\\.bam$", "", name)
|
||||||
|
tools::file_path_sans_ext(name)
|
||||||
|
}
|
||||||
|
|
||||||
|
bambu_resolve_inputs <- function(
|
||||||
|
argv,
|
||||||
|
bamfile_list_ctor = Rsamtools::BamFileList,
|
||||||
|
list_files_fn = base::list.files
|
||||||
|
) {
|
||||||
|
sample_df <- NULL
|
||||||
|
if (!workflow_glue_r_arg_missing(argv$sample_sheet)) {
|
||||||
|
sample_df <- workflow_glue_r_read_csv(argv$sample_sheet)
|
||||||
|
if (!"alias" %in% names(sample_df)) {
|
||||||
|
stop("Sample sheet must contain an 'alias' column.", call. = FALSE)
|
||||||
|
}
|
||||||
|
duplicate_sample_aliases <- unique(sample_df$alias[duplicated(sample_df$alias)])
|
||||||
|
if (length(duplicate_sample_aliases) > 0) {
|
||||||
|
stop(
|
||||||
|
sprintf(
|
||||||
|
"Sample sheet aliases must be unique; duplicated aliases: %s",
|
||||||
|
paste(duplicate_sample_aliases, collapse = ", ")
|
||||||
|
),
|
||||||
|
call. = FALSE
|
||||||
|
)
|
||||||
|
}
|
||||||
|
}
|
||||||
|
|
||||||
|
if (!workflow_glue_r_arg_missing(argv$bam_dir)) {
|
||||||
|
bam_paths <- sort(list_files_fn(argv$bam_dir, pattern = "\\.bam$", full.names = TRUE))
|
||||||
|
if (length(bam_paths) < 1) {
|
||||||
|
stop("No BAM files were found in bam_dir.", call. = FALSE)
|
||||||
|
}
|
||||||
|
aliases <- unname(vapply(bam_paths, bambu_strip_alias, character(1)))
|
||||||
|
} else {
|
||||||
|
bam_paths <- argv$bam_path
|
||||||
|
aliases <- if (!workflow_glue_r_arg_missing(argv$sample_alias)) {
|
||||||
|
argv$sample_alias
|
||||||
|
} else {
|
||||||
|
bambu_strip_alias(argv$bam_path)
|
||||||
|
}
|
||||||
|
}
|
||||||
|
|
||||||
|
duplicate_bam_aliases <- unique(aliases[duplicated(aliases)])
|
||||||
|
if (length(duplicate_bam_aliases) > 0) {
|
||||||
|
stop(
|
||||||
|
sprintf(
|
||||||
|
"BAM aliases must be unique; duplicated aliases: %s",
|
||||||
|
paste(duplicate_bam_aliases, collapse = ", ")
|
||||||
|
),
|
||||||
|
call. = FALSE
|
||||||
|
)
|
||||||
|
}
|
||||||
|
|
||||||
|
if (!is.null(sample_df)) {
|
||||||
|
missing_aliases <- setdiff(aliases, sample_df$alias)
|
||||||
|
if (length(missing_aliases) > 0) {
|
||||||
|
stop(
|
||||||
|
sprintf(
|
||||||
|
"Sample sheet is missing alias rows for BAM files: %s",
|
||||||
|
paste(missing_aliases, collapse = ", ")
|
||||||
|
),
|
||||||
|
call. = FALSE
|
||||||
|
)
|
||||||
|
}
|
||||||
|
sample_df <- sample_df[match(aliases, sample_df$alias), , drop = FALSE]
|
||||||
|
} else {
|
||||||
|
sample_df <- data.frame(alias = aliases, stringsAsFactors = FALSE)
|
||||||
|
}
|
||||||
|
|
||||||
|
reads <- if (length(bam_paths) == 1) {
|
||||||
|
bam_paths
|
||||||
|
} else {
|
||||||
|
bamfile_list_ctor(bam_paths, yieldSize = 1000000)
|
||||||
|
}
|
||||||
|
|
||||||
|
list(
|
||||||
|
bam_paths = bam_paths,
|
||||||
|
aliases = aliases,
|
||||||
|
sample_df = sample_df,
|
||||||
|
reads = reads
|
||||||
|
)
|
||||||
|
}
|
||||||
|
|
||||||
|
bambu_discovery_enabled <- function(argv) {
|
||||||
|
identical(argv$transcriptome_mode, "discover")
|
||||||
|
}
|
||||||
|
|
||||||
|
bambu_resolve_ndr <- function(argv, default_ndr = 0.1) {
|
||||||
|
if (workflow_glue_r_arg_missing(argv$ndr)) {
|
||||||
|
default_ndr
|
||||||
|
} else {
|
||||||
|
as.numeric(argv$ndr)
|
||||||
|
}
|
||||||
|
}
|
||||||
|
|
||||||
|
bambu_build_args <- function(argv, reads, annotation_obj) {
|
||||||
|
bambu_args <- list(
|
||||||
|
reads = reads,
|
||||||
|
annotations = annotation_obj,
|
||||||
|
genome = argv$genome,
|
||||||
|
ncore = as.integer(argv$threads),
|
||||||
|
discovery = bambu_discovery_enabled(argv)
|
||||||
|
)
|
||||||
|
|
||||||
|
if (bambu_discovery_enabled(argv)) {
|
||||||
|
bambu_args$NDR <- bambu_resolve_ndr(argv)
|
||||||
|
}
|
||||||
|
|
||||||
|
bambu_args
|
||||||
|
}
|
||||||
|
|
||||||
|
bambu_filter_transcripts <- function(se) {
|
||||||
|
counts_mat <- SummarizedExperiment::assays(se)$counts
|
||||||
|
full_length_mat <- SummarizedExperiment::assays(se)$fullLengthCounts
|
||||||
|
|
||||||
|
gene_ids <- SummarizedExperiment::rowData(se)$GENEID
|
||||||
|
qc_stats <- list(
|
||||||
|
total_transcripts_before_filter = nrow(se),
|
||||||
|
total_genes_before_filter = length(unique(gene_ids)),
|
||||||
|
samples = ncol(se)
|
||||||
|
)
|
||||||
|
|
||||||
|
if (is.null(full_length_mat)) {
|
||||||
|
keep_idx <- rowSums(counts_mat) > 0
|
||||||
|
} else {
|
||||||
|
keep_idx <- rowSums(full_length_mat) > 0
|
||||||
|
}
|
||||||
|
if (!any(keep_idx)) {
|
||||||
|
keep_idx <- rowSums(counts_mat) >= 0
|
||||||
|
}
|
||||||
|
|
||||||
|
qc_stats$transcripts_filtered <- sum(!keep_idx)
|
||||||
|
se <- se[keep_idx, ]
|
||||||
|
qc_stats$total_transcripts_after_filter <- nrow(se)
|
||||||
|
qc_stats$total_genes_after_filter <- length(unique(SummarizedExperiment::rowData(se)$GENEID))
|
||||||
|
|
||||||
|
list(se = se, qc_stats = qc_stats)
|
||||||
|
}
|
||||||
|
|
||||||
|
bambu_matrix_to_df <- function(se_obj, assay_name, id_col, meta_df) {
|
||||||
|
assay_df <- as.data.frame(SummarizedExperiment::assays(se_obj)[[assay_name]])
|
||||||
|
assay_df[[id_col]] <- rownames(se_obj)
|
||||||
|
assay_df <- assay_df[, c(id_col, setdiff(names(assay_df), id_col)), drop = FALSE]
|
||||||
|
merge(meta_df, assay_df, by.x = id_col, by.y = id_col, all.y = TRUE, sort = FALSE)
|
||||||
|
}
|
||||||
|
|
||||||
|
bambu_extract_gtf_attribute <- function(attr_field, key) {
|
||||||
|
match <- regexec(sprintf('%s "([^"]*)";', key), attr_field, perl = TRUE)
|
||||||
|
captures <- regmatches(attr_field, match)[[1]]
|
||||||
|
if (length(captures) < 2) {
|
||||||
|
return(NULL)
|
||||||
|
}
|
||||||
|
captures[2]
|
||||||
|
}
|
||||||
|
|
||||||
|
bambu_normalise_gtf_attribute_value <- function(value) {
|
||||||
|
if (is.null(value)) {
|
||||||
|
return(NULL)
|
||||||
|
}
|
||||||
|
value <- gsub('[";]', "", value)
|
||||||
|
value <- trimws(gsub("\\s+", " ", value))
|
||||||
|
if (!nzchar(value)) {
|
||||||
|
return(NULL)
|
||||||
|
}
|
||||||
|
value
|
||||||
|
}
|
||||||
|
|
||||||
|
bambu_sanitise_gtf_file <- function(path) {
|
||||||
|
lines <- readLines(path, warn = FALSE)
|
||||||
|
cleaned_lines <- vapply(lines, function(line) {
|
||||||
|
if (!nzchar(line) || startsWith(line, "#")) {
|
||||||
|
return(line)
|
||||||
|
}
|
||||||
|
|
||||||
|
fields <- strsplit(line, "\t", fixed = TRUE)[[1]]
|
||||||
|
if (length(fields) < 9) {
|
||||||
|
return(line)
|
||||||
|
}
|
||||||
|
|
||||||
|
attr_field <- fields[9]
|
||||||
|
transcript_id <- bambu_normalise_gtf_attribute_value(
|
||||||
|
bambu_extract_gtf_attribute(attr_field, "transcript_id")
|
||||||
|
)
|
||||||
|
gene_id <- bambu_normalise_gtf_attribute_value(
|
||||||
|
bambu_extract_gtf_attribute(attr_field, "gene_id")
|
||||||
|
)
|
||||||
|
|
||||||
|
if (!is.null(gene_id) && identical(gene_id, "transcript_id")) {
|
||||||
|
warning(
|
||||||
|
sprintf(
|
||||||
|
"Replaced malformed gene_id 'transcript_id' with transcript_id '%s'.",
|
||||||
|
transcript_id
|
||||||
|
),
|
||||||
|
call. = FALSE
|
||||||
|
)
|
||||||
|
gene_id <- transcript_id
|
||||||
|
}
|
||||||
|
if (is.null(gene_id)) {
|
||||||
|
gene_id <- transcript_id
|
||||||
|
}
|
||||||
|
|
||||||
|
if (!is.null(gene_id)) {
|
||||||
|
attr_field <- sub(
|
||||||
|
'gene_id "([^"]*)";',
|
||||||
|
sprintf('gene_id "%s";', gene_id),
|
||||||
|
attr_field,
|
||||||
|
perl = TRUE
|
||||||
|
)
|
||||||
|
}
|
||||||
|
if (!is.null(transcript_id)) {
|
||||||
|
attr_field <- sub(
|
||||||
|
'transcript_id "([^"]*)";',
|
||||||
|
sprintf('transcript_id "%s";', transcript_id),
|
||||||
|
attr_field,
|
||||||
|
perl = TRUE
|
||||||
|
)
|
||||||
|
}
|
||||||
|
|
||||||
|
fields[9] <- attr_field
|
||||||
|
paste(fields, collapse = "\t")
|
||||||
|
}, character(1))
|
||||||
|
|
||||||
|
writeLines(cleaned_lines, path)
|
||||||
|
}
|
||||||
|
|
||||||
|
bambu_format_count <- function(value) {
|
||||||
|
if (length(value) == 0 || all(is.na(value))) {
|
||||||
|
return("NA")
|
||||||
|
}
|
||||||
|
format(
|
||||||
|
round(as.numeric(value), 0),
|
||||||
|
scientific = FALSE,
|
||||||
|
trim = TRUE,
|
||||||
|
big.mark = ","
|
||||||
|
)
|
||||||
|
}
|
||||||
|
|
||||||
|
bambu_write_outputs <- function(se, gene_se, sample_df, argv, qc_stats, write_gtf_fn = bambu::writeToGTF) {
|
||||||
|
write_gtf_fn(
|
||||||
|
SummarizedExperiment::rowRanges(se),
|
||||||
|
file = file.path(argv$out_dir, "transcripts.gtf")
|
||||||
|
)
|
||||||
|
|
||||||
|
saveRDS(se, file.path(argv$out_dir, "bambu_transcripts.rds"))
|
||||||
|
saveRDS(gene_se, file.path(argv$out_dir, "bambu_genes.rds"))
|
||||||
|
utils::write.csv(
|
||||||
|
sample_df,
|
||||||
|
file.path(argv$out_dir, "samples.csv"),
|
||||||
|
row.names = FALSE,
|
||||||
|
quote = FALSE
|
||||||
|
)
|
||||||
|
|
||||||
|
tx_meta <- as.data.frame(SummarizedExperiment::rowData(se))
|
||||||
|
if (!"TXNAME" %in% names(tx_meta)) {
|
||||||
|
tx_meta$TXNAME <- rownames(se)
|
||||||
|
}
|
||||||
|
if (!"GENEID" %in% names(tx_meta)) {
|
||||||
|
tx_meta$GENEID <- NA_character_
|
||||||
|
}
|
||||||
|
|
||||||
|
gene_meta <- as.data.frame(SummarizedExperiment::rowData(gene_se))
|
||||||
|
if (!"GENEID" %in% names(gene_meta)) {
|
||||||
|
gene_meta$GENEID <- rownames(gene_se)
|
||||||
|
}
|
||||||
|
|
||||||
|
tx_meta <- workflow_glue_r_normalise_tsv_df(tx_meta)
|
||||||
|
gene_meta <- workflow_glue_r_normalise_tsv_df(gene_meta)
|
||||||
|
|
||||||
|
utils::write.table(
|
||||||
|
tx_meta,
|
||||||
|
file = file.path(argv$out_dir, "transcript_metadata.tsv"),
|
||||||
|
sep = "\t",
|
||||||
|
quote = FALSE,
|
||||||
|
row.names = FALSE
|
||||||
|
)
|
||||||
|
utils::write.table(
|
||||||
|
gene_meta,
|
||||||
|
file = file.path(argv$out_dir, "gene_metadata.tsv"),
|
||||||
|
sep = "\t",
|
||||||
|
quote = FALSE,
|
||||||
|
row.names = FALSE
|
||||||
|
)
|
||||||
|
|
||||||
|
tx_counts <- bambu_matrix_to_df(se, "counts", "TXNAME", tx_meta)
|
||||||
|
tx_cpm <- bambu_matrix_to_df(se, "CPM", "TXNAME", tx_meta)
|
||||||
|
gene_counts <- bambu_matrix_to_df(gene_se, "counts", "GENEID", gene_meta)
|
||||||
|
gene_cpm <- bambu_matrix_to_df(gene_se, "CPM", "GENEID", gene_meta)
|
||||||
|
|
||||||
|
utils::write.table(
|
||||||
|
tx_counts,
|
||||||
|
file = file.path(argv$out_dir, "transcript_counts.tsv"),
|
||||||
|
sep = "\t",
|
||||||
|
quote = FALSE,
|
||||||
|
row.names = FALSE
|
||||||
|
)
|
||||||
|
utils::write.table(
|
||||||
|
tx_cpm,
|
||||||
|
file = file.path(argv$out_dir, "transcript_cpm.tsv"),
|
||||||
|
sep = "\t",
|
||||||
|
quote = FALSE,
|
||||||
|
row.names = FALSE
|
||||||
|
)
|
||||||
|
utils::write.table(
|
||||||
|
gene_counts,
|
||||||
|
file = file.path(argv$out_dir, "gene_counts.tsv"),
|
||||||
|
sep = "\t",
|
||||||
|
quote = FALSE,
|
||||||
|
row.names = FALSE
|
||||||
|
)
|
||||||
|
utils::write.table(
|
||||||
|
gene_cpm,
|
||||||
|
file = file.path(argv$out_dir, "gene_cpm.tsv"),
|
||||||
|
sep = "\t",
|
||||||
|
quote = FALSE,
|
||||||
|
row.names = FALSE
|
||||||
|
)
|
||||||
|
|
||||||
|
bambu_sanitise_gtf_file(file.path(argv$out_dir, "transcripts.gtf"))
|
||||||
|
|
||||||
|
qc_stats$transcriptome_mode <- argv$transcriptome_mode
|
||||||
|
qc_stats$ndr_used <- if (bambu_discovery_enabled(argv)) {
|
||||||
|
bambu_resolve_ndr(argv)
|
||||||
|
} else {
|
||||||
|
"N/A"
|
||||||
|
}
|
||||||
|
qc_stats$timestamp <- format(Sys.time(), "%Y-%m-%d %H:%M:%S")
|
||||||
|
|
||||||
|
jsonlite::write_json(
|
||||||
|
qc_stats,
|
||||||
|
file.path(argv$out_dir, "bambu_qc_stats.json"),
|
||||||
|
pretty = TRUE,
|
||||||
|
auto_unbox = TRUE
|
||||||
|
)
|
||||||
|
|
||||||
|
qc_summary <- c(
|
||||||
|
"Bambu Quantification QC Summary",
|
||||||
|
"================================",
|
||||||
|
"",
|
||||||
|
sprintf("Timestamp: %s", qc_stats$timestamp),
|
||||||
|
sprintf("Mode: %s", argv$transcriptome_mode),
|
||||||
|
if (bambu_discovery_enabled(argv)) sprintf("NDR: %.3f", bambu_resolve_ndr(argv)) else NULL,
|
||||||
|
"",
|
||||||
|
"Sample Statistics:",
|
||||||
|
sprintf(" Samples analyzed: %s", bambu_format_count(qc_stats$samples)),
|
||||||
|
sprintf(" Median library size: %s reads", bambu_format_count(qc_stats$median_library_size)),
|
||||||
|
sprintf(
|
||||||
|
" Library size range: %s - %s reads",
|
||||||
|
bambu_format_count(qc_stats$min_library_size),
|
||||||
|
bambu_format_count(qc_stats$max_library_size)
|
||||||
|
),
|
||||||
|
if (!is.null(qc_stats$library_size_warning)) sprintf(" WARNING: %s", qc_stats$library_size_warning) else NULL,
|
||||||
|
"",
|
||||||
|
"Transcript Discovery:",
|
||||||
|
sprintf(" Transcripts before filtering: %s", bambu_format_count(qc_stats$total_transcripts_before_filter)),
|
||||||
|
sprintf(" Transcripts after filtering: %s", bambu_format_count(qc_stats$total_transcripts_after_filter)),
|
||||||
|
sprintf(" Transcripts removed: %s", bambu_format_count(qc_stats$transcripts_filtered)),
|
||||||
|
sprintf(
|
||||||
|
" Median transcripts detected per sample: %s",
|
||||||
|
bambu_format_count(qc_stats$median_transcripts_detected)
|
||||||
|
),
|
||||||
|
"",
|
||||||
|
"Gene-Level Summary:",
|
||||||
|
sprintf(" Unique genes (before filter): %s", bambu_format_count(qc_stats$total_genes_before_filter)),
|
||||||
|
sprintf(" Unique genes (after filter): %s", bambu_format_count(qc_stats$total_genes_after_filter)),
|
||||||
|
""
|
||||||
|
)
|
||||||
|
|
||||||
|
writeLines(qc_summary, file.path(argv$out_dir, "bambu_qc_summary.txt"))
|
||||||
|
writeLines(capture.output(sessionInfo()), file.path(argv$out_dir, "session_info.txt"))
|
||||||
|
}
|
||||||
|
|
||||||
|
main_run_bambu <- function(
|
||||||
|
argv,
|
||||||
|
analysis_fn = bambu::bambu,
|
||||||
|
prepare_annotations_fn = bambu::prepareAnnotations,
|
||||||
|
gene_expression_fn = bambu::transcriptToGeneExpression,
|
||||||
|
write_gtf_fn = bambu::writeToGTF,
|
||||||
|
bamfile_list_ctor = Rsamtools::BamFileList,
|
||||||
|
list_files_fn = base::list.files
|
||||||
|
) {
|
||||||
|
set.seed(42)
|
||||||
|
suppressPackageStartupMessages({
|
||||||
|
library(GenomicRanges)
|
||||||
|
library(Rsamtools)
|
||||||
|
})
|
||||||
|
|
||||||
|
bambu_validate_args(argv)
|
||||||
|
dir.create(argv$out_dir, showWarnings = FALSE, recursive = TRUE)
|
||||||
|
|
||||||
|
inputs <- bambu_resolve_inputs(
|
||||||
|
argv,
|
||||||
|
bamfile_list_ctor = bamfile_list_ctor,
|
||||||
|
list_files_fn = list_files_fn
|
||||||
|
)
|
||||||
|
annotation_obj <- prepare_annotations_fn(argv$annotation)
|
||||||
|
ndr_value <- bambu_resolve_ndr(argv)
|
||||||
|
|
||||||
|
if (!workflow_glue_r_arg_missing(argv$ndr)) {
|
||||||
|
message(sprintf("Using user-specified NDR = %.3f", ndr_value))
|
||||||
|
} else {
|
||||||
|
message(sprintf("Using default NDR = %.3f", ndr_value))
|
||||||
|
}
|
||||||
|
|
||||||
|
if (bambu_discovery_enabled(argv)) {
|
||||||
|
message("Novel Discovery Rate (NDR) controls transcript discovery stringency:")
|
||||||
|
message(" Lower NDR (e.g., 0.05) = fewer false positive transcripts, may miss real ones")
|
||||||
|
message(" Higher NDR (e.g., 0.2) = more sensitive discovery, more false positives")
|
||||||
|
message(sprintf(" Current NDR = %.3f balances precision and recall", ndr_value))
|
||||||
|
}
|
||||||
|
|
||||||
|
message("Running bambu...")
|
||||||
|
se <- do.call(analysis_fn, bambu_build_args(argv, inputs$reads, annotation_obj))
|
||||||
|
message("Bambu completed successfully")
|
||||||
|
colnames(se) <- inputs$aliases
|
||||||
|
|
||||||
|
filtered <- bambu_filter_transcripts(se)
|
||||||
|
se <- filtered$se
|
||||||
|
qc_stats <- filtered$qc_stats
|
||||||
|
message(
|
||||||
|
sprintf(
|
||||||
|
"Filtering: keeping %d / %d transcripts",
|
||||||
|
qc_stats$total_transcripts_after_filter,
|
||||||
|
qc_stats$total_transcripts_before_filter
|
||||||
|
)
|
||||||
|
)
|
||||||
|
|
||||||
|
lib_sizes <- colSums(SummarizedExperiment::assays(se)$counts)
|
||||||
|
qc_stats$library_sizes <- as.list(lib_sizes)
|
||||||
|
qc_stats$min_library_size <- min(lib_sizes)
|
||||||
|
qc_stats$max_library_size <- max(lib_sizes)
|
||||||
|
qc_stats$median_library_size <- stats::median(lib_sizes)
|
||||||
|
|
||||||
|
if (length(lib_sizes) > 1) {
|
||||||
|
lib_size_ratio <- max(lib_sizes) / min(lib_sizes)
|
||||||
|
qc_stats$library_size_ratio <- lib_size_ratio
|
||||||
|
if (lib_size_ratio > 3) {
|
||||||
|
warning(
|
||||||
|
sprintf(
|
||||||
|
paste0(
|
||||||
|
"Large library size variation detected (%.1fx difference).\n",
|
||||||
|
" Min: %d, Max: %d reads.\n",
|
||||||
|
" CPM normalization may not be appropriate for such variation."
|
||||||
|
),
|
||||||
|
lib_size_ratio,
|
||||||
|
min(lib_sizes),
|
||||||
|
max(lib_sizes)
|
||||||
|
)
|
||||||
|
)
|
||||||
|
qc_stats$library_size_warning <- sprintf("%.1fx variation (>3x threshold)", lib_size_ratio)
|
||||||
|
}
|
||||||
|
}
|
||||||
|
|
||||||
|
detected_per_sample <- colSums(SummarizedExperiment::assays(se)$counts > 0)
|
||||||
|
qc_stats$transcripts_detected_per_sample <- as.list(detected_per_sample)
|
||||||
|
qc_stats$median_transcripts_detected <- stats::median(detected_per_sample)
|
||||||
|
|
||||||
|
gene_se <- gene_expression_fn(se)
|
||||||
|
colnames(gene_se) <- inputs$aliases
|
||||||
|
|
||||||
|
bambu_write_outputs(
|
||||||
|
se,
|
||||||
|
gene_se,
|
||||||
|
inputs$sample_df,
|
||||||
|
argv,
|
||||||
|
qc_stats,
|
||||||
|
write_gtf_fn = write_gtf_fn
|
||||||
|
)
|
||||||
|
|
||||||
|
invisible(
|
||||||
|
list(
|
||||||
|
se = se,
|
||||||
|
gene_se = gene_se,
|
||||||
|
sample_df = inputs$sample_df,
|
||||||
|
qc_stats = qc_stats
|
||||||
|
)
|
||||||
|
)
|
||||||
|
}
|
||||||
|
|
||||||
|
run_bambu_cli <- function(argv = commandArgs(trailingOnly = TRUE)) {
|
||||||
|
parsed <- argparser::parse_args(bambu_arg_parser(), argv = argv)
|
||||||
|
main_run_bambu(parsed)
|
||||||
|
}
|
||||||
72
bin/workflow_glue_r/R/cli.R
Normal file
72
bin/workflow_glue_r/R/cli.R
Normal file
@ -0,0 +1,72 @@
|
|||||||
|
workflow_glue_r_components <- function(env = globalenv()) {
|
||||||
|
parser_suffix <- "_arg_parser"
|
||||||
|
parser_names <- grep(
|
||||||
|
paste0(parser_suffix, "$"),
|
||||||
|
ls(envir = env, all.names = TRUE),
|
||||||
|
value = TRUE
|
||||||
|
)
|
||||||
|
|
||||||
|
components <- list()
|
||||||
|
for (parser_name in parser_names) {
|
||||||
|
component <- sub(paste0(parser_suffix, "$"), "", parser_name)
|
||||||
|
cli_name <- paste0("run_", component, "_cli")
|
||||||
|
parser <- get(parser_name, envir = env)
|
||||||
|
runner <- if (exists(cli_name, envir = env, mode = "function")) {
|
||||||
|
get(cli_name, envir = env)
|
||||||
|
} else {
|
||||||
|
NULL
|
||||||
|
}
|
||||||
|
|
||||||
|
if (is.function(parser) && is.function(runner)) {
|
||||||
|
components[[component]] <- list(
|
||||||
|
name = component,
|
||||||
|
parser_name = parser_name,
|
||||||
|
runner_name = cli_name,
|
||||||
|
parser = parser,
|
||||||
|
runner = runner
|
||||||
|
)
|
||||||
|
}
|
||||||
|
}
|
||||||
|
|
||||||
|
components[sort(names(components))]
|
||||||
|
}
|
||||||
|
|
||||||
|
workflow_glue_r_usage <- function(components = workflow_glue_r_components()) {
|
||||||
|
component_names <- names(components)
|
||||||
|
lines <- c(
|
||||||
|
"Usage: supeRglue <command> [options]",
|
||||||
|
"",
|
||||||
|
"Commands:",
|
||||||
|
if (length(component_names) > 0) {
|
||||||
|
paste0(" ", component_names)
|
||||||
|
} else {
|
||||||
|
" <none found>"
|
||||||
|
},
|
||||||
|
"",
|
||||||
|
"Use 'supeRglue <command> --help' for command-specific options."
|
||||||
|
)
|
||||||
|
paste(lines, collapse = "\n")
|
||||||
|
}
|
||||||
|
|
||||||
|
workflow_glue_r_cli <- function(argv = commandArgs(trailingOnly = TRUE), env = globalenv()) {
|
||||||
|
components <- workflow_glue_r_components(env = env)
|
||||||
|
|
||||||
|
if (length(argv) < 1 || argv[[1]] %in% c("-h", "--help", "help")) {
|
||||||
|
cat(workflow_glue_r_usage(components), "\n")
|
||||||
|
return(invisible(0L))
|
||||||
|
}
|
||||||
|
|
||||||
|
command <- argv[[1]]
|
||||||
|
if (!command %in% names(components)) {
|
||||||
|
stop(
|
||||||
|
sprintf(
|
||||||
|
"Unknown supeRglue command '%s'. Available commands: %s",
|
||||||
|
command,
|
||||||
|
paste(names(components), collapse = ", ")
|
||||||
|
),
|
||||||
|
call. = FALSE
|
||||||
|
)
|
||||||
|
}
|
||||||
|
|
||||||
|
components[[command]]$runner(argv[-1])
|
||||||
|
}
|
||||||
94
bin/workflow_glue_r/R/common.R
Normal file
94
bin/workflow_glue_r/R/common.R
Normal file
@ -0,0 +1,94 @@
|
|||||||
|
workflow_glue_r_arg_missing <- function(value) {
|
||||||
|
if (is.null(value) || length(value) == 0 || all(is.na(value))) {
|
||||||
|
return(TRUE)
|
||||||
|
}
|
||||||
|
if (is.character(value)) {
|
||||||
|
return(all(!nzchar(value)))
|
||||||
|
}
|
||||||
|
FALSE
|
||||||
|
}
|
||||||
|
|
||||||
|
workflow_glue_r_require_args <- function(argv, required_args) {
|
||||||
|
missing_args <- required_args[vapply(required_args, function(arg_name) {
|
||||||
|
workflow_glue_r_arg_missing(argv[[arg_name]])
|
||||||
|
}, logical(1))]
|
||||||
|
|
||||||
|
if (length(missing_args) > 0) {
|
||||||
|
stop(
|
||||||
|
sprintf(
|
||||||
|
"Missing required arguments: %s",
|
||||||
|
paste(sprintf("--%s", missing_args), collapse = ", ")
|
||||||
|
),
|
||||||
|
call. = FALSE
|
||||||
|
)
|
||||||
|
}
|
||||||
|
}
|
||||||
|
|
||||||
|
workflow_glue_r_parse_csv_list <- function(value) {
|
||||||
|
if (workflow_glue_r_arg_missing(value)) {
|
||||||
|
return(character(0))
|
||||||
|
}
|
||||||
|
|
||||||
|
values <- trimws(strsplit(value, ",", fixed = TRUE)[[1]])
|
||||||
|
values[nzchar(values)]
|
||||||
|
}
|
||||||
|
|
||||||
|
workflow_glue_r_is_r_formula_name <- function(name) {
|
||||||
|
is.character(name) &&
|
||||||
|
length(name) == 1 &&
|
||||||
|
grepl("^[A-Za-z][A-Za-z0-9_.]*$", name)
|
||||||
|
}
|
||||||
|
|
||||||
|
workflow_glue_r_validate_r_formula_names <- function(names, label = "Column") {
|
||||||
|
invalid <- names[!vapply(names, workflow_glue_r_is_r_formula_name, logical(1))]
|
||||||
|
if (length(invalid) > 0) {
|
||||||
|
stop(
|
||||||
|
sprintf(
|
||||||
|
"%s names must be safe for R formulas. Invalid names: %s. Names must start with a letter and contain only letters, numbers, underscores, and dots.",
|
||||||
|
label,
|
||||||
|
paste(invalid, collapse = ", ")
|
||||||
|
),
|
||||||
|
call. = FALSE
|
||||||
|
)
|
||||||
|
}
|
||||||
|
invisible(names)
|
||||||
|
}
|
||||||
|
|
||||||
|
workflow_glue_r_read_csv <- function(path) {
|
||||||
|
utils::read.csv(path, check.names = FALSE, stringsAsFactors = FALSE)
|
||||||
|
}
|
||||||
|
|
||||||
|
workflow_glue_r_normalise_tsv_value <- function(value) {
|
||||||
|
if (length(value) == 0 || all(is.na(value))) {
|
||||||
|
return(NA_character_)
|
||||||
|
}
|
||||||
|
if (is.list(value)) {
|
||||||
|
value <- unlist(value, recursive = TRUE, use.names = FALSE)
|
||||||
|
}
|
||||||
|
if (length(value) == 0 || all(is.na(value))) {
|
||||||
|
return(NA_character_)
|
||||||
|
}
|
||||||
|
paste(as.character(value), collapse = ";")
|
||||||
|
}
|
||||||
|
|
||||||
|
workflow_glue_r_normalise_tsv_df <- function(df) {
|
||||||
|
as.data.frame(
|
||||||
|
lapply(df, function(column) {
|
||||||
|
if (is.list(column)) {
|
||||||
|
vapply(column, workflow_glue_r_normalise_tsv_value, character(1))
|
||||||
|
} else {
|
||||||
|
column
|
||||||
|
}
|
||||||
|
}),
|
||||||
|
stringsAsFactors = FALSE,
|
||||||
|
check.names = FALSE
|
||||||
|
)
|
||||||
|
}
|
||||||
|
|
||||||
|
bambu_normalise_tsv_df <- workflow_glue_r_normalise_tsv_df
|
||||||
|
|
||||||
|
workflow_glue_r_empty_tsv <- function(columns) {
|
||||||
|
out <- as.data.frame(matrix(nrow = 0, ncol = length(columns)))
|
||||||
|
names(out) <- columns
|
||||||
|
out
|
||||||
|
}
|
||||||
781
bin/workflow_glue_r/R/de_analysis.R
Normal file
781
bin/workflow_glue_r/R/de_analysis.R
Normal file
@ -0,0 +1,781 @@
|
|||||||
|
de_analysis_arg_parser <- function() {
|
||||||
|
parser <- argparser::arg_parser("Run DESeq2 and DEXSeq on bambu output.")
|
||||||
|
parser <- argparser::add_argument(parser, "--transcript_rds", help = "bambu transcript RDS.")
|
||||||
|
parser <- argparser::add_argument(parser, "--gene_rds", help = "bambu gene RDS.")
|
||||||
|
parser <- argparser::add_argument(parser, "--sample_sheet", help = "Sample sheet CSV.")
|
||||||
|
parser <- argparser::add_argument(
|
||||||
|
parser,
|
||||||
|
"--condition_column",
|
||||||
|
help = "Primary condition column.",
|
||||||
|
default = "condition"
|
||||||
|
)
|
||||||
|
parser <- argparser::add_argument(parser, "--covariates", help = "Comma-separated nuisance covariates.")
|
||||||
|
parser <- argparser::add_argument(parser, "--reference_level", help = "Reference level for the condition column.")
|
||||||
|
argparser::add_argument(parser, "--out_dir", help = "Output directory.", default = "de_analysis")
|
||||||
|
}
|
||||||
|
|
||||||
|
de_parse_covariates <- function(value) {
|
||||||
|
workflow_glue_r_parse_csv_list(value)
|
||||||
|
}
|
||||||
|
|
||||||
|
de_validate_inputs <- function(tx_se, gene_se, sample_df, argv) {
|
||||||
|
workflow_glue_r_require_args(argv, c("transcript_rds", "gene_rds", "sample_sheet"))
|
||||||
|
|
||||||
|
covariates <- de_parse_covariates(argv$covariates)
|
||||||
|
workflow_glue_r_validate_r_formula_names(
|
||||||
|
c(argv$condition_column, covariates),
|
||||||
|
label = "Design column"
|
||||||
|
)
|
||||||
|
|
||||||
|
if (!"alias" %in% names(sample_df)) {
|
||||||
|
stop("Sample sheet must contain an 'alias' column.", call. = FALSE)
|
||||||
|
}
|
||||||
|
duplicate_sample_aliases <- unique(sample_df$alias[duplicated(sample_df$alias)])
|
||||||
|
if (length(duplicate_sample_aliases) > 0) {
|
||||||
|
stop(
|
||||||
|
sprintf(
|
||||||
|
"Sample sheet aliases must be unique; duplicated aliases: %s",
|
||||||
|
paste(duplicate_sample_aliases, collapse = ", ")
|
||||||
|
),
|
||||||
|
call. = FALSE
|
||||||
|
)
|
||||||
|
}
|
||||||
|
if (!(argv$condition_column %in% names(sample_df))) {
|
||||||
|
stop(
|
||||||
|
sprintf("Sample sheet must contain the '%s' column.", argv$condition_column),
|
||||||
|
call. = FALSE
|
||||||
|
)
|
||||||
|
}
|
||||||
|
|
||||||
|
missing_covariates <- setdiff(covariates, names(sample_df))
|
||||||
|
if (length(missing_covariates) > 0) {
|
||||||
|
stop(
|
||||||
|
sprintf(
|
||||||
|
"Missing covariate columns: %s",
|
||||||
|
paste(missing_covariates, collapse = ", ")
|
||||||
|
),
|
||||||
|
call. = FALSE
|
||||||
|
)
|
||||||
|
}
|
||||||
|
|
||||||
|
if (any(duplicated(colnames(tx_se)))) {
|
||||||
|
stop("Transcript RDS sample names must be unique.", call. = FALSE)
|
||||||
|
}
|
||||||
|
if (any(duplicated(colnames(gene_se)))) {
|
||||||
|
stop("Gene RDS sample names must be unique.", call. = FALSE)
|
||||||
|
}
|
||||||
|
if (!setequal(colnames(tx_se), colnames(gene_se))) {
|
||||||
|
stop("Transcript and gene RDS sample names must match.", call. = FALSE)
|
||||||
|
}
|
||||||
|
|
||||||
|
tx_counts <- SummarizedExperiment::assays(tx_se)$counts
|
||||||
|
gene_counts <- SummarizedExperiment::assays(gene_se)$counts
|
||||||
|
if (is.null(tx_counts)) {
|
||||||
|
stop("Transcript RDS must contain a 'counts' assay.", call. = FALSE)
|
||||||
|
}
|
||||||
|
if (is.null(gene_counts)) {
|
||||||
|
stop("Gene RDS must contain a 'counts' assay.", call. = FALSE)
|
||||||
|
}
|
||||||
|
if (anyNA(tx_counts) || anyNA(gene_counts)) {
|
||||||
|
stop("Count matrices must not contain NA values.", call. = FALSE)
|
||||||
|
}
|
||||||
|
zero_count_samples <- unique(c(
|
||||||
|
colnames(tx_counts)[colSums(tx_counts) == 0],
|
||||||
|
colnames(gene_counts)[colSums(gene_counts) == 0]
|
||||||
|
))
|
||||||
|
if (length(zero_count_samples) > 0) {
|
||||||
|
stop(
|
||||||
|
sprintf(
|
||||||
|
"Count matrices contain samples with zero total counts: %s",
|
||||||
|
paste(zero_count_samples, collapse = ", ")
|
||||||
|
),
|
||||||
|
call. = FALSE
|
||||||
|
)
|
||||||
|
}
|
||||||
|
|
||||||
|
sample_df <- sample_df[match(colnames(tx_se), sample_df$alias), , drop = FALSE]
|
||||||
|
if (any(is.na(sample_df$alias))) {
|
||||||
|
stop("Sample sheet aliases do not match the bambu output sample names.", call. = FALSE)
|
||||||
|
}
|
||||||
|
|
||||||
|
condition_values <- unique(sample_df[[argv$condition_column]])
|
||||||
|
if (length(condition_values) < 2) {
|
||||||
|
stop("Differential analysis requires at least two condition levels.", call. = FALSE)
|
||||||
|
}
|
||||||
|
|
||||||
|
reference_level <- argv$reference_level
|
||||||
|
if (workflow_glue_r_arg_missing(reference_level)) {
|
||||||
|
if ("control" %in% condition_values) {
|
||||||
|
reference_level <- "control"
|
||||||
|
} else {
|
||||||
|
stop(
|
||||||
|
"Provide --reference_level when the condition column does not contain 'control'.",
|
||||||
|
call. = FALSE
|
||||||
|
)
|
||||||
|
}
|
||||||
|
}
|
||||||
|
if (!(reference_level %in% condition_values)) {
|
||||||
|
stop("The requested reference level is not present in the condition column.", call. = FALSE)
|
||||||
|
}
|
||||||
|
|
||||||
|
sample_df[[argv$condition_column]] <- factor(sample_df[[argv$condition_column]])
|
||||||
|
for (covariate in covariates) {
|
||||||
|
sample_df[[covariate]] <- factor(sample_df[[covariate]])
|
||||||
|
}
|
||||||
|
|
||||||
|
list(
|
||||||
|
tx_se = tx_se,
|
||||||
|
gene_se = gene_se,
|
||||||
|
sample_df = sample_df,
|
||||||
|
covariates = covariates,
|
||||||
|
condition_values = condition_values,
|
||||||
|
reference_level = reference_level
|
||||||
|
)
|
||||||
|
}
|
||||||
|
|
||||||
|
de_build_contrast_name <- function(condition_column, target_level, reference_level) {
|
||||||
|
sprintf("%s_%s_vs_%s", condition_column, target_level, reference_level)
|
||||||
|
}
|
||||||
|
|
||||||
|
de_set_dispersions <- function(object, value) {
|
||||||
|
setter <- get("dispersions<-", envir = asNamespace("DESeq2"))
|
||||||
|
setter(object, value = value)
|
||||||
|
}
|
||||||
|
|
||||||
|
de_run_deseq_with_fallback <- function(dds, contrast_name, out_dir) {
|
||||||
|
tryCatch(
|
||||||
|
DESeq2::DESeq(dds, quiet = TRUE),
|
||||||
|
error = function(err) {
|
||||||
|
if (!grepl(
|
||||||
|
"all gene-wise dispersion estimates are within 2 orders of magnitude",
|
||||||
|
conditionMessage(err),
|
||||||
|
fixed = TRUE
|
||||||
|
)) {
|
||||||
|
stop(err)
|
||||||
|
}
|
||||||
|
|
||||||
|
warning(
|
||||||
|
"STATISTICAL POWER REDUCED: DESeq2 dispersion estimation failed for ",
|
||||||
|
contrast_name,
|
||||||
|
".\n",
|
||||||
|
"This usually indicates:\n",
|
||||||
|
" 1. Too few replicates (recommend n>=3 per group)\n",
|
||||||
|
" 2. High biological variability\n",
|
||||||
|
" 3. Poor data quality\n",
|
||||||
|
"Falling back to gene-wise dispersion (no information sharing).\n",
|
||||||
|
"Results will have reduced power and wider confidence intervals."
|
||||||
|
)
|
||||||
|
|
||||||
|
dds <- DESeq2::estimateSizeFactors(dds)
|
||||||
|
dds <- DESeq2::estimateDispersionsGeneEst(dds)
|
||||||
|
dds <- de_set_dispersions(dds, S4Vectors::mcols(dds)$dispGeneEst)
|
||||||
|
|
||||||
|
diag_content <- c(
|
||||||
|
"DESeq2 Dispersion Estimation Fallback Applied",
|
||||||
|
"==============================================",
|
||||||
|
"",
|
||||||
|
sprintf("Timestamp: %s", format(Sys.time(), "%Y-%m-%d %H:%M:%S")),
|
||||||
|
sprintf("Contrast: %s", contrast_name),
|
||||||
|
sprintf("Samples: %d", ncol(dds)),
|
||||||
|
sprintf("Genes tested: %d", nrow(dds)),
|
||||||
|
sprintf(
|
||||||
|
"Dispersion range: %.3f to %.3f",
|
||||||
|
min(DESeq2::dispersions(dds)),
|
||||||
|
max(DESeq2::dispersions(dds))
|
||||||
|
),
|
||||||
|
"",
|
||||||
|
"WHAT HAPPENED:",
|
||||||
|
" Curve fitting failed. Using gene-wise dispersion estimates.",
|
||||||
|
"",
|
||||||
|
"IMPLICATIONS:",
|
||||||
|
" - No information sharing across genes",
|
||||||
|
" - Reduced statistical power",
|
||||||
|
" - Wider confidence intervals",
|
||||||
|
" - More conservative results (fewer discoveries)",
|
||||||
|
"",
|
||||||
|
"LIKELY CAUSES:",
|
||||||
|
" 1. Too few replicates (recommend n>=3 per group)",
|
||||||
|
" 2. High biological variability",
|
||||||
|
" 3. Poor data quality or outlier samples",
|
||||||
|
"",
|
||||||
|
"RECOMMENDATIONS:",
|
||||||
|
" - Add more biological replicates if possible",
|
||||||
|
" - Check sample quality metrics",
|
||||||
|
" - Consider filtering low-count genes more stringently"
|
||||||
|
)
|
||||||
|
|
||||||
|
diag_file <- file.path(
|
||||||
|
out_dir,
|
||||||
|
sprintf(
|
||||||
|
"DESeq2_dispersion_fallback_%s.txt",
|
||||||
|
gsub("[^A-Za-z0-9_-]", "_", contrast_name)
|
||||||
|
)
|
||||||
|
)
|
||||||
|
writeLines(diag_content, diag_file)
|
||||||
|
|
||||||
|
DESeq2::nbinomWaldTest(dds)
|
||||||
|
}
|
||||||
|
)
|
||||||
|
}
|
||||||
|
|
||||||
|
de_estimate_dispersions_with_fallback <- function(object, context_label, allow_gene_est = TRUE) {
|
||||||
|
tryCatch(
|
||||||
|
DESeq2::estimateDispersions(object),
|
||||||
|
error = function(err) {
|
||||||
|
if (!grepl(
|
||||||
|
"all gene-wise dispersion estimates are within 2 orders of magnitude",
|
||||||
|
conditionMessage(err),
|
||||||
|
fixed = TRUE
|
||||||
|
)) {
|
||||||
|
stop(err)
|
||||||
|
}
|
||||||
|
|
||||||
|
message(context_label, " dispersion fitting failed; retrying with fitType='local'.")
|
||||||
|
tryCatch(
|
||||||
|
DESeq2::estimateDispersions(object, fitType = "local"),
|
||||||
|
error = function(local_err) {
|
||||||
|
if (!grepl(
|
||||||
|
"all gene-wise dispersion estimates are within 2 orders of magnitude",
|
||||||
|
conditionMessage(local_err),
|
||||||
|
fixed = TRUE
|
||||||
|
)) {
|
||||||
|
stop(local_err)
|
||||||
|
}
|
||||||
|
|
||||||
|
message(context_label, " local-fit dispersion retry failed; retrying with fitType='mean'.")
|
||||||
|
tryCatch(
|
||||||
|
DESeq2::estimateDispersions(object, fitType = "mean"),
|
||||||
|
error = function(mean_err) {
|
||||||
|
if (!grepl(
|
||||||
|
"all gene-wise dispersion estimates are within 2 orders of magnitude",
|
||||||
|
conditionMessage(mean_err),
|
||||||
|
fixed = TRUE
|
||||||
|
)) {
|
||||||
|
stop(mean_err)
|
||||||
|
}
|
||||||
|
if (!allow_gene_est) {
|
||||||
|
stop(mean_err)
|
||||||
|
}
|
||||||
|
|
||||||
|
message(
|
||||||
|
context_label,
|
||||||
|
" mean-fit dispersion retry failed; falling back to gene-wise dispersion estimates."
|
||||||
|
)
|
||||||
|
object <- DESeq2::estimateDispersionsGeneEst(object)
|
||||||
|
object <- de_set_dispersions(object, S4Vectors::mcols(object)$dispGeneEst)
|
||||||
|
object
|
||||||
|
}
|
||||||
|
)
|
||||||
|
}
|
||||||
|
)
|
||||||
|
}
|
||||||
|
)
|
||||||
|
}
|
||||||
|
|
||||||
|
de_is_recoverable_dexseq_error <- function(message_text) {
|
||||||
|
grepl(
|
||||||
|
"all gene-wise dispersion estimates are within 2 orders of magnitude",
|
||||||
|
message_text,
|
||||||
|
fixed = TRUE
|
||||||
|
) || grepl(
|
||||||
|
"model matrix is not full rank",
|
||||||
|
message_text,
|
||||||
|
fixed = TRUE
|
||||||
|
) || grepl(
|
||||||
|
"replacement has 1 row, data has 0",
|
||||||
|
message_text,
|
||||||
|
fixed = TRUE
|
||||||
|
)
|
||||||
|
}
|
||||||
|
|
||||||
|
de_write_placeholder_pdf <- function(path, label) {
|
||||||
|
grDevices::pdf(path)
|
||||||
|
graphics::plot.new()
|
||||||
|
graphics::text(0.5, 0.5, label, cex = 0.9)
|
||||||
|
grDevices::dev.off()
|
||||||
|
}
|
||||||
|
|
||||||
|
de_run_deseq2_result <- function(
|
||||||
|
count_mat,
|
||||||
|
coldata,
|
||||||
|
target_level,
|
||||||
|
reference_level,
|
||||||
|
condition_column,
|
||||||
|
covariates,
|
||||||
|
out_dir,
|
||||||
|
contrast_name
|
||||||
|
) {
|
||||||
|
design_terms <- c(covariates, condition_column)
|
||||||
|
design_formula <- stats::as.formula(paste("~", paste(design_terms, collapse = " + ")))
|
||||||
|
dds <- DESeq2::DESeqDataSetFromMatrix(
|
||||||
|
countData = round(count_mat),
|
||||||
|
colData = coldata,
|
||||||
|
design = design_formula
|
||||||
|
)
|
||||||
|
dds <- de_run_deseq_with_fallback(dds, contrast_name, out_dir)
|
||||||
|
result <- DESeq2::results(
|
||||||
|
dds,
|
||||||
|
contrast = c(condition_column, target_level, reference_level),
|
||||||
|
independentFiltering = TRUE
|
||||||
|
)
|
||||||
|
list(dds = dds, result = result)
|
||||||
|
}
|
||||||
|
|
||||||
|
de_run_dexseq_result <- function(
|
||||||
|
tx_counts,
|
||||||
|
tx_meta,
|
||||||
|
coldata,
|
||||||
|
condition_column,
|
||||||
|
covariates
|
||||||
|
) {
|
||||||
|
coldata$sample <- factor(coldata$alias)
|
||||||
|
coldata[[condition_column]] <- factor(coldata[[condition_column]])
|
||||||
|
for (covariate in covariates) {
|
||||||
|
coldata[[covariate]] <- factor(coldata[[covariate]])
|
||||||
|
}
|
||||||
|
|
||||||
|
run_inner <- function(active_covariates) {
|
||||||
|
covariate_exon_terms <- if (length(active_covariates) > 0) {
|
||||||
|
paste0(active_covariates, ":exon")
|
||||||
|
} else {
|
||||||
|
character(0)
|
||||||
|
}
|
||||||
|
design_terms <- c("sample", "exon", covariate_exon_terms, paste0(condition_column, ":exon"))
|
||||||
|
reduced_terms <- c("sample", "exon", covariate_exon_terms)
|
||||||
|
full_formula <- stats::as.formula(paste("~", paste(design_terms, collapse = " + ")))
|
||||||
|
reduced_formula <- stats::as.formula(paste("~", paste(reduced_terms, collapse = " + ")))
|
||||||
|
|
||||||
|
tryCatch({
|
||||||
|
dxd <- DEXSeq::DEXSeqDataSet(
|
||||||
|
countData = round(tx_counts),
|
||||||
|
sampleData = as.data.frame(coldata),
|
||||||
|
design = full_formula,
|
||||||
|
featureID = tx_meta$TXNAME,
|
||||||
|
groupID = tx_meta$GENEID
|
||||||
|
)
|
||||||
|
dxd <- DESeq2::estimateSizeFactors(dxd)
|
||||||
|
dxd <- de_estimate_dispersions_with_fallback(dxd, "DEXSeq", allow_gene_est = TRUE)
|
||||||
|
dxd <- DEXSeq::testForDEU(dxd, reducedModel = reduced_formula)
|
||||||
|
dxd <- DEXSeq::estimateExonFoldChanges(dxd, fitExpToVar = condition_column)
|
||||||
|
dxr <- DEXSeq::DEXSeqResults(dxd, independentFiltering = FALSE)
|
||||||
|
list(dxd = dxd, dxr = dxr)
|
||||||
|
}, error = function(err) {
|
||||||
|
if (length(active_covariates) == 0 || !grepl(
|
||||||
|
"model matrix is not full rank",
|
||||||
|
conditionMessage(err),
|
||||||
|
fixed = TRUE
|
||||||
|
)) {
|
||||||
|
stop(err)
|
||||||
|
}
|
||||||
|
|
||||||
|
dropped_covariate <- tail(active_covariates, 1)
|
||||||
|
kept_covariates <- head(active_covariates, -1)
|
||||||
|
message(
|
||||||
|
"DEXSeq design was not full rank with covariate '",
|
||||||
|
dropped_covariate,
|
||||||
|
"'; retrying without it."
|
||||||
|
)
|
||||||
|
run_inner(kept_covariates)
|
||||||
|
})
|
||||||
|
}
|
||||||
|
|
||||||
|
run_inner(covariates)
|
||||||
|
}
|
||||||
|
|
||||||
|
main_run_de_analysis <- function(
|
||||||
|
argv,
|
||||||
|
deseq_runner = de_run_deseq2_result,
|
||||||
|
dexseq_runner = de_run_dexseq_result,
|
||||||
|
pdf_fn = grDevices::pdf,
|
||||||
|
dev_off_fn = grDevices::dev.off,
|
||||||
|
plot_ma_fn = DESeq2::plotMA,
|
||||||
|
plot_disp_fn = DESeq2::plotDispEsts,
|
||||||
|
per_gene_q_fn = DEXSeq::perGeneQValue,
|
||||||
|
placeholder_pdf_fn = de_write_placeholder_pdf
|
||||||
|
) {
|
||||||
|
set.seed(42)
|
||||||
|
dir.create(argv$out_dir, showWarnings = FALSE, recursive = TRUE)
|
||||||
|
|
||||||
|
tx_se <- readRDS(argv$transcript_rds)
|
||||||
|
gene_se <- readRDS(argv$gene_rds)
|
||||||
|
sample_df <- workflow_glue_r_read_csv(argv$sample_sheet)
|
||||||
|
validated <- de_validate_inputs(tx_se, gene_se, sample_df, argv)
|
||||||
|
sample_df <- validated$sample_df
|
||||||
|
covariates <- validated$covariates
|
||||||
|
condition_values <- validated$condition_values
|
||||||
|
reference_level <- validated$reference_level
|
||||||
|
|
||||||
|
tx_meta <- as.data.frame(SummarizedExperiment::rowData(tx_se))
|
||||||
|
if (!"TXNAME" %in% names(tx_meta)) {
|
||||||
|
tx_meta$TXNAME <- rownames(tx_se)
|
||||||
|
}
|
||||||
|
if (!"GENEID" %in% names(tx_meta)) {
|
||||||
|
stop("Transcript rowData must contain GENEID for DEXSeq.", call. = FALSE)
|
||||||
|
}
|
||||||
|
|
||||||
|
gene_meta <- as.data.frame(SummarizedExperiment::rowData(gene_se))
|
||||||
|
if (!"GENEID" %in% names(gene_meta)) {
|
||||||
|
gene_meta$GENEID <- rownames(gene_se)
|
||||||
|
}
|
||||||
|
|
||||||
|
targets <- setdiff(as.character(condition_values), reference_level)
|
||||||
|
|
||||||
|
de_qc_stats <- list(
|
||||||
|
timestamp = format(Sys.time(), "%Y-%m-%d %H:%M:%S"),
|
||||||
|
total_samples = nrow(sample_df),
|
||||||
|
condition_column = argv$condition_column,
|
||||||
|
reference_level = reference_level,
|
||||||
|
covariates = if (length(covariates) > 0) covariates else "none",
|
||||||
|
num_contrasts = length(targets),
|
||||||
|
contrasts = list()
|
||||||
|
)
|
||||||
|
|
||||||
|
n_per_group <- table(sample_df[[argv$condition_column]])
|
||||||
|
de_qc_stats$samples_per_group <- as.list(n_per_group)
|
||||||
|
|
||||||
|
sample_size_warnings <- character(0)
|
||||||
|
if (any(n_per_group < 3)) {
|
||||||
|
warning(
|
||||||
|
"WARNING: Some condition groups have fewer than 3 replicates.\n",
|
||||||
|
"Recommended minimum for DGE: n=3 per group\n",
|
||||||
|
"Current sample sizes: ",
|
||||||
|
paste(names(n_per_group), "=", n_per_group, collapse = ", "),
|
||||||
|
"\nResults may have reduced statistical power."
|
||||||
|
)
|
||||||
|
sample_size_warnings <- c(sample_size_warnings, "Some groups have n<3 (recommended minimum)")
|
||||||
|
}
|
||||||
|
if (any(n_per_group < 2)) {
|
||||||
|
stop(
|
||||||
|
"ERROR: Some condition groups have fewer than 2 replicates. Cannot perform statistical testing.",
|
||||||
|
call. = FALSE
|
||||||
|
)
|
||||||
|
}
|
||||||
|
de_qc_stats$sample_size_warnings <- if (length(sample_size_warnings) > 0) {
|
||||||
|
sample_size_warnings
|
||||||
|
} else {
|
||||||
|
"none"
|
||||||
|
}
|
||||||
|
|
||||||
|
if (length(targets) > 1) {
|
||||||
|
fwer <- (1 - (1 - 0.05)^length(targets)) * 100
|
||||||
|
mt_warning <- sprintf(
|
||||||
|
"Multiple contrasts tested (%d). Per-contrast FDR < 0.05 yields family-wise error rate of ~%.1f%%",
|
||||||
|
length(targets),
|
||||||
|
fwer
|
||||||
|
)
|
||||||
|
message("WARNING: ", mt_warning)
|
||||||
|
de_qc_stats$multiple_testing_note <- mt_warning
|
||||||
|
mt_content <- c(
|
||||||
|
"Multiple Testing Across Contrasts",
|
||||||
|
"==================================",
|
||||||
|
"",
|
||||||
|
sprintf("Timestamp: %s", format(Sys.time(), "%Y-%m-%d %H:%M:%S")),
|
||||||
|
sprintf("Number of contrasts tested: %d", length(targets)),
|
||||||
|
sprintf("Contrasts: %s", paste(sprintf("%s vs %s", targets, reference_level), collapse = ", ")),
|
||||||
|
"",
|
||||||
|
"PER-CONTRAST FDR THRESHOLD: 0.05",
|
||||||
|
sprintf("FAMILY-WISE ERROR RATE: ~%.1f%%", fwer),
|
||||||
|
"",
|
||||||
|
"WHAT THIS MEANS:",
|
||||||
|
" Each contrast uses FDR < 0.05 independently.",
|
||||||
|
" When testing multiple contrasts, the overall false positive rate increases.",
|
||||||
|
sprintf(" Expected: %.1f%% chance of at least one false positive across all contrasts", fwer),
|
||||||
|
"",
|
||||||
|
"RECOMMENDATIONS:",
|
||||||
|
sprintf(" 1. Bonferroni correction: 0.05 / %d = %.4f", length(targets), 0.05 / length(targets)),
|
||||||
|
" 2. Focus on pre-specified contrasts of interest",
|
||||||
|
" 3. Treat results as exploratory and validate key findings",
|
||||||
|
" 4. Consider using hierarchical testing procedures",
|
||||||
|
""
|
||||||
|
)
|
||||||
|
writeLines(mt_content, file.path(argv$out_dir, "MULTIPLE_TESTING_WARNING.txt"))
|
||||||
|
}
|
||||||
|
|
||||||
|
for (target_level in targets) {
|
||||||
|
contrast_name <- de_build_contrast_name(argv$condition_column, target_level, reference_level)
|
||||||
|
contrast_dir <- file.path(argv$out_dir, contrast_name)
|
||||||
|
dir.create(contrast_dir, showWarnings = FALSE, recursive = TRUE)
|
||||||
|
|
||||||
|
keep_samples <- sample_df[[argv$condition_column]] %in% c(reference_level, target_level)
|
||||||
|
contrast_samples <- droplevels(sample_df[keep_samples, , drop = FALSE])
|
||||||
|
contrast_samples[[argv$condition_column]] <- stats::relevel(
|
||||||
|
factor(contrast_samples[[argv$condition_column]]),
|
||||||
|
ref = reference_level
|
||||||
|
)
|
||||||
|
|
||||||
|
contrast_qc <- list(
|
||||||
|
name = contrast_name,
|
||||||
|
target_level = target_level,
|
||||||
|
reference_level = reference_level,
|
||||||
|
n_samples = nrow(contrast_samples),
|
||||||
|
n_target = sum(contrast_samples[[argv$condition_column]] == target_level),
|
||||||
|
n_reference = sum(contrast_samples[[argv$condition_column]] == reference_level)
|
||||||
|
)
|
||||||
|
|
||||||
|
if (nrow(contrast_samples) < 6) {
|
||||||
|
contrast_qc$dtu_power_warning <- sprintf(
|
||||||
|
"DTU analysis may be underpowered (n=%d, recommend n>=6 with >=3 per group)",
|
||||||
|
nrow(contrast_samples)
|
||||||
|
)
|
||||||
|
warning(contrast_qc$dtu_power_warning)
|
||||||
|
}
|
||||||
|
|
||||||
|
gene_counts <- SummarizedExperiment::assays(gene_se)$counts[, contrast_samples$alias, drop = FALSE]
|
||||||
|
tx_counts <- SummarizedExperiment::assays(tx_se)$counts[, contrast_samples$alias, drop = FALSE]
|
||||||
|
contrast_qc$genes_tested <- nrow(gene_counts)
|
||||||
|
contrast_qc$transcripts_tested <- nrow(tx_counts)
|
||||||
|
|
||||||
|
dge_run <- deseq_runner(
|
||||||
|
gene_counts,
|
||||||
|
contrast_samples,
|
||||||
|
target_level,
|
||||||
|
reference_level,
|
||||||
|
argv$condition_column,
|
||||||
|
covariates,
|
||||||
|
argv$out_dir,
|
||||||
|
contrast_name
|
||||||
|
)
|
||||||
|
dge_res <- as.data.frame(dge_run$result)
|
||||||
|
dge_res$GENEID <- rownames(dge_res)
|
||||||
|
dge_res <- merge(gene_meta, dge_res, by = "GENEID", all.y = TRUE, sort = FALSE)
|
||||||
|
dge_res <- workflow_glue_r_normalise_tsv_df(dge_res)
|
||||||
|
|
||||||
|
contrast_qc$dge_total_genes <- nrow(dge_res)
|
||||||
|
contrast_qc$dge_significant_fdr05 <- sum(dge_res$padj < 0.05, na.rm = TRUE)
|
||||||
|
contrast_qc$dge_significant_fdr01 <- sum(dge_res$padj < 0.01, na.rm = TRUE)
|
||||||
|
contrast_qc$dge_upregulated <- sum(
|
||||||
|
dge_res$padj < 0.05 & dge_res$log2FoldChange > 0,
|
||||||
|
na.rm = TRUE
|
||||||
|
)
|
||||||
|
contrast_qc$dge_downregulated <- sum(
|
||||||
|
dge_res$padj < 0.05 & dge_res$log2FoldChange < 0,
|
||||||
|
na.rm = TRUE
|
||||||
|
)
|
||||||
|
|
||||||
|
utils::write.table(
|
||||||
|
dge_res[order(dge_res$padj), ],
|
||||||
|
file = file.path(contrast_dir, "results_dge.tsv"),
|
||||||
|
sep = "\t",
|
||||||
|
quote = FALSE,
|
||||||
|
row.names = FALSE
|
||||||
|
)
|
||||||
|
|
||||||
|
pdf_fn(file.path(contrast_dir, "results_dge.pdf"))
|
||||||
|
plot_ma_fn(dge_run$result)
|
||||||
|
dev_off_fn()
|
||||||
|
|
||||||
|
dex_res <- tryCatch(
|
||||||
|
dexseq_runner(
|
||||||
|
tx_counts,
|
||||||
|
tx_meta,
|
||||||
|
contrast_samples,
|
||||||
|
argv$condition_column,
|
||||||
|
covariates
|
||||||
|
),
|
||||||
|
error = function(err) {
|
||||||
|
message_text <- conditionMessage(err)
|
||||||
|
if (!de_is_recoverable_dexseq_error(message_text)) {
|
||||||
|
stop(err)
|
||||||
|
}
|
||||||
|
|
||||||
|
warning(
|
||||||
|
"DEXSeq failed for contrast ",
|
||||||
|
target_level,
|
||||||
|
" vs ",
|
||||||
|
reference_level,
|
||||||
|
"\nError: ",
|
||||||
|
message_text
|
||||||
|
)
|
||||||
|
|
||||||
|
failure_content <- c(
|
||||||
|
"DTU Analysis Failed",
|
||||||
|
"===================",
|
||||||
|
"",
|
||||||
|
sprintf("Timestamp: %s", format(Sys.time(), "%Y-%m-%d %H:%M:%S")),
|
||||||
|
sprintf("Contrast: %s vs %s", target_level, reference_level),
|
||||||
|
sprintf(
|
||||||
|
"Samples: %d (%d %s, %d %s)",
|
||||||
|
nrow(contrast_samples),
|
||||||
|
sum(contrast_samples[[argv$condition_column]] == target_level),
|
||||||
|
target_level,
|
||||||
|
sum(contrast_samples[[argv$condition_column]] == reference_level),
|
||||||
|
reference_level
|
||||||
|
),
|
||||||
|
sprintf("Transcripts: %d", nrow(tx_counts)),
|
||||||
|
"",
|
||||||
|
"ERROR MESSAGE:",
|
||||||
|
sprintf(" %s", message_text),
|
||||||
|
"",
|
||||||
|
"DTU RESULTS CANNOT BE INTERPRETED",
|
||||||
|
""
|
||||||
|
)
|
||||||
|
writeLines(failure_content, file.path(contrast_dir, "DTU_ANALYSIS_FAILED.txt"))
|
||||||
|
NULL
|
||||||
|
}
|
||||||
|
)
|
||||||
|
|
||||||
|
if (is.null(dex_res)) {
|
||||||
|
dex_df <- workflow_glue_r_empty_tsv(c(
|
||||||
|
"featureID",
|
||||||
|
"groupID",
|
||||||
|
"log2fold",
|
||||||
|
"pvalue",
|
||||||
|
"padj",
|
||||||
|
"exonBaseMean"
|
||||||
|
))
|
||||||
|
tx_dtu <- dex_df
|
||||||
|
gene_dtu <- workflow_glue_r_empty_tsv(c("GENEID", "qval"))
|
||||||
|
placeholder_pdf_fn(
|
||||||
|
file.path(contrast_dir, "results_dtu.pdf"),
|
||||||
|
"DEXSeq did not converge for this contrast.\nSee DTU_ANALYSIS_FAILED.txt for details."
|
||||||
|
)
|
||||||
|
contrast_qc$dtu_status <- "FAILED"
|
||||||
|
contrast_qc$dtu_significant_transcripts <- 0
|
||||||
|
contrast_qc$dtu_significant_genes <- 0
|
||||||
|
} else {
|
||||||
|
dex_df <- as.data.frame(dex_res$dxr)
|
||||||
|
dex_df <- workflow_glue_r_normalise_tsv_df(dex_df)
|
||||||
|
tx_dtu <- dex_df[, intersect(
|
||||||
|
c("featureID", "groupID", "log2fold", "pvalue", "padj", "exonBaseMean"),
|
||||||
|
names(dex_df)
|
||||||
|
), drop = FALSE]
|
||||||
|
tx_dtu <- workflow_glue_r_normalise_tsv_df(tx_dtu)
|
||||||
|
|
||||||
|
gene_q <- per_gene_q_fn(dex_res$dxr)
|
||||||
|
gene_dtu <- data.frame(
|
||||||
|
GENEID = names(gene_q),
|
||||||
|
qval = unname(gene_q),
|
||||||
|
row.names = NULL
|
||||||
|
)
|
||||||
|
|
||||||
|
contrast_qc$dtu_status <- "SUCCESS"
|
||||||
|
contrast_qc$dtu_significant_transcripts <- sum(tx_dtu$padj < 0.05, na.rm = TRUE)
|
||||||
|
contrast_qc$dtu_significant_genes <- sum(gene_dtu$qval < 0.05, na.rm = TRUE)
|
||||||
|
|
||||||
|
pdf_fn(file.path(contrast_dir, "results_dtu.pdf"))
|
||||||
|
plot_ma_fn(dex_res$dxr, cex = 0.8, alpha = 0.05)
|
||||||
|
plot_disp_fn(dex_res$dxd)
|
||||||
|
dev_off_fn()
|
||||||
|
}
|
||||||
|
|
||||||
|
utils::write.table(
|
||||||
|
dex_df,
|
||||||
|
file = file.path(contrast_dir, "results_dexseq.tsv"),
|
||||||
|
sep = "\t",
|
||||||
|
quote = FALSE,
|
||||||
|
row.names = FALSE
|
||||||
|
)
|
||||||
|
utils::write.table(
|
||||||
|
tx_dtu[order(tx_dtu$padj), ],
|
||||||
|
file = file.path(contrast_dir, "results_dtu_transcript.tsv"),
|
||||||
|
sep = "\t",
|
||||||
|
quote = FALSE,
|
||||||
|
row.names = FALSE
|
||||||
|
)
|
||||||
|
utils::write.table(
|
||||||
|
gene_dtu[order(gene_dtu$qval), ],
|
||||||
|
file = file.path(contrast_dir, "results_dtu_gene.tsv"),
|
||||||
|
sep = "\t",
|
||||||
|
quote = FALSE,
|
||||||
|
row.names = FALSE
|
||||||
|
)
|
||||||
|
utils::write.table(
|
||||||
|
contrast_samples,
|
||||||
|
file = file.path(contrast_dir, "samples_used.tsv"),
|
||||||
|
sep = "\t",
|
||||||
|
quote = FALSE,
|
||||||
|
row.names = FALSE
|
||||||
|
)
|
||||||
|
|
||||||
|
contrast_qc_summary <- c(
|
||||||
|
sprintf("Contrast QC Summary: %s", contrast_name),
|
||||||
|
paste(rep("=", 50), collapse = ""),
|
||||||
|
"",
|
||||||
|
"Sample Information:",
|
||||||
|
sprintf(" Target level (%s): %d samples", target_level, contrast_qc$n_target),
|
||||||
|
sprintf(" Reference level (%s): %d samples", reference_level, contrast_qc$n_reference),
|
||||||
|
sprintf(" Total samples: %d", contrast_qc$n_samples),
|
||||||
|
"",
|
||||||
|
"DGE Results:",
|
||||||
|
sprintf(" Genes tested: %d", contrast_qc$genes_tested),
|
||||||
|
sprintf(" Significant (FDR < 0.05): %d", contrast_qc$dge_significant_fdr05),
|
||||||
|
sprintf(" Significant (FDR < 0.01): %d", contrast_qc$dge_significant_fdr01),
|
||||||
|
sprintf(" Upregulated: %d", contrast_qc$dge_upregulated),
|
||||||
|
sprintf(" Downregulated: %d", contrast_qc$dge_downregulated),
|
||||||
|
"",
|
||||||
|
"DTU Results:",
|
||||||
|
sprintf(" Status: %s", contrast_qc$dtu_status),
|
||||||
|
sprintf(" Transcripts tested: %d", contrast_qc$transcripts_tested),
|
||||||
|
if (contrast_qc$dtu_status == "SUCCESS") {
|
||||||
|
c(
|
||||||
|
sprintf(
|
||||||
|
" Significant transcripts (FDR < 0.05): %d",
|
||||||
|
contrast_qc$dtu_significant_transcripts
|
||||||
|
),
|
||||||
|
sprintf(" Genes with DTU (q < 0.05): %d", contrast_qc$dtu_significant_genes)
|
||||||
|
)
|
||||||
|
} else {
|
||||||
|
" See DTU_ANALYSIS_FAILED.txt for details"
|
||||||
|
},
|
||||||
|
if (!is.null(contrast_qc$dtu_power_warning)) paste0(" WARNING: ", contrast_qc$dtu_power_warning) else NULL,
|
||||||
|
""
|
||||||
|
)
|
||||||
|
writeLines(contrast_qc_summary, file.path(contrast_dir, "contrast_qc_summary.txt"))
|
||||||
|
de_qc_stats$contrasts[[contrast_name]] <- contrast_qc
|
||||||
|
}
|
||||||
|
|
||||||
|
jsonlite::write_json(
|
||||||
|
de_qc_stats,
|
||||||
|
file.path(argv$out_dir, "de_qc_stats.json"),
|
||||||
|
pretty = TRUE,
|
||||||
|
auto_unbox = TRUE
|
||||||
|
)
|
||||||
|
|
||||||
|
overall_summary <- c(
|
||||||
|
"Differential Expression/Usage Analysis Summary",
|
||||||
|
paste(rep("=", 50), collapse = ""),
|
||||||
|
"",
|
||||||
|
sprintf("Timestamp: %s", de_qc_stats$timestamp),
|
||||||
|
sprintf("Total samples: %d", de_qc_stats$total_samples),
|
||||||
|
sprintf("Condition column: %s", de_qc_stats$condition_column),
|
||||||
|
sprintf("Reference level: %s", de_qc_stats$reference_level),
|
||||||
|
sprintf("Covariates: %s", paste(de_qc_stats$covariates, collapse = ", ")),
|
||||||
|
"",
|
||||||
|
"Sample Sizes:",
|
||||||
|
sapply(names(de_qc_stats$samples_per_group), function(grp) {
|
||||||
|
sprintf(" %s: %d samples", grp, de_qc_stats$samples_per_group[[grp]])
|
||||||
|
}),
|
||||||
|
if (de_qc_stats$sample_size_warnings != "none") paste0(" WARNING: ", de_qc_stats$sample_size_warnings) else NULL,
|
||||||
|
"",
|
||||||
|
sprintf("Number of contrasts tested: %d", de_qc_stats$num_contrasts),
|
||||||
|
if (!is.null(de_qc_stats$multiple_testing_note)) paste0(" NOTE: ", de_qc_stats$multiple_testing_note) else NULL,
|
||||||
|
"",
|
||||||
|
"Per-Contrast Results:",
|
||||||
|
sapply(names(de_qc_stats$contrasts), function(cname) {
|
||||||
|
cqc <- de_qc_stats$contrasts[[cname]]
|
||||||
|
c(
|
||||||
|
"",
|
||||||
|
sprintf(" %s:", cname),
|
||||||
|
sprintf(" Samples: %d (%d vs %d)", cqc$n_samples, cqc$n_target, cqc$n_reference),
|
||||||
|
sprintf(" DGE significant: %d genes (FDR<0.05)", cqc$dge_significant_fdr05),
|
||||||
|
sprintf(" DTU status: %s", cqc$dtu_status),
|
||||||
|
if (cqc$dtu_status == "SUCCESS") sprintf(" DTU significant: %d genes", cqc$dtu_significant_genes) else NULL
|
||||||
|
)
|
||||||
|
}),
|
||||||
|
"",
|
||||||
|
"For detailed per-contrast statistics, see:",
|
||||||
|
" - <contrast>/contrast_qc_summary.txt",
|
||||||
|
" - <contrast>/results_dge.tsv",
|
||||||
|
" - <contrast>/results_dtu_gene.tsv",
|
||||||
|
""
|
||||||
|
)
|
||||||
|
writeLines(unlist(overall_summary), file.path(argv$out_dir, "de_overall_summary.txt"))
|
||||||
|
writeLines(capture.output(sessionInfo()), file.path(argv$out_dir, "session_info.txt"))
|
||||||
|
|
||||||
|
invisible(list(qc = de_qc_stats))
|
||||||
|
}
|
||||||
|
|
||||||
|
run_de_analysis_cli <- function(argv = commandArgs(trailingOnly = TRUE)) {
|
||||||
|
parsed <- argparser::parse_args(de_analysis_arg_parser(), argv = argv)
|
||||||
|
main_run_de_analysis(parsed)
|
||||||
|
}
|
||||||
16
bin/workflow_glue_r/load.R
Normal file
16
bin/workflow_glue_r/load.R
Normal file
@ -0,0 +1,16 @@
|
|||||||
|
workflow_glue_r_load <- function(pkg_dir = NULL, env = globalenv()) {
|
||||||
|
if (is.null(pkg_dir)) {
|
||||||
|
stop("pkg_dir must be provided when loading workflow_glue_r.", call. = FALSE)
|
||||||
|
}
|
||||||
|
|
||||||
|
r_dir <- file.path(pkg_dir, "R")
|
||||||
|
if (!dir.exists(r_dir)) {
|
||||||
|
stop(sprintf("R source directory not found: %s", r_dir), call. = FALSE)
|
||||||
|
}
|
||||||
|
|
||||||
|
for (path in sort(list.files(r_dir, pattern = "\\.[Rr]$", full.names = TRUE))) {
|
||||||
|
sys.source(path, envir = env)
|
||||||
|
}
|
||||||
|
|
||||||
|
invisible(pkg_dir)
|
||||||
|
}
|
||||||
36
bin/workflow_glue_r/tests/testthat.R
Normal file
36
bin/workflow_glue_r/tests/testthat.R
Normal file
@ -0,0 +1,36 @@
|
|||||||
|
args <- commandArgs(trailingOnly = FALSE)
|
||||||
|
file_arg <- grep("^--file=", args, value = TRUE)
|
||||||
|
if (length(file_arg) < 1) {
|
||||||
|
stop("Unable to determine testthat.R path", call. = FALSE)
|
||||||
|
}
|
||||||
|
|
||||||
|
test_dir <- dirname(normalizePath(sub("^--file=", "", file_arg[[1]])))
|
||||||
|
pkg_dir <- normalizePath(file.path(test_dir, ".."))
|
||||||
|
repo_root <- normalizePath(file.path(pkg_dir, "..", ".."))
|
||||||
|
|
||||||
|
source(file.path(pkg_dir, "load.R"))
|
||||||
|
workflow_glue_r_load(pkg_dir)
|
||||||
|
|
||||||
|
Sys.setenv(
|
||||||
|
WORKFLOW_GLUE_R_PACKAGE_DIR = pkg_dir,
|
||||||
|
WORKFLOW_GLUE_R_REPO_ROOT = repo_root,
|
||||||
|
TEST_DATA = Sys.getenv("TEST_DATA", unset = file.path(repo_root, "test_data"))
|
||||||
|
)
|
||||||
|
|
||||||
|
testthat_dir <- file.path(pkg_dir, "tests", "testthat")
|
||||||
|
test_files <- if (dir.exists(testthat_dir)) {
|
||||||
|
list.files(testthat_dir, pattern = "\\.[Rr]$", full.names = TRUE)
|
||||||
|
} else {
|
||||||
|
character(0)
|
||||||
|
}
|
||||||
|
|
||||||
|
if (length(test_files) < 1) {
|
||||||
|
message("No workflow-local R test files found; skipping.")
|
||||||
|
quit(save = "no", status = 0)
|
||||||
|
}
|
||||||
|
|
||||||
|
testthat::test_dir(
|
||||||
|
testthat_dir,
|
||||||
|
reporter = "summary",
|
||||||
|
stop_on_failure = TRUE
|
||||||
|
)
|
||||||
@ -60,7 +60,7 @@ params {
|
|||||||
"--sample_sheet 'wf-transcriptomes-demo/sample_sheet.csv'",
|
"--sample_sheet 'wf-transcriptomes-demo/sample_sheet.csv'",
|
||||||
]
|
]
|
||||||
common_sha = "sha21d552f9910c575766e5d465fcb7b52fefda4b79"
|
common_sha = "sha21d552f9910c575766e5d465fcb7b52fefda4b79"
|
||||||
container_sha = "shae31de64635251eca32c3f5619cb144a27f8233c6"
|
container_sha = "shaff0012055c9e1e71caf5b7857d717d17a3465a77"
|
||||||
pychopper_sha = "shaaaf20a5a0e76f9e18bad21af639a6b69e4a31a2f"
|
pychopper_sha = "shaaaf20a5a0e76f9e18bad21af639a6b69e4a31a2f"
|
||||||
sqanti_sha = "sha5bd775836492699e2537ebf846098eb117191d87"
|
sqanti_sha = "sha5bd775836492699e2537ebf846098eb117191d87"
|
||||||
agent = null
|
agent = null
|
||||||
|
|||||||
@ -37,7 +37,7 @@ process runDifferentialAnalysis {
|
|||||||
String covariates_arg = params.covariates ? "--covariates '${params.covariates}'" : ""
|
String covariates_arg = params.covariates ? "--covariates '${params.covariates}'" : ""
|
||||||
String reference_arg = params.reference_level ? "--reference_level '${params.reference_level}'" : ""
|
String reference_arg = params.reference_level ? "--reference_level '${params.reference_level}'" : ""
|
||||||
"""
|
"""
|
||||||
run_de_analysis.R \
|
supeRglue de_analysis \
|
||||||
--transcript_rds "${transcript_rds}" \
|
--transcript_rds "${transcript_rds}" \
|
||||||
--gene_rds "${gene_rds}" \
|
--gene_rds "${gene_rds}" \
|
||||||
--sample_sheet "${sample_sheet}" \
|
--sample_sheet "${sample_sheet}" \
|
||||||
|
|||||||
@ -179,7 +179,7 @@ process runJointBambu {
|
|||||||
String sample_sheet_arg = sample_sheet.name == OPTIONAL_FILE.name ? "" : "--sample_sheet ${sample_sheet}"
|
String sample_sheet_arg = sample_sheet.name == OPTIONAL_FILE.name ? "" : "--sample_sheet ${sample_sheet}"
|
||||||
String ndr_arg = params.ndr != null ? "--ndr ${params.ndr}" : ""
|
String ndr_arg = params.ndr != null ? "--ndr ${params.ndr}" : ""
|
||||||
"""
|
"""
|
||||||
run_bambu.R \
|
supeRglue bambu \
|
||||||
--bam_dir bams \
|
--bam_dir bams \
|
||||||
${sample_sheet_arg} \
|
${sample_sheet_arg} \
|
||||||
--annotation "${annotation}" \
|
--annotation "${annotation}" \
|
||||||
@ -210,7 +210,7 @@ process runPerSampleBambu {
|
|||||||
tuple val(meta), path("${meta.alias}/transcript_metadata.tsv"), emit: transcript_metadata
|
tuple val(meta), path("${meta.alias}/transcript_metadata.tsv"), emit: transcript_metadata
|
||||||
script:
|
script:
|
||||||
"""
|
"""
|
||||||
run_bambu.R \
|
supeRglue bambu \
|
||||||
--bam_path "${bam}" \
|
--bam_path "${bam}" \
|
||||||
--sample_alias "${meta.alias}" \
|
--sample_alias "${meta.alias}" \
|
||||||
--annotation "${annotation}" \
|
--annotation "${annotation}" \
|
||||||
|
|||||||
Loading…
Reference in New Issue
Block a user