report: added mash figure (ssiue #59)

workflow_report/Snakefile

@@ -47,6 +47,7 @@ rule render_report:
         diamond=os.path.join(RESULTS_DIR, "annotation/diamond/summary_qs_len_ratio.tsv"),
         prodigal=os.path.join(RESULTS_DIR, "annotation/prodigal/summary_gene_counts.txt"),
         cdhit=os.path.join(RESULTS_DIR, "analysis/cdhit/summary.tsv"),
+        mash=os.path.join(RESULTS_DIR, "analysis/mash/contigs.dist"),
     output:
         html=os.path.join(RESULTS_DIR, "report/report.html"),
         pdf_mappability=os.path.join(RESULTS_DIR, "report/fig_mappability.pdf"),
@@ -57,6 +58,7 @@ rule render_report:
         pdf_diamond=os.path.join(RESULTS_DIR, "report/fig_diamond.pdf"),
         pdf_prodigal=os.path.join(RESULTS_DIR, "report/fig_prodigal.pdf"),
         pdf_cdhit=os.path.join(RESULTS_DIR, "report/fig_cdhit.pdf"),
+        pdf_mash=os.path.join(RESULTS_DIR, "report/fig_mash.pdf"),
     log:
         out="logs/report.out.log",
         err="logs/report.err.log"

workflow_report/scripts/report.R

@@ -60,6 +60,21 @@ if(!is.null(TABS$mappability$metat)){
         mappability_theme
 }
 
+##############################
+# Mash
+
+TABS$mash <- read_mash_dist(snakemake@input$mash)
+mash_htree <- hclust(d=as.dist(TABS$mash), method="complete")
+FIGS$mash <- pheatmap(
+    TABS$mash,
+    color=viridis(20),
+    border_color="white",
+    scale="none",
+    cluster_rows=mash_htree,
+    cluster_cols=mash_htree,
+    silent=TRUE
+)
+
 ##############################
 # QAUST: assembly quality
 

workflow_report/scripts/report.Rmd

@@ -36,6 +36,13 @@ if("metat" %in% names(FIGS$mappability)){
 }
 ```
 
+# Assembly similarity
+
+```{r figures-mash, echo=FALSE, fig.width=10, fig.height=10, fig.cap='Assembly similarity (Mash)'}
+grid.newpage()
+grid.draw(FIGS$mash)
+```
+
 # Assembly quality
 
 ```{r tables-quast, echo=FALSE}

workflow_report/scripts/report_render.R

@@ -10,17 +10,21 @@ sink(file=file(snakemake@log$err, open="wt"), type="message")
 # XXX
 
 ## IMPORT
-suppressMessages(library(testit))
+suppressMessages(library(testit)) # assertions
 # suppressMessages(library(tools))
 # suppressMessages(library(gtools))
 # suppressMessages(library(rjson))
-suppressMessages(library(rmarkdown))
+suppressMessages(library(rmarkdown)) #Rmarkdown
 # suppressMessages(library(knitr))
 # suppressMessages(library(DT))
 # suppressMessages(require(grid))
-suppressMessages(library(ggplot2))
-suppressMessages(library(reshape2))
-suppressMessages(library(scales))
+suppressMessages(library(ggplot2)) # plotting
+suppressMessages(library(reshape2)) # reshaping dataframes
+suppressMessages(library(scales)) # plot scales
+suppressMessages(library(pheatmap)) # heatmaps
+suppressMessages(require(grid)) # for heatmaps
+suppressMessages(library(viridis)) # color palette
+
 # custom
 source(snakemake@params$utils)
 
@@ -45,6 +49,12 @@ pdf(snakemake@output$pdf_mappability, width=7, height=5)
 for(pp in FIGS$mappability){ print(pp) }
 dev.off()
 
+# Mash
+pdf(snakemake@output$pdf_mash, width=7, height=7)
+grid.newpage()
+grid.draw(FIGS$mash)
+dev.off()
+
 # QUAST
 pdf(snakemake@output$pdf_quast, width=7, height=10)
 print(FIGS$quast)

workflow_report/scripts/utils.R

@@ -142,6 +142,22 @@ read_cdhit <- function(fname){
     return(df)
 }
 
+read_mash_dist <- function(fname){
+    proc_name <- function(x){
+        x <- dirname(x)
+        x <- basename(x)
+        return(x)
+    }
+    dm <- read.csv(file=fname, sep='\t', header=TRUE, row.names=1, check.names=FALSE)
+    colnames(dm) <- sapply(colnames(dm), proc_name)
+    rownames(dm) <- sapply(rownames(dm), proc_name)
+    testit::assert(all(colnames(dm) %in% names(ASM_TOOL_NAMES)))
+    colnames(dm) <- ASM_TOOL_NAMES[colnames(dm)]
+    rownames(dm) <- ASM_TOOL_NAMES[rownames(dm)]
+    return(dm)
+}
+
+
 ##############################
 # PLOTS