Merge pull request #595 from maxplanck-ie/develop

Develop

.ci_stuff/test_dag.sh

@@ -54,6 +54,9 @@ mkdir -p allelic_input
 mkdir -p allelic_input/Ngenome
 touch allelic_input/file.vcf.gz allelic_input/snpfile.txt
 
+# Ensure an empty snakePipes config doesn't muck anything up
+snakePipes config
+
 # DNA mapping
 WC=`DNA-mapping -i PE_input -o output .ci_stuff/organism.yaml --snakemakeOptions " --dryrun --conda-prefix /tmp" | tee >(cat 1>&2) | grep -v "Conda environment" | wc -l`
 if [ ${PIPESTATUS[0]} -ne 0 ] || [ $WC -ne 886 ]; then exit 1 ; fi
@@ -65,6 +68,8 @@ WC=`DNA-mapping -i PE_input -o output .ci_stuff/organism.yaml --snakemakeOptions
 if [ ${PIPESTATUS[0]} -ne 0 ] || [ $WC -ne 970 ]; then exit 1 ; fi
 WC=`DNA-mapping -i PE_input -o output .ci_stuff/organism.yaml --snakemakeOptions " --dryrun --conda-prefix /tmp" --trim --mapq 20 --UMIDedup --properPairs | tee >(cat 1>&2) | grep -v "Conda environment" | wc -l`
 if [ ${PIPESTATUS[0]} -ne 0 ] || [ $WC -ne 1010 ]; then exit 1 ; fi
+WC=`DNA-mapping -i PE_input -o output .ci_stuff/organism.yaml --snakemakeOptions " --dryrun --conda-prefix /tmp" --DAG --trim --mapq 20 --UMIDedup --properPairs | tee >(cat 1>&2) | grep -v "Conda environment" | wc -l`
+if [ ${PIPESTATUS[0]} -ne 0 ] || [ $WC -ne 1010 ]; then exit 1 ; fi
 WC=`DNA-mapping -i SE_input -o output .ci_stuff/organism.yaml --snakemakeOptions " --dryrun --conda-prefix /tmp" | tee >(cat 1>&2) | grep -v "Conda environment" | wc -l`
 if [ ${PIPESTATUS[0]} -ne 0 ] || [ $WC -ne 786 ]; then exit 1 ; fi
 WC=`DNA-mapping -i SE_input -o output .ci_stuff/organism.yaml --snakemakeOptions " --dryrun --conda-prefix /tmp" --trim --mapq 20 --dedup --properPairs | tee >(cat 1>&2) | grep -v "Conda environment" | wc -l`
@@ -140,6 +145,8 @@ WC=`WGBS -i PE_input -o output --sampleSheet .ci_stuff/test_sampleSheet.tsv --sn
 if [ ${PIPESTATUS[0]} -ne 0 ] || [ $WC -ne 970 ]; then exit 1 ; fi
 WC=`WGBS -i BAM_input/filtered_bam -o output --sampleSheet .ci_stuff/test_sampleSheet.tsv --fromBAM --snakemakeOptions " --dryrun --conda-prefix /tmp" --GCbias .ci_stuff/organism.yaml | tee >(cat 1>&2) | grep -v "Conda environment" | wc -l`
 if [ ${PIPESTATUS[0]} -ne 0 ] || [ $WC -ne 690 ]; then exit 1 ; fi
+WC=`WGBS -i BAM_input/filtered_bam -o output --sampleSheet .ci_stuff/test_sampleSheet.tsv --fromBAM --fastqc --snakemakeOptions " --dryrun --conda-prefix /tmp" --GCbias .ci_stuff/organism.yaml | tee >(cat 1>&2) | grep -v "Conda environment" | wc -l`
+if [ ${PIPESTATUS[0]} -ne 0 ] || [ $WC -ne 690 ]; then exit 1 ; fi
 
 # ATAC-seq
 WC=`ATAC-seq -d BAM_input --sampleSheet .ci_stuff/test_sampleSheet.tsv --snakemakeOptions " --dryrun --conda-prefix /tmp" .ci_stuff/organism.yaml | tee >(cat 1>&2) | grep -v "Conda environment" | wc -l`
