Merge pull request #1096 from maxplanck-ie/develop

* QC in the ChIPseq workflow is now performed with the ChIPQC R package
* added allelic-whatshap mode to mRNA seq
* fixes #1048
* fixes #1085
* fixes #1083
* fixes #1082
* fixes #1063
* fixes #1058
* fixes #1024
* fixes minor issues with mRNAseq allelic-whatshap mode

Author: katsikora

.ci_stuff/test_dag.sh

@@ -1,4 +1,4 @@
-{% set version = "3.1.0" %}
+{% set version = "3.2.0" %}
 
 package:
   name: snakepipes
@@ -24,7 +24,7 @@ requirements:
     - snakemake-executor-plugin-cluster-generic >=1.0.9
     - pandas
     - thefuzz
-    - pyyaml >=5.1
+    - ruamel.yaml
 
 test:
   commands:

conda-recipe/meta.yaml

@@ -1,6 +1,22 @@
 snakePipes News
 ===============
 
+snakePipes 3.2.0
+________________
+
+
+* QC in the ChIPseq workflow is now performed with the ChIPQC R package
+* added allelic-whatshap mode to mRNA seq
+* fixes #1048
+* fixes #1085
+* fixes #1083
+* fixes #1082
+* fixes #1063
+* fixes #1058
+* fixes #1024
+* fixes minor issues with mRNAseq allelic-whatshap mode
+
+
 snakePipes 3.1.0
 ________________
 

docs/content/News.rst

@@ -163,41 +163,28 @@ Understanding the outputs
 The ChIPseq pipeline will generate additional output as follows::
 
     .
+    ├── AnnotatedResults_MACS2_diffChIP_k4me3
+    ├── Annotation
+    ├── bamCoverage
+    ├── CSAW_MACS2_diffChIP_k4me3
     ├── deepTools_ChIP
     │   ├── bamCompare
-    │   │   ├── sample1.filtered.log2ratio.over_SRR6761502.bw
-    │   │   ├── sample1.filtered.subtract.SRR6761502.bw
-    │   │   ├── sample2.filtered.log2ratio.over_SRR6761502.bw
-    │   │   └── sample2.filtered.subtract.SRR6761502.bw
     │   └── plotFingerprint
-    │       ├── plotFingerprint.metrics.txt
-    │       └── plotFingerprint.png
+    ├── deepTools_qc
+    │   ├── bamPEFragmentSize
+    │   ├── estimateReadFiltering
+    │   ├── multiBamSummary
+    │   ├── plotCorrelation
+    │   ├── plotCoverage
+    │   └── plotPCA
+    ├── filtered_bam
     ├── histoneHMM
-    │   ├── sample2.filtered.histoneHMM-em-posterior.txt.gz
-    │   ├── sample2.filtered.histoneHMM-regions.gff.gz
-    │   ├── sample2.filtered.histoneHMM-regions.gff.gz.tbi
-    │   ├── sample2.filtered.histoneHMM.txt.gz
-    │   ├── sample2.filtered.histoneHMM-zinba-emfit.pdf
-    │   ├── sample2.filtered.histoneHMM-zinba-params-em.RData
-    │   └── sample2.filtered.histoneHMM-zinba-params-em.txt
-    ├── Genrich
-    │   └── sample2.narrowPeak
-    └── MACS2
-        ├── sample1.filtered.BAM_peaks.narrowPeak
-        ├── sample1.filtered.BAM_peaks.qc.txt
-        ├── sample1.filtered.BAM_peaks.xls
-        ├── sample1.filtered.BAMPE_peaks.narrowPeak
-        ├── sample1.filtered.BAMPE_peaks.xls
-        ├── sample1.filtered.BAMPE_summits.bed
-        ├── sample1.filtered.BAM_summits.bed
-        ├── sample2.filtered.BAM_peaks.broadPeak
-        ├── sample2.filtered.BAM_peaks.gappedPeak
-        ├── sample2.filtered.BAM_peaks.qc.txt
-        ├── sample2.filtered.BAM_peaks.xls
-        ├── sample2.filtered.BAMPE_peaks.broadPeak
-        ├── sample2.filtered.BAMPE_peaks.gappedPeak
-        └── sample2.filtered.BAMPE_peaks.xls
-    
+    ├── histoneHMM_chipqc
+    ├── logs
+    ├── MACS2
+    ├── MACS2_chipqc
+    ├── QC_report
+    ├── Sambamba
 
 
 Following up on the DNAmapping module results (see :doc:`DNAmapping`), the workflow produces the following output directories :
