Merge pull request #82 from qbic-pipelines/dev

Release 1.3.2

Author: ggabernet

.github/workflows/ci.yml

@@ -31,13 +31,13 @@ jobs:
 
       - name: Build new docker image
         if: env.GIT_DIFF
-        run: docker build --no-cache . -t qbicpipelines/rnadeseq:1.3.1
+        run: docker build --no-cache . -t qbicpipelines/rnadeseq:1.3.2
 
       - name: Pull docker image
         if: ${{ !env.GIT_DIFF }}
         run: |
           docker pull qbicpipelines/rnadeseq:dev
-          docker tag qbicpipelines/rnadeseq:dev qbicpipelines/rnadeseq:1.3.1
+          docker tag qbicpipelines/rnadeseq:dev qbicpipelines/rnadeseq:1.3.2
 
       - name: Install Nextflow
         run: |
@@ -69,13 +69,13 @@ jobs:
       
       - name: Build new docker image
         if: env.GIT_DIFF
-        run: docker build --no-cache . -t qbicpipelines/rnadeseq:1.3.1
+        run: docker build --no-cache . -t qbicpipelines/rnadeseq:1.3.2
 
       - name: Pull docker image
         if: ${{ !env.GIT_DIFF }}
         run: |
           docker pull qbicpipelines/rnadeseq:dev
-          docker tag qbicpipelines/rnadeseq:dev qbicpipelines/rnadeseq:1.3.1
+          docker tag qbicpipelines/rnadeseq:dev qbicpipelines/rnadeseq:1.3.2
 
       - name: Install Nextflow
         run: |

CHANGELOG.md

@@ -1,5 +1,17 @@
 # qbic-pipelines/rnadeseq: Changelog
 
+## 1.3.2 - Almond Blossoms hotfix II
+
+### Added
+
+- Bump versions to 1.3.2
+- write deseq2 table to file
+
+### Fixed
+- Contrast names in report plots
+- Convert species name to lower case also in report
+- LogFC is also reported in the report and set in volcano plots
+
 ## 1.3.1 - Almond Blossoms hotfix
 
 ### Added

Dockerfile

@@ -9,10 +9,10 @@ RUN apt-get update -qq && \
     apt-get install -y zip procps ghostscript
 
 # Add conda installation dir to PATH
-ENV PATH /opt/conda/envs/qbic-pipelines-rnadeseq-1.3.1/bin:$PATH
+ENV PATH /opt/conda/envs/qbic-pipelines-rnadeseq-1.3.2/bin:$PATH
 
 # Dump the details of the installed packates to a file for posterity
-RUN conda env export --name qbic-pipelines-rnadeseq-1.3.1 > qbic-pipelines-rnadeseq-1.3.1.yml
+RUN conda env export --name qbic-pipelines-rnadeseq-1.3.2 > qbic-pipelines-rnadeseq-1.3.2.yml
 
 # Instruct R processes to use these empty files instead of clashing with a local config
 RUN touch .Rprofile

assets/RNAseq_report.Rmd

@@ -22,6 +22,7 @@ params:
   path_quote: ""
   organism: ""
   batch_effect: ""
+  log_FC: ""
   revision: ""
 # Author: Silvia Morini, Gisela Gabernet, Simon Heumos
 ---
@@ -54,7 +55,8 @@ description = as.character(df$Description)
 library(yaml)
 report_options <- read_yaml(params$path_report_options)
 design_deseq2 <- paste(readLines(params$path_design), collapse=" ")
-#if(batch_effect=="Yes"){batch_effect=TRUE}else{batch_effect=FALSE}
+log_FC_text <- as.character(params$log_FC)
+log_FC_num <- as.numeric(params$log_FC)
 ```
 
 ---
@@ -355,7 +357,8 @@ for (i in 1:length(contrnames)){
 ## List of differentially expressed genes
 
 The analysis of the differential gene expression was performed using `DESeq2`. 
-Genes were considered differentially expressed (DE) when the adjusted p-value was lower than 0.05 (padj < 0.05). 
+Genes were considered differentially expressed (DE) when the adjusted p-value was lower than 0.05 (padj < 0.05).
+Genes were further considered differentially expressed (DE) if the log2 Fold Change in expression between the two considered groups was above the threshold of `r log_FC_text` in absolute terms.
 The adjusted p-value is calculated in the `DESeq2` package with the Benjamini-Hochberg method. P-value adjustment 
 helps reduce the number of false postives (not real differentially expressed genes). 
 
@@ -389,7 +392,7 @@ kable(df, escape=F,row.names =F) %>%
 
 Volcano plots display the DE genes in Log2 Fold Change values (x axis) 
 against their adjusted p-value (y axis) in a form of a -log10. Only significantly DE genes are shown (padj < 0.05). 
-It is recommended to select differentially expressed genes with a Log2FC > 1 or < -1 (blue lines show that threshold).
+Genes were further considered DE with a Log2FC > `r log_FC_text` or < -`r log_FC_text` (the blue lines show this threshold).
 Here, volcano plots for all considered contrasts are shown. Hover over the dots to show the gene names.
 
 ```{r, echo=F, message =F}