@@ -168,11 +175,15 @@ if [ ${PIPESTATUS[0]} -ne 0 ] || [ $WC -ne 508 ]; then exit 1 ; fi
 # preprocessing
 WC=`preprocessing -i PE_input -o output --snakemakeOptions " --dryrun --conda-prefix /tmp"  --fastqc --optDedupDist 2500 | tee >(cat 1>&2) | grep -v "Conda environment" | wc -l`
 if [ ${PIPESTATUS[0]} -ne 0 ] || [ $WC -ne 488 ]; then exit 1 ; fi
+WC=`preprocessing -i PE_input -o output --snakemakeOptions " --dryrun --conda-prefix /tmp"  --DAG --fastqc --optDedupDist 2500 | tee >(cat 1>&2) | grep -v "Conda environment" | wc -l`
+if [ ${PIPESTATUS[0]} -ne 0 ] || [ $WC -ne 488 ]; then exit 1 ; fi
 
 # createIndices
 WC=`createIndices -o output --snakemakeOptions " --dryrun --conda-prefix /tmp" --genome ftp://ftp.ensembl.org/pub/release-93/fasta/mus_musculus/dna/Mus_musculus.GRCm38.dna_sm.primary_assembly.fa.gz --gtf ftp://ftp.ensembl.org/pub/release-93/gtf/mus_musculus/Mus_musculus.GRCm38.93.gtf.gz blah | tee >(cat 1>&2) | grep -v "Conda environment" | wc -l`
 if [ ${PIPESTATUS[0]} -ne 0 ] || [ $WC -ne 143 ]; then exit 1 ; fi
 WC=`createIndices -o output --snakemakeOptions " --dryrun --conda-prefix /tmp" --genome ftp://ftp.ensembl.org/pub/release-93/fasta/mus_musculus/dna/Mus_musculus.GRCm38.dna_sm.primary_assembly.fa.gz --gtf ftp://ftp.ensembl.org/pub/release-93/gtf/mus_musculus/Mus_musculus.GRCm38.93.gtf.gz --rmskURL http://hgdownload.soe.ucsc.edu/goldenPath/dm6/database/rmsk.txt.gz blah | tee >(cat 1>&2) | grep -v "Conda environment" | wc -l`
 if [ ${PIPESTATUS[0]} -ne 0 ] || [ $WC -ne 150 ]; then exit 1 ; fi
+WC=`createIndices -o output --snakemakeOptions " --dryrun --conda-prefix /tmp" --DAG --genome ftp://ftp.ensembl.org/pub/release-93/fasta/mus_musculus/dna/Mus_musculus.GRCm38.dna_sm.primary_assembly.fa.gz --gtf ftp://ftp.ensembl.org/pub/release-93/gtf/mus_musculus/Mus_musculus.GRCm38.93.gtf.gz --rmskURL http://hgdownload.soe.ucsc.edu/goldenPath/dm6/database/rmsk.txt.gz blah | tee >(cat 1>&2) | grep -v "Conda environment" | wc -l`
+if [ ${PIPESTATUS[0]} -ne 0 ] || [ $WC -ne 150 ]; then exit 1 ; fi
 
 rm -rf SE_input PE_input BAM_input output allelic_input /tmp/genes.gtf /tmp/genome.fa /tmp/genome.fa.fai

.readthedocs.yaml

@@ -1,5 +1,7 @@
-conda:
-  file: docs/environment.yaml
+#conda:
+#  file: docs/environment.yaml
 python:
   version: 3.7
   pip_install: true
+  install:
+      - requirements: docs/requirements.txt

bin/snakePipes

@@ -100,6 +100,11 @@ def parse_arguments():
                                    "(Default: %(default)s)",
                               default=defaults['tempDir'])
 
+    configParser.add_argument('--noToolsVersion',
+                              dest='toolsVersion',
+                              help="By default, tool versions are printed to a workflow-specific file. Specifying this disables that behavior.",
+                              action="store_false")
+
     email = configParser.add_argument_group(title="Email/SMTP options", description="These options are only used if/when --emailAddress is used.")
     email.add_argument('--smtpServer',
                        help="SMTP server address. (Default: %(default)s)",
@@ -302,9 +307,13 @@ def updateConfig(args):
          'onlySSL': args.onlySSL,
          'emailSender': args.emailSender,
          'smtpUsername': args.smtpUsername,
-         'smtpPassword': args.smtpPassword
+         'smtpPassword': args.smtpPassword,
+         'toolsVersion': args.toolsVersion
         }
     baseDir = os.path.dirname(snakePipes.__file__)