@@ -206,12 +193,16 @@ Following up on the DNAmapping module results (see :doc:`DNAmapping`), the workf
 
 * **Genrich**: This folder contains the output of `Genrich <https://github.com/jsh58/Genrich>`__. This will only exist IF you specified ``--peakCaller Genrich`` and you have samples with non-broad peaks. The output is in narrowPeak format, like that from MACS2.
 
-* **MACS2**: This folder contains the output of `MACS2 <https://github.com/taoliu/MACS>`__ on the ChIP samples, MACS2 would perform either a **narrow** or **broad** peak calling on the samples, as indicated by the ChIP sample configuration file (see :ref:`ChIPconfig`). The outputs files would contain the respective tags (**narrowPeak** or **broadPeak**). This folder will only exist if you have non-broad marks and use MACS2 for peak calling
+* **MACS2**: This folder contains the output of `MACS2 <https://github.com/taoliu/MACS>`__ on the ChIP samples, MACS2 would perform either a **narrow** or **broad** peak calling on the samples, as indicated by the ChIP sample configuration file (see :ref:`ChIPconfig`). The outputs files would contain the respective tags (**narrowPeak** or **broadPeak**). This folder will only exist if you have non-broad marks and use MACS2 for peak calling (default).
+
+* **MACS2_chipqc**: This folder contains the output of `ChIPQC <https://bioconductor.org/packages/release/bioc/vignettes/ChIPQC/inst/doc/ChIPQC.pdf>`__ analysis of the peaks called by MACS2. If you used a different peak caller, a chipqc output folder with the peak caller in its name will be listed.
 
 * **histoneHMM**: This folder contains the output of `histoneHMM <https://github.com/matthiasheinig/histoneHMM>`__. This folder will only exist if you have broad marks.
 
-* **CSAW_sampleSheet**: This folder is created optionally, if you provide a sample sheet for differential binding analysis. (see :ref:`diffBinding`) CSAW will be run using peaks called by the chosen peak caller, and the output folder will be named accordingly.
-* **AnnotatedResults_sampleSheet**: This folder is created optionally, if you provide a sample sheet for differential binding analysis. (see :ref:`diffBinding`). Differentially bound regions annotated with distance to nearest gene are stored here.
+* **histoneHMM_chipqc**: This folder contains the output of `ChIPQC <https://bioconductor.org/packages/release/bioc/vignettes/ChIPQC/inst/doc/ChIPQC.pdf>`__ analysis of the peaks called by histoneHMM. This folder will only exist if you have broad marks.
+
+* **CSAW_peakCaller_sampleSheet**: This folder is created optionally, if you provide a sample sheet for differential binding analysis. (see :ref:`diffBinding`) CSAW will be run using peaks called by the chosen peak caller, and the output folder will be named accordingly.
+* **AnnotatedResults_peakCaller_sampleSheet**: This folder is created optionally, if you provide a sample sheet for differential binding analysis. (see :ref:`diffBinding`). Differentially bound regions annotated with distance to nearest gene are stored here.
 
 .. note:: Although in case of broad marks, we also perform the MACS2 `broadpeak` analysis (output available as ``MACS2/<sample>.filtered.BAM_peaks.broadPeak``), we would recommend using the histoneHMM outputs in these cases, since histoneHMM produces better results than MACS2 for broad peaks.
 

