From 19b54ba2ffb6b72ad125e7f900c8973bc9c903f7 Mon Sep 17 00:00:00 2001 From: Neil Horner Date: Wed, 14 May 2025 13:55:01 +0000 Subject: [PATCH] Fix hyphens in sample alias CW-5484 --- CHANGELOG.md | 1 + bin/de_analysis.R | 22 +++++++++++++---- subworkflows/differential_expression.nf | 33 ++++++++++++++++++------- 3 files changed, 42 insertions(+), 14 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 67c73b7..9a0f992 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -17,6 +17,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 - dacite.exceptions.WrongTypeError during report generation when barcode is null. - Sequence summary read length N50 incorrectly displayed minimum read length, it now correctly shows the N50. - Sequence summary component alignment and coverage plots failed to plot under some conditions. +- Error in `deAnalysis` process - `mode(counts) %in% "numeric" is not TRUE` - caused by hyphens in sample sheet aliases. ## [v1.7.0] ### Changed diff --git a/bin/de_analysis.R b/bin/de_analysis.R index 83f13b3..8888b87 100755 --- a/bin/de_analysis.R +++ b/bin/de_analysis.R @@ -18,6 +18,14 @@ suppressMessages(library("DRIMSeq")) suppressMessages(library("GenomicFeatures")) suppressMessages(library("edgeR")) +# Some functions, including dmDSdata, converts '.' in sample IDs to '-'. +# Make the output DF match the sample IDs in the sample sheet +# start_col is the index of the first sample column in the data frame +rename_sample_columns <- function(df, sample_ids, start_col) { + colnames(df)[start_col:ncol(df)] <- sample_ids + return(df) +} + # Create output directories if (!dir.exists(argv$de_out_dir)){ dir.create(argv$de_out_dir, recursive=TRUE) @@ -27,13 +35,15 @@ if (!dir.exists(argv$merged_out_dir)){ } cat("Loading counts, conditions and parameters.\n") -cts <- as.matrix(read.csv(argv$all_counts, sep="\t", row.names="Reference", stringsAsFactors=FALSE)) - +cts <- as.matrix(read.csv(argv$all_counts, sep="\t", row.names="Reference", stringsAsFactors=FALSE, check.names=FALSE)) # Set up sample data frame: #changed this to sample_id -coldata <- read.csv(argv$sample_sheet, row.names="alias", sep=",", stringsAsFactors=TRUE) +coldata <- read.csv(argv$sample_sheet, row.names="alias", sep=",", stringsAsFactors=TRUE, check.names=FALSE) coldata$sample_id <- rownames(coldata) +# Reorder the input counts columns to match the sample sheet order. +# This ensures we don't missassign column names, when renaming output DF columns. +cts <- cts[, coldata$sample_id, drop = FALSE] # check if control condition exists, sets as reference if(!"control" %in% coldata$condition) stop("sample_sheet.csv does not contain 'control' @@ -84,7 +94,6 @@ txdf$ntx<- tab[match(txdf$GENEID, names(tab))] cts <- cts[rownames(cts) %in% txdf$TXNAME, ] # FIXME: filter for transcripts which are in the annotation. Why they are not all there? - # Reorder transcript/gene database to match input counts: txdf <- txdf[match(rownames(cts), txdf$TXNAME), ] rownames(txdf) <- NULL @@ -92,6 +101,7 @@ rownames(txdf) <- NULL # Create counts data frame: counts<-data.frame(gene_id=txdf$GENEID, feature_id=txdf$TXNAME, cts) +counts <- rename_sample_columns(counts, coldata$sample_id, 3) # output unfiltered version of the counts table now we have paired transcripts with gene ids write.table(counts, file=file.path(argv$de_out_dir, "unfiltered_transcript_counts_with_genes.tsv"), sep="\t", row.names = FALSE, quote=FALSE) @@ -107,17 +117,18 @@ cat("Building model matrix.\n") design <- model.matrix(~condition, data=DRIMSeq::samples(d)) - suppressMessages(library("dplyr")) # Sum transcript counts into gene counts: cat("Sum transcript counts into gene counts.\n") trs_cts <- counts(d) +trs_cts <- rename_sample_columns(trs_cts, coldata$sample_id, 3) write.table(trs_cts, file=file.path(argv$merged_out_dir, "filtered_transcript_counts_with_genes.tsv"), sep="\t", row.names = FALSE, quote=FALSE) gene_cts <- trs_cts_unfiltered %>% dplyr::select(c(1, 3:ncol(trs_cts))) %>% group_by(gene_id) %>% summarise_all(tibble::lst(sum)) %>% data.frame() rownames(gene_cts) <- gene_cts$gene_id gene_cts$gene_id <- NULL +gene_cts <- rename_sample_columns(gene_cts, coldata$sample_id, 1) write.table(gene_cts, file=file.path(argv$merged_out_dir, "all_gene_counts.tsv"), sep="\t", quote=FALSE) # Output count per million of the gene counts using edgeR CPM @@ -126,6 +137,7 @@ cpm_gene_counts <- cpm(gene_cts) cpm_gene_counts <- cbind(var_name = rownames(cpm_gene_counts), cpm_gene_counts) rownames(cpm_gene_counts) <- NULL colnames(cpm_gene_counts)[1] <- "gene_id" +cpm_gene_counts <- rename_sample_columns(cpm_gene_counts, coldata$sample_id, 2) write.table(cpm_gene_counts, file=file.path(argv$de_out_dir, "cpm_gene_counts.tsv"), sep="\t", quote=FALSE, row.names = FALSE) # Differential gene expression using edgeR: diff --git a/subworkflows/differential_expression.nf b/subworkflows/differential_expression.nf index 5e1c97b..981b1e7 100644 --- a/subworkflows/differential_expression.nf +++ b/subworkflows/differential_expression.nf @@ -80,15 +80,30 @@ process deAnalysis { path "de_analysis/cpm_gene_counts.tsv", emit: cpm """ de_analysis.R \ - --annotation annotation.gtf \ - --min_samps_gene_expr $params.min_samps_gene_expr \ - --min_samps_feature_expr $params.min_samps_feature_expr \ - --min_gene_expr $params.min_gene_expr \ - --min_feature_expr $params.min_feature_expr \ - --sample_sheet sample_sheet.csv \ - --all_counts all_counts.tsv \ - --de_out_dir de_analysis \ - --merged_out_dir merged + --annotation annotation.gtf \ + --min_samps_gene_expr $params.min_samps_gene_expr \ + --min_samps_feature_expr $params.min_samps_feature_expr \ + --min_gene_expr $params.min_gene_expr \ + --min_feature_expr $params.min_feature_expr \ + --sample_sheet sample_sheet.csv \ + --all_counts all_counts.tsv \ + --de_out_dir de_analysis \ + --merged_out_dir merged + + # Check that the original aliases in the input TSV have not been mangled by R's read.csv or other functions + head -1 all_counts.tsv | cut -f2- | tr '\t' '\n' > expected_colnames + + head -1 de_analysis/cpm_gene_counts.tsv | cut -f2- | tr '\t' '\n' > cpm_gene_counts_colnames + head -1 merged/all_gene_counts.tsv | tr '\t' '\n' > merged_counts_colnames + head -1 merged/filtered_transcript_counts_with_genes.tsv | cut -f3- | tr '\t' '\n' > merged_filtered_colnames + + # Check for mismatches in sample column names + for file in cpm_gene_counts_colnames merged_counts_colnames merged_filtered_colnames; do + if ! diff -q \$file expected_colnames > /dev/null; then + echo "Column names in \$file do not match expected aliases." + exit 70 + fi + done """ }