Merge branch 'sample_col_check_CW-6532' into 'dev'

Fix and simplify column name checking

Closes CW-6532

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

bin/de_analysis.R

@@ -21,8 +21,8 @@ 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
+rename_sample_columns <- function(df, alias, start_col) {
+  colnames(df)[start_col:ncol(df)] <- alias
   return(df)
 }
 
@@ -37,13 +37,14 @@ 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, 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, check.names=FALSE)
+coldata$alias <- rownames(coldata)
+#dmDSdata looks up sample_id in coldata, we need this to be sample alias
+coldata$sample_id <- coldata$alias
 
-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]
+cts <- cts[, coldata$alias, drop = FALSE]
 # check if control condition exists, sets as reference 
 if(!"control" %in% coldata$condition)
   stop("sample_sheet.csv does not contain 'control' 
@@ -101,10 +102,11 @@ 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)
+counts <- rename_sample_columns(counts, coldata$alias, 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)
 
+
 cat("Filtering counts using DRIMSeq.\n")
 
 d <- dmDSdata(counts=counts, samples=coldata)
@@ -122,13 +124,13 @@ 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)
+trs_cts <- rename_sample_columns(trs_cts, coldata$alias, 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)
+gene_cts <- rename_sample_columns(gene_cts, coldata$alias, 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
@@ -137,7 +139,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)
+cpm_gene_counts <- rename_sample_columns(cpm_gene_counts, coldata$alias, 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

@@ -78,6 +78,7 @@ process deAnalysis {
         path "de_analysis/results_dtu_stageR.tsv", emit: dtu_stageR
         path "de_analysis/results_dtu.pdf", emit: dtu_pdf
         path "de_analysis/cpm_gene_counts.tsv", emit: cpm
+    script:
     """
     de_analysis.R \
         --annotation annotation.gtf \
@@ -91,11 +92,13 @@ process deAnalysis {
         --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
+    # Sample column order should be in the same order as the input sample sheet
+    alias_col=\$(awk -v RS=',' '/alias/{print NR; exit}' "sample_sheet.csv")
+    cut -d',' -f3 sample_sheet.csv | tail -n +2 | paste -sd '\t' - | sed 's/\t*\$//' > 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
+    head -1 de_analysis/cpm_gene_counts.tsv | cut -f2-  > cpm_gene_counts_colnames
+    head -1 merged/all_gene_counts.tsv  > merged_counts_colnames
+    head -1 merged/filtered_transcript_counts_with_genes.tsv | cut -f3-  > merged_filtered_colnames
     
     # Check for mismatches in sample column names
     for file in cpm_gene_counts_colnames merged_counts_colnames merged_filtered_colnames; do