docs/content/workflows/ChIPseq.rst

@@ -9,7 +9,7 @@ What it does
 The snakePipes mRNAseq workflow allows users to process their single or paired-end
 mRNAseq fastq files upto the point of gene/transcript-counts and differential expression.
 It also allows full allele-specific mRNAseq analysis (up to allele-specific
-differential expression) using the *allelic-mapping* mode.
+differential expression) using the *allelic-mapping* or *allelic-whatshap* mode.
 
 .. image:: ../images/RNAseq_pipeline.png
 
@@ -94,6 +94,8 @@ There is a configuration file in ``snakePipes/workflows/mRNAseq/defaults.yaml``:
     # for allele-spcific mapping
     SNPFile:
     NMaskedIndex:
+    # for allelic-whatshap
+    pvcf:
     #### Flag to control the pipeline entry point
     fromBAM: False
     bamExt: '.bam'
@@ -202,6 +204,15 @@ Allele-specific, gene-level differential expression analysis is then performed u
 
 **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. 
 
+
+"allelic-whatshap"
+~~~~~~~~~~~~~~~~~~
+
+**allelic-whatshap** mode applies a standard alignment to a nonmasked genome with STAR, followed by allele-specific splitting
+of mapped files with whatshap, requiring a phased vcf file as input ( ``--phased-vcf`` ). The file must be bzip-compressed and tabix-indexed as well. Gene-level quantification is performed for each allele using **featureCounts**.
+Allele-specific, gene-level differential expression analysis is then performed using **DESeq2**.
+
+
 "alignment-free"
 ~~~~~~~~~~~~~~~~
 

docs/content/workflows/mRNAseq.rst

@@ -6,7 +6,7 @@ build-backend = "setuptools.build_meta"
 name = "snakePipes"
 description = 'Snakemake workflows and wrappers for NGS data processing from the MPI-IE'
 readme = "README.md"
-version = "3.1.0"
+version = "3.2.0"
 keywords = [
     "DNAmapping",
     "ChIPSeq",
@@ -35,7 +35,8 @@ dependencies = [
     "snakemake >= 8",
     "pandas",
     "thefuzz",
-    "pyyaml >= 5.1",
+#    "pyyaml >= 5.1",
+    "ruamel.yaml",
     "snakemake-executor-plugin-cluster-generic >= 1.0.9",
     "graphviz"
 ]

pyproject.toml

@@ -5,7 +5,8 @@
 import subprocess
 import os
 import re
-import yaml
+#import yaml
+from ruamel.yaml import YAML
 import glob
 import sys
 import shutil
@@ -43,7 +44,10 @@ def set_env_yamls():
             'CONDA_SAMBAMBA_ENV': 'envs/sambamba.yaml',
             'CONDA_pysam_ENV': 'envs/pysam.yaml',
             'CONDA_SEACR_ENV': 'envs/chip_seacr.yaml',
-            'CONDA_FQLINT_ENV': 'envs/fqlint.yaml'
+            'CONDA_FQLINT_ENV': 'envs/fqlint.yaml',
+            'CONDA_WHATSHAP_ENV': 'envs/whatshap.yaml',
+            'CONDA_PICARD_ENV': 'envs/picard.yaml',
+            'CONDA_CHIPQC_ENV': 'envs/chipqc.yaml'
             }
 
 
@@ -86,8 +90,10 @@ def namesOKinR(sampleNames):
 
 
 def load_configfile(configFiles, verbose, info='Config'):
+    yaml=YAML(typ='safe')
     with open(configFiles, "r") as f:
-        config = yaml.load(f, Loader=yaml.FullLoader)
+        #config = yaml.load(f, Loader=yaml.FullLoader)
+        config = yaml.load(f)
 
     config = sanity_dict_clean(config)
 
@@ -100,9 +106,15 @@ def load_configfile(configFiles, verbose, info='Config'):
     return config
 
 
