Merge pull request #21 from icbi-lab/Add_volcano

Add module volcano plot

bin/VolcanoPlots_script.R

@@ -0,0 +1,80 @@
+#!/usr/bin/env Rscript
+'
+Usage:
+  VolcanoPlots_script.R --de_res=<de_res> --goi=<goi> --prefix=<prefix> [options]
+
+Mandatory arguments:
+  --de_res=<de_res>           TopTable from DESeq2 in TSV format
+  --goi=<goi>                 Genes of interest in txt format
+  --prefix=<prefix>           Prefix for output filenames
+
+Optional arguments:
+  --pCutoff=<pCutoff>         Cut-off for statistical significance [default: 0.05]
+  --FCcutoff=<FCcutoff>       Cut-off for absolute log2 fold-change [default: 2]
+  --results_dir=<dir>         Output directory [default: ./]
+' -> doc
+
+library(conflicted)
+library(docopt)
+arguments <- docopt(doc, version = "0.1")
+print(arguments)
+
+library(readr)
+library(dplyr)
+library(EnhancedVolcano)
+
+
+# Load parameters
+de_res <- read_tsv(arguments$de_res)
+goi <- read_lines(arguments$goi)
+pCutoff <- as.numeric(arguments$pCutoff)
+FCcutoff <- as.numeric(arguments$FCcutoff)
+prefix <- arguments$prefix
+results_dir <- arguments$results_dir
+
+
+# Make volcano plots using "EnhancedVolcano" package
+
+message(paste0("Drawing volcano plots..."))
+
+p <- EnhancedVolcano(
+  toptable = de_res,
+  lab = NA,
+  x = "log2FoldChange",
+  y = "pvalue",
+  pCutoff = pCutoff,
+  FCcutoff = FCcutoff,
+  title = paste0(prefix, "_volcano_plot"),
+  caption = paste0("fold change cutoff: ", pCutoff, ", p-value cutoff: ", FCcutoff)
+)
+
+ggsave(file.path(results_dir, paste0(prefix, "_volcano_plot.pdf")), plot = p, width = 297, height = 210, units = "mm")
+
+p <- EnhancedVolcano(
+  toptable = de_res,
+  lab = NA,
+  x = "log2FoldChange",
+  y = "padj",
+  pCutoff = pCutoff,
+  FCcutoff = FCcutoff,
+  title = paste0(prefix, "_volcano_plot"),
+  caption = paste0("fold change cutoff: ", pCutoff, ", adj.p-value cutoff:: ", FCcutoff)
+)
+
+ggsave(file.path(results_dir, paste0(prefix, "_volcano_padj.pdf")), plot = p, width = 297, height = 210, units = "mm")
+
+p <- EnhancedVolcano(
+  toptable = de_res,
+  lab = de_res$gene_name,
+  selectLab = goi,
+  x = "log2FoldChange",
+  y = "padj",
+  pCutoff = pCutoff,
+  FCcutoff = FCcutoff,
+  drawConnectors = TRUE,
+  title = paste0(prefix, "_volcano_plot_genes_of_interest"),
+  caption = paste0("fold change cutoff: ", pCutoff, ", adj.p-value cutoff:: ", FCcutoff)
+)
+
+ggsave(file.path(results_dir, paste0(prefix, "_volcano_padj_GoI.pdf")), plot = p, width = 297, height = 210, units = "mm")
+

main.nf

@@ -4,6 +4,7 @@ nextflow.enable.dsl = 2
 
 include { dummyCheckInput } from "./modules/dummy_check_input"
 include { DESeq2 } from "./modules/DESeq2"
+include { VolcanoPlot } from "./modules/VolcanoPlot"
 include { CLUSTERPROFILER_ORA } from "./modules/clusterprofiler_ora"
 
 workflow {
@@ -16,10 +17,13 @@ workflow {
     ch_ora_pathway_dbs = Channel.from(params.ora_pathway_dbs)
     prefix = params.prefix
     de_fdr_cutoff = params.fdr
+    pCutoff = params.pCutoff
+    FCcutoff = params.FCcutoff
 
     // start workflow
     dummyCheckInput(samplesheet, gene_expression, genes_of_interest)
     DESeq2(gene_expression, samplesheet, prefix)
+    VolcanoPlot(DESeq2.out.de_res, genes_of_interest, pCutoff, FCcutoff,prefix)
 
     if (!params.skip_ora) {
         CLUSTERPROFILER_ORA(DESeq2.out.de_res, ch_ora_pathway_dbs, prefix, de_fdr_cutoff)

modules/DESeq2.nf

@@ -1,16 +1,5 @@
 nextflow.enable.dsl=2
 
-def s = """\
-        _________________________________________________
-        Nextflow DESeq2:
-
-        ▶ Inpu count data matrix: ${params.gene_expression}
-        ▶ Input sample annotation data: ${params.samplesheet}
-        ▶ result output directory: ${params.resDir}
-        _________________________________________________
-        """
-println s.stripIndent()
-
 out_dir = file(params.resDir)
 mode = params.publish_dir_mode
 

modules/VolcanoPlot.nf

@@ -0,0 +1,33 @@
+nextflow.enable.dsl=2
+
+out_dir = file(params.resDir)
+mode = params.publish_dir_mode
+
+process VolcanoPlot {
+    //Packages dependencies
+    conda "conda-forge::r-base=4.1.2 conda-forge::r-docopt=0.7.1 conda-forge::r-conflicted=1.1.0 conda-forge::r-readr=2.1.1 conda-forge::dplyr=1.0.7 bioconda::bioconductor-enhancedvolcano=1.12.0"
+    publishDir "${out_dir}", mode: "$mode"
+
+    input:
+        path(de_res)
+	      path(goi)
+	      val(pCutoff)
+	      val(FCcutoff)
+	      val(prefix)
+
+    output:
+        path("${prefix}_volcano_plot.pdf"), emit: volcano
+        path("${prefix}_volcano_padj.pdf"), emit: volcano_padj
+        path("${prefix}_volcano_padj_GoI.pdf"), emit: volcano_GoI
+        path("*.pdf"), emit: plots, optional: true
+
+	script:
+	"""
+    VolcanoPlots_script.R \\
+    --de_res=${de_res} \\
+    --goi=${goi} \\
+    --prefix=${prefix} \\
+    --pCutoff=${pCutoff} \\
+    --FCcutoff=${FCcutoff}
+	"""
+}

nextflow.config

@@ -20,6 +20,8 @@ params {
   fdr = 0.1
   c1 = "grpA"
   c2 = "grpB"
+  pCutoff = 0.05
+  FCcutoff= 2
   readPaths = false
 
   ora_pathway_dbs = ["KEGG", "Reactome", "WikiPathway", "GO_BP", "GO_MF"]