Added optional parameter --nsubgenes

CHANGELOG.md

@@ -6,6 +6,7 @@
 
 - Bump versions to 1.4.0dev
 - Add parameter "--min_DE_genes"
+- Add parameter "--nsubgenes"
 
 ## 1.3.2 - Almond Blossoms hotfix II
 

assets/RNAseq_report.Rmd

@@ -24,6 +24,7 @@ params:
   organism: ""
   batch_effect: ""
   log_FC: ""
+  nsub_genes: ""
   revision: ""
 # Author: Silvia Morini, Gisela Gabernet, Simon Heumos
 ---
@@ -58,6 +59,7 @@ report_options <- read_yaml(params$path_report_options)
 design_deseq2 <- paste(readLines(params$path_design), collapse=" ")
 log_FC_text <- as.character(params$log_FC)
 log_FC_num <- as.numeric(params$log_FC)
+nsub_genes_text <-as.character(params$nsub_genes)
 ```
 
 ---
@@ -305,6 +307,7 @@ The raw count table and normalized count tables are available [here](./different
 The differential expression analysis is performed using the raw gene count table. 
 For PCA analysis and heatmap plotting, the regularized logarithm (rlog) normalized gene counts were employed.
 The vst (variant stabilizing transformation) normalized counts are also provided.
+The number of subset genes for the transformation is `r nsub_genes_text`.
 
 ## Principal component analysis (PCA) 
 A PCA plot of the rlog normalized gene expression visualizes the clustering of samples according to their gene expression. 

bin/DESeq2.R

@@ -50,6 +50,7 @@ option_list = list(
   make_option(c("-p", "--contrasts_pairs"), type="character", default=NULL, help="path to contrasts pairs file", metavar="character"),
   make_option(c("-l", "--genelist"), type="character", default=NULL, help="path to gene list file", metavar="character"),
   make_option(c("-t", "--logFCthreshold"), type="integer", default=0, help="Log 2 Fold Change threshold for DE genes", metavar="character"),
+  make_option(c("-n", "--nsubgenes"), type="integer", default=1000, help="subset number of genes for vst", metavar="character"),
   make_option(c("-b", "--batchEffect"), default=FALSE, action="store_true", help="Whether to consider batch effects in the DESeq2 analysis", metavar="character")
 )
 
@@ -375,7 +376,7 @@ write.table(DE_genes_final_table, "differential_gene_expression/final_gene_table
 # rlog transformation
 rld <- rlog(cds, blind=FALSE)
 # vst transformation
-vsd <- vst(cds, blind=FALSE)
+vsd <- vst(cds, blind=FALSE, nsub = opt$nsubgenes)
 
 # write normalized values to a file
 rld_names <- merge(x=gene_names, y=assay(rld), by.x = "Ensembl_ID", by.y="row.names")

bin/Execute_report.R

@@ -16,6 +16,7 @@ option_list = list(
   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("-n", "--nsub_genes"), type="integer", default=NULL, help="subset number of genes for vst."),
   make_option(c("-x", "--revision"), type="character", default=NULL, help="rnadeseq workflow revision", metavar="character"),
   make_option(c("-p", "--min_DEG_pathway"), type="integer", default=NULL, help="min. number of genes DE in a pathway for this pathway to be considered enriched.", metavar="integer")
 )
@@ -44,4 +45,5 @@ rmarkdown::render(opt$report, output_file = opt$output, knit_root_dir = wd, outp
                                 organism = opt$organism,
                                 batch_effect = opt$batch_effect,
                                 log_FC = opt$log_FC,
+                                nsub_genes = opt$nsub_genes,
                                 revision = opt$revision))

main.nf

@@ -39,6 +39,7 @@ def helpMessage() {
       --relevel                     Tsv indicating list of factors (conditions in the metadata table) and the new level on which to relevel the factor. Check contrasts docs.
       --logFCthreshold              Threshold (int) to apply to Log 2 Fold Change to consider a gene as differentially expressed.
       --genelist                    Txt file with list of genes (one per line) of which to plot heatmaps for normalized counts across all samples.
+      --nsubgenes                   Number (int) of genes to subset in the vst step of deseq. Reduce default if encounter errors with smaller datasets
       --batch_effect                Turn on this flag if you wish to consider batch effects. You need to add the batch effect to the linear model too!                
       --quote                       Signed copy of the offer.
       --kegg_blacklist              Txt file with list of pathways (one per line) that should be discarded for the KEGG pathway plotting (e.g. because the xml file in KEGG contains errors).
@@ -267,7 +268,8 @@ process DESeq2 {
     def batch_effect_opt = params.batch_effect ? "--batchEffect" : ''
     """
     DESeq2.R --counts $gene_counts --metadata $metadata --design $model \
-    --logFCthreshold $params.logFCthreshold $relevel_opt $contrast_mat_opt \
+    --logFCthreshold $params.logFCthreshold --nsubgenes $params.nsubgenes \
+    $relevel_opt $contrast_mat_opt \
     $contrast_list_opt $contrast_pairs_opt $gene_list_opt $batch_effect_opt
     zip -r differential_gene_expression.zip differential_gene_expression
     """
@@ -346,6 +348,7 @@ process Report {
     $genelistopt \
     --organism $params.species \
     --log_FC $params.logFCthreshold \
+    --nsub_genes $params.nsubgenes \
     $batchopt \
     --min_DEG_pathway $params.min_DEG_pathway
     zip -r report.zip RNAseq_report.html differential_gene_expression/ QC/ pathway_analysis/ $quoteopt

nextflow.config

@@ -23,6 +23,7 @@ params {
   kegg_blacklist = 'NO_FILE3'
   quote = 'NO_FILE4'
   min_DEG_pathway = 1
+  nsubgenes = 1000
 
   // Boilerplate options
   name = false