@@ -417,7 +420,7 @@ Here, volcano plots for all considered contrasts are shown. Hover over the dots
     # label: gene_name
     # x: log2FoldChange
     # y: -log10(padj)
-    table_list <- append(table_list, list(DE_genes), 0)
+    table_list <- append(table_list, list(DE_genes))
   }
   DE_all <- ldply(table_list, rbind)
   DE_all$logpval <- -log10(DE_all$padj)
@@ -426,8 +429,8 @@ Here, volcano plots for all considered contrasts are shown. Hover over the dots
   pg <- ggplot(DE_all, aes(x=log2FoldChange, y=logpval, text=paste("Gene: ", gene_name, "<br>", "Log2FC: ", formatC(log2FoldChange, digits=2)))) +
         geom_point(aes(alpha=0.3)) +
         geom_hline(yintercept = 16, linetype= "dashed", size = 0.2, color = "grey") + 
-        geom_vline(xintercept = 1, size = 0.2, color = "blue") +
-        geom_vline(xintercept = -1, size = 0.2, color = "blue") +
+        geom_vline(xintercept = log_FC_num, size = 0.2, color = "blue") +
+        geom_vline(xintercept = -log_FC_num, size = 0.2, color = "blue") +
         scale_y_continuous(limits = c(1, 18)) +
         scale_x_continuous(limits = c(-5,5), breaks = c(-5:5)) +
         facet_wrap(~contrast, ncol=1, scales = "free_x", shrink = FALSE) +
@@ -511,15 +514,15 @@ The table below provides more detail on all enriched pathways."))
       next
     }
     q = as.character(DE_genes$Ensembl_ID)
-    q_list <- append(q_list, list(q), 0)
+    q_list <- append(q_list, list(q))
     q_names <- append(q_names, fname)
   }
   datasources <- c("KEGG","REAC")
   names(q_list) <- q_names
   
   #gost query
   gostres <- gost(query=q_list,
-                  organism=params$organism,
+                  organism=tolower(params$organism),
                   significant=T,
                   correction_method="fdr",
                   sources=datasources,

bin/DESeq2.R

@@ -177,6 +177,9 @@ cds <- DESeq(cds,  parallel = FALSE)
 # SizeFactors(cds) as indicator of library sequencing depth
 write.table(sizeFactors(cds),paste("differential_gene_expression/gene_counts_tables/sizeFactor_libraries.tsv",sep=""), append = FALSE, quote = FALSE, sep = "\t",eol = "\n", na = "NA", dec = ".", row.names = T,  col.names = F, qmethod = c("escape", "double"))
 
+# Write cds assay table to file
+write.table(counts(cds, normalized=T), paste("differential_gene_expression/gene_counts_tables/deseq2_table.tsv", sep=""), append=F, quote = F, sep = "\t", eol = "\n", na = "NA", dec = ".", row.names = T, col.names = T, qmethod = c("escape", "double"))
+
 # Write raw counts to file
 count_table_names <- merge(x=gene_names, y=count.table, by.x = "Ensembl_ID", by.y="row.names")
 write.table(count_table_names, paste("differential_gene_expression/gene_counts_tables/raw_gene_counts.tsv",sep=""), append = FALSE, quote = FALSE, sep = "\t",eol = "\n", na = "NA", dec = ".", row.names = F, qmethod = c("escape", "double"))

bin/Execute_report.R

@@ -15,6 +15,7 @@ option_list = list(
   make_option(c("-q", "--quote"), type="character", default=NULL, help="path to the signed quote PDF file", metavar="character"),
   make_option(c("-g", "--organism"), type="character", default=NULL, help="Organism, e.g. Hsapiens."),
   make_option(c("-b", "--batch_effect"), action="store_true", default=FALSE, help="Batch effect correction."),
+  make_option(c("-f", "--log_FC"), type="double", default=NULL, help="Log Fold Change threshold to consider a gene DE."),
   make_option(c("-x", "--revision"), type="character", default=NULL, help="rnadeseq workflow revision", metavar="character")
 )
 
@@ -40,4 +41,5 @@ rmarkdown::render(opt$report, output_file = opt$output, knit_root_dir = wd, outp
                                 path_quote = opt$quote,
                                 organism = opt$organism,
                                 batch_effect = opt$batch_effect,
+                                log_FC = opt$log_FC,
                                 revision = opt$revision))

bin/pathway_analysis.R

@@ -117,12 +117,13 @@ theme_set(theme_classic())
 
 # For each contrast do pathway analysis
 
+# TODO debugging summary table
 # Defining summary variables
