Merge branch 'alias_hyphens_CW-5484' into 'dev'

Fix hyphens in sample alias CW-5484

See merge request epi2melabs/workflows/wf-transcriptomes!212

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

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:

subworkflows/differential_expression.nf

@@ -89,6 +89,21 @@ process deAnalysis {
         --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
     """
 }