-def write_configfile(configFile, config):
+def write_configfile(configFile, config, trafo):
+    yaml=YAML(typ='safe')
+    yaml.default_flow_style = False
     with open(configFile, 'w') as f:
-        yaml.dump(config, f, default_flow_style=False)
+        #yaml.dump(config, f, default_flow_style=False)
+        if trafo:
+            yaml.dump(config, f, transform=trafo)
+        else:
+            yaml.dump(config, f)
 
 
 # returns all key-value pairs that are different from dict1 to dict2
@@ -630,7 +642,7 @@ def commonYAMLandLogs(baseDir, workflowDir, defaults, args, callingScript):
     # save to configs.yaml in outdir
     config = defaults
     config.update(vars(args))  # This allows modifications of args after handling a user config file to still make it to the YAML given to snakemake!
-    write_configfile(os.path.join(args.outdir, '{}.config.yaml'.format(workflowName)), config)
+    write_configfile(os.path.join(args.outdir, '{}.config.yaml'.format(workflowName)), config, trafo=None)
 
     # merge cluster config files: 1) global one, 2) workflow specific one, 3) user provided one
     cfg = load_configfile(os.path.join(baseDir, "shared", "defaults.yaml"), False, "defaults")
@@ -717,7 +729,7 @@ def print_DAG(args, snakemake_cmd, callingScript, defaults):
         config['verbose'] = False
         write_configfile(
             os.path.join(args.outdir,
-                         '{}.config.yaml'.format(workflowName)), config)
+                         '{}.config.yaml'.format(workflowName)), config, trafo=None)
 
         DAGproc = subprocess.Popen(
             snakemake_cmd + " --rulegraph -q ",
@@ -732,7 +744,7 @@ def print_DAG(args, snakemake_cmd, callingScript, defaults):
         config['verbose'] = oldVerbose
         write_configfile(
             os.path.join(args.outdir, '{}.config.yaml'.format(workflowName)),
-            config)
+            config, trafo=None)
 
 
 def logAndExport(args, workflowName):
@@ -772,6 +784,8 @@ def runAndCleanup(args, cmd, logfile_name):
     for _l in p.stdout:
         sys.stdout.write(_l.strip() + '\n')
         f.write(_l.strip() + '\n')
+        sys.stdout.flush()
+        f.flush()
     p.wait()
 
     # Exit with an error if snakemake encountered an error
@@ -792,6 +806,9 @@ def runAndCleanup(args, cmd, logfile_name):
     if args.emailAddress:
        sendEmail(args, 0)
 