-contrast <- c()
-number_DE_genes <- c()
-number_enriched_pathways <- c()
-list_DE_genes <- list()
-list_enriched_pathways <- list()
+# contrast <- c()
+# number_DE_genes <- c()
+# number_enriched_pathways <- c()
+# list_DE_genes <- list()
+# list_enriched_pathways <- list()
 
 for (file in contrast_files){
   #Reading DE genes list
@@ -183,12 +184,13 @@ for (file in contrast_files){
               file = paste0(outdir, "/", fname, "/", fname, "_pathway_enrichment_results.tsv"), 
               sep="\t", quote = F, col.names = T, row.names = F)
 
+  # TODO debugging summary table
   # Collecting summary variables
-  contrast <- append(contrast, fname)
-  number_DE_genes <- append(number_DE_genes, length(DE_genes$Ensembl_ID))
-  number_enriched_pathways <- append(number_enriched_pathways, summary(as.factor(pathway_gostres_table$source)))
-  list_DE_genes <- append(list_DE_genes, DE_genes$Ensembl_ID)
-  list_enriched_pathways <- append(list_enriched_pathways, pathway_gostres_table$term_id)
+  # contrast <- append(contrast, fname)
+  # number_DE_genes <- append(number_DE_genes, length(DE_genes$Ensembl_ID))
+  # number_enriched_pathways <- append(number_enriched_pathways, summary(as.factor(pathway_gostres_table$source)))
+  # list_DE_genes <- append(list_DE_genes, DE_genes$Ensembl_ID)
+  # list_enriched_pathways <- append(list_enriched_pathways, pathway_gostres_table$term_id)
 
   # Printing summary variables
   print("------------------------------------")
@@ -289,10 +291,10 @@ for (file in contrast_files){
 }
 
 # Writing DE gene summary table
-df_summary <- data.frame(Contrast = contrast, Number_DE_genes = number_DE_genes, Number_enriched_pathways = number_enriched_pathways, DE_genes_list = list_DE_genes, Enriched_pathways_list = list_enriched_pathways)
-write.table(df_summary, 
-        file = paste0(outdir, "/", fname, "/", fname, "_pathway_enrichment_summary.tsv"), 
-        sep="\t", quote = F, col.names = T, row.names = F)
+# df_summary <- data.frame(Contrast = contrast, Number_DE_genes = number_DE_genes, Number_enriched_pathways = number_enriched_pathways, DE_genes_list = list_DE_genes, Enriched_pathways_list = list_enriched_pathways)
+# write.table(df_summary, 
+#         file = paste0(outdir, "/", fname, "/", fname, "_pathway_enrichment_summary.tsv"), 
+#         sep="\t", quote = F, col.names = T, row.names = F)
 
 
 # Plotting heatmap for provided gene list

environment.yml

@@ -1,7 +1,7 @@
 # You can use this file to create a conda environment for this pipeline:
 #   conda env create -f environment.yml
 # use this to find packages: https://anaconda.org/
-name: qbic-pipelines-rnadeseq-1.3.1
+name: qbic-pipelines-rnadeseq-1.3.2
 channels:
   - conda-forge
   - bioconda

main.nf

@@ -333,9 +333,17 @@ process Report {
     mkdir QC
     mv MultiQC/multiqc_plots/ MultiQC/multiqc_data/ MultiQC/multiqc_report.html QC/
     Execute_report.R --report '$baseDir/assets/RNAseq_report.Rmd' \
-    --output 'RNAseq_report.html' --proj_summary $proj_summary \
-    --versions $softwareversions --model $model --report_options $report_options --revision $workflow.revision \
-    --contrasts $contrnames $genelistopt --organism $params.species $batchopt
+    --output 'RNAseq_report.html' \
+    --proj_summary $proj_summary \
+    --versions $softwareversions \
+    --model $model \
+    --report_options $report_options \
+    --revision $workflow.revision \
+    --contrasts $contrnames \
+    $genelistopt \
+    --organism $params.species \
+    --log_FC $params.logFCthreshold \
+    $batchopt
     zip -r report.zip RNAseq_report.html differential_gene_expression/ QC/ pathway_analysis/ $quoteopt
     """
 }

nextflow.config

@@ -46,7 +46,7 @@ params {
 
 // Container slug. Stable releases should specify release tag!
 // Developmental code should specify :dev
-process.container = 'qbicpipelines/rnadeseq:1.3.1'
+process.container = 'qbicpipelines/rnadeseq:1.3.2'
 
 // Load base.config by default for all pipelines
 includeConfig 'conf/base.config'
@@ -102,7 +102,7 @@ manifest {
   description = 'Downstream differential gene expression analysis with DESeq2'
   mainScript = 'main.nf'
   nextflowVersion = '>=19.04'
-  version = '1.3.1'
+  version = '1.3.2'
 }
 
 // Function to ensure that resource requirements don't go beyond