Merge pull request #1005 from maxplanck-ie/develop

Develop

.ci_stuff/test_dag.sh

@@ -101,6 +101,30 @@ touch allelic_BAM_input/allelic_bams/sample1.genome1.sorted.bam \
       allelic_BAM_input/allelic_bams/sample5.genome2.sorted.bam.bai \
       allelic_BAM_input/allelic_bams/sample6.genome1.sorted.bam.bai \
       allelic_BAM_input/allelic_bams/sample6.genome2.sorted.bam.bai \
+      allelic_BAM_input/allelic_bams/sample1.allele_flagged.sorted.bam \
+      allelic_BAM_input/allelic_bams/sample1.unassigned.sorted.bam \
+      allelic_BAM_input/allelic_bams/sample2.allele_flagged.sorted.bam \
+      allelic_BAM_input/allelic_bams/sample2.unassigned.sorted.bam \
+      allelic_BAM_input/allelic_bams/sample3.allele_flagged.sorted.bam \
+      allelic_BAM_input/allelic_bams/sample3.unassigned.sorted.bam \
+      allelic_BAM_input/allelic_bams/sample4.allele_flagged.sorted.bam \
+      allelic_BAM_input/allelic_bams/sample4.unassigned.sorted.bam \
+      allelic_BAM_input/allelic_bams/sample5.allele_flagged.sorted.bam \
+      allelic_BAM_input/allelic_bams/sample5.unassigned.sorted.bam \
+      allelic_BAM_input/allelic_bams/sample6.allele_flagged.sorted.bam \
+      allelic_BAM_input/allelic_bams/sample6.unassigned.sorted.bam \
+      allelic_BAM_input/allelic_bams/sample1.allele_flagged.sorted.bam.bai \
+      allelic_BAM_input/allelic_bams/sample1.unassigned.sorted.bam.bai \
+      allelic_BAM_input/allelic_bams/sample2.allele_flagged.sorted.bam.bai \
+      allelic_BAM_input/allelic_bams/sample2.unassigned.sorted.bam.bai \
+      allelic_BAM_input/allelic_bams/sample3.allele_flagged.sorted.bam.bai \
+      allelic_BAM_input/allelic_bams/sample3.unassigned.sorted.bam.bai \
+      allelic_BAM_input/allelic_bams/sample4.allele_flagged.sorted.bam.bai \
+      allelic_BAM_input/allelic_bams/sample4.unassigned.sorted.bam.bai \
+      allelic_BAM_input/allelic_bams/sample5.allele_flagged.sorted.bam.bai \
+      allelic_BAM_input/allelic_bams/sample5.unassigned.sorted.bam.bai \
+      allelic_BAM_input/allelic_bams/sample6.allele_flagged.sorted.bam.bai \
+      allelic_BAM_input/allelic_bams/sample6.unassigned.sorted.bam.bai \
       allelic_BAM_input/deepTools_qc/bamPEFragmentSize/fragmentSize.metric.tsv \
       allelic_BAM_input/filtered_bam/sample1.filtered.bam \
       allelic_BAM_input/filtered_bam/sample2.filtered.bam \
@@ -197,18 +221,18 @@ if [ ${PIPESTATUS[0]} -ne 0 ] || [ $WC -ne 628 ]; then exit 1 ; fi
 WC=`ChIP-seq -d BAM_input --sampleSheet .ci_stuff/test_sampleSheet.tsv --snakemakeOptions " --dryrun --conda-prefix /tmp" --bigWigType log2ratio .ci_stuff/organism.yaml .ci_stuff/ChIP.sample_config.yaml | tee >(cat 1>&2) | grep -v "conda installation\|Conda environment" | sed '/^\s*$/d' | wc -l`
 if [ ${PIPESTATUS[0]} -ne 0 ] || [ $WC -ne 562 ]; then exit 1 ; fi
 WC=`ChIP-seq -d BAM_input --sampleSheet .ci_stuff/test_sampleSheet.tsv --snakemakeOptions " --dryrun --conda-prefix /tmp" --peakCaller SEACR .ci_stuff/organism.yaml .ci_stuff/ChIP.sample_config.yaml | tee >(cat 1>&2) | grep -v "conda installation\|Conda environment" | sed '/^\s*$/d' | wc -l`