+    # Load, update and rewrite the default dictionary
+    currentDict = cof.load_configfile(os.path.join(baseDir, "shared", "defaults.yaml"), False)
+    cof.merge_dicts(currentDict, d)
     cof.write_configfile(os.path.join(baseDir, "shared", "defaults.yaml"), d)
 
 

conda-recipe/meta.yaml

@@ -1,6 +1,6 @@
 package:
   name: snakepipes
-  version: 2.0.0
+  version: 2.0.2
 
 source:
   path: ../

docs/content/News.rst

@@ -1,6 +1,21 @@
 snakePipes News
 ===============
 
+snakePipes 2.0.2
+----------------
+
+ * DAG print is now moved to _after_ workflow run execution such that any error messages from e.g. input file evaluation do not interfere with the DAG and are visible to the user.
+ * Fixed fastqc for --forBAM .
+ * Fixed DESeq2 report failure with just 1 DEG.
+ * Updated links to test data and commands on zenodo in the docs.
+ * SampleSheet check now explicitly checks for tab-delimited header.
+ * Fixed metilene groups, as well methylation density plots in WGBS.
+
+snakePipes 2.0.1
+----------------
+
+ * Fixed a bug in `snakePipes config` that caused the `toolsVersion` variable to be removed from `defaults.yaml`. This is likely related to issue #579.
+
 snakePipes 2.0.0
 ----------------
 

docs/content/setting_up.rst

@@ -290,10 +290,11 @@ Test data
 
 Test data for the various workflows is available at the following locations:
 
- - `DNA mapping <https://zenodo.org/record/1346303>`__
+ - `DNA mapping <https://zenodo.org/record/3707259>`__
  - `ChIP-seq <https://zenodo.org/record/2624281>`__
- - `ATAC-seq <https://zenodo.org/record/2624323>`__
- - `RNA-seq <https://zenodo.org/record/2624408>`__
- - `HiC <https://zenodo.org/record/2624479>`__
- - `WGBS <https://zenodo.org/record/2624498>`__
- - `scRNA-seq <https://zenodo.org/record/2624518>`__
+ - `ATAC-seq <https://zenodo.org/record/3707666>`__
+ - `mRNA-seq <https://zenodo.org/record/3707602>`__
+ - `noncoding-RNA-seq <https://zenodo.org/deposit/3707749>`__
+ - `HiC <https://zenodo.org/record/3707714>`__
+ - `WGBS <https://zenodo.org/record/3707727>`__
+ - `scRNA-seq <https://zenodo.org/record/3707747>`__

docs/index.rst

@@ -36,9 +36,15 @@ Quick start
 
     conda create -n snakePipes -c mpi-ie -c bioconda -c conda-forge snakePipes
 
+* You can update snakePipes to the latest version available on conda with:
+
+.. code:: bash
+
+    conda update -n snakePipes -c mpi-ie -c bioconda -c conda-forge snakePipes
+
 * Download genome fasta and annotations for an your organism, and build indexes, Check in :ref:`createIndices`
 
-* Download example fastq files for the human genome `here <https://zenodo.org/record/1346303>`_
+* Download example fastq files for the human genome `here <https://zenodo.org/record/3707259>`_
 
 * Execute the DNA-mapping pipeline using the example **command.sh** in the test data directory.
 

snakePipes/__init__.py

@@ -1 +1 @@
-__version__ = '2.0.0'
+__version__ = '2.0.2'

snakePipes/common_functions.py

@@ -320,11 +320,11 @@ def check_sample_info_header(sampleSheet_file):
     if not os.path.isfile(sampleSheet_file):
         sys.exit("ERROR: Cannot find sample info file! (--sampleSheet {})\n".format(sampleSheet_file))
     sampleSheet_file = os.path.abspath(sampleSheet_file)
-    ret = open(sampleSheet_file).read().split("\n")[0].split()
+    ret = open(sampleSheet_file).read().split("\n")[0].split("\t")
     if "name" in ret and "condition" in ret:
-        sys.stderr.write("Sample sheet found and format is ok!\n")
+        sys.stderr.write("Sample sheet found and header is ok!\n")
     else:
-        sys.exit("ERROR: Please use 'name' and 'condition' as column headers in sample info file! ({})\n".format(sampleSheet_file))
+        sys.exit("ERROR: Please use 'name' and 'condition' as column headers in sample info file! Please use a tab-delimited file! ({})\n".format(sampleSheet_file))
     return sampleSheet_file
 
 
@@ -506,16 +506,6 @@ def commonYAMLandLogs(baseDir, workflowDir, defaults, args, callingScript):
                                tempDir=cfg["tempDir"],
                                configFile=os.path.join(args.outdir, '{}.config.yaml'.format(workflowName))).split()
 
-    # Produce the DAG if desired
-    if args.createDAG:
-        oldVerbose = config['verbose']
-        config['verbose'] = False
-        write_configfile(os.path.join(args.outdir, '{}.config.yaml'.format(workflowName)), config)
-        DAGproc = subprocess.Popen(" ".join(snakemake_cmd + ["--rulegraph"]), stdout=subprocess.PIPE, shell=True)
-        subprocess.check_call("dot -Tpdf -o{}/{}_pipeline.pdf".format(args.outdir, workflowName), stdin=DAGproc.stdout, shell=True)
-        config['verbose'] = oldVerbose
-        write_configfile(os.path.join(args.outdir, '{}.config.yaml'.format(workflowName)), config)
-
     if args.verbose:
         snakemake_cmd.append("--printshellcmds")
 
@@ -526,6 +516,20 @@ def commonYAMLandLogs(baseDir, workflowDir, defaults, args, callingScript):
     return " ".join(snakemake_cmd)
 
 
