Merge branch 'cw-7221' into 'dev'
[CW-7221] Fix various pathological cases See merge request epi2melabs/workflows/wf-transcriptomes!266
This commit is contained in:
Chris Wright
commit
e9b68d6222
@ -65,7 +65,7 @@ docker-run:
|
||||
parallel:
|
||||
matrix:
|
||||
- MATRIX_NAME: [
|
||||
"int_discover_dna", "int_fixed_rna",
|
||||
"int_discover_dna", "int_fixed_rna", "int_de_control_vs_control",
|
||||
"smoke_discover", "smoke_fixed", "smoke_direct_rna", "smoke_de",
|
||||
]
|
||||
rules:
|
||||
@ -77,7 +77,7 @@ docker-run:
|
||||
when: never
|
||||
|
||||
# Integration: larger discover-mode run on representative cDNA test bundle.
|
||||
- if: $MATRIX_NAME == "discover"
|
||||
- if: $MATRIX_NAME == "int_discover_dna"
|
||||
variables:
|
||||
NF_BEFORE_SCRIPT: "mkdir -p ${CI_PROJECT_NAME}/data/ && wget -nv -O ${CI_PROJECT_NAME}/data/wf-isoforms_test_data.tar.gz https://ont-exd-int-s3-euwst1-epi2me-labs.s3.amazonaws.com/wf-isoforms/wf-isoforms_test_data.tar.gz && tar -xzvf ${CI_PROJECT_NAME}/data/wf-isoforms_test_data.tar.gz -C ${CI_PROJECT_NAME}/data/ && wget -nv -O ${CI_PROJECT_NAME}/data/demo.nextflow.config https://ont-exd-int-s3-euwst1-epi2me-labs.s3.amazonaws.com/wf-isoforms/demo.nextflow.config;"
|
||||
NF_WORKFLOW_OPTS: "--fastq ${CI_PROJECT_NAME}/data/ERR6053095_chr20.fastq \
|
||||
@ -103,6 +103,37 @@ docker-run:
|
||||
gz_faidx,merge_transcriptomes,assemble_transcripts,decompress_annotation,decompress_ref,
|
||||
build_minimap_index,validate_ref_annotation,get_transcriptome,merge_gff_bundles,run_gffcompare,build_minimap_index,split_bam,decompress_transcriptome,preprocess_ref_transcriptome
|
||||
|
||||
# Integration: demo-data control-vs-control to test pathological DE/DTU failure modes and graceful failure handling.
|
||||
# See also unit tests in test_de_analysis.R
|
||||
- if: $MATRIX_NAME == "int_de_control_vs_control"
|
||||
variables:
|
||||
NF_BEFORE_SCRIPT: "\
|
||||
mkdir -p ${CI_PROJECT_NAME}/data/ \
|
||||
&& echo 'Downlading demo data bundle and config' \
|
||||
&& wget -nv -O ${CI_PROJECT_NAME}/data/differential_expression.tar.gz https://ont-exd-int-s3-euwst1-epi2me-labs.s3.amazonaws.com/wf-isoforms/differential_expression.tar.gz \
|
||||
&& tar -xzvf ${CI_PROJECT_NAME}/data/differential_expression.tar.gz -C ${CI_PROJECT_NAME}/data/ \
|
||||
&& wget -nv -O ${CI_PROJECT_NAME}/data/demo.nextflow.config https://ont-exd-int-s3-euwst1-epi2me-labs.s3.amazonaws.com/wf-isoforms/demo.nextflow.config \
|
||||
&& echo 'Simulating control-vs-control by copying each sample to a replicate with identical reads' \
|
||||
&& rm -rf ${CI_PROJECT_NAME}/data/differential_expression/differential_expression_fastq/barcode04 ${CI_PROJECT_NAME}/data/differential_expression/differential_expression_fastq/barcode05 ${CI_PROJECT_NAME}/data/differential_expression/differential_expression_fastq/barcode06 \
|
||||
&& cp -R ${CI_PROJECT_NAME}/data/differential_expression/differential_expression_fastq/barcode01 ${CI_PROJECT_NAME}/data/differential_expression/differential_expression_fastq/barcode04 \
|
||||
&& cp -R ${CI_PROJECT_NAME}/data/differential_expression/differential_expression_fastq/barcode02 ${CI_PROJECT_NAME}/data/differential_expression/differential_expression_fastq/barcode05 \
|
||||
&& cp -R ${CI_PROJECT_NAME}/data/differential_expression/differential_expression_fastq/barcode03 ${CI_PROJECT_NAME}/data/differential_expression/differential_expression_fastq/barcode06 \
|
||||
&& echo 'Finished data preparation step' \
|
||||
"
|
||||
NF_WORKFLOW_OPTS: "--fastq ${CI_PROJECT_NAME}/data/differential_expression/differential_expression_fastq \
|
||||
--de_analysis \
|
||||
--ref_genome ${CI_PROJECT_NAME}/data/differential_expression/hg38_chr20.fa \
|
||||
--ref_annotation ${CI_PROJECT_NAME}/data/differential_expression/gencode.v22.annotation.chr20.gtf \
|
||||
--direct_rna \
|
||||
--sample_sheet test_data/sample_sheet.csv \
|
||||
-c ${CI_PROJECT_NAME}/data/demo.nextflow.config"
|
||||
AFTER_NEXTFLOW_CMD: >
|
||||
test -f ${CI_PROJECT_NAME}/de_analysis/de_qc_stats.json &&
|
||||
test -f ${CI_PROJECT_NAME}/de_analysis/condition_treated_vs_control/results_dge.tsv &&
|
||||
test -f ${CI_PROJECT_NAME}/de_analysis/condition_treated_vs_control/results_dtu_transcript.tsv &&
|
||||
jq -e '.contrasts["condition_treated_vs_control"].dge_status | IN("SUCCESS", "FAILED")' ${CI_PROJECT_NAME}/de_analysis/de_qc_stats.json >/dev/null &&
|
||||
jq -e '.contrasts["condition_treated_vs_control"].dtu_status | IN("SUCCESS", "FAILED")' ${CI_PROJECT_NAME}/de_analysis/de_qc_stats.json >/dev/null
|
||||
|
||||
# Smoke: quick discover-mode sanity check for core cohort and per-sample outputs.
|
||||
- if: $MATRIX_NAME == "smoke_discover"
|
||||
variables:
|
||||
|
||||
@ -177,6 +177,7 @@ def _collect_de_method_rows(de_qc):
|
||||
deseq2_gene_wise = []
|
||||
dexseq_gene_wise = []
|
||||
dexseq_covariate_drops = []
|
||||
failed_dge = []
|
||||
|
||||
for contrast_name, contrast_data in de_qc.get("contrasts", {}).items():
|
||||
fallback = contrast_data.get("deseq2_dispersion_fallback") or {}
|
||||
@ -201,6 +202,9 @@ def _collect_de_method_rows(de_qc):
|
||||
if dropped_covariates:
|
||||
dexseq_covariate_drops.append((contrast_name, dropped_covariates))
|
||||
|
||||
if contrast_data.get("dge_status") == "FAILED":
|
||||
failed_dge.append(contrast_name)
|
||||
|
||||
rows.append(
|
||||
{
|
||||
"Contrast": contrast_name,
|
||||
@ -215,11 +219,12 @@ def _collect_de_method_rows(de_qc):
|
||||
"DEXSeq covariates dropped": (
|
||||
", ".join(dropped_covariates) if dropped_covariates else "none"
|
||||
),
|
||||
"DGE status": contrast_data.get("dge_status", "N/A"),
|
||||
"DTU status": contrast_data.get("dtu_status", "N/A"),
|
||||
}
|
||||
)
|
||||
|
||||
return rows, deseq2_gene_wise, dexseq_gene_wise, dexseq_covariate_drops
|
||||
return rows, deseq2_gene_wise, dexseq_gene_wise, dexseq_covariate_drops, failed_dge
|
||||
|
||||
|
||||
def main(args):
|
||||
@ -529,6 +534,7 @@ def main(args):
|
||||
deseq2_gene_wise,
|
||||
dexseq_gene_wise,
|
||||
dexseq_covariate_drops,
|
||||
failed_dge,
|
||||
) = _collect_de_method_rows(de_qc)
|
||||
|
||||
if sample_size_warnings:
|
||||
@ -578,6 +584,17 @@ def main(args):
|
||||
)
|
||||
has_warnings = True
|
||||
|
||||
if failed_dge:
|
||||
_create_warning_banner(
|
||||
"DGE Analysis Failed: "
|
||||
f"{len(failed_dge)} contrast(s) could not complete "
|
||||
"DGE testing. "
|
||||
f"Affected: {', '.join(failed_dge)}. "
|
||||
"See DGE_ANALYSIS_FAILED.txt files for details.",
|
||||
level="danger",
|
||||
)
|
||||
has_warnings = True
|
||||
|
||||
# Check for failed DTU analyses
|
||||
failed_dtu = []
|
||||
for contrast_name, contrast_data in de_qc.get(
|
||||
@ -681,16 +698,25 @@ def main(args):
|
||||
f"vs "
|
||||
f"{contrast_data.get('n_reference', 0)}"
|
||||
),
|
||||
"DGE Status": contrast_data.get(
|
||||
"dge_status", "N/A"
|
||||
),
|
||||
"DGE Significant (FDR<0.05)": (
|
||||
contrast_data.get(
|
||||
"dge_significant_fdr05", 0
|
||||
)
|
||||
if contrast_data.get("dge_status") == "SUCCESS"
|
||||
else "N/A"
|
||||
),
|
||||
"DGE Up": contrast_data.get(
|
||||
"dge_upregulated", 0
|
||||
"DGE Up": (
|
||||
contrast_data.get("dge_upregulated", 0)
|
||||
if contrast_data.get("dge_status") == "SUCCESS"
|
||||
else "N/A"
|
||||
),
|
||||
"DGE Down": contrast_data.get(
|
||||
"dge_downregulated", 0
|
||||
"DGE Down": (
|
||||
contrast_data.get("dge_downregulated", 0)
|
||||
if contrast_data.get("dge_status") == "SUCCESS"
|
||||
else "N/A"
|
||||
),
|
||||
"DTU Status": contrast_data.get(
|
||||
"dtu_status", "N/A"
|
||||
@ -751,6 +777,15 @@ def main(args):
|
||||
),
|
||||
}
|
||||
)
|
||||
if failed_dge:
|
||||
warnings_data.append(
|
||||
{
|
||||
"Warning Type": "DGE Failure",
|
||||
"Details": (
|
||||
f"{len(failed_dge)} contrasts failed"
|
||||
),
|
||||
}
|
||||
)
|
||||
|
||||
warnings_df = pd.DataFrame(warnings_data)
|
||||
DataTable.from_pandas(
|
||||
@ -766,7 +801,19 @@ def main(args):
|
||||
# Check for contrast-specific warnings
|
||||
if de_qc and contrast in de_qc.get("contrasts", {}):
|
||||
contrast_data = de_qc["contrasts"][contrast]
|
||||
if contrast_data.get("dtu_power_warning"):
|
||||
if contrast_data.get("dge_status") == "FAILED":
|
||||
_create_warning_banner(
|
||||
"DGE analysis failed for this contrast. "
|
||||
f"See {contrast}/"
|
||||
"DGE_ANALYSIS_FAILED.txt for detailed "
|
||||
"explanation.",
|
||||
level="danger",
|
||||
)
|
||||
p(
|
||||
"Empty results indicate analysis failure, "
|
||||
"not 'no DGE detected'."
|
||||
)
|
||||
elif contrast_data.get("dtu_power_warning"):
|
||||
with div(
|
||||
style=(
|
||||
"padding: 10px; margin-bottom: 10px; "
|
||||
@ -801,6 +848,9 @@ def main(args):
|
||||
"explanation.",
|
||||
level="danger",
|
||||
)
|
||||
failure_hint = contrast_data.get("dtu_failure_hint")
|
||||
if failure_hint:
|
||||
p(f"Probable cause: {failure_hint}")
|
||||
p(
|
||||
"Empty results indicate analysis failure, "
|
||||
"not 'no DTU detected'."
|
||||
|
||||
@ -201,16 +201,50 @@ de_run_deseq_with_fallback <- function(
|
||||
)
|
||||
|
||||
de_out <- tryCatch(
|
||||
DESeq2::DESeq(dds, quiet = TRUE, sfType = sf_type),
|
||||
list(
|
||||
dds = DESeq2::DESeq(dds, quiet = TRUE, sfType = sf_type),
|
||||
deseq2_dispersion_fallback = fallback_info
|
||||
),
|
||||
error = function(err) {
|
||||
if (!grepl(
|
||||
"all gene-wise dispersion estimates are within 2 orders of magnitude",
|
||||
conditionMessage(err),
|
||||
msg <- conditionMessage(err)
|
||||
zero_geometric_means <- grepl(
|
||||
"every gene contains at least one zero",
|
||||
msg,
|
||||
fixed = TRUE
|
||||
)) {
|
||||
)
|
||||
|
||||
# Size factor estimation failure: all genes have at least one zero.
|
||||
# de_choose_size_factor_type should prevent this proactively, but as
|
||||
# a safety net switch to poscounts so the gene-wise fallback below
|
||||
# can call estimateSizeFactors without failing again.
|
||||
fallback_sf_type <- if (zero_geometric_means) "poscounts" else sf_type
|
||||
|
||||
recoverable <- grepl(
|
||||
"all gene-wise dispersion estimates are within 2 orders of magnitude",
|
||||
msg,
|
||||
fixed = TRUE
|
||||
) || grepl(
|
||||
"newsplit: out of vertex space",
|
||||
msg,
|
||||
fixed = TRUE
|
||||
) || grepl(
|
||||
"every gene contains at least one zero",
|
||||
msg,
|
||||
fixed = TRUE
|
||||
)
|
||||
if (!recoverable) {
|
||||
stop(err)
|
||||
}
|
||||
|
||||
if (zero_geometric_means) {
|
||||
warning(
|
||||
"STATISTICAL POWER REDUCED: DESeq2 default size-factor estimation failed for ",
|
||||
contrast_name,
|
||||
" because every gene contains at least one zero.\n",
|
||||
"Switching to sfType='poscounts' and using gene-wise dispersion estimates.\n",
|
||||
"Results will have reduced power and wider confidence intervals."
|
||||
)
|
||||
} else {
|
||||
warning(
|
||||
"STATISTICAL POWER REDUCED: DESeq2 dispersion estimation failed for ",
|
||||
contrast_name,
|
||||
@ -222,8 +256,9 @@ de_run_deseq_with_fallback <- function(
|
||||
"Falling back to gene-wise dispersion (no information sharing).\n",
|
||||
"Results will have reduced power and wider confidence intervals."
|
||||
)
|
||||
}
|
||||
|
||||
dds <- DESeq2::estimateSizeFactors(dds, type = sf_type)
|
||||
dds <- DESeq2::estimateSizeFactors(dds, type = fallback_sf_type)
|
||||
dds <- DESeq2::estimateDispersionsGeneEst(dds)
|
||||
dispersions_setter <- get("dispersions<-", envir = asNamespace("DESeq2"))
|
||||
dds <- dispersions_setter(dds, value = S4Vectors::mcols(dds)$dispGeneEst)
|
||||
@ -241,17 +276,35 @@ de_run_deseq_with_fallback <- function(
|
||||
}
|
||||
|
||||
diag_content <- c(
|
||||
"DESeq2 Dispersion Estimation Fallback Applied",
|
||||
"==============================================",
|
||||
if (zero_geometric_means) {
|
||||
"DESeq2 Poscounts and Gene-wise Fallback Applied"
|
||||
} else {
|
||||
"DESeq2 Dispersion Estimation Fallback Applied"
|
||||
},
|
||||
if (zero_geometric_means) {
|
||||
"============================================="
|
||||
} else {
|
||||
"=============================================="
|
||||
},
|
||||
"",
|
||||
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("Size factor method: %s", fallback_sf_type),
|
||||
dispersion_range,
|
||||
"",
|
||||
"WHAT HAPPENED:",
|
||||
" Curve fitting failed. Using gene-wise dispersion estimates.",
|
||||
if (zero_geometric_means) {
|
||||
" Every gene contained at least one zero, so DESeq2's default size-factor estimator was undefined."
|
||||
} else {
|
||||
" Curve fitting failed. Using gene-wise dispersion estimates."
|
||||
},
|
||||
if (zero_geometric_means) {
|
||||
" The analysis switched to sfType='poscounts' and used gene-wise dispersion estimates."
|
||||
} else {
|
||||
NULL
|
||||
},
|
||||
"",
|
||||
"IMPLICATIONS:",
|
||||
" - No information sharing across genes",
|
||||
@ -260,15 +313,35 @@ de_run_deseq_with_fallback <- function(
|
||||
" - More conservative results (fewer discoveries)",
|
||||
"",
|
||||
"LIKELY CAUSES:",
|
||||
if (zero_geometric_means) {
|
||||
c(
|
||||
" 1. Sparse count matrices with many zeros",
|
||||
" 2. Low-expression features dominating the design",
|
||||
" 3. Small or weakly separated groups"
|
||||
)
|
||||
} else {
|
||||
c(
|
||||
" 1. Too few replicates (recommend n>=3 per group)",
|
||||
" 2. High biological variability",
|
||||
" 3. Poor data quality or outlier samples",
|
||||
" 3. Poor data quality or outlier samples"
|
||||
)
|
||||
},
|
||||
"",
|
||||
"RECOMMENDATIONS:",
|
||||
if (zero_geometric_means) {
|
||||
c(
|
||||
" - Check whether the contrast is dominated by zeros or nearly identical samples",
|
||||
" - Add more biological replicates if possible",
|
||||
" - Consider filtering very low-count genes more stringently"
|
||||
)
|
||||
} else {
|
||||
c(
|
||||
" - Add more biological replicates if possible",
|
||||
" - Check sample quality metrics",
|
||||
" - Consider filtering low-count genes more stringently"
|
||||
)
|
||||
}
|
||||
)
|
||||
|
||||
diag_file <- file.path(
|
||||
out_dir,
|
||||
@ -278,18 +351,19 @@ de_run_deseq_with_fallback <- function(
|
||||
)
|
||||
)
|
||||
writeLines(diag_content, diag_file)
|
||||
fallback_info <<- list(
|
||||
list(
|
||||
dds = DESeq2::nbinomWaldTest(dds),
|
||||
deseq2_dispersion_fallback = list(
|
||||
applied = TRUE,
|
||||
method_used = "gene-wise",
|
||||
reason = conditionMessage(err),
|
||||
diagnostic_file = basename(diag_file),
|
||||
size_factor_type = sf_type
|
||||
size_factor_type = fallback_sf_type
|
||||
)
|
||||
)
|
||||
|
||||
DESeq2::nbinomWaldTest(dds)
|
||||
}
|
||||
)
|
||||
list(dds = de_out, deseq2_dispersion_fallback = fallback_info)
|
||||
de_out
|
||||
}
|
||||
|
||||
de_estimate_dispersions_with_fallback <- function(
|
||||
@ -297,6 +371,15 @@ de_estimate_dispersions_with_fallback <- function(
|
||||
context_label,
|
||||
allow_gene_est = TRUE
|
||||
) {
|
||||
de_dispersion_recoverable <- function(msg) {
|
||||
grepl(
|
||||
"all gene-wise dispersion estimates are within 2 orders of magnitude",
|
||||
msg, fixed = TRUE
|
||||
) || grepl(
|
||||
"newsplit: out of vertex space",
|
||||
msg, fixed = TRUE
|
||||
)
|
||||
}
|
||||
tryCatch(
|
||||
list(
|
||||
object = DESeq2::estimateDispersions(object),
|
||||
@ -304,11 +387,7 @@ de_estimate_dispersions_with_fallback <- function(
|
||||
fallback_applied = FALSE
|
||||
),
|
||||
error = function(err) {
|
||||
if (!grepl(
|
||||
"all gene-wise dispersion estimates are within 2 orders of magnitude",
|
||||
conditionMessage(err),
|
||||
fixed = TRUE
|
||||
)) {
|
||||
if (!de_dispersion_recoverable(conditionMessage(err))) {
|
||||
stop(err)
|
||||
}
|
||||
|
||||
@ -323,11 +402,7 @@ de_estimate_dispersions_with_fallback <- function(
|
||||
fallback_applied = TRUE
|
||||
),
|
||||
error = function(local_err) {
|
||||
if (!grepl(
|
||||
"all gene-wise dispersion estimates are within 2 orders of magnitude",
|
||||
conditionMessage(local_err),
|
||||
fixed = TRUE
|
||||
)) {
|
||||
if (!de_dispersion_recoverable(conditionMessage(local_err))) {
|
||||
stop(local_err)
|
||||
}
|
||||
|
||||
@ -342,11 +417,7 @@ de_estimate_dispersions_with_fallback <- function(
|
||||
fallback_applied = TRUE
|
||||
),
|
||||
error = function(mean_err) {
|
||||
if (!grepl(
|
||||
"all gene-wise dispersion estimates are within 2 orders of magnitude",
|
||||
conditionMessage(mean_err),
|
||||
fixed = TRUE
|
||||
)) {
|
||||
if (!de_dispersion_recoverable(conditionMessage(mean_err))) {
|
||||
stop(mean_err)
|
||||
}
|
||||
if (!allow_gene_est) {
|
||||
@ -389,9 +460,41 @@ de_is_recoverable_dexseq_error <- function(message_text) {
|
||||
"replacement has 1 row, data has 0",
|
||||
message_text,
|
||||
fixed = TRUE
|
||||
) || grepl(
|
||||
"newsplit: out of vertex space",
|
||||
message_text,
|
||||
fixed = TRUE
|
||||
) || grepl(
|
||||
"every gene contains at least one zero",
|
||||
message_text,
|
||||
fixed = TRUE
|
||||
)
|
||||
}
|
||||
|
||||
de_dexseq_failure_hint <- function(message_text) {
|
||||
if (grepl("nb of cols in 'assay'", message_text, fixed = TRUE)) {
|
||||
return(
|
||||
paste(
|
||||
"Probable cause: degenerate or near-identical sample groups can",
|
||||
"drive DEXSeq into an internal object-shape mismatch during model",
|
||||
"fitting. This is often seen when the contrast has little real",
|
||||
"separation or too few informative non-zero features."
|
||||
)
|
||||
)
|
||||
}
|
||||
if (grepl("'x' must be an array of at least two dimensions", message_text, fixed = TRUE)) {
|
||||
return(
|
||||
paste(
|
||||
"Probable cause: degenerate or near-identical sample groups can",
|
||||
"leave DEXSeq with too little informative structure for downstream",
|
||||
"gene-level DTU summarisation. This is often seen when the contrast",
|
||||
"has little real separation or too few informative non-zero features."
|
||||
)
|
||||
)
|
||||
}
|
||||
NULL
|
||||
}
|
||||
|
||||
de_write_placeholder_pdf <- function(path, label) {
|
||||
grDevices::pdf(path)
|
||||
graphics::plot.new()
|
||||
@ -399,6 +502,50 @@ de_write_placeholder_pdf <- function(path, label) {
|
||||
grDevices::dev.off()
|
||||
}
|
||||
|
||||
de_write_dtu_failure_outputs <- function(
|
||||
contrast_dir,
|
||||
target_level,
|
||||
reference_level,
|
||||
contrast_samples,
|
||||
condition_column,
|
||||
tx_counts,
|
||||
message_text,
|
||||
failure_hint
|
||||
) {
|
||||
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[[condition_column]] == target_level),
|
||||
target_level,
|
||||
sum(contrast_samples[[condition_column]] == reference_level),
|
||||
reference_level
|
||||
),
|
||||
sprintf("Transcripts: %d", nrow(tx_counts)),
|
||||
"",
|
||||
"ERROR MESSAGE:",
|
||||
sprintf(" %s", message_text),
|
||||
if (!is.null(failure_hint)) {
|
||||
c(
|
||||
"",
|
||||
"PROBABLE CAUSE:",
|
||||
sprintf(" %s", failure_hint)
|
||||
)
|
||||
} else {
|
||||
NULL
|
||||
},
|
||||
"",
|
||||
"DTU RESULTS CANNOT BE INTERPRETED",
|
||||
""
|
||||
)
|
||||
writeLines(failure_content, file.path(contrast_dir, "DTU_ANALYSIS_FAILED.txt"))
|
||||
}
|
||||
|
||||
de_run_deseq2_result <- function(
|
||||
count_mat,
|
||||
coldata,
|
||||
@ -453,9 +600,7 @@ de_run_dexseq_result <- function(
|
||||
for (covariate in covariates) {
|
||||
coldata[[covariate]] <- factor(coldata[[covariate]])
|
||||
}
|
||||
dropped_covariates <- character(0)
|
||||
|
||||
run_inner <- function(active_covariates) {
|
||||
run_inner <- function(active_covariates, dropped_covariates = character(0)) {
|
||||
covariate_exon_terms <- if (length(active_covariates) > 0) {
|
||||
paste0(active_covariates, ":exon")
|
||||
} else {
|
||||
@ -498,7 +643,8 @@ de_run_dexseq_result <- function(
|
||||
dxd = dxd,
|
||||
dxr = dxr,
|
||||
dexseq_dispersion_method = dispersion_method,
|
||||
dexseq_size_factor_type = dexseq_sf_type
|
||||
dexseq_size_factor_type = dexseq_sf_type,
|
||||
dexseq_covariates_dropped = dropped_covariates
|
||||
)
|
||||
}, error = function(err) {
|
||||
if (length(active_covariates) == 0 || !grepl(
|
||||
@ -511,21 +657,23 @@ de_run_dexseq_result <- function(
|
||||
|
||||
dropped_covariate <- tail(active_covariates, 1)
|
||||
kept_covariates <- head(active_covariates, -1)
|
||||
dropped_covariates <<- c(dropped_covariates, dropped_covariate)
|
||||
message(
|
||||
"DEXSeq design was not full rank with covariate '",
|
||||
dropped_covariate,
|
||||
"'; retrying without it."
|
||||
)
|
||||
run_inner(kept_covariates)
|
||||
run_inner(
|
||||
kept_covariates,
|
||||
c(dropped_covariates, dropped_covariate)
|
||||
)
|
||||
})
|
||||
}
|
||||
|
||||
result <- run_inner(covariates)
|
||||
result$dexseq_covariates_dropped <- dropped_covariates
|
||||
result
|
||||
run_inner(covariates)
|
||||
}
|
||||
|
||||
de_dge_columns <- c("GENEID", "baseMean", "log2FoldChange", "lfcSE", "stat", "pvalue", "padj")
|
||||
|
||||
de_dtu_transcript_columns <- c(
|
||||
"featureID",
|
||||
"groupID",
|
||||
@ -557,6 +705,81 @@ de_extract_dtu_transcript_table <- function(dex_df, contrast_name) {
|
||||
workflow_glue_r_normalise_tsv_df(tx_dtu)
|
||||
}
|
||||
|
||||
de_postprocess_dexseq_result <- function(
|
||||
dex_res,
|
||||
contrast_dir,
|
||||
target_level,
|
||||
reference_level,
|
||||
contrast_samples,
|
||||
condition_column,
|
||||
tx_counts,
|
||||
contrast_name,
|
||||
per_gene_qvalue_fn = DEXSeq::perGeneQValue,
|
||||
plot_writer = function(dex_res, contrast_dir) {
|
||||
grDevices::pdf(file.path(contrast_dir, "results_dtu.pdf"))
|
||||
DESeq2::plotMA(dex_res$dxr, cex = 0.8, alpha = 0.05)
|
||||
DESeq2::plotDispEsts(dex_res$dxd)
|
||||
grDevices::dev.off()
|
||||
}
|
||||
) {
|
||||
tryCatch(
|
||||
{
|
||||
dex_df <- as.data.frame(dex_res$dxr)
|
||||
dex_df <- workflow_glue_r_normalise_tsv_df(dex_df)
|
||||
tx_dtu <- de_extract_dtu_transcript_table(dex_df, contrast_name)
|
||||
gene_q <- per_gene_qvalue_fn(dex_res$dxr)
|
||||
gene_dtu <- data.frame(
|
||||
GENEID = names(gene_q),
|
||||
qval = unname(gene_q),
|
||||
row.names = NULL
|
||||
)
|
||||
plot_writer(dex_res, contrast_dir)
|
||||
list(
|
||||
ok = TRUE,
|
||||
dex_df = dex_df,
|
||||
tx_dtu = tx_dtu,
|
||||
gene_dtu = gene_dtu,
|
||||
failure_hint = NULL
|
||||
)
|
||||
},
|
||||
error = function(err) {
|
||||
message_text <- conditionMessage(err)
|
||||
failure_hint <- de_dexseq_failure_hint(message_text)
|
||||
|
||||
warning(
|
||||
"DEXSeq post-processing failed for contrast ",
|
||||
target_level,
|
||||
" vs ",
|
||||
reference_level,
|
||||
if (!is.null(failure_hint)) {
|
||||
paste0("\nProbable cause: ", failure_hint)
|
||||
} else {
|
||||
""
|
||||
},
|
||||
"\nError: ",
|
||||
message_text
|
||||
)
|
||||
de_write_dtu_failure_outputs(
|
||||
contrast_dir = contrast_dir,
|
||||
target_level = target_level,
|
||||
reference_level = reference_level,
|
||||
contrast_samples = contrast_samples,
|
||||
condition_column = condition_column,
|
||||
tx_counts = tx_counts,
|
||||
message_text = message_text,
|
||||
failure_hint = failure_hint
|
||||
)
|
||||
list(
|
||||
ok = FALSE,
|
||||
dex_df = NULL,
|
||||
tx_dtu = NULL,
|
||||
gene_dtu = NULL,
|
||||
failure_hint = failure_hint
|
||||
)
|
||||
}
|
||||
)
|
||||
}
|
||||
|
||||
main_run_de_analysis <- function(args) {
|
||||
set.seed(42)
|
||||
dir.create(args$out_dir, showWarnings = FALSE, recursive = TRUE)
|
||||
@ -693,7 +916,9 @@ main_run_de_analysis <- function(args) {
|
||||
),
|
||||
dexseq_size_factor_method = "ratio",
|
||||
dexseq_dispersion_method = "parametric",
|
||||
dexseq_covariates_dropped = list()
|
||||
dexseq_covariates_dropped = list(),
|
||||
dtu_failure_hint = NULL,
|
||||
dge_status = "SUCCESS"
|
||||
)
|
||||
|
||||
if (nrow(contrast_samples) < 6) {
|
||||
@ -709,7 +934,8 @@ main_run_de_analysis <- function(args) {
|
||||
contrast_qc$genes_tested <- nrow(gene_counts)
|
||||
contrast_qc$transcripts_tested <- nrow(tx_counts)
|
||||
|
||||
dge_run <- de_run_deseq2_result(
|
||||
dge_run <- tryCatch(
|
||||
de_run_deseq2_result(
|
||||
gene_counts,
|
||||
contrast_samples,
|
||||
target_level,
|
||||
@ -718,7 +944,53 @@ main_run_de_analysis <- function(args) {
|
||||
covariates,
|
||||
args$out_dir,
|
||||
contrast_name
|
||||
),
|
||||
error = function(err) {
|
||||
warning(
|
||||
"DESeq2 DGE failed for contrast ",
|
||||
target_level, " vs ", reference_level,
|
||||
"\nError: ", conditionMessage(err)
|
||||
)
|
||||
failure_content <- c(
|
||||
"DGE 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[[args$condition_column]] == target_level),
|
||||
target_level,
|
||||
sum(contrast_samples[[args$condition_column]] == reference_level),
|
||||
reference_level
|
||||
),
|
||||
sprintf("Genes: %d", nrow(gene_counts)),
|
||||
"",
|
||||
"ERROR MESSAGE:",
|
||||
sprintf(" %s", conditionMessage(err)),
|
||||
"",
|
||||
"DGE RESULTS CANNOT BE INTERPRETED",
|
||||
""
|
||||
)
|
||||
writeLines(failure_content, file.path(contrast_dir, "DGE_ANALYSIS_FAILED.txt"))
|
||||
NULL
|
||||
}
|
||||
)
|
||||
|
||||
if (is.null(dge_run)) {
|
||||
dge_res <- workflow_glue_r_empty_tsv(de_dge_columns)
|
||||
contrast_qc$dge_status <- "FAILED"
|
||||
contrast_qc$dge_total_genes <- 0
|
||||
contrast_qc$dge_significant_fdr05 <- 0
|
||||
contrast_qc$dge_significant_fdr01 <- 0
|
||||
contrast_qc$dge_upregulated <- 0
|
||||
contrast_qc$dge_downregulated <- 0
|
||||
de_write_placeholder_pdf(
|
||||
file.path(contrast_dir, "results_dge.pdf"),
|
||||
"DESeq2 did not converge for this contrast.\nSee DGE_ANALYSIS_FAILED.txt for details."
|
||||
)
|
||||
} else {
|
||||
if (!is.null(dge_run$deseq2_dispersion_fallback)) {
|
||||
fallback <- dge_run$deseq2_dispersion_fallback
|
||||
fallback_applied <- isTRUE(fallback$applied)
|
||||
@ -740,7 +1012,7 @@ main_run_de_analysis <- function(args) {
|
||||
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_status <- "SUCCESS"
|
||||
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)
|
||||
@ -752,6 +1024,10 @@ main_run_de_analysis <- function(args) {
|
||||
dge_res$padj < 0.05 & dge_res$log2FoldChange < 0,
|
||||
na.rm = TRUE
|
||||
)
|
||||
grDevices::pdf(file.path(contrast_dir, "results_dge.pdf"))
|
||||
DESeq2::plotMA(dge_run$result)
|
||||
grDevices::dev.off()
|
||||
}
|
||||
|
||||
utils::write.table(
|
||||
dge_res[order(dge_res$padj), ],
|
||||
@ -761,60 +1037,63 @@ main_run_de_analysis <- function(args) {
|
||||
row.names = FALSE
|
||||
)
|
||||
|
||||
grDevices::pdf(file.path(contrast_dir, "results_dge.pdf"))
|
||||
DESeq2::plotMA(dge_run$result)
|
||||
grDevices::dev.off()
|
||||
|
||||
dex_res <- tryCatch(
|
||||
de_run_dexseq_result(
|
||||
dex_run <- tryCatch(
|
||||
list(
|
||||
ok = TRUE,
|
||||
result = de_run_dexseq_result(
|
||||
tx_counts,
|
||||
tx_meta,
|
||||
contrast_samples,
|
||||
args$condition_column,
|
||||
covariates
|
||||
),
|
||||
failure_hint = NULL
|
||||
),
|
||||
error = function(err) {
|
||||
message_text <- conditionMessage(err)
|
||||
if (!de_is_recoverable_dexseq_error(message_text)) {
|
||||
stop(err)
|
||||
}
|
||||
is_known_failure <- de_is_recoverable_dexseq_error(message_text)
|
||||
failure_hint <- de_dexseq_failure_hint(message_text)
|
||||
|
||||
warning(
|
||||
"DEXSeq failed for contrast ",
|
||||
target_level,
|
||||
" vs ",
|
||||
reference_level,
|
||||
if (!is_known_failure) {
|
||||
"\nThis error did not match a known statistical fallback pattern."
|
||||
} else {
|
||||
""
|
||||
},
|
||||
if (!is.null(failure_hint)) {
|
||||
paste0("\nProbable cause: ", failure_hint)
|
||||
} else {
|
||||
""
|
||||
},
|
||||
"\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[[args$condition_column]] == target_level),
|
||||
target_level,
|
||||
sum(contrast_samples[[args$condition_column]] == reference_level),
|
||||
reference_level
|
||||
),
|
||||
sprintf("Transcripts: %d", nrow(tx_counts)),
|
||||
"",
|
||||
"ERROR MESSAGE:",
|
||||
sprintf(" %s", message_text),
|
||||
"",
|
||||
"DTU RESULTS CANNOT BE INTERPRETED",
|
||||
""
|
||||
de_write_dtu_failure_outputs(
|
||||
contrast_dir = contrast_dir,
|
||||
target_level = target_level,
|
||||
reference_level = reference_level,
|
||||
contrast_samples = contrast_samples,
|
||||
condition_column = args$condition_column,
|
||||
tx_counts = tx_counts,
|
||||
message_text = message_text,
|
||||
failure_hint = failure_hint
|
||||
)
|
||||
writeLines(failure_content, file.path(contrast_dir, "DTU_ANALYSIS_FAILED.txt"))
|
||||
NULL
|
||||
list(ok = FALSE, result = NULL, failure_hint = failure_hint)
|
||||
}
|
||||
)
|
||||
|
||||
if (!isTRUE(dex_run$ok)) {
|
||||
dex_res <- NULL
|
||||
dtu_failure_hint <- dex_run$failure_hint
|
||||
} else {
|
||||
dex_res <- dex_run$result
|
||||
dtu_failure_hint <- NULL
|
||||
}
|
||||
|
||||
if (is.null(dex_res)) {
|
||||
dex_df <- workflow_glue_r_empty_tsv(de_dtu_transcript_columns)
|
||||
tx_dtu <- dex_df
|
||||
@ -824,6 +1103,7 @@ main_run_de_analysis <- function(args) {
|
||||
"DEXSeq did not converge for this contrast.\nSee DTU_ANALYSIS_FAILED.txt for details."
|
||||
)
|
||||
contrast_qc$dtu_status <- "FAILED"
|
||||
contrast_qc$dtu_failure_hint <- dtu_failure_hint
|
||||
contrast_qc$dtu_significant_transcripts <- 0
|
||||
contrast_qc$dtu_significant_genes <- 0
|
||||
} else {
|
||||
@ -836,28 +1116,39 @@ main_run_de_analysis <- function(args) {
|
||||
if (!is.null(dex_res$dexseq_covariates_dropped)) {
|
||||
contrast_qc$dexseq_covariates_dropped <- as.list(dex_res$dexseq_covariates_dropped)
|
||||
}
|
||||
dex_df <- as.data.frame(dex_res$dxr)
|
||||
dex_df <- workflow_glue_r_normalise_tsv_df(dex_df)
|
||||
tx_dtu <- de_extract_dtu_transcript_table(
|
||||
dex_df,
|
||||
paste(target_level, reference_level, sep = "_")
|
||||
|
||||
postprocess_run <- de_postprocess_dexseq_result(
|
||||
dex_res = dex_res,
|
||||
contrast_dir = contrast_dir,
|
||||
target_level = target_level,
|
||||
reference_level = reference_level,
|
||||
contrast_samples = contrast_samples,
|
||||
condition_column = args$condition_column,
|
||||
tx_counts = tx_counts,
|
||||
contrast_name = paste(target_level, reference_level, sep = "_")
|
||||
)
|
||||
|
||||
gene_q <- DEXSeq::perGeneQValue(dex_res$dxr)
|
||||
gene_dtu <- data.frame(
|
||||
GENEID = names(gene_q),
|
||||
qval = unname(gene_q),
|
||||
row.names = NULL
|
||||
if (!isTRUE(postprocess_run$ok)) {
|
||||
dex_df <- workflow_glue_r_empty_tsv(de_dtu_transcript_columns)
|
||||
tx_dtu <- dex_df
|
||||
gene_dtu <- workflow_glue_r_empty_tsv(c("GENEID", "qval"))
|
||||
de_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_failure_hint <- postprocess_run$failure_hint
|
||||
contrast_qc$dtu_significant_transcripts <- 0
|
||||
contrast_qc$dtu_significant_genes <- 0
|
||||
} else {
|
||||
dex_df <- postprocess_run$dex_df
|
||||
tx_dtu <- postprocess_run$tx_dtu
|
||||
gene_dtu <- postprocess_run$gene_dtu
|
||||
|
||||
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)
|
||||
|
||||
grDevices::pdf(file.path(contrast_dir, "results_dtu.pdf"))
|
||||
DESeq2::plotMA(dex_res$dxr, cex = 0.8, alpha = 0.05)
|
||||
DESeq2::plotDispEsts(dex_res$dxd)
|
||||
grDevices::dev.off()
|
||||
}
|
||||
}
|
||||
|
||||
utils::write.table(
|
||||
@ -899,12 +1190,25 @@ main_run_de_analysis <- function(args) {
|
||||
sprintf(" Total samples: %d", contrast_qc$n_samples),
|
||||
"",
|
||||
"DGE Results:",
|
||||
sprintf(" Status: %s", contrast_qc$dge_status),
|
||||
sprintf(" Genes tested: %d", contrast_qc$genes_tested),
|
||||
sprintf(" Size factor method: %s", contrast_qc$deseq2_size_factor_method),
|
||||
if (contrast_qc$dge_status == "SUCCESS") {
|
||||
c(
|
||||
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),
|
||||
sprintf(" Downregulated: %d", contrast_qc$dge_downregulated)
|
||||
)
|
||||
} else {
|
||||
c(
|
||||
" Significant (FDR < 0.05): N/A",
|
||||
" Significant (FDR < 0.01): N/A",
|
||||
" Upregulated: N/A",
|
||||
" Downregulated: N/A",
|
||||
" See DGE_ANALYSIS_FAILED.txt for details"
|
||||
)
|
||||
},
|
||||
"",
|
||||
"DTU Results:",
|
||||
sprintf(" Status: %s", contrast_qc$dtu_status),
|
||||
@ -919,7 +1223,14 @@ main_run_de_analysis <- function(args) {
|
||||
sprintf(" Genes with DTU (q < 0.05): %d", contrast_qc$dtu_significant_genes)
|
||||
)
|
||||
} else {
|
||||
c(
|
||||
if (!is.null(contrast_qc$dtu_failure_hint)) {
|
||||
paste0(" Probable cause: ", contrast_qc$dtu_failure_hint)
|
||||
} else {
|
||||
NULL
|
||||
},
|
||||
" See DTU_ANALYSIS_FAILED.txt for details"
|
||||
)
|
||||
},
|
||||
if (!is.null(contrast_qc$dtu_power_warning)) paste0(" WARNING: ", contrast_qc$dtu_power_warning) else NULL,
|
||||
""
|
||||
@ -928,6 +1239,15 @@ main_run_de_analysis <- function(args) {
|
||||
de_qc_stats$contrasts[[contrast_name]] <- contrast_qc
|
||||
}
|
||||
|
||||
failed_dge_contrasts <- names(Filter(
|
||||
function(cqc) identical(cqc$dge_status, "FAILED"),
|
||||
de_qc_stats$contrasts
|
||||
))
|
||||
failed_dtu_contrasts <- names(Filter(
|
||||
function(cqc) identical(cqc$dtu_status, "FAILED"),
|
||||
de_qc_stats$contrasts
|
||||
))
|
||||
|
||||
deseq2_dispersion_fallbacks <- names(Filter(
|
||||
function(cqc) isTRUE(cqc$deseq2_dispersion_fallback$applied),
|
||||
de_qc_stats$contrasts
|
||||
@ -965,6 +1285,10 @@ main_run_de_analysis <- function(args) {
|
||||
integer(1)
|
||||
))
|
||||
de_qc_stats$analysis_fallbacks <- list(
|
||||
failed_dge_contrasts = length(failed_dge_contrasts),
|
||||
failed_dge_contrast_names = as.list(failed_dge_contrasts),
|
||||
failed_dtu_contrasts = length(failed_dtu_contrasts),
|
||||
failed_dtu_contrast_names = as.list(failed_dtu_contrasts),
|
||||
deseq2_dispersion_fallback_contrasts = length(deseq2_dispersion_fallbacks),
|
||||
deseq2_dispersion_fallback_contrast_names = as.list(deseq2_dispersion_fallbacks),
|
||||
deseq2_gene_wise_contrasts = length(deseq2_gene_wise),
|
||||
@ -1015,8 +1339,14 @@ main_run_de_analysis <- function(args) {
|
||||
"",
|
||||
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(" DGE status: %s", cqc$dge_status),
|
||||
if (cqc$dge_status == "SUCCESS") {
|
||||
sprintf(" DGE significant: %d genes (FDR<0.05)", cqc$dge_significant_fdr05)
|
||||
} else {
|
||||
" DGE significant: N/A (analysis failed)"
|
||||
},
|
||||
sprintf(" DESeq2 size factors: %s", cqc$deseq2_size_factor_method),
|
||||
if (cqc$dge_status != "SUCCESS") " See DGE_ANALYSIS_FAILED.txt for details" else NULL,
|
||||
sprintf(" DTU status: %s", cqc$dtu_status),
|
||||
sprintf(" DEXSeq size factors: %s", cqc$dexseq_size_factor_method),
|
||||
if (cqc$dtu_status == "SUCCESS") sprintf(" DTU significant: %d genes", cqc$dtu_significant_genes) else NULL
|
||||
|
||||
@ -308,3 +308,86 @@ expect_bam_fixture_built <- function(reference, reads, out_dir, alias = "sampleA
|
||||
|
||||
bam_path
|
||||
}
|
||||
|
||||
#' Create a degenerate DE fixture where both condition groups have identical
|
||||
#' zero-heavy counts, simulating the control-vs-control case from the bug reports.
|
||||
#'
|
||||
#' Properties:
|
||||
#' - Every transcript has at least one zero (forces DEXSeq poscounts path)
|
||||
#' - Every gene has at least one zero (forces DESeq2 poscounts path via
|
||||
#' de_choose_size_factor_type)
|
||||
#' - Both condition groups have identical counts (no differential signal,
|
||||
#' stresses dispersion estimation)
|
||||
#'
|
||||
#' @param out_dir Directory to write the output files.
|
||||
#' @return A list with paths to transcript RDS, gene RDS, and sample sheet CSV.
|
||||
#' @export
|
||||
write_degenerate_de_fixture_bundle <- function(out_dir) {
|
||||
sample_names <- c(
|
||||
"ctrl_rep1", "ctrl_rep2", "ctrl_rep3",
|
||||
"case_rep1", "case_rep2", "case_rep3"
|
||||
)
|
||||
|
||||
# Identical counts across both groups; zeros in every transcript.
|
||||
tx_counts <- matrix(
|
||||
c(
|
||||
0, 0, 4, 0, 0, 4, # tx1: zero in rep1 and rep2 of both groups
|
||||
3, 4, 3, 3, 4, 3, # tx2: all nonzero but no group difference
|
||||
2, 2, 0, 2, 2, 0, # tx3: zero in rep3
|
||||
0, 1, 1, 0, 1, 1 # tx4: zero in rep1
|
||||
),
|
||||
nrow = 4, byrow = TRUE,
|
||||
dimnames = list(c("tx1", "tx2", "tx3", "tx4"), sample_names)
|
||||
)
|
||||
tx_cpm <- t(t(tx_counts) / pmax(colSums(tx_counts), 1)) * 1e6
|
||||
|
||||
row_ranges <- GenomicRanges::GRanges(
|
||||
seqnames = rep("chr1", 4),
|
||||
ranges = IRanges::IRanges(start = c(1, 101, 201, 301), width = 50),
|
||||
strand = rep("+", 4),
|
||||
TXNAME = c("tx1", "tx2", "tx3", "tx4"),
|
||||
GENEID = c("gene1", "gene1", "gene2", "gene2"),
|
||||
eqClassById = IRanges::CharacterList(list(c("1", "2"), "3", "4", "5"))
|
||||
)
|
||||
|
||||
tx_se <- SummarizedExperiment::SummarizedExperiment(
|
||||
assays = list(counts = tx_counts, CPM = tx_cpm),
|
||||
rowRanges = row_ranges
|
||||
)
|
||||
|
||||
# Gene counts: gene1 = tx1+tx2, gene2 = tx3+tx4.
|
||||
# Each gene has at least one zero per group (triggers poscounts for DGE).
|
||||
gene_counts <- rbind(
|
||||
gene1 = tx_counts["tx1", ] + tx_counts["tx2", ],
|
||||
gene2 = tx_counts["tx3", ] + tx_counts["tx4", ]
|
||||
)
|
||||
gene_cpm <- t(t(gene_counts) / pmax(colSums(gene_counts), 1)) * 1e6
|
||||
|
||||
gene_se <- SummarizedExperiment::SummarizedExperiment(
|
||||
assays = list(counts = gene_counts, CPM = gene_cpm),
|
||||
rowData = S4Vectors::DataFrame(GENEID = rownames(gene_counts))
|
||||
)
|
||||
|
||||
tx_path <- file.path(out_dir, "transcripts.rds")
|
||||
gene_path <- file.path(out_dir, "genes.rds")
|
||||
saveRDS(tx_se, tx_path)
|
||||
saveRDS(gene_se, gene_path)
|
||||
|
||||
sample_sheet <- file.path(out_dir, "sample_sheet.csv")
|
||||
utils::write.csv(
|
||||
data.frame(
|
||||
alias = sample_names,
|
||||
condition = rep(c("control", "case"), each = 3),
|
||||
stringsAsFactors = FALSE
|
||||
),
|
||||
sample_sheet,
|
||||
row.names = FALSE,
|
||||
quote = FALSE
|
||||
)
|
||||
|
||||
list(
|
||||
transcript_rds = tx_path,
|
||||
gene_rds = gene_path,
|
||||
sample_sheet = sample_sheet
|
||||
)
|
||||
}
|
||||
|
||||
@ -535,9 +535,88 @@ testthat::test_that("de_is_recoverable_dexseq_error recognises expected messages
|
||||
))
|
||||
testthat::expect_true(de_is_recoverable_dexseq_error("model matrix is not full rank"))
|
||||
testthat::expect_true(de_is_recoverable_dexseq_error("replacement has 1 row, data has 0"))
|
||||
# Internal object-shape mismatches should still be reported as failures,
|
||||
# but they are not known statistical fallback patterns.
|
||||
testthat::expect_false(de_is_recoverable_dexseq_error(
|
||||
"nb of cols in 'assay' (7) must equal nb of rows in 'colData' (12)"
|
||||
))
|
||||
testthat::expect_true(de_is_recoverable_dexseq_error(
|
||||
"newsplit: out of vertex space"
|
||||
))
|
||||
testthat::expect_true(de_is_recoverable_dexseq_error(
|
||||
"every gene contains at least one zero, cannot compute log geometric means"
|
||||
))
|
||||
testthat::expect_false(de_is_recoverable_dexseq_error("random unrelated failure"))
|
||||
})
|
||||
|
||||
testthat::test_that("de_dexseq_failure_hint identifies degenerate-data shape mismatch", {
|
||||
hint <- de_dexseq_failure_hint(
|
||||
"nb of cols in 'assay' (7) must equal nb of rows in 'colData' (12)"
|
||||
)
|
||||
testthat::expect_match(hint, "degenerate or near-identical sample groups")
|
||||
row_sums_hint <- de_dexseq_failure_hint(
|
||||
"Error in rowSums(numerator) : 'x' must be an array of at least two dimensions"
|
||||
)
|
||||
testthat::expect_match(row_sums_hint, "gene-level DTU summarisation")
|
||||
testthat::expect_null(de_dexseq_failure_hint("random unrelated failure"))
|
||||
})
|
||||
|
||||
testthat::test_that("de_postprocess_dexseq_result degrades gracefully on gene-level DTU failure", {
|
||||
contrast_dir <- tempfile("dexseq-postprocess-")
|
||||
dir.create(contrast_dir)
|
||||
|
||||
dex_res <- list(
|
||||
dxr = data.frame(
|
||||
featureID = c("tx1", "tx2"),
|
||||
groupID = c("gene1", "gene1"),
|
||||
log2fold_case_control = c(1, -1),
|
||||
pvalue = c(0.01, 0.2),
|
||||
padj = c(0.02, 0.3),
|
||||
exonBaseMean = c(10, 20),
|
||||
check.names = FALSE
|
||||
),
|
||||
dxd = structure(list(), class = "mock_dxd")
|
||||
)
|
||||
contrast_samples <- data.frame(
|
||||
alias = c("s1", "s2"),
|
||||
condition = c("case", "control"),
|
||||
stringsAsFactors = FALSE
|
||||
)
|
||||
tx_counts <- matrix(
|
||||
c(10, 0, 4, 8),
|
||||
nrow = 2,
|
||||
dimnames = list(c("tx1", "tx2"), c("s1", "s2"))
|
||||
)
|
||||
|
||||
result <- testthat::expect_warning(
|
||||
de_postprocess_dexseq_result(
|
||||
dex_res = dex_res,
|
||||
contrast_dir = contrast_dir,
|
||||
target_level = "case",
|
||||
reference_level = "control",
|
||||
contrast_samples = contrast_samples,
|
||||
condition_column = "condition",
|
||||
tx_counts = tx_counts,
|
||||
contrast_name = "case_control",
|
||||
per_gene_qvalue_fn = function(...) {
|
||||
stop("'x' must be an array of at least two dimensions", call. = FALSE)
|
||||
},
|
||||
plot_writer = function(...) NULL
|
||||
),
|
||||
"DEXSeq post-processing failed"
|
||||
)
|
||||
|
||||
testthat::expect_false(result$ok)
|
||||
testthat::expect_match(
|
||||
result$failure_hint,
|
||||
"gene-level DTU summarisation"
|
||||
)
|
||||
testthat::expect_true(file.exists(file.path(contrast_dir, "DTU_ANALYSIS_FAILED.txt")))
|
||||
failure_text <- readLines(file.path(contrast_dir, "DTU_ANALYSIS_FAILED.txt"))
|
||||
testthat::expect_true(any(grepl("DTU Analysis Failed", failure_text, fixed = TRUE)))
|
||||
testthat::expect_true(any(grepl("at least two dimensions", failure_text, fixed = TRUE)))
|
||||
})
|
||||
|
||||
testthat::test_that("de_run_dexseq_result records rank-deficiency covariate drops", {
|
||||
testthat::skip_if_not_installed("DESeq2")
|
||||
testthat::skip_if_not_installed("DEXSeq")
|
||||
@ -779,7 +858,89 @@ testthat::test_that("CLI integration produces expected outputs", {
|
||||
contrast_qc <- de_qc$contrasts[["condition_treated_vs_control"]]
|
||||
testthat::expect_true("deseq2_size_factor_method" %in% names(contrast_qc))
|
||||
testthat::expect_true("deseq2_dispersion_fallback" %in% names(contrast_qc))
|
||||
testthat::expect_true("dge_status" %in% names(contrast_qc))
|
||||
testthat::expect_true("dexseq_size_factor_method" %in% names(contrast_qc))
|
||||
testthat::expect_true("dexseq_dispersion_method" %in% names(contrast_qc))
|
||||
testthat::expect_true("dexseq_covariates_dropped" %in% names(contrast_qc))
|
||||
})
|
||||
|
||||
###
|
||||
# Degenerate input (control-vs-control / zero-heavy counts)
|
||||
#
|
||||
# Exercises the graceful-degradation paths added for three cases:
|
||||
# - newsplit/lfproc OOM in DESeq2 local regression
|
||||
# - DEXSeq colData dimension mismatch on degenerate input
|
||||
# - size-factor estimation failure when every gene has a zero
|
||||
#
|
||||
# The process must run to completion and produce placeholder outputs + a
|
||||
# machine-readable status in de_qc_stats.json rather than crashing.
|
||||
|
||||
testthat::test_that("degenerate control-vs-control data completes gracefully", {
|
||||
testthat::skip_if_not_installed("DESeq2")
|
||||
testthat::skip_if_not_installed("DEXSeq")
|
||||
|
||||
fixture_dir <- tempfile("de-degenerate-")
|
||||
dir.create(fixture_dir)
|
||||
bundle <- write_degenerate_de_fixture_bundle(fixture_dir)
|
||||
|
||||
argv <- c(
|
||||
bundle,
|
||||
list(
|
||||
condition_column = "condition",
|
||||
covariates = NULL,
|
||||
reference_level = "control",
|
||||
out_dir = file.path(fixture_dir, "out")
|
||||
)
|
||||
)
|
||||
argv <- workflow_glue_r_normalise_args(argv, de_analysis_arg_spec())
|
||||
|
||||
# Must complete without throwing — any DESeq2/DEXSeq failure should be
|
||||
# caught and degraded to placeholder outputs rather than crashing the process.
|
||||
result <- tryCatch(
|
||||
suppressWarnings(suppressMessages(main_run_de_analysis(argv))),
|
||||
error = function(err) err
|
||||
)
|
||||
testthat::expect_false(
|
||||
inherits(result, "error"),
|
||||
info = if (inherits(result, "error")) conditionMessage(result) else ""
|
||||
)
|
||||
|
||||
# QC JSON must exist with per-contrast status fields.
|
||||
qc_path <- file.path(argv$out_dir, "de_qc_stats.json")
|
||||
testthat::expect_true(file.exists(qc_path))
|
||||
de_qc <- jsonlite::read_json(qc_path, simplifyVector = TRUE)
|
||||
|
||||
contrast_name <- "condition_case_vs_control"
|
||||
contrast_qc <- de_qc$contrasts[[contrast_name]]
|
||||
testthat::expect_false(is.null(contrast_qc))
|
||||
testthat::expect_true(contrast_qc$dge_status %in% c("SUCCESS", "FAILED"))
|
||||
testthat::expect_true(contrast_qc$dtu_status %in% c("SUCCESS", "FAILED"))
|
||||
|
||||
# DGE results file must always exist (empty placeholder on failure).
|
||||
contrast_dir <- file.path(argv$out_dir, contrast_name)
|
||||
dge_tsv <- file.path(contrast_dir, "results_dge.tsv")
|
||||
testthat::expect_true(file.exists(dge_tsv))
|
||||
|
||||
if (identical(contrast_qc$dge_status, "FAILED")) {
|
||||
# Failure log and placeholder PDF written.
|
||||
testthat::expect_true(file.exists(file.path(contrast_dir, "DGE_ANALYSIS_FAILED.txt")))
|
||||
testthat::expect_true(file.exists(file.path(contrast_dir, "results_dge.pdf")))
|
||||
# Empty TSV must have the expected column headers.
|
||||
dge_df <- utils::read.delim(dge_tsv, check.names = FALSE)
|
||||
testthat::expect_equal(nrow(dge_df), 0)
|
||||
testthat::expect_true(all(c("GENEID", "log2FoldChange", "padj") %in% names(dge_df)))
|
||||
}
|
||||
|
||||
# DTU outputs must always exist.
|
||||
testthat::expect_true(file.exists(file.path(contrast_dir, "results_dtu_transcript.tsv")))
|
||||
testthat::expect_true(file.exists(file.path(contrast_dir, "results_dtu_gene.tsv")))
|
||||
|
||||
if (identical(contrast_qc$dtu_status, "FAILED")) {
|
||||
testthat::expect_true(file.exists(file.path(contrast_dir, "DTU_ANALYSIS_FAILED.txt")))
|
||||
}
|
||||
|
||||
# analysis_fallbacks must record failed counts.
|
||||
fallbacks <- de_qc$analysis_fallbacks
|
||||
testthat::expect_true(!is.null(fallbacks$failed_dge_contrasts))
|
||||
testthat::expect_true(!is.null(fallbacks$failed_dtu_contrasts))
|
||||
})
|
||||
|
||||
Loading…
Reference in New Issue
Block a user