-if [ ${PIPESTATUS[0]} -ne 0 ] || [ $WC -ne 759 ]; then exit 1 ; fi
+if [ ${PIPESTATUS[0]} -ne 0 ] || [ $WC -ne 815 ]; then exit 1 ; fi
 WC=`ChIP-seq -d BAM_input --sampleSheet .ci_stuff/test_sampleSheet.tsv --snakemakeOptions " --dryrun --conda-prefix /tmp" --peakCaller SEACR --useSpikeInForNorm .ci_stuff/spikein_organism.yaml .ci_stuff/ChIP.sample_config.yaml | tee >(cat 1>&2) | grep -v "conda installation\|Conda environment" | sed '/^\s*$/d' | wc -l`
-if [ ${PIPESTATUS[0]} -ne 0 ] || [ $WC -ne 1203 ]; then exit 1 ; fi
+if [ ${PIPESTATUS[0]} -ne 0 ] || [ $WC -ne 1259 ]; then exit 1 ; fi
 #noInput
 WC=`ChIP-seq -d BAM_input --sampleSheet .ci_stuff/test_sampleSheet.tsv --snakemakeOptions " --dryrun --conda-prefix /tmp" .ci_stuff/organism.yaml .ci_stuff/ChIP.sample_noControl_config.yaml | tee >(cat 1>&2) | grep -v "conda installation\|Conda environment" | sed '/^\s*$/d' | wc -l`
 if [ ${PIPESTATUS[0]} -ne 0 ] || [ $WC -ne 403 ]; then exit 1 ; fi
 WC=`ChIP-seq -d BAM_input --sampleSheet .ci_stuff/test_sampleSheet.tsv --snakemakeOptions " --dryrun --conda-prefix /tmp" --peakCaller Genrich .ci_stuff/organism.yaml .ci_stuff/ChIP.sample_noControl_config.yaml | tee >(cat 1>&2) | grep -v "conda installation\|Conda environment" | sed '/^\s*$/d' | wc -l`
 if [ ${PIPESTATUS[0]} -ne 0 ] || [ $WC -ne 349 ]; then exit 1 ; fi
 WC=`ChIP-seq -d BAM_input --sampleSheet .ci_stuff/test_sampleSheet.tsv --snakemakeOptions " --dryrun --conda-prefix /tmp" --peakCaller SEACR .ci_stuff/organism.yaml .ci_stuff/ChIP.sample_noControl_config.yaml | tee >(cat 1>&2) | grep -v "conda installation\|Conda environment" | sed '/^\s*$/d' | wc -l`
-if [ ${PIPESTATUS[0]} -ne 0 ] || [ $WC -ne 469 ]; then exit 1 ; fi
+if [ ${PIPESTATUS[0]} -ne 0 ] || [ $WC -ne 525 ]; then exit 1 ; fi
 WC=`ChIP-seq -d BAM_input --sampleSheet .ci_stuff/test_sampleSheet.tsv --snakemakeOptions " --dryrun --conda-prefix /tmp" --peakCaller SEACR --useSpikeInForNorm .ci_stuff/spikein_organism.yaml .ci_stuff/ChIP.sample_noControl_config.yaml | tee >(cat 1>&2) | grep -v "conda installation\|Conda environment" | sed '/^\s*$/d' | wc -l`
-if [ ${PIPESTATUS[0]} -ne 0 ] || [ $WC -ne 696 ]; then exit 1 ; fi
+if [ ${PIPESTATUS[0]} -ne 0 ] || [ $WC -ne 752 ]; then exit 1 ; fi
 # fromBAM
 WC=`ChIP-seq -d outdir --sampleSheet .ci_stuff/test_sampleSheet.tsv --snakemakeOptions " --dryrun --conda-prefix /tmp" --fromBAM BAM_input/filtered_bam/ .ci_stuff/organism.yaml .ci_stuff/ChIP.sample_config.yaml | tee >(cat 1>&2) | grep -v "conda installation\|Conda environment" | sed '/^\s*$/d' | wc -l`
 if [ ${PIPESTATUS[0]} -ne 0 ] || [ $WC -ne 1187 ]; then exit 1 ; fi
