Fix non-conformable arrays bug

CHANGELOG.md

@@ -9,6 +9,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 
 - [#229](https://github.com/qbic-pipelines/rnadeseq/pull/229) Added param for clustering (or not) the heatmaps
 - [#225](https://github.com/qbic-pipelines/rnadeseq/pull/225) Added param for pathway analysis datasources
+- [#221](https://github.com/qbic-pipelines/rnadeseq/pull/221) Added padj to volcano hovertext
 
 ### Changed
 
@@ -18,6 +19,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 
 - [#229](https://github.com/qbic-pipelines/rnadeseq/pull/229) Fixed cutoff enrichment plot labels, fixed wrong plotMA function being called (also fixed this changelog)
 - [#225](https://github.com/qbic-pipelines/rnadeseq/pull/225) Fixed too many devices error from tryCatch around normalized heatmaps
+- [#221](https://github.com/qbic-pipelines/rnadeseq/pull/221) Fixed non-conformable arrays bug, fix wrong volcano colors when no DE genes
 
 ## 2.2 Avenue of Poplars
 

assets/RNAseq_report.Rmd

@@ -698,6 +698,7 @@ if (run_rlog){
     # rlog transformation
     rld <- rlog(cds, blind=FALSE)
     rld_names <- merge(x=gene_names, y=assay(rld), by.x = "Ensembl_ID", by.y="row.names")
+    rld_names <- rld_names[order(rld_names$Ensembl_ID),]
     write.table(rld_names, "differential_gene_expression/gene_counts_tables/rlog_transformed_gene_counts.tsv", append = FALSE, quote = FALSE, sep = "\t",eol = "\n", na = "NA", dec = ".", row.names = F,  qmethod = c("escape", "double"))
 
     # save table to another variable if pathway analysis
@@ -711,6 +712,7 @@ if (run_rlog){
     # vst transformation
     vsd <- vst(cds, blind=FALSE, nsub = params$nsub_genes)
     vsd_names <- merge(x=gene_names, y=assay(vsd), by.x = "Ensembl_ID", by.y="row.names")
+    vsd_names <- vsd_names[order(vsd_names$Ensembl_ID),]
     write.table(vsd_names, "differential_gene_expression/gene_counts_tables/vst_transformed_gene_counts.tsv", append = FALSE, quote = FALSE, sep = "\t",eol = "\n", na = "NA", dec = ".", row.names = F, qmethod = c("escape", "double"))
 
     # save table to another variable if pathway analysis
@@ -1254,9 +1256,20 @@ A PCA of the batch-effect corrected data is shown below and can be found [here](
 if (run_rlog) {
     assay(rld) <- limma::removeBatchEffect(assay(rld), rld$batch)
     rld_corrected <- rld
+    out_rld <- as.data.frame(assay(rld_corrected))
+    out_rld <- out_rld[order(rownames(out_rld)),]
+    out_rld <- cbind(Ensembl_ID = rownames(out_rld), gene_name = rld_names$gene_name, out_rld[,!names(out_rld) %in% c("gene_name")])
+
+    write.table(out_rld, file="differential_gene_expression/gene_counts_tables/rlog_transformed_gene_counts_batchcorrected.tsv", quote=F, sep="\t", row.names=F)
 } else {
     assay(vsd) <- limma::removeBatchEffect(assay(vsd), vsd$batch)
     vsd_corrected <- vsd
+
+    out_vsd <- as.data.frame(assay(vsd_corrected))
+    out_vsd <- out_vsd[order(rownames(out_vsd)),]
+    out_vsd <- cbind(Ensembl_ID = rownames(out_vsd), gene_name = vsd_names$gene_name, out_vsd[,!names(out_vsd) %in% c("gene_name")])
+
+    write.table(out_vsd, file="differential_gene_expression/gene_counts_tables/vst_transformed_gene_counts_batchcorrected.tsv", quote=F, sep="\t", row.names=F)
 }
 pcaData2 <- plotPCA(if (run_rlog) rld else vsd, intgroup=c("combfactor"), ntop = dim(if (run_rlog) rld else vsd)[1], returnData=TRUE)
 percentVar <- round(100*attr(pcaData, "percentVar"))
@@ -1579,6 +1592,11 @@ write(contrast_names, file="contrast_names.txt", sep="\t")
 # Remove identical columns of DE_genes_df
 DE_genes_df$DE_genes_df <- NULL
 idx <- duplicated(t(DE_genes_df))
+
+# If any padj or logFC columns were marked for removal (TRUE), undo that by setting FALSE as we need these cols
+idx[grepl("padj", as.character(names(idx)))] <- FALSE
+idx[grepl("log2FoldChange", as.character(names(idx)))] <- FALSE
+
 DE_genes_df <- DE_genes_df[, !idx]
 DE_genes_df$Ensembl_ID <- row.names(DE_genes_df)
 DE_genes_df <- DE_genes_df[,c(dim(DE_genes_df)[2],1:dim(DE_genes_df)[2]-1)]
@@ -1694,17 +1712,21 @@ for (file in allgenes_files){
 }
 DE_all <- ldply(table_list, rbind)
 DE_all$logpval <- -log10(DE_all$padj)
+DE_all$truelogpval <- DE_all$logpval # Save true value for hovertext
 DE_all$logpval[DE_all$logpval > 16] <- 17
 DE_all <- na.omit(DE_all)
+DE_all <- DE_all %>% mutate(color = ifelse(abs(DE_all$log2FoldChange) >= params$logFC_threshold & DE_all$logpval >= -log10(params$adj_pval_threshold), "red", "black")) # Create column which assigns each entry a plot color so that coloring is consistent even when no DE genes are found
 log2FoldChange_min <- min(DE_all$log2FoldChange)
 log2FoldChange_max <- max(DE_all$log2FoldChange)
 
 # get quotient of log2FoldChange range and 10 to determine if xticks of plot need to be scaled discretely (otherwise, for large FC ranges, there will be too many xticks and they will overlap)
 log2FoldChange_quotient <- length(seq(log2FoldChange_min, log2FoldChange_max, 10))
 
-pg <- ggplot(DE_all, aes(x=log2FoldChange, y=logpval, text=paste("Gene: ", gene_name, "<br>", "Log2FC: ", formatC(log2FoldChange, digits=2)))) +
+pg <- ggplot(DE_all, aes(x=log2FoldChange, y=logpval, text=paste("Gene: ", gene_name,
+                                                                    "<br>", "Log2FC: ", formatC(log2FoldChange, digits=2),
+                                                                    "<br>", "-Log10padj: ", formatC(truelogpval, digits=2)))) +
         #this goes after alpha=0.5:
-        geom_jitter(alpha=0.5, width = 0.2, aes(color=ifelse(abs(DE_all$log2FoldChange) >= params$logFC_threshold &  DE_all$logpval >= -log10(params$adj_pval_threshold), "Differentially expressed genes", "Non-differentially expressed genes"))) +
+        geom_jitter(alpha=0.5, width = 0.2, aes(color=ifelse(DE_all$color == "red", "Differentially expressed genes", "Non-differentially expressed genes"))) +
         geom_hline(yintercept = 16, linetype= "dashed", size = 0.2, color = "grey") +
         geom_hline(yintercept = -log10(params$adj_pval_threshold), size = 0.2, color = "grey") +
         geom_vline(xintercept = -params$logFC_threshold, size = 0.2, color = "grey") +

lib/NfcoreTemplate.groovy

@@ -3,6 +3,7 @@
 //
 
 import org.yaml.snakeyaml.Yaml
+import groovy.json.JsonOutput
 
 class NfcoreTemplate {
 
@@ -222,6 +223,21 @@ class NfcoreTemplate {
         }
     }
 
+    //
+    // Dump pipeline parameters in a json file
+    //
+    public static void dump_parameters(workflow, params) {
+        def output_d = new File("${params.outdir}/pipeline_info/")
+        if (!output_d.exists()) {
+            output_d.mkdirs()
+        }
+
+        def timestamp  = new java.util.Date().format( 'yyyy-MM-dd_HH-mm-ss')
+        def output_pf  = new File(output_d, "params_${timestamp}.json")
+        def jsonStr    = JsonOutput.toJson(params)
+        output_pf.text = JsonOutput.prettyPrint(jsonStr)
+    }
+
     //
     // Print pipeline summary on completion
     //

tests/test_full.yml

@@ -82,4 +82,8 @@
       md5sum: 6594fdb8c0eb591fa131982a8be00877
     - path: results_test/pathway_analysis/DE_contrast_condition_treatment_Treated_vs_Control/enrichment_plots/REAC_pathway_enrichment_plot.pdf
     - path: results_test/pathway_analysis/DE_contrast_condition_treatment_Treated_vs_Control/enrichment_plots/REAC_pathway_enrichment_plot.png
+    - path: results_test/pathway_analysis/heatmap_gene_list/Heatmap_normalized_counts_gene_list.pdf
+    - path: results_test/pathway_analysis/DE_contrast_condition_treatment_Treated_vs_Control/gost_pathway_venn_diagram.pdf
+    - path: results_test/pathway_analysis/DE_contrast_condition_treatment_Treated_vs_Control/gost_pathway_venn_diagram.png
+    - path: results_test/pathway_analysis/DE_contrast_condition_treatment_Treated_vs_Control/gost_pathway_venn_diagram.svg
     - path: results_test/RNAseq_report.html