+def print_DAG(args, snakemake_cmd, callingScript, defaults):
+    if args.createDAG:
+        config = defaults
+        config.update(vars(args))
+        workflowName = os.path.basename(callingScript)
+        oldVerbose = config['verbose']
+        config['verbose'] = False
+        write_configfile(os.path.join(args.outdir, '{}.config.yaml'.format(workflowName)), config)
+        DAGproc = subprocess.Popen(snakemake_cmd + " --rulegraph ", stdout=subprocess.PIPE, shell=True)
+        subprocess.check_call("dot -Tpdf -o{}/{}_pipeline.pdf".format(args.outdir, workflowName), stdin=DAGproc.stdout, shell=True)
+        config['verbose'] = oldVerbose
+        write_configfile(os.path.join(args.outdir, '{}.config.yaml'.format(workflowName)), config)
+
+
 def logAndExport(args, workflowName):
     """
     Set up logging

snakePipes/shared/rscripts/DESeq2Report.Rmd

@@ -157,7 +157,7 @@ Most of the exploration below is based on the annotated DESeq2 output file gener
     ## Extract data for heatmap
     # order by fold change (by abs foldch if only few top genes requested)
     #print("Preparing data: heatmap")
-    if (nrow(df.filt) > 0) {
+    if (nrow(df.filt) > 1) {
         d <- data.frame(id = rownames(df.filt), padj = df.filt$padj)
         if (length(rownames(d)) < heatmap_topN ) {
             heatmap_topN <- nrow(d)
@@ -178,7 +178,7 @@ Most of the exploration below is based on the annotated DESeq2 output file gener
         ## scaling is not so useful for already normalized data, or?
         #heatmap_data <- scale(heatmap_data, center = TRUE, scale = TRUE)
     } else {
-        print("No DE genes detected!")
+        if(nrow(df.filt) == 0){ message("No DE genes detected!") }
     }
 ```
 
@@ -397,7 +397,7 @@ This plot shows a heatmap of DESeq2-normlized counts for the **top 20** differen
 
 ```{r 'heatmap2', fig.show=TRUE}
 # 7.1 Heatmap topN genes z-score
-if (nrow(df.filt) > 0) {
+if (nrow(df.filt) > 1) {
     pheatmap(heatmap_data,
              cluster_rows = TRUE,
              clustering_method = "average",
@@ -407,7 +407,7 @@ if (nrow(df.filt) > 0) {
              color = colorRampPalette(rev(RColorBrewer::brewer.pal(9,"RdBu")))(255),
              main = sprintf("Heatmap : Top %d DE genes (by p-value) color: z-score  ", heatmap_topN))
 } else {
-  message("No DE genes detected")
+  if(nrow(df.filt) == 0){ message("No DE genes detected!") }
 }
 ```
 

snakePipes/shared/rscripts/WGBS_dmrseq.Rmd

@@ -35,6 +35,9 @@ infiles = sprintf("MethylDackel/%s_CpG.bedGraph", ss$name)
 
 g1 = which(ss$condition == groups[1])
 g2 = which(ss$condition == groups[2])
+
+message(sprintf("Samples belonging to group 1 (%s) are %s",groups[1],paste(ss$name[g1],collapse=" ")))
+message(sprintf("Samples belonging to group 2 (%s) are %s",groups[2],paste(ss$name[g2],collapse=" ")))
 ```
 
 # Overview
@@ -67,7 +70,7 @@ Using a minimum methylation difference of `r minMethDiff` and FDR of `r FDR` the
 
 ```{r warning=FALSE, message=FALSE}
 if (length(regions) > 0) {
-  g = ggplot(as.data.frame(regions), aes(x=beta)) + geom_histogram() + labs(x="Methylation Difference (per-DMR)")
+  g = ggplot(as.data.frame(regions,stringsAsFactors=FALSE), aes(x=beta)) + geom_histogram() + labs(x="Methylation Difference (per-DMR)")
   g = g + geom_vline(xintercept=minMethDiff) + geom_vline(xintercept=-1*minMethDiff)
   g
 } else {
@@ -79,7 +82,7 @@ Similarly, the test statistic, which is used to compute the p-value, is shown be
 
 ```{r warning=FALSE, message=FALSE}
 if (length(regions) > 0) {
-  g = ggplot(as.data.frame(regions), aes(x=stat)) + geom_histogram() + labs(x="Methylation Difference Statistic")
+  g = ggplot(as.data.frame(regions,stringsAsFactors=FALSE), aes(x=stat)) + geom_histogram() + labs(x="Methylation Difference Statistic")
   g
 } else { 
   message('No DMRs found.')
@@ -90,7 +93,7 @@ It can be useful to plot the unadjusted p-values for the DMRs to help in diagnos
 
 ```{r warning=FALSE, message=FALSE}
 if (length(regions) > 0) {
-  g = ggplot(as.data.frame(regions), aes(x=pval)) + geom_histogram(bins=100) + labs(x="Unadjusted p-value")
+  g = ggplot(as.data.frame(regions,stringsAsFactors=FALSE), aes(x=pval)) + geom_histogram() + labs(x="Unadjusted p-value")
   print(g)
 } else {
   message('No DMRs found.')

snakePipes/shared/rscripts/WGBS_metileneQC.Rmd

@@ -18,7 +18,7 @@ gtfFile = snakemake@params[["genes_gtf"]]
 FDR = snakemake@params[["FDR"]]
 minMethDiff = snakemake@params[["minMethDiff"]]
 # Only a single output file is allows for Rmd files
-outputTxtFile = sprintf("%s/DMRs.annotated.txt", snakemake@params[["outdir"]])
+outputTxtFile = sprintf("%s/DMRs.FDR%s.annotated.txt", snakemake@params[["outdir"]], snakemake@params[["FDR"]])
 outputDiffMethyl = sprintf("%s/MethylationDensity.png", snakemake@params[["outdir"]])
 outputQValue = sprintf("%s/QValueDistribution.png", snakemake@params[["outdir"]])
 outputVolcano = sprintf("%s/Volcano.png", snakemake@params[["outdir"]])
@@ -34,7 +34,7 @@ This report summarizes the differentially methylated region (DMR) calling with m
 ```{r echo=FALSE, warning=FALSE, message=FALSE}
 # Read in the input data
 d = read.delim(inputFile)
-colnames(d) = c("CHROM", "START", "END", "qvalue", "MeanDiff", "NumCpGs", "pMWU", "p2DKS", "meanA", "meanB")
+colnames(d)[1:(length(colnames(d))-2)] = c("CHROM", "START", "END", "qvalue", "MeanDiff", "NumCpGs", "pMWU", "p2DKS")
 gr = GRanges(seqnames=d$CHROM,
              ranges=IRanges(start=d$START, end=d$END),
              mcols=d[,c(4:10)])
@@ -59,7 +59,7 @@ d = t(sapply(split(og, queryHits(og)), function(x) {
       distances = paste0(elementMetadata(x)$distance, collapse=",")
       return(as.character(c(idx, gnames, distances)))
 }))
-d = as.data.frame(d)
+d = as.data.frame(d,stringsAsFactors=FALSE)
 colnames(d) = c("idx", "NearestGene", "DistanceToNearestGene")
 d$idx = as.numeric(d$idx)
 gr$NearestGene = as.character("")
@@ -84,17 +84,15 @@ Averge methylation values for candidate DMRs in each of the groups are plotted b
 
 ```{r echo=FALSE, warning=FALSE, message=FALSE}
 # Select the mean* columns, strip "mean" and make long
-d = as.data.frame(elementMetadata(gr)) %>% select(starts_with("mean", ignore.case=FALSE)) %>% 
-  rename_all(funs(gsub("^mean", "", .))) %>% gather("Group", "MeanMethylation")
+d = as.data.frame(elementMetadata(gr,stringsAsFactor=FALSE)) %>% select(starts_with("mean", ignore.case=FALSE)) %>% 
+  rename_all(funs(gsub("^mean_", "", .))) %>% gather("Group", "MeanMethylation")
 
-g = ggplot(d, aes(x=MeanMethylation, fill=Group, color=Group)) + geom_density(alpha=0.3)
+g = ggplot(d, aes(x=MeanMethylation, fill=Group, color=Group, group=Group)) + geom_density(alpha=0.3)
 g = g + theme(text=element_text(size=16),
               axis.text = element_text(size=12),
               axis.title = element_text(size=14))
-g = g + xlab("Mean methylation ratio")
-g = g + scale_fill_manual(values=c("grey28", "red", "darkblue", "darkgreen"))
-g = g + scale_colour_manual(values=c("grey28", "red", "darkblue", "darkgreen"))
-g = g + xlim(0,1)
+g = g + xlab("Mean methylation ratio") + xlim(0,100)
+g = g + scale_fill_manual(values=c("grey28", "red", "darkblue", "darkgreen"),aesthetics = c("colour", "fill"))
 ggsave(outputDiffMethyl, plot=g)
 g
 ```

snakePipes/shared/rules/FastQC.snakefile

@@ -17,7 +17,7 @@ else:
     if pairedEnd:
         rule FastQC:
             input:
-                "FASTQ/{sample}{read}.fastq.gz"
+                "EXTERNAL_BAM/{sample}."+bamExt if fromBAM else "FASTQ/{sample}{read}.fastq.gz"
             output:
                 "FastQC/{sample}{read}_fastqc.html"
             log:
@@ -32,7 +32,7 @@ else:
     else:
         rule FastQC_singleEnd:
             input:
-                "FASTQ/{sample}"+reads[0]+".fastq.gz"
+                "EXTERNAL_BAM/{sample}."+bamExt if fromBAM else "FASTQ/{sample}"+reads[0]+".fastq.gz"
             output:
                 "FastQC/{sample}"+reads[0]+"_fastqc.html"
             params:

snakePipes/shared/rules/WGBS.snakefile

@@ -406,9 +406,9 @@ rule run_metilene:
     benchmark: '{}/.benchmark/run_metilene.benchmark'.format(get_outdir("metilene", minCoverage))
     conda: CONDA_WGBS_ENV
     shell: """
-        echo -e "chrom\tstart\tend\tq-value\tmean methylation difference\tnCpGs\tp (MWU)\tp (2D KS)\tmean_{params.groups[0]}\tmean_{params.groups[1]}" > {output}
-        metilene --groupA {params.groups[0]} \
-                 --groupB {params.groups[1]} \
+        echo -e "chrom\tstart\tend\tq-value\tmean methylation difference\tnCpGs\tp (MWU)\tp (2D KS)\tmean_{params.groups[1]}\tmean_{params.groups[0]}" > {output}
+        metilene --groupA {params.groups[1]} \
+                 --groupB {params.groups[0]} \
                  --maxdist {params.maxDist} \
                  --mincpgs {params.minCpGs} \
                  --minMethDiff {params.minMethDiff} \