@@ -314,6 +338,10 @@ WC=`mRNA-seq -m allelic-mapping,deepTools_qc -i PE_input -o output --sampleSheet
 if [ ${PIPESTATUS[0]} -ne 0 ] || [ $WC -ne 2539 ]; then exit 1 ; fi
 WC=`mRNA-seq -m allelic-mapping,deepTools_qc,alignment-free -i PE_input -o output --sampleSheet .ci_stuff/test_sampleSheet.tsv --snakemakeOptions " --dryrun --conda-prefix /tmp" --VCFfile allelic_input/file.vcf.gz --strains strain1 .ci_stuff/organism.yaml | tee >(cat 1>&2) | grep -v "conda installation\|Conda environment" | sed '/^\s*$/d' | wc -l`
 if [ ${PIPESTATUS[0]} -ne 0 ] || [ $WC -ne 3294 ]; then exit 1 ; fi
+WC=`mRNA-seq -m allelic-counting -i allelic_BAM_input/allelic_bams --fromBAM --bamExt '.sorted.bam' -o output --sampleSheet .ci_stuff/test_sampleSheet.tsv --snakemakeOptions " --dryrun --conda-prefix /tmp"  .ci_stuff/organism.yaml | tee >(cat 1>&2) | grep -v "conda installation\|Conda environment" | sed '/^\s*$/d' | wc -l`
+if [ ${PIPESTATUS[0]} -ne 0 ] || [ $WC -ne 659 ]; then exit 1 ; fi
+
+
 
 WC=`noncoding-RNA-seq -i PE_input -o output --snakemakeOptions " --dryrun --conda-prefix /tmp" .ci_stuff/organism.yaml | tee >(cat 1>&2) | grep -v "conda installation\|Conda environment" | sed '/^\s*$/d' | wc -l`
 if [ ${PIPESTATUS[0]} -ne 0 ] || [ $WC -ne 1370 ]; then exit 1 ; fi

.github/workflows/linux.yml

@@ -127,4 +127,4 @@ jobs:
       run: |
         micromamba activate snakePipes_CI
         snakePipes config --tempDir /tmp
-        snakePipes createEnvs --autodetectCondaEnvDir --force --only ${{matrix.envs}}
+        snakePipes createEnvs --autodetectCondaEnvDir --only ${{matrix.envs}}

.github/workflows/osx.yml

@@ -52,5 +52,5 @@ jobs:
       - name: createEnvsOSX
         run: |
           micromamba activate snakePipes_CI
-          snakePipes createEnvs --force --autodetectCondaEnvDir --only ${{matrix.envs}}
+          snakePipes createEnvs --autodetectCondaEnvDir --only ${{matrix.envs}}
      

bin/snakePipes

@@ -54,11 +54,11 @@ def parse_arguments():
         "possible environments are: {}".format(cof.set_env_yamls().keys()),
     )
 
-    createEnvsParser.add_argument(
-        "--force",
-        action="store_true",
-        help="Force creation of conda environments, even if they apparently exist.",
-    )
+#    createEnvsParser.add_argument(
+#        "--force",
+#        action="store_true",
+#        help="Force creation of conda environments, even if they apparently exist.",
+#    )
 
     createEnvsParser.add_argument(
         "--info",
@@ -334,15 +334,14 @@ def createCondaEnvs(args):
             "mamba",
             "env",
             "create",
-            "--force",
             "--file",
             os.path.join(baseDir, "shared/rules", env),
         ]
         cmd += ["--prefix", os.path.join(condaDirUse, h)]
 
         # Don't actually create the env if either --info is set or it already exists and --force is NOT set
         if not args.info:
-            if not os.path.exists(os.path.join(condaDirUse, h)) or args.force:
+            if not os.path.exists(os.path.join(condaDirUse, h)):
                 try:
                     os.makedirs(os.path.join(condaDirUse, h), exist_ok=True)
                     subprocess.check_call(cmd)

conda-recipe/meta.yaml

@@ -1,6 +1,6 @@
 package:
   name: snakepipes
-  version: 2.8.1
+  version: 2.8.2
 
 source:
   path: ../

docs/content/News.rst

@@ -1,5 +1,22 @@
 snakePipes News
 ===============
+
+snakePipes 2.8.2
+________________
+
+* added SEACR peaks qc
+* added external bed functionality for differential binding in ChIP-seq and ATAC-seq workflows
+* added "allelic-counting" mode to mRNA-seq, allowing to count reads and run DGE from allelic bam files split e.g. by whatshap
+* added support for custom model formula to mRNA-seq workflow
+* fixed copyfile command for sampleSheet
+* removed deprecated --force argument from mamba commands
+* fixes #998
+* fixes #997
+* fixes #996
+* fixes #994
+* fixes #1000
+* fixes #1001
+
 snakePipes 2.8.1
 ----------------
 * Boosted versions on shared_env, as python 3.7 and multiqc no longer work together.

docs/content/workflows/mRNA-seq.rst

@@ -144,6 +144,9 @@ Like the other workflows, differential expression can be performed using the ``-
 
 .. note:: The first entry defines which group of samples are control. This way, the order of comparison and likewise the sign of values can be changed. The DE analysis might fail if your sample names begin with a number. So watch out for that!
 
+
+Optionally, the user may submit their desired model formula (without the leading tilda ``~``)  with ``--formula``.
+
 Differential Splicing
 ---------------------
 
@@ -194,6 +197,11 @@ Allele-specific, gene-level differential expression analysis is then performed u
 
 .. note:: **allelic-mapping** mode is mutually exclusive with **mapping** mode
 
+"allelic-counting"
+~~~~~~~~~~~~~~~~~~
+
+**allelic-counting** mode requires the user to input, per sample, 4 bam files, corresponding to haplotype1, haplotype2, unassigned and haplotagged , e.g. as generated by whatshap. The respective suffixes ".genome1", ".genome2", ".unassigned", ".alelle_flagged" are required to be followed by the bam extention ".sorted.bam". This mode is mutually exclusive with "deepTools_qc". Only the allelic version of deepTools qc will be run, by default. Allelic version of featureCounts will be run by default. If sample sheet is provided, allelic DESeq2- or allelic Salmon-based differential gene expression analysis will be run. 
+
 "alignment-free"
 ~~~~~~~~~~~~~~~~
 

snakePipes/__init__.py

@@ -1 +1 @@
-__version__ = '2.8.1'
+__version__ = '2.8.2'

snakePipes/common_functions.py

@@ -173,6 +173,20 @@ def get_sample_names_bam(infiles, bamExt):
     return sorted(list(set(s)))
 
 
+def get_sample_names_suffix_bam(infiles, bamExt):
+    """
+    Get sample names without file extensions
+    """
+    bamSuff = [x + bamExt for x in [".genome1", ".genome2", ".unassigned", ".allele_flagged"]]
+    s = []
+    for x in infiles:
+        for y in bamSuff:
+            if y in os.path.basename(x):
+                x = os.path.basename(x).replace(y, "")
+                s.append(x)
+    return sorted(list(set(s)))
+
+
 def is_paired(infiles, ext, reads):
     """
     Check for paired-end input files
@@ -848,7 +862,7 @@ def copySampleSheet(sampleSheet, wdir):
     if os.path.isfile(sampleSheet) and os.path.exists(wdir):
         bname = os.path.basename(sampleSheet)
         try:
-            shutil.copy(sampleSheet, os.path.join(wdir, bname))
+            shutil.copyfile(sampleSheet, os.path.join(wdir, bname))
         except Exception as err:
             print("Unexpected error:\n{}".format(err))
             raise

snakePipes/parserCommon.py

@@ -168,7 +168,8 @@ def snpArguments(defaults):
 
     snpargs.add_argument("--NMaskedIndex",
                          default='',
-                         help="N-masked index of the reference genome (default: 'None')")
+                         help="N-masked index of the reference genome (default: 'None'). "
+                         "Note that this should point to a file (i.e. 'Genome' for STAR indices, genome.1.bt2 for bowtie2 indices).")
 
     return parser
 

snakePipes/shared/defaults.yaml

@@ -7,7 +7,7 @@
 # permitted here.
 ################################################################################
 #
-condaEnvDir: '/package/mamba/envs/snakepipesenvs/2.8.1'
+condaEnvDir: '/package/mamba/envs/snakepipesenvs/2.8.2'
 snakemakeOptions: ''
 organismsDir: 'shared/organisms'
 clusterConfig: 'shared/cluster.yaml'

snakePipes/shared/rscripts/CSAW.R

@@ -15,6 +15,9 @@ outdir<-snakemake@params[["outdir"]]
 yaml_path<-snakemake@params[["yaml_path"]]
 useSpikeInForNorm<-snakemake@params[["useSpikeInForNorm"]]
 scale_factors<-snakemake@params[["scale_factors"]]
+external_bed<-as.logical(snakemake@params[["externalBed"]])
+
+
 
 bam_pfx<-ifelse(useSpikeInForNorm,"_host",".filtered")
 bam_folder<-ifelse(useSpikeInForNorm,"split_bam","filtered_bam")
@@ -42,6 +45,7 @@ message(paste("FDR:", fdr, "\n"))
 message(paste("LFC:", lfc, "\n"))
 message(paste("paired-end? :", pairedEnd, "\n"))
 message(paste("allele-specific? :", allelic_info, "\n"))
+message(paste("External bed? ;", external_bed, "\n"))
 
 ## sampleInfo (setup of the experiment)
 sampleInfo <- read.table(sampleInfoFilePath, header = TRUE, colClasses = c("character", "character"))
@@ -98,35 +102,43 @@ if (!is.null(sampleInfo$UseRegions)) {
     fnames<-sampleInfo$name
 }
 
-if(snakemake@params[['peakCaller']] == "MACS2") {
-    allpeaks <- lapply(fnames, function(x) {
-        narrow <- paste0("../MACS2/",x,bam_pfx,".BAM_peaks.narrowPeak") #bam_pfx
-        if(snakemake@params[["pipeline"]] %in% "ATAC-seq"){
-            narrow <- paste0("../MACS2/",x,".filtered.short.BAM_peaks.narrowPeak")
-        }
-        broad <- paste0("../MACS2/",x,bam_pfx,".BAM_peaks.broadPeak") #bam_pfx
-        # first look for narrowpeak then braod peak
-        if(file.exists(narrow)) {
-            bed <- read.delim(narrow, header = FALSE)
-        } else if (file.exists(broad)) {
-            bed <- read.delim(broad, header = FALSE)
-        } else {
-            stop("MACS2 output doesn't exist. Neither ", narrow, " , nor ", broad)
-        }
-
-        bed.gr <- GRanges(seqnames = bed$V1, ranges = IRanges(start = bed$V2, end = bed$V3), name = bed$V4)
-        return(bed.gr)
-    })
-} else {
-    allpeaks = lapply(snakemake@input[['peaks']], function(x) {
-        bed = read.delim(paste0("../", x), header=FALSE)
+if (! external_bed) {
+    if(snakemake@params[['peakCaller']] == "MACS2") {
+        allpeaks <- lapply(fnames, function(x) {
+            narrow <- paste0("../MACS2/",x,bam_pfx,".BAM_peaks.narrowPeak") #bam_pfx
+            if(snakemake@params[["pipeline"]] %in% "ATAC-seq"){
+                narrow <- paste0("../MACS2/",x,".filtered.short.BAM_peaks.narrowPeak")
+            }
+            broad <- paste0("../MACS2/",x,bam_pfx,".BAM_peaks.broadPeak") #bam_pfx
+             #first look for narrowpeak then braod peak
+            if(file.exists(narrow)) {
+                bed <- read.delim(narrow, header = FALSE)
+            } else if (file.exists(broad)) {
+                bed <- read.delim(broad, header = FALSE)
+            } else {
+                stop("MACS2 output doesn't exist. Neither ", narrow, " , nor ", broad)
+            }
+
+            bed.gr <- GRanges(seqnames = bed$V1, ranges = IRanges(start = bed$V2, end = bed$V3), name = bed$V4)
+            return(bed.gr)
+        })
+    } else {
+        allpeaks = lapply(snakemake@input[['peaks']], function(x) {
+            bed = read.delim(paste0("../", x), header=FALSE)
+            bed.gr = GRanges(seqnames = bed$V1, ranges = IRanges(start = bed$V2, end = bed$V3), name = bed$V4)
+            return(bed.gr)
+        })
+    # merge
+    allpeaks <- Reduce(function(x,y) GenomicRanges::union(x,y), allpeaks)
+    }
+    } else {
+        bed = read.delim(snakemake@input[['peaks']],header=FALSE)
         bed.gr = GRanges(seqnames = bed$V1, ranges = IRanges(start = bed$V2, end = bed$V3), name = bed$V4)
-        return(bed.gr)
-    })
+        allpeaks <- bed.gr
 }
 
 # merge
-allpeaks <- Reduce(function(x,y) GenomicRanges::union(x,y), allpeaks)
+#allpeaks <- Reduce(function(x,y) GenomicRanges::union(x,y), allpeaks)
 
 ## keep only these peaks for testing DB
 message(paste0("Filtering windows using MACS2 output : ", length(allpeaks) , " regions used (Union of peaks)"))
@@ -165,4 +177,3 @@ sink("CSAW.session_info.txt")
 sessionInfo()
 sink()
 
-print("DONE..!")

snakePipes/shared/rscripts/DESeq2.R

@@ -9,28 +9,47 @@
 # args 5 : path to DE_functions
 # args 6 : T/F whether or not the workflow is allele-sepecific
 # args 7 : tx2gene file for salmon --> DESeq mode
+# args 9:  model formula
 
 .libPaths(R.home("library"))
 
-args = commandArgs(TRUE)
+#args = commandArgs(TRUE)
 
 
 ## re invented RNaseq workflow
-sampleInfoFilePath <- args[1]
-countFilePath <- args[2]
-fdr <- as.numeric(args[3])
-geneNamesFilePath <- args[4]
-importfunc <- args[5]
-allelic_info <- as.logical(args[6])
+#sampleInfoFilePath <- args[1]
+#countFilePath <- args[2]
+#fdr <- as.numeric(args[3])
+#geneNamesFilePath <- args[4]
+#importfunc <- args[5]
+#allelic_info <- as.logical(args[6])
 ## if output is from salmon then tx2gene file should be present
-tx2gene_file <- args[7]
+#tx2gene_file <- args[7]
+
+sampleInfoFilePath <- snakemake@params[["sampleSheet"]]
+countFilePath <- snakemake@params[["counts_table"]]
+fdr <- as.numeric(snakemake@params[["fdr"]])
+geneNamesFilePath <- snakemake@params[["symbol_file"]]
+importfunc <- snakemake@params[["importfunc"]]
+allelic_info <- as.logical(snakemake@params[["allele_info"]])
+tx2gene_file <- snakemake@params[["tx2gene_file"]]
+rmdTemplate <- snakemake@params[["rmdTemplate"]]
+formulaInput <- snakemake@params[["formula"]]
+wdir <- snakemake@params[["outdir"]]
+
+setwd(wdir)
+
+
 if(file.exists(tx2gene_file)) {
   tximport <- TRUE
 } else {
   tximport <- FALSE
 }
 
-rmdTemplate <- args[8]
+#rmdTemplate <- args[8]
+
+#formulaInput <- args[9]
+
 topN <- 50
 ## include functions
 suppressPackageStartupMessages(library(ggplot2))
@@ -49,6 +68,7 @@ cat(paste("Working dir:", getwd(), "\n"))
 cat(paste("Sample info CSV:", sampleInfoFilePath, "\n"))
 cat(paste("Count file:", countFilePath, "\n"))
 cat(paste("FDR:", fdr, "\n"))
+cat(paste("Custom formula:",formulaInput,"\n"))
 cat(paste("Gene names:", geneNamesFilePath, "\n"))
 cat(paste("Number of top N genes:", topN, "\n"))
 cat(paste("Salmon --> DESeq2 : ", tximport, "\n"))
@@ -96,7 +116,7 @@ if(length(unique(sampleInfo$condition))>1){
     if(tximport & allelic_info){
         message("Detected allelic Salmon counts. Skipping DESeq_basic.")
     }else{
-        seqout <- DESeq_basic(countdata, coldata = sampleInfo, fdr = fdr, alleleSpecific = allelic_info, from_salmon = tximport)
+        seqout <- DESeq_basic(countdata, coldata = sampleInfo, fdr = fdr, alleleSpecific = allelic_info, from_salmon = tximport, customFormula = formulaInput)
 
         DESeq_writeOutput(DEseqout = seqout,
                 fdr = fdr, outprefix = "DEseq_basic",