+def tr(s):
+    return s.replace('null', 'None')
+
 
 def predict_chip_dict(wdir, input_pattern_str, bamExt, fromBAM=None):
     """
@@ -854,14 +871,14 @@ def predict_chip_dict(wdir, input_pattern_str, bamExt, fromBAM=None):
             print("No control sample found!")
 
         chip_dict_pred["chip_dict"][i] = {}
-        chip_dict_pred["chip_dict"][i]['control'] = tmp
+        chip_dict_pred["chip_dict"][i]['Control'] = tmp if tmp !=  "" else None
         if re.match(".*(H3K4me1|H3K36me3|H3K9me3|H3K27me3).*", i, re.IGNORECASE):
-            chip_dict_pred["chip_dict"][i]['broad'] = True
+            chip_dict_pred["chip_dict"][i]['Broad'] = True
         else:
-            chip_dict_pred["chip_dict"][i]['broad'] = False
+            chip_dict_pred["chip_dict"][i]['Broad'] = False
 
     outfile = os.path.join(wdir, "chip_seq_sample_config.PREDICTED.yaml")
-    write_configfile(outfile, chip_dict_pred)
+    write_configfile(outfile, chip_dict_pred,trafo=tr)
     print("---------------------------------------------------------------------------------------")
     print("ChIPseq sample configuration is written to file ", outfile)
     print("Please check and modify this file - this is just a guess! Then run the workflow with it.")

snakePipes/common_functions.py

@@ -3,7 +3,7 @@
 # Note that due to limitations in yaml.dump, only very basic structures are
 # permitted here.
 ################################################################################
-snakemakeOptions: ''
+snakemakeOptions: ' --keep-going '
 organismsDir: 'shared/organisms'
 snakemakeProfile: 'shared/profiles/local'
 tempDir: /scratch/local

snakePipes/shared/defaults.yaml

@@ -91,6 +91,5 @@ ignoreForNormalization: chrX chrY chrM GL000008.2 GL000009.2 GL000194.1 GL000195
   KI270748.1 KI270749.1 KI270750.1 KI270751.1 KI270752.1 KI270753.1 KI270754.1 KI270755.1
   KI270756.1 KI270757.1
 known_splicesites: /data/repository/organisms/GRCh38_gencode_40/gencode/release-40/HISAT2/genome.ss
-rmsk_file: ''
 star_index: /data/repository/organisms/GRCh38_gencode_40/Indices/STAR_2.7.10
 rmsk_file: /data/repository/organisms/GRCh38_gencode_40/repeatMasker/genome.fa.tbl

snakePipes/shared/organisms/GRCh38_gencode40.yaml

@@ -60,6 +60,8 @@ set-resources:
     mem: 2G
   CollectInsertSizeMetrics:
     mem: 1G
+  conversionRate:
+    mem: 4G
   correct_matrix:
     mem: 7G
   CpG_report:
@@ -205,4 +207,4 @@ set-resources:
   velo_to_sce:
     mem: 30G
   velocyto:
-    mem: 20G
+    mem: 20G

snakePipes/shared/profiles/local/config.yaml

@@ -6,7 +6,7 @@ cluster-generic-submit-cmd:
     --cpus-per-task={threads} \
     --mem-per-cpu={resources.mem} \
     --job-name=snakePipes_{rule}-{wildcards} \
-    --output=logs/{rule}/{rule}-{wildcards}-%j.out \
+    --output=logs/{rule}/{rule}-{wildcards}-%j-%N.out \
     --time={resources.time} \
     --parsable"
 jobs: 20
@@ -72,6 +72,8 @@ set-resources:
     mem: 2G
   CollectInsertSizeMetrics:
     mem: 1G
+  conversionRate:
+    mem: 4G
   correct_matrix:
     mem: 7G
   CpG_report:

snakePipes/shared/profiles/snakepipes_genericprofile/config.yaml

@@ -118,9 +118,13 @@ if (! external_bed) {
         })
     } 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)
+            peakfile<-paste0("../", x)
+            if(file.exists(peakfile) & file.info(peakfile)$size > 0){
+                bed = read.delim(peakfile, header=FALSE)
+                bed.gr = GRanges(seqnames = bed$V1, ranges = IRanges(start = bed$V2, end = bed$V3), name = bed$V4)
+                }else{message(paste0("Skipping peakfile ",peakfile))
+                      bed.gr=GRanges(c(seqnames=NULL,ranges=NULL,strand=NULL,name=NULL))}
+                return(bed.gr)
         })
     }
     # merge

snakePipes/shared/rscripts/CSAW.R

@@ -343,6 +343,8 @@ writeOutput_chip <- function(chipResultObject, outfile_prefix, fdrcutoff,lfccuto
     full_res[,2]<-full_res[,2]-1
     full_res[,2]<-format(full_res[,2], scientific = FALSE,trim=TRUE)
     full_res[,3]<-format(full_res[,3], scientific = FALSE,trim=TRUE)
+    #write full results
+    write.table(full_res,file="Full.results.bed",row.names=FALSE,col.names=FALSE,sep="\t",quote=FALSE)
     ##filter full result for FDR and LFC and write to output
     full_res.filt<-subset(full_res,(FDR<=fdrcutoff)&(abs(best.logFC)>=lfccutoff))
     if(nrow(full_res.filt)>0){