diff --git a/example/templates/hg19/README b/example/templates/dev/hg19/README
similarity index 100%
rename from example/templates/hg19/README
rename to example/templates/dev/hg19/README
diff --git a/example/templates/dev/hg19/config.yaml b/example/templates/dev/hg19/config.yaml
new file mode 100644
index 0000000000000000000000000000000000000000..6ac31ea5e381f1a13a695958b4b15f719a1f68c0
--- /dev/null
+++ b/example/templates/dev/hg19/config.yaml
@@ -0,0 +1,239 @@
+################################################
+################################################
+# CONFIGURATION FILE FOR THE ATAC-SEQ PIPELINE #
+################################################
+################################################
+
+# This format allows comments and is therefore easier to work with than the json format we used before.
+# Quotation marks are optional for strings. Make sure to put ": " (that is, hyphen space) as separator
+
+
+##################
+# SECTION output #
+##################
+output:
+
+  # STRING. Absolute path to the output directory. Will be created if not yet present.
+  outdir: "/FULL/PATH/TO/OUTPUT/SAME/AS/IN/CLUSTER.json"
+
+  # BOOLEAN. “true” or “false”. Default “true”. Should the pipeline use MACS2 to call peaks? If yes, peaks will be called in 3 different flavors (ENCODE, stringent, non-stringent). See the section “par_peakCalling” in the configuration file for details.
+  doPeakCalling: true
+
+  # BOOLEAN. "true" or "false". Default "false". Do peak calling also, in addition to the "stringent" and "non-stringent" flavour, according to ENCODE guidelines (as of 2018, we so far rarely used them)
+  encodePeaks: false
+
+  # BOOLEAN. “true” or “false”. Default “false”. Should, in addition to treating all input files separately, the pipeline merge all replicate files and do all downstream analysis in addition? If set to true, for each individual as specified in the sample table, all samples belonging to this individual ID will be merged to one file after the post-processing. This makes only sense if you have more than 1 replicate per individual.
+  alsoMergeReplicates: true
+
+  # BOOLEAN. "true" or "false". Default "false". Should peaks be annotated to find out which genomic regions they overlap with? Still being tested, set to false if it causes errors
+  annotatePeaks: true
+
+  # BOOLEAN. “true” or “false”. Default “false”. Should GC bias be assessed and corrected for the final BAM files in addition to the non-corrected ones? If set to true, the GC bias will be assessed and corrected using deepTools. Additionally, all downstream steps of the pipeline will be done for both GC-corrected and original files (including peak calling, PCA, coverage, etc). This greatly increases running time of th pipeline and essentially doubles the number of files. We recommend setting this to false unless the GC bias is important for downstream analyses
+  # !!! Important: Leave at false for now, I have to fix an issue in the Snakefile that currently causes an error !!
+  correctGCBias: false
+
+  # BOOLEAN. “true” or “false”. Default “false”. Should base qualities be recalibrated using GATK to detect and correct systematic errors in base quality scores? This step can be time-consuming and needs the following other parameters: GATK_jar, knownSNPs, knownIndels. I recommend turning this off for now.
+  doBaseRecalibration: false
+
+  # BOOLEAN. “true” or “false”. Default “true”. Should the pipeline produce coverage files and diagnostic plots? If set to true, coverage files for the final BAM files in bigwig and bedgraph.gz format will be produced as well as a coverage plot using deepTools.
+  generateCoverageFiles: true
+
+  # BOOLEAN. “true” or “false”. Default “true”. Not implemented yet.
+  doIDR: false
+
+
+###################
+# SECTION general #
+###################
+general:
+
+  # INTEGER. Default 12. Maximum number of cores per rule. For local computation, the minimum of this value and the --cores parameter will define the number of CPUs per rule, while in a cluster setting, the minimum of this value and the number of cores on the node the jobs runs is used.
+  maxCoresPerRule: 12
+
+###################
+# SECTION samples #
+###################
+samples:
+
+  # STRING. No default. Absolute path to the sample summary file. See section 2.4 for details.
+  summaryFile: "/FULL/PATH/TO/samples.csv"
+
+  # BOOLEAN. “true” or “false”. Default “true”. Paired-end data? Single-end ATAC-Seq data is not yet supported with this pipeline. If set to "false", the Snakemake pipeline will abort in the beginning.
+  pairedEnd: true
+
+  # STRING. "ATACseq" or "ChIPseq". Only data type specific steps are executed. Currently, if set to ChIP-Seq, ATAC-seq specific steps like RSS adjustment are not executed.
+  dataType: "ATACseq"
+
+###########################
+# SECTION additionalInput #
+###########################
+additionalInput:
+
+  # STRING. Absolute path to the scripts directory. Either "/g/scb2/zaugg/carnold/Projects/AtacSeq/src/Snakemake/R" or "/g/scb2/zaugg/carnold/Projects/AtacSeq/src/Snakemake/dev/R"
+  scriptsDir: "/g/scb2/zaugg/carnold/Projects/AtacSeq/src/Snakemake/dev/R"
+
+  # STRING. Absolute path to the adapters file for Trimmomatic in fasta format. Default “/g/scb2/zaugg/zaugg_shared/Programs/Trimmomatic-0.33/adapters/NexteraPE-PE.fa”. There is usually no need to change this unless for your experiment, this adapter file is not suited.
+  trimmomatic_adapters: "/g/scb2/zaugg/zaugg_shared/Programs/Trimmomatic-0.33/adapters/NexteraPE-PE.fa"
+
+  # STRING. Absolute path to a BED file that contains the genomic regions that should be filtered from the peaks. The default depends on the genome assembly, see the templates for details. Only needed if doPeakCalling is set to true.
+  blacklistRegions: "/g/scb2/zaugg/zaugg_shared/annotations/hg19/blacklisted/hg19-blacklist.v2.bed"
+
+  # STRING. Absolute path to a database of known polymorphic sites (SNPs and Indels, respectively). This is only needed if (1) doBaseRecalibration is set to true and (2) the genome is either hg19 or hg38 and ignored otherwise. The default depends on the genome assembly, see the templates for details. Supported formats from GATK: BCF2, BEAGLE, BED, BEDTABLE, EXAMPLEBINARY, GELITEXT, RAWHAPMAP, REFSEQ, SAMPILEUP, SAMREAD, TABLE, VCF, VCF3.
+  knownSNPs: "/g/scb2/zaugg/zaugg_shared/annotations/hg19/GATK_bundle/dbsnp_138.hg19.vcf.gz"
+  knownIndels: "/g/scb2/zaugg/zaugg_shared/annotations/hg19/GATK_bundle/Mills_and_1000G_gold_standard.indels.hg19.sites.vcf.gz"
+
+  # STRING. The default depends on the genome assembly, see the templates for details. Absolute path to the reference genome in fasta and 2bit format, respectively, both of which have to correspond to the same genome assembly version as used for the alignment as well as the database of polymorphic sites (knownSNPs and knownIndels, if applicable).
+  refGenome_fasta: "/g/scb2/zaugg/zaugg_shared/annotations/hg19/GATK_bundle/ucsc.hg19.onlyRefChr.fasta"
+  refGenome_2bit: "/g/scb2/zaugg/zaugg_shared/annotations/hg19/GATK_bundle/ucsc.hg19.2bit"
+
+  # STRING. The default depends on the genome assembly, see the templates for details. Absolute path to an genome annotation file in GTF format.
+  annotationGTF: "/g/scb2/zaugg/zaugg_shared/annotations/hg19/Gencode_v19/gencode.v19.annotation.gtf"
+
+
+#######################
+# SECTION executables #
+#######################
+
+# If run with Singularity, this section can be ignored
+
+executables:
+
+  # STRING. Default “java”. Name of the executable for Java. Java version must be at least 1.8! Only needed if doBaseRecalibration is set to true
+  java: "java"
+
+  # STRING. Default "/g/scb2/zaugg/carnold/Projects/AtacSeq/src/Snakemake/tools/GenomeAnalysisTK.jar". Absolute path to a JAR file for the GATK suite. Only needed if doBaseRecalibration is set to true
+  GATK_jar: "/g/scb2/zaugg/carnold/Projects/AtacSeq/src/Snakemake/tools/GenomeAnalysisTK.jar"
+
+  # STRING. Default "/g/scb2/zaugg/zaugg_shared/Programs/Picardtools/picardOld.jar". Absolute path to a JAR file for the Picard suite.
+  PICARD_jar: "/g/scb2/zaugg/carnold/Projects/AtacSeq/src/Snakemake/tools/picardOld.jar"
+
+########################
+# SECTION par_trimming #
+########################
+par_trimming:
+
+  # STRING. Default "1:30:4:5:true". ILLUMINACLIP value. See trimmomatic manual
+  trimmomatic_ILLUMINACLIP: "1:30:4:1:true"
+
+  # INTEGER. Default 3. TRAILING value. See trimmomatic manual
+  trimmomatic_trailing: 3
+
+  # INTEGER. Default 20. MINLEN value. See trimmomatic manual
+  trimmomatic_minlen: 20
+
+  # STRING. Default "phred33". Phred type. See trimmomatic manual. The "-" is added automatically by the Snakemake pipeline.
+  #trimmomatic_phredType: "phred33"
+
+#####################
+# SECTION par_align #
+#####################
+par_align:
+
+  # STRING. Default "--very-sensitive". Sensitivity. Leave empty for the default sensitivity. See bowtie2 manual.
+  bowtie2_sensitivity: "--very-sensitive"
+
+  # INTEGER. Default 2000. Value for parameter X. See bowtie2 manual.
+  bowtie2_maxFragmentLength: 2000
+
+  # STRING. Default "/g/scb/zaugg/zaugg_shared/annotations/hg19/referenceGenome/Bowtie2/hg19". Absolute path to the reference genome + the prefix of the file with the index. In the default case, we use those files within /g/scb/zaugg/zaugg_shared/annotations/hg19/referenceGenome/Bowtie2 that start with hg19. See bowtie2 manual for more details.
+  bowtie2_refGenome: "/g/scb2/zaugg/zaugg_shared/annotations/hg19/referenceGenome/Bowtie2/hg19"
+
+  # STRING. Default “hg19”. Reference genome assembly version. Must match the one used by the alignment program.
+  assemblyVersion: "hg19"
+
+
+#########################
+# SECTION par_postalign #
+#########################
+par_postalign:
+
+  # INTEGER. Default 10. Minimum MAPQ score. Reads with a lower MAPQ quality will be removed during the post-alignment processing.
+  minMAPQscore: 10
+
+  # STRING. Default “LENIENT”. Value of the VALIDATION_STRINGENCY from SortSam (Picard tools). See the manual for details.
+  ValidationStringencySortSam: "LENIENT"
+
+  # STRING. Default “SILENT”. Value of the VALIDATION_STRINGENCY from MarkDuplicates (Picard tools). See the manual for details.
+  ValidationStringencyMarkDuplicates: "SILENT"
+
+  # STRING. Defauled “ID”. Used for filtering reads. Relates to the one letter abbreviations for CIGAR strings such as I for insertion and D for deletion.  Specify all the one letter abbreviations in the CIGAR string of a read here that should be filtered. "ID" keeps a read only if the CIGAR string does not contain the letters "I" and "D" (e.g., only M for example)
+  CIGAR: "ID"
+
+  # INTEGER. Default 4 and -5, respectively. Adjustment of the read start positions on the forward and reverse strand and should be a positive or negative number, respectively. See the Buenrostro paper for details. FOR DATA TYPE ATAC-SEQ ONLY
+  adjustRSS_forward: 4
+  adjustRSS_reverse: 5
+
+#######################
+# SECTION par_scripts #
+#######################
+
+# NOTE: There is usually no need to modify the parameters here.
+
+par_scripts:
+
+  # INTEGER. Default 4000. The region size in bp that specifies what is considered within a TSS. The STATS script does a TSS enrichment test to test whether or not ATAC-Seq reads are primarily located within annotated TSS as opposed to outside of TSS regions. A value of 4000 means the region from -2kb up to +2kb of annotated TSS.
+  STATS_script_withinThr: 4000
+
+  # INTEGER. Default 1000. The size of the region adjacent to the within TSS region that is considered outside of a TSS. A value of 1000 therefore denotes the 1kb region up- and downstream of the within TSS region (from -3 to -2kb upstream and from +2 to +3 kb downstream of annotated TSS.)
+  STATS_script_outsideThr: 1000
+
+  # STRING. Default "protein_coding". Gene type to keep / do the analayses for. Allowed are gene types as specified by GENCODE. The default is "protein_coding"
+  STATS_script_geneTypesToKeep: "protein_coding"
+
+  # INTEGER. Default 600. Fragment length cutoff. All reads with a fragment length less than this value will be filtered for the purpose of this script.
+  FL_distr_script_cutoff: 600
+
+
+
+###########################
+# SECTION par_peakCalling #
+###########################
+par_peakCalling:
+
+  # STRING. Default "--nolambda –nomodel". Peak calling model for non-stringent peak calling.
+  modelNonStringent: "--nolambda --nomodel"
+
+  # STRING. Default "—nomodel”. Peak calling model for stringent peak calling.
+  modelStringent: "--nomodel"
+
+  # NUMERIC [0, 1]. Default 0.01. Minimum q-value threshold for stringent peak calling.
+  modelStringent_minQValue: 0.01
+
+  # NUMERIC [0, 1]. Default 0.1. Minimum q-value threshold for non-stringent peak calling.
+  modelNonStringent_minQValue: 0.1
+
+  # INTEGER. Default 10000. Value for slocal parameter for non-stringent peak calling.
+  modelNonStringent_slocal: 10000
+
+  # the enxt 3 options are only relevant if encodePeaks=true. ignored otherwise
+
+  # NUMERIC [0, 1]. Default 0.1. p-value threshold for ENCODE peak calling.
+  Encode_pValThreshold: 0.1
+
+  # STRING. Default "--nomodel --shift -75 --extsize 150 --broad --keep-dup all". Model for Encode peak calling (broad and gapped mode)
+  Encode_modelBroadAndGapped: "--nomodel --shift -75 --extsize 150 --broad --keep-dup all"
+
+  # STRING. Default "--nomodel --shift -75 --extsize 150 -B --SPMR --keep-dup all –call-summits". Model for Encode peak calling (narrow mode)
+  Encode_modelNarrow: "--nomodel --shift -75 --extsize 150 -B --SPMR --keep-dup all --call-summits"
+
+
+#########################
+# SECTION par_deepTools #
+#########################
+par_deepTools:
+
+  # INTEGER. The default depends on the genome assembly, see the templates for details. Length of the “mappable” genome in bp as defined by deepTools (see https://deeptools.readthedocs.io/en/develop/content/feature/effectiveGenomeSize.html).
+  #effectiveGenomeSize: 2750000000
+  readLength: ""
+
+  # STRING. Default “normalizeTo1x”. Either "normalizeTo1x NUMBER" OR "normalizeUsingRPKM" (note the missing (!) leading "--"), where NUMBER denotes the number of base pairs that are mappable. Beware of the mapping dependence on the read length for this number: The reported numbers on the websites are for 30bp reads, and we now have much longer reads usually. See Koehler et al. (2011) for numbers1.
+  bamCoverage_normalizationCoverage: "RPGC"
+
+  # INTEGER. Default 10. Size of the bins, in bases
+  bamCoverage_binSize: 10
+
+  # BOOLEAN. Default true. If the data is pair ended a read will by default be extended to match its pair.
+  bamCoverage_extendReads: true
+
+  # STRING. Default "--extendReads  --centerReads". Additional options that are supported by bamCoverage. Note that the "--" or "-" has to be present here.
+  bamCoverage_otherOptions: "--centerReads"
diff --git a/example/templates/hg38/README b/example/templates/dev/hg38/README
similarity index 100%
rename from example/templates/hg38/README
rename to example/templates/dev/hg38/README
diff --git a/example/templates/dev/hg38/config.yaml b/example/templates/dev/hg38/config.yaml
new file mode 100644
index 0000000000000000000000000000000000000000..a25eaf5864bcbd9a85e5091641d4ef022c1dff63
--- /dev/null
+++ b/example/templates/dev/hg38/config.yaml
@@ -0,0 +1,240 @@
+################################################
+################################################
+# CONFIGURATION FILE FOR THE ATAC-SEQ PIPELINE #
+################################################
+################################################
+
+# This format allows comments and is therefore easier to work with than the json format we used before.
+# Quotation marks are optional for strings. Make sure to put ": " (that is, hyphen space) as separator
+
+
+##################
+# SECTION output #
+##################
+output:
+
+  # STRING. Absolute path to the output directory. Will be created if not yet present.
+  outdir: "/FULL/PATH/TO/OUTPUT/SAME/AS/IN/CLUSTER.json"
+
+  # BOOLEAN. “true” or “false”. Default “true”. Should the pipeline use MACS2 to call peaks? If yes, peaks will be called in 3 different flavors (ENCODE, stringent, non-stringent). See the section “par_peakCalling” in the configuration file for details.
+  doPeakCalling: true
+
+  # BOOLEAN. "true" or "false". Default "false". Do peak calling also, in addition to the "stringent" and "non-stringent" flavour, according to ENCODE guidelines (as of 2018, we so far rarely used them)
+  encodePeaks: false
+
+  # BOOLEAN. “true” or “false”. Default “false”. Should, in addition to treating all input files separately, the pipeline merge all replicate files and do all downstream analysis in addition? If set to true, for each individual as specified in the sample table, all samples belonging to this individual ID will be merged to one file after the post-processing. This makes only sense if you have more than 1 replicate per individual.
+  alsoMergeReplicates: true
+
+  # BOOLEAN. "true" or "false". Default "false". Should peaks be annotated to find out which genomic regions they overlap with? Still being tested, set to false if it causes errors
+  annotatePeaks: true
+
+  # BOOLEAN. “true” or “false”. Default “false”. Should GC bias be assessed and corrected for the final BAM files in addition to the non-corrected ones? If set to true, the GC bias will be assessed and corrected using deepTools. Additionally, all downstream steps of the pipeline will be done for both GC-corrected and original files (including peak calling, PCA, coverage, etc). This greatly increases running time of th pipeline and essentially doubles the number of files. We recommend setting this to false unless the GC bias is important for downstream analyses
+  # !!! Important: Leave at false for now, I have to fix an issue in the Snakefile that currently causes an error !!
+  correctGCBias: false
+
+  # BOOLEAN. “true” or “false”. Default “false”. Should base qualities be recalibrated using GATK to detect and correct systematic errors in base quality scores? This step can be time-consuming and needs the following other parameters: GATK_jar, knownSNPs, knownIndels. I recommend turning this off for now.
+  doBaseRecalibration: false
+
+  # BOOLEAN. “true” or “false”. Default “true”. Should the pipeline produce coverage files and diagnostic plots? If set to true, coverage files for the final BAM files in bigwig and bedgraph.gz format will be produced as well as a coverage plot using deepTools.
+  generateCoverageFiles: true
+
+  # BOOLEAN. “true” or “false”. Default “true”. Not implemented yet.
+  doIDR: false
+
+
+###################
+# SECTION general #
+###################
+general:
+
+  # INTEGER. Default 12. Maximum number of cores per rule. For local computation, the minimum of this value and the --cores parameter will define the number of CPUs per rule, while in a cluster setting, the minimum of this value and the number of cores on the node the jobs runs is used.
+  maxCoresPerRule: 12
+
+###################
+# SECTION samples #
+###################
+samples:
+
+  # STRING. No default. Absolute path to the sample summary file. See section 2.4 for details.
+  summaryFile: "/FULL/PATH/TO/samples.csv"
+
+  # BOOLEAN. “true” or “false”. Default “true”. Paired-end data? Single-end ATAC-Seq data is not yet supported with this pipeline. If set to "false", the Snakemake pipeline will abort in the beginning.
+  pairedEnd: true
+
+  # STRING. "ATACseq" or "ChIPseq". Only data type specific steps are executed. Currently, if set to ChIP-Seq, ATAC-seq specific steps like RSS adjustment are not executed.
+  dataType: "ATACseq"
+
+###########################
+# SECTION additionalInput #
+###########################
+additionalInput:
+
+  # STRING. Absolute path to the scripts directory. Either "/g/scb2/zaugg/carnold/Projects/AtacSeq/src/Snakemake/R" or "/g/scb2/zaugg/carnold/Projects/AtacSeq/src/Snakemake/dev/R"
+  scriptsDir: "/g/scb2/zaugg/carnold/Projects/AtacSeq/src/Snakemake/dev/R"
+
+  # STRING. Absolute path to the adapters file for Trimmomatic in fasta format. Default “/g/scb2/zaugg/zaugg_shared/Programs/Trimmomatic-0.33/adapters/NexteraPE-PE.fa”. There is usually no need to change this unless for your experiment, this adapter file is not suited.
+  trimmomatic_adapters: "/g/scb2/zaugg/zaugg_shared/Programs/Trimmomatic-0.33/adapters/NexteraPE-PE.fa"
+
+  # STRING. Absolute path to a BED file that contains the genomic regions that should be filtered from the peaks. The default depends on the genome assembly, see the templates for details. Only needed if doPeakCalling is set to true.
+  blacklistRegions: "/g/scb2/zaugg/zaugg_shared/annotations/hg38/blacklisted/hg38-blacklist.v2.bed"
+
+  # STRING. Absolute path to a database of known polymorphic sites (SNPs and Indels, respectively). This is only needed if (1) doBaseRecalibration is set to true and (2) the genome is either hg19 or hg38 and ignored otherwise. The default depends on the genome assembly, see the templates for details. Supported formats from GATK: BCF2, BEAGLE, BED, BEDTABLE, EXAMPLEBINARY, GELITEXT, RAWHAPMAP, REFSEQ, SAMPILEUP, SAMREAD, TABLE, VCF, VCF3.
+  knownSNPs: "/g/scb2/zaugg/zaugg_shared/annotations/hg38/GATK/Homo_sapiens_assembly38.dbsnp138.vcf.gz"
+  knownIndels: "/g/scb2/zaugg/zaugg_shared/annotations/hg38/GATK/Mills_and_1000G_gold_standard.indels.hg38.vcf.gz"
+
+  # STRING. The default depends on the genome assembly, see the templates for details. Absolute path to the reference genome in fasta and 2bit format, respectively, both of which have to correspond to the same genome assembly version as used for the alignment as well as the database of polymorphic sites (knownSNPs and knownIndels, if applicable).
+  refGenome_fasta: "/g/scb2/zaugg/zaugg_shared/annotations/hg38/GATK/Homo_sapiens_assembly38.fasta"
+  refGenome_2bit: "/g/scb2/zaugg/zaugg_shared/annotations/hg38/GATK/Homo_sapiens_assembly38.2bit"
+
+  # STRING. The default depends on the genome assembly, see the templates for details. Absolute path to an genome annotation file in GTF format.
+  annotationGTF: "/g/scb2/zaugg/zaugg_shared/annotations/hg38/Gencode_v34/gencode.v34.annotation.gtf"
+
+
+
+#######################
+# SECTION executables #
+#######################
+
+# If run with Singularity, this section can be ignored
+
+executables:
+
+  # STRING. Default “java”. Name of the executable for Java. Java version must be at least 1.8! Only needed if doBaseRecalibration is set to true
+  java: "java"
+
+  # STRING. Default "/g/scb2/zaugg/carnold/Projects/AtacSeq/src/Snakemake/tools/GenomeAnalysisTK.jar". Absolute path to a JAR file for the GATK suite. Only needed if doBaseRecalibration is set to true
+  GATK_jar: "/g/scb2/zaugg/carnold/Projects/AtacSeq/src/Snakemake/tools/GenomeAnalysisTK.jar"
+
+  # STRING. Default "/g/scb2/zaugg/zaugg_shared/Programs/Picardtools/picardOld.jar". Absolute path to a JAR file for the Picard suite.
+  PICARD_jar: "/g/scb2/zaugg/carnold/Projects/AtacSeq/src/Snakemake/tools/picardOld.jar"
+
+########################
+# SECTION par_trimming #
+########################
+par_trimming:
+
+  # STRING. Default "1:30:4:5:true". ILLUMINACLIP value. See trimmomatic manual
+  trimmomatic_ILLUMINACLIP: "1:30:4:1:true"
+
+  # INTEGER. Default 3. TRAILING value. See trimmomatic manual
+  trimmomatic_trailing: 3
+
+  # INTEGER. Default 20. MINLEN value. See trimmomatic manual
+  trimmomatic_minlen: 20
+
+  # STRING. Default "phred33". Phred type. See trimmomatic manual. The "-" is added automatically by the Snakemake pipeline.
+  #trimmomatic_phredType: "phred33"
+
+#####################
+# SECTION par_align #
+#####################
+par_align:
+
+  # STRING. Default "--very-sensitive". Sensitivity. Leave empty for the default sensitivity. See bowtie2 manual.
+  bowtie2_sensitivity: "--very-sensitive"
+
+  # INTEGER. Default 2000. Value for parameter X. See bowtie2 manual.
+  bowtie2_maxFragmentLength: 2000
+
+  # STRING. Default "/g/scb/zaugg/zaugg_shared/annotations/hg19/referenceGenome/Bowtie2/hg19". Absolute path to the reference genome + the prefix of the file with the index. In the default case, we use those files within /g/scb/zaugg/zaugg_shared/annotations/hg19/referenceGenome/Bowtie2 that start with hg19. See bowtie2 manual for more details.
+  bowtie2_refGenome: "/g/scb2/zaugg/zaugg_shared/annotations/hg38/referenceGenome/Bowtie2/hg38"
+
+  # STRING. Default “hg38”. Reference genome assembly version. Must match the one used by the alignment program.
+  assemblyVersion: "hg38"
+
+
+#########################
+# SECTION par_postalign #
+#########################
+par_postalign:
+
+  # INTEGER. Default 10. Minimum MAPQ score. Reads with a lower MAPQ quality will be removed during the post-alignment processing.
+  minMAPQscore: 10
+
+  # STRING. Default “LENIENT”. Value of the VALIDATION_STRINGENCY from SortSam (Picard tools). See the manual for details.
+  ValidationStringencySortSam: "LENIENT"
+
+  # STRING. Default “SILENT”. Value of the VALIDATION_STRINGENCY from MarkDuplicates (Picard tools). See the manual for details.
+  ValidationStringencyMarkDuplicates: "SILENT"
+
+  # STRING. Defauled “ID”. Used for filtering reads. Relates to the one letter abbreviations for CIGAR strings such as I for insertion and D for deletion.  Specify all the one letter abbreviations in the CIGAR string of a read here that should be filtered. "ID" keeps a read only if the CIGAR string does not contain the letters "I" and "D" (e.g., only M for example)
+  CIGAR: "ID"
+
+  # INTEGER. Default 4 and -5, respectively. Adjustment of the read start positions on the forward and reverse strand and should be a positive or negative number, respectively. See the Buenrostro paper for details. FOR DATA TYPE ATAC-SEQ ONLY
+  adjustRSS_forward: 4
+  adjustRSS_reverse: 5
+
+#######################
+# SECTION par_scripts #
+#######################
+
+# NOTE: There is usually no need to modify the parameters here.
+
+par_scripts:
+
+  # INTEGER. Default 4000. The region size in bp that specifies what is considered within a TSS. The STATS script does a TSS enrichment test to test whether or not ATAC-Seq reads are primarily located within annotated TSS as opposed to outside of TSS regions. A value of 4000 means the region from -2kb up to +2kb of annotated TSS.
+  STATS_script_withinThr: 4000
+
+  # INTEGER. Default 1000. The size of the region adjacent to the within TSS region that is considered outside of a TSS. A value of 1000 therefore denotes the 1kb region up- and downstream of the within TSS region (from -3 to -2kb upstream and from +2 to +3 kb downstream of annotated TSS.)
+  STATS_script_outsideThr: 1000
+
+  # STRING. Default "protein_coding". Gene type to keep / do the analayses for. Allowed are gene types as specified by GENCODE. The default is "protein_coding"
+  STATS_script_geneTypesToKeep: "protein_coding"
+
+  # INTEGER. Default 600. Fragment length cutoff. All reads with a fragment length less than this value will be filtered for the purpose of this script.
+  FL_distr_script_cutoff: 600
+
+
+
+###########################
+# SECTION par_peakCalling #
+###########################
+par_peakCalling:
+
+  # STRING. Default "--nolambda –nomodel". Peak calling model for non-stringent peak calling.
+  modelNonStringent: "--nolambda --nomodel"
+
+  # STRING. Default "—nomodel”. Peak calling model for stringent peak calling.
+  modelStringent: "--nomodel"
+
+  # NUMERIC [0, 1]. Default 0.01. Minimum q-value threshold for stringent peak calling.
+  modelStringent_minQValue: 0.01
+
+  # NUMERIC [0, 1]. Default 0.1. Minimum q-value threshold for non-stringent peak calling.
+  modelNonStringent_minQValue: 0.1
+
+  # INTEGER. Default 10000. Value for slocal parameter for non-stringent peak calling.
+  modelNonStringent_slocal: 10000
+
+  # the enxt 3 options are only relevant if encodePeaks=true. ignored otherwise
+
+  # NUMERIC [0, 1]. Default 0.1. p-value threshold for ENCODE peak calling.
+  Encode_pValThreshold: 0.1
+
+  # STRING. Default "--nomodel --shift -75 --extsize 150 --broad --keep-dup all". Model for Encode peak calling (broad and gapped mode)
+  Encode_modelBroadAndGapped: "--nomodel --shift -75 --extsize 150 --broad --keep-dup all"
+
+  # STRING. Default "--nomodel --shift -75 --extsize 150 -B --SPMR --keep-dup all –call-summits". Model for Encode peak calling (narrow mode)
+  Encode_modelNarrow: "--nomodel --shift -75 --extsize 150 -B --SPMR --keep-dup all --call-summits"
+
+
+#########################
+# SECTION par_deepTools #
+#########################
+par_deepTools:
+
+  # INTEGER. The default depends on the genome assembly, see the templates for details. Length of the “mappable” genome in bp as defined by deepTools (see https://deeptools.readthedocs.io/en/develop/content/feature/effectiveGenomeSize.html).
+  #effectiveGenomeSize: 2750000000
+  readLength: ""
+
+  # STRING. Default “normalizeTo1x”. Either "normalizeTo1x NUMBER" OR "normalizeUsingRPKM" (note the missing (!) leading "--"), where NUMBER denotes the number of base pairs that are mappable. Beware of the mapping dependence on the read length for this number: The reported numbers on the websites are for 30bp reads, and we now have much longer reads usually. See Koehler et al. (2011) for numbers1.
+  bamCoverage_normalizationCoverage: "RPGC"
+
+  # INTEGER. Default 10. Size of the bins, in bases
+  bamCoverage_binSize: 10
+
+  # BOOLEAN. Default true. If the data is pair ended a read will by default be extended to match its pair.
+  bamCoverage_extendReads: true
+
+  # STRING. Default "--extendReads  --centerReads". Additional options that are supported by bamCoverage. Note that the "--" or "-" has to be present here.
+  bamCoverage_otherOptions: "--centerReads"
diff --git a/example/templates/mm10/README b/example/templates/dev/mm10/README
similarity index 100%
rename from example/templates/mm10/README
rename to example/templates/dev/mm10/README
diff --git a/example/templates/dev/mm10/config.yaml b/example/templates/dev/mm10/config.yaml
new file mode 100644
index 0000000000000000000000000000000000000000..ba21d851ed32b3b9275157b781bf5acf36b1939a
--- /dev/null
+++ b/example/templates/dev/mm10/config.yaml
@@ -0,0 +1,238 @@
+
+  ################################################
+  ################################################
+  # CONFIGURATION FILE FOR THE ATAC-SEQ PIPELINE #
+  ################################################
+  ################################################
+
+  # This format allows comments and is therefore easier to work with than the json format we used before.
+  # Quotation marks are optional for strings. Make sure to put ": " (that is, hyphen space) as separator
+
+
+  ##################
+  # SECTION output #
+  ##################
+  output:
+
+    # STRING. Absolute path to the output directory. Will be created if not yet present.
+    outdir: "/FULL/PATH/TO/OUTPUT/SAME/AS/IN/CLUSTER.json"
+
+    # BOOLEAN. “true” or “false”. Default “true”. Should the pipeline use MACS2 to call peaks? If yes, peaks will be called in 3 different flavors (ENCODE, stringent, non-stringent). See the section “par_peakCalling” in the configuration file for details.
+    doPeakCalling: true
+
+    # BOOLEAN. "true" or "false". Default "false". Do peak calling also, in addition to the "stringent" and "non-stringent" flavour, according to ENCODE guidelines (as of 2018, we so far rarely used them)
+    encodePeaks: false
+
+    # BOOLEAN. “true” or “false”. Default “false”. Should, in addition to treating all input files separately, the pipeline merge all replicate files and do all downstream analysis in addition? If set to true, for each individual as specified in the sample table, all samples belonging to this individual ID will be merged to one file after the post-processing. This makes only sense if you have more than 1 replicate per individual.
+    alsoMergeReplicates: true
+
+    # BOOLEAN. "true" or "false". Default "false". Should peaks be annotated to find out which genomic regions they overlap with? Still being tested, set to false if it causes errors
+    annotatePeaks: true
+
+    # BOOLEAN. “true” or “false”. Default “false”. Should GC bias be assessed and corrected for the final BAM files in addition to the non-corrected ones? If set to true, the GC bias will be assessed and corrected using deepTools. Additionally, all downstream steps of the pipeline will be done for both GC-corrected and original files (including peak calling, PCA, coverage, etc). This greatly increases running time of th pipeline and essentially doubles the number of files. We recommend setting this to false unless the GC bias is important for downstream analyses
+    # !!! Important: Leave at false for now, I have to fix an issue in the Snakefile that currently causes an error !!
+    correctGCBias: false
+
+    # BOOLEAN. “true” or “false”. Default “false”. Should base qualities be recalibrated using GATK to detect and correct systematic errors in base quality scores? This step can be time-consuming and needs the following other parameters: GATK_jar, knownSNPs, knownIndels. I recommend turning this off for now.
+    doBaseRecalibration: false
+
+    # BOOLEAN. “true” or “false”. Default “true”. Should the pipeline produce coverage files and diagnostic plots? If set to true, coverage files for the final BAM files in bigwig and bedgraph.gz format will be produced as well as a coverage plot using deepTools.
+    generateCoverageFiles: true
+
+    # BOOLEAN. “true” or “false”. Default “true”. Not implemented yet.
+    doIDR: false
+
+
+  ###################
+  # SECTION general #
+  ###################
+  general:
+
+    # INTEGER. Default 12. Maximum number of cores per rule. For local computation, the minimum of this value and the --cores parameter will define the number of CPUs per rule, while in a cluster setting, the minimum of this value and the number of cores on the node the jobs runs is used.
+    maxCoresPerRule: 12
+
+  ###################
+  # SECTION samples #
+  ###################
+  samples:
+
+    # STRING. No default. Absolute path to the sample summary file. See section 2.4 for details.
+  summaryFile: "/FULL/PATH/TO/samples.csv"
+
+    # BOOLEAN. “true” or “false”. Default “true”. Paired-end data? Single-end ATAC-Seq data is not yet supported with this pipeline. If set to "false", the Snakemake pipeline will abort in the beginning.
+    pairedEnd: true
+
+    # STRING. "ATACseq" or "ChIPseq". Only data type specific steps are executed. Currently, if set to ChIP-Seq, ATAC-seq specific steps like RSS adjustment are not executed.
+    dataType: "ATACseq"
+
+  ###########################
+  # SECTION additionalInput #
+  ###########################
+  additionalInput:
+
+    # STRING. Absolute path to the scripts directory. Either "/g/scb2/zaugg/carnold/Projects/AtacSeq/src/Snakemake/R" or "/g/scb2/zaugg/carnold/Projects/AtacSeq/src/Snakemake/dev/R"
+    scriptsDir: "/g/scb2/zaugg/carnold/Projects/AtacSeq/src/Snakemake/dev/R"
+
+    # STRING. Absolute path to the adapters file for Trimmomatic in fasta format. Default “/g/scb2/zaugg/zaugg_shared/Programs/Trimmomatic-0.33/adapters/NexteraPE-PE.fa”. There is usually no need to change this unless for your experiment, this adapter file is not suited.
+    trimmomatic_adapters: "/g/scb2/zaugg/zaugg_shared/Programs/Trimmomatic-0.33/adapters/NexteraPE-PE.fa"
+
+  # STRING. Absolute path to a BED file that contains the genomic regions that should be filtered from the peaks. The default depends on the genome assembly, see the templates for details. Only needed if doPeakCalling is set to true.
+  blacklistRegions: "/g/scb2/zaugg/zaugg_shared/annotations/mm10/blacklisted/mm10-blacklist.v2.bed"
+
+
+  # STRING. The default depends on the genome assembly, see the templates for details. Absolute path to the reference genome in fasta and 2bit format, respectively, both of which have to correspond to the same genome assembly version as used for the alignment as well as the database of polymorphic sites (knownSNPs and knownIndels, if applicable).
+  refGenome_fasta: "/g/scb2/zaugg/zaugg_shared/annotations/mm10/reference/mm10.fa"
+  refGenome_2bit: "/g/scb2/zaugg/zaugg_shared/annotations/mm10/reference/mm10.2bit"
+
+  # STRING. The default depends on the genome assembly, see the templates for details. Absolute path to an genome annotation file in GTF format.
+  annotationGTF: "/g/scb2/zaugg/zaugg_shared/annotations/mm10/Gencode_M16/gencode.vM16.annotation.gtf"
+
+
+
+  #######################
+  # SECTION executables #
+  #######################
+
+  # If run with Singularity, this section can be ignored
+
+  executables:
+
+    # STRING. Default “java”. Name of the executable for Java. Java version must be at least 1.8! Only needed if doBaseRecalibration is set to true
+    java: "java"
+
+    # STRING. Default "/g/scb2/zaugg/carnold/Projects/AtacSeq/src/Snakemake/tools/GenomeAnalysisTK.jar". Absolute path to a JAR file for the GATK suite. Only needed if doBaseRecalibration is set to true
+    GATK_jar: "/g/scb2/zaugg/carnold/Projects/AtacSeq/src/Snakemake/tools/GenomeAnalysisTK.jar"
+
+    # STRING. Default "/g/scb2/zaugg/zaugg_shared/Programs/Picardtools/picardOld.jar". Absolute path to a JAR file for the Picard suite.
+    PICARD_jar: "/g/scb2/zaugg/carnold/Projects/AtacSeq/src/Snakemake/tools/picardOld.jar"
+
+  ########################
+  # SECTION par_trimming #
+  ########################
+  par_trimming:
+
+    # STRING. Default "1:30:4:5:true". ILLUMINACLIP value. See trimmomatic manual
+    trimmomatic_ILLUMINACLIP: "1:30:4:1:true"
+
+    # INTEGER. Default 3. TRAILING value. See trimmomatic manual
+    trimmomatic_trailing: 3
+
+    # INTEGER. Default 20. MINLEN value. See trimmomatic manual
+    trimmomatic_minlen: 20
+
+    # STRING. Default "phred33". Phred type. See trimmomatic manual. The "-" is added automatically by the Snakemake pipeline.
+    #trimmomatic_phredType: "phred33"
+
+  #####################
+  # SECTION par_align #
+  #####################
+  par_align:
+
+    # STRING. Default "--very-sensitive". Sensitivity. Leave empty for the default sensitivity. See bowtie2 manual.
+    bowtie2_sensitivity: "--very-sensitive"
+
+    # INTEGER. Default 2000. Value for parameter X. See bowtie2 manual.
+    bowtie2_maxFragmentLength: 2000
+
+    # STRING. Default "/g/scb/zaugg/zaugg_shared/annotations/hg19/referenceGenome/Bowtie2/hg19". Absolute path to the reference genome + the prefix of the file with the index. In the default case, we use those files within /g/scb/zaugg/zaugg_shared/annotations/hg19/referenceGenome/Bowtie2 that start with hg19. See bowtie2 manual for more details.
+    bowtie2_refGenome: "/g/scb2/zaugg/zaugg_shared/annotations/mm10/Bowtie2_new_useThis/mm10"
+
+    # STRING. Default “mm10”. Reference genome assembly version. Must match the one used by the alignment program.
+    assemblyVersion: "mm10"
+
+
+  #########################
+  # SECTION par_postalign #
+  #########################
+  par_postalign:
+
+    # INTEGER. Default 10. Minimum MAPQ score. Reads with a lower MAPQ quality will be removed during the post-alignment processing.
+    minMAPQscore: 10
+
+    # STRING. Default “LENIENT”. Value of the VALIDATION_STRINGENCY from SortSam (Picard tools). See the manual for details.
+    ValidationStringencySortSam: "LENIENT"
+
+    # STRING. Default “SILENT”. Value of the VALIDATION_STRINGENCY from MarkDuplicates (Picard tools). See the manual for details.
+    ValidationStringencyMarkDuplicates: "SILENT"
+
+    # STRING. Defauled “ID”. Used for filtering reads. Relates to the one letter abbreviations for CIGAR strings such as I for insertion and D for deletion.  Specify all the one letter abbreviations in the CIGAR string of a read here that should be filtered. "ID" keeps a read only if the CIGAR string does not contain the letters "I" and "D" (e.g., only M for example)
+    CIGAR: "ID"
+
+    # INTEGER. Default 4 and -5, respectively. Adjustment of the read start positions on the forward and reverse strand and should be a positive or negative number, respectively. See the Buenrostro paper for details. FOR DATA TYPE ATAC-SEQ ONLY
+    adjustRSS_forward: 4
+    adjustRSS_reverse: 5
+
+  #######################
+  # SECTION par_scripts #
+  #######################
+
+  # NOTE: There is usually no need to modify the parameters here.
+
+  par_scripts:
+
+    # INTEGER. Default 4000. The region size in bp that specifies what is considered within a TSS. The STATS script does a TSS enrichment test to test whether or not ATAC-Seq reads are primarily located within annotated TSS as opposed to outside of TSS regions. A value of 4000 means the region from -2kb up to +2kb of annotated TSS.
+    STATS_script_withinThr: 4000
+
+    # INTEGER. Default 1000. The size of the region adjacent to the within TSS region that is considered outside of a TSS. A value of 1000 therefore denotes the 1kb region up- and downstream of the within TSS region (from -3 to -2kb upstream and from +2 to +3 kb downstream of annotated TSS.)
+    STATS_script_outsideThr: 1000
+
+    # STRING. Default "protein_coding". Gene type to keep / do the analayses for. Allowed are gene types as specified by GENCODE. The default is "protein_coding"
+    STATS_script_geneTypesToKeep: "protein_coding"
+
+    # INTEGER. Default 600. Fragment length cutoff. All reads with a fragment length less than this value will be filtered for the purpose of this script.
+    FL_distr_script_cutoff: 600
+
+
+
+  ###########################
+  # SECTION par_peakCalling #
+  ###########################
+  par_peakCalling:
+
+    # STRING. Default "--nolambda –nomodel". Peak calling model for non-stringent peak calling.
+    modelNonStringent: "--nolambda --nomodel"
+
+    # STRING. Default "—nomodel”. Peak calling model for stringent peak calling.
+    modelStringent: "--nomodel"
+
+    # NUMERIC [0, 1]. Default 0.01. Minimum q-value threshold for stringent peak calling.
+    modelStringent_minQValue: 0.01
+
+    # NUMERIC [0, 1]. Default 0.1. Minimum q-value threshold for non-stringent peak calling.
+    modelNonStringent_minQValue: 0.1
+
+    # INTEGER. Default 10000. Value for slocal parameter for non-stringent peak calling.
+    modelNonStringent_slocal: 10000
+
+    # the enxt 3 options are only relevant if encodePeaks=true. ignored otherwise
+
+    # NUMERIC [0, 1]. Default 0.1. p-value threshold for ENCODE peak calling.
+    Encode_pValThreshold: 0.1
+
+    # STRING. Default "--nomodel --shift -75 --extsize 150 --broad --keep-dup all". Model for Encode peak calling (broad and gapped mode)
+    Encode_modelBroadAndGapped: "--nomodel --shift -75 --extsize 150 --broad --keep-dup all"
+
+    # STRING. Default "--nomodel --shift -75 --extsize 150 -B --SPMR --keep-dup all –call-summits". Model for Encode peak calling (narrow mode)
+    Encode_modelNarrow: "--nomodel --shift -75 --extsize 150 -B --SPMR --keep-dup all --call-summits"
+
+
+  #########################
+  # SECTION par_deepTools #
+  #########################
+  par_deepTools:
+
+    # INTEGER. The default depends on the genome assembly, see the templates for details. Length of the “mappable” genome in bp as defined by deepTools (see https://deeptools.readthedocs.io/en/develop/content/feature/effectiveGenomeSize.html).
+    #effectiveGenomeSize: 2750000000
+    readLength: ""
+
+    # STRING. Default “normalizeTo1x”. Either "normalizeTo1x NUMBER" OR "normalizeUsingRPKM" (note the missing (!) leading "--"), where NUMBER denotes the number of base pairs that are mappable. Beware of the mapping dependence on the read length for this number: The reported numbers on the websites are for 30bp reads, and we now have much longer reads usually. See Koehler et al. (2011) for numbers1.
+    bamCoverage_normalizationCoverage: "RPGC"
+
+    # INTEGER. Default 10. Size of the bins, in bases
+    bamCoverage_binSize: 10
+
+    # BOOLEAN. Default true. If the data is pair ended a read will by default be extended to match its pair.
+    bamCoverage_extendReads: true
+
+    # STRING. Default "--extendReads  --centerReads". Additional options that are supported by bamCoverage. Note that the "--" or "-" has to be present here.
+    bamCoverage_otherOptions: "--centerReads"
diff --git a/example/templates/stable/hg19/README b/example/templates/stable/hg19/README
new file mode 100644
index 0000000000000000000000000000000000000000..7127c06b8941dbcc978579039164c66b8ae1fc29
--- /dev/null
+++ b/example/templates/stable/hg19/README
@@ -0,0 +1,5 @@
+The config file is now in a new and more suitable format: yaml instead of the previous json.
+
+For parameter explanations, see the yaml file and the comments above each parameter.
+
+Feel free to add as many comments to the yaml file as you want
diff --git a/example/templates/hg19/config.yaml b/example/templates/stable/hg19/config.yaml
similarity index 100%
rename from example/templates/hg19/config.yaml
rename to example/templates/stable/hg19/config.yaml
diff --git a/example/templates/stable/hg38/README b/example/templates/stable/hg38/README
new file mode 100644
index 0000000000000000000000000000000000000000..7127c06b8941dbcc978579039164c66b8ae1fc29
--- /dev/null
+++ b/example/templates/stable/hg38/README
@@ -0,0 +1,5 @@
+The config file is now in a new and more suitable format: yaml instead of the previous json.
+
+For parameter explanations, see the yaml file and the comments above each parameter.
+
+Feel free to add as many comments to the yaml file as you want
diff --git a/example/templates/hg38/config.yaml b/example/templates/stable/hg38/config.yaml
similarity index 100%
rename from example/templates/hg38/config.yaml
rename to example/templates/stable/hg38/config.yaml
diff --git a/example/templates/stable/mm10/README b/example/templates/stable/mm10/README
new file mode 100644
index 0000000000000000000000000000000000000000..7127c06b8941dbcc978579039164c66b8ae1fc29
--- /dev/null
+++ b/example/templates/stable/mm10/README
@@ -0,0 +1,5 @@
+The config file is now in a new and more suitable format: yaml instead of the previous json.
+
+For parameter explanations, see the yaml file and the comments above each parameter.
+
+Feel free to add as many comments to the yaml file as you want
diff --git a/example/templates/mm10/config.yaml b/example/templates/stable/mm10/config.yaml
similarity index 100%
rename from example/templates/mm10/config.yaml
rename to example/templates/stable/mm10/config.yaml
diff --git a/src/Snakemake/dev/R/aut_stats.R b/src/Snakemake/dev/R/aut_stats.R
index 24ab6be75013a2cc83afcfae549633fe883c6b77..62b9f1d170855691d47aa4d520ad8e690f5bd5c7 100755
--- a/src/Snakemake/dev/R/aut_stats.R
+++ b/src/Snakemake/dev/R/aut_stats.R
@@ -7,7 +7,7 @@ source("/g/scb2/zaugg/carnold/Projects/AtacSeq/src/Snakemake/R/functions.R")
 
 
 initFunctionsScript(packagesReq = NULL, minRVersion = "3.1.0", warningsLevel = 1, disableScientificNotation = TRUE)
-checkAndLoadPackages(c("checkmate", "futile.logger", "tidyverse", "reshape2", "tools", "rlist", "scales", "grDevices", "gridExtra", "Rsamtools", "GenomicRanges"), verbose = TRUE)
+checkAndLoadPackages(c("checkmate", "futile.logger", "tidyverse", "reshape2", "tools", "rlist", "scales", "grDevices", "gridExtra", "GenomicRanges"), verbose = TRUE)
     
 
 assertClass(snakemake, "Snakemake")
diff --git a/src/Snakemake/dev/Snakefile b/src/Snakemake/dev/Snakefile
index 6bcc2cc5c61669b3066df75fa81a07734b9495ef..58b44462675cf8366b6bf3bd732dfc06372659a0 100755
--- a/src/Snakemake/dev/Snakefile
+++ b/src/Snakemake/dev/Snakefile
@@ -24,7 +24,7 @@ import warnings
 
 # Enforce a minimum Snakemake version because of various features
 #min_version("5.9.1")
-min_version("6.0.0") 
+min_version("6.0.0")
 
 __author__  = "Christian Arnold"
 __license__ = "MIT"
@@ -174,12 +174,9 @@ genomeSizeDict = {
     'mm10': { "":2652783500, "50":2308125349, "75":2407883318, "100":2467481108, "150":2494787188,"200":2520869189 }
 }
 
-# genomeSizeDict = {
-#     'shit': { "":2913022398, "50":2701495761, "75":2747877777, "100":2805636331, "150":2862010578, "200":2887553303 },
-#     'hg19': { "":2864785220, "50":2685511504, "75":2736124973, "100":2776919808, "150":2827437033, "200":2855464000 }
-# }
 
-    # 'mm10': { "":2652783500, "50":2308125349, "75":2407883318, "100":2467481108, "150":2494787188,"200":2520869189 }
+# TODO: Check whether config["par_deepTools"]["readLength"] is among the supported options
+
 effectiveGenomeSize = genomeSizeDict[config["par_align"]["assemblyVersion"]][config["par_deepTools"]["readLength"]]
 
 
@@ -211,6 +208,8 @@ if pairedEnd:
 
 else:
     macs2_pairedEndMode = ""
+
+    # TODO: Currently not used at all. When single-end reads, this is not used at all
     duplicates_keepReadsFlag = 16 # reverse strand
 
 
@@ -386,7 +385,7 @@ mergeReplicates = config["output"]["alsoMergeReplicates"]
 doIDR           = config["output"]["doIDR"]
 correctGCBias   = config["output"]["correctGCBias"]
 coverageFiles   = config["output"]["generateCoverageFiles"]
-encodePeaks     = config["output"]["encodePeaks"] 
+encodePeaks     = config["output"]["encodePeaks"]
 
 
 
@@ -468,6 +467,14 @@ if doPeakCalling:
                             ext = ["pdf", "csv"])
 
     allResultFiles.append(individualPeaksNonStringent)
+    #
+    # # FRIP score
+    frip_StringentnonStringent = expand('{dir}/{dir2}/{sample}{GCBias}.final.{analysisType}.{peakType}_frip.pdf', dir = REPORTS_dir_frip, dir2 = ["stringent", "nonStringent"] , sample = allSamplesUnique, GCBias = GCBiasStr, analysisType = ["stringent", "nonStringent"], peakType = ["narrow"])
+
+    print(frip_StringentnonStringent)
+    # exit(1)
+    # TODO: Causes error in filterPeaks, something messed up with the wildcards
+    #allResultFiles.append(frip_StringentnonStringent)
 
     if annotatePeaks:
         allResultFiles.append(annotation_individualPeaksEncode)
@@ -475,7 +482,7 @@ if doPeakCalling:
         allResultFiles.append(annotation_individualPeaksNonStringent)
 
     # Consensus peaks]
-    
+
 
     if nSamplesUnique <=10:
         rangeOverlap = list(range(1, nSamplesUnique + 1, 1))
@@ -496,19 +503,19 @@ if doPeakCalling:
         #consensusPeaks_stringent_frip = expand("{dir1}/{dir2}/consensusPeaks{GCBias}.minOverlap{overlap}_frip.pdf", dir1 = REPORTS_dir_frip, dir2 = "stringent", GCBias = GCBiasStr, overlap = rangeOverlap)
         #consensusPeaks_nonStringent_frip = expand("{dir1}/{dir2}/consensusPeaks{GCBias}.minOverlap{overlap}_frip.pdf", dir1 = REPORTS_dir_frip, dir2 = "nonStringent", GCBias = GCBiasStr, overlap = rangeOverlap)
 
-        
+
         annotation_consensusPeaks_stringent = expand("{dir}/consensusPeaks{GCBias}.minOverlap{overlap}.{ext}", dir = REPORTS_dir_anno_STRINGENT, GCBias = GCBiasStr, overlap = rangeOverlap, ext = ["pdf", "csv"])
         annotation_consensusPeaks_nonStringent = expand("{dir}/consensusPeaks{GCBias}.minOverlap{overlap}.{ext}", dir = REPORTS_dir_anno_NONSTRINGENT, GCBias = GCBiasStr, overlap = rangeOverlap, ext = ["pdf", "csv"])
-        
-        annotation_consensusPeaks_ENCODE = expand("{dir}/consensusPeaks{GCBias}{peakType}.minOverlap{overlap}_annotation.{ext}", dir = REPORTS_dir_anno_ENCODE, GCBias = GCBiasStr, peakType = [".broadPeak", ".gappedPeak", ".narrowPeak"], overlap = rangeOverlap, ext = ["pdf", "csv"]) 
+
+        annotation_consensusPeaks_ENCODE = expand("{dir}/consensusPeaks{GCBias}{peakType}.minOverlap{overlap}_annotation.{ext}", dir = REPORTS_dir_anno_ENCODE, GCBias = GCBiasStr, peakType = [".broadPeak", ".gappedPeak", ".narrowPeak"], overlap = rangeOverlap, ext = ["pdf", "csv"])
         if encodePeaks:
             consensusPeaks_ENCODE = expand("{dir}/consensusPeaks{GCBias}{peakType}.minOverlap{overlap}.bed.gz", dir = PEAK_ENCODE_dir, GCBias = GCBiasStr, peakType = [".broadPeak", ".gappedPeak", ".narrowPeak"], overlap = rangeOverlap)
 
             #consensusPeaks_ENCODE_frip = expand("{dir1}/{dir2}/consensusPeaks{GCBias}{peakType}.minOverlap{overlap}_frip.pdf", dir1 = REPORTS_dir_frip, dir2 = "Encode", GCBias = GCBiasStr, peakType = [".broadPeak", ".gappedPeak", ".narrowPeak"], overlap = rangeOverlap)
             #allResultFiles.append(consensusPeaks_ENCODE_frip)
 
-            PCA_peaks = expand("{dir}/{basename}_PCAPlot_peaks{GCBias}_{peakType}_minOverlap{overlap}.pdf", 
-                dir = REPORTS_dir_PCA,basename = ["allSamples"], GCBias = GCBiasStr, overlap = rangeOverlap, 
+            PCA_peaks = expand("{dir}/{basename}_PCAPlot_peaks{GCBias}_{peakType}_minOverlap{overlap}.pdf",
+                dir = REPORTS_dir_PCA,basename = ["allSamples"], GCBias = GCBiasStr, overlap = rangeOverlap,
                 peakType = ["stringent", "nonStringent", "EncodeGapped", "EncodeBroad", "EncodeNarrow"])
             #allResultFiles.append(consensusPeaks_ENCODE)
         else:
@@ -517,7 +524,7 @@ if doPeakCalling:
                 dir = REPORTS_dir_PCA,basename = ["allSamples"], GCBias = GCBiasStr, overlap = rangeOverlap,
                 peakType = ["stringent", "nonStringent"])
 
-        
+
 
         allResultFiles.append(consensusPeaks_stringent)
         allResultFiles.append(consensusPeaks_nonStringent)
@@ -834,7 +841,7 @@ if pairedEnd:
             """
 else:
     rule trimming_SE:
-        input: 
+        input:
             forwardStrand = rules.link_inputFiles.output.forwardStrand,
             report = rules.fastqc_BT_SE.output # Not really needed, but force execution here
         output:
@@ -1434,7 +1441,7 @@ if dataType == "ATACseq":
 #         bam     = postalign_AlignementSieve_outputName,
 #         stats   = postalign_AlignementSieve_outputName + ".stats",
 #         csv     =  postalign_AlignementSieve_outputName + ".csv.gz"
-#     log: 
+#     log:
 #     message: "{ruleDisplayMessage}Adjust read start sites for {input.bam:q} using alignmentSieve"
 #     threads: threadsMax
 #     conda: "condaEnvironments/deepTools.yaml"
@@ -1453,7 +1460,7 @@ if dataType == "ATACseq":
 #             samtools flagstat {output.bam:q} > {output.stats:q} &&
 #             samtools view {output.bam:q} | cut -f3,5,9 | gzip  -f > {output.csv:q}
 #         """
-        
+
 # --shift {params.forwardShift} -{params.reverseShift} {params.reverseShift} -{params.forwardShift} \
 
 Picard_sortFinal_outputName = ADJRSS_dir + '/{sample}.cleaned4.BQrecal.rmChrM.markedDup.rmDup.minMapQ' + str(par_minMAPQscore) + '.adjRSS.cleaned1.s'
@@ -1462,7 +1469,7 @@ Picard_sortFinal_outputName = ADJRSS_dir + '/{sample}.cleaned4.BQrecal.rmChrM.ma
 # VALIDATION_STRINGENCY=SILENT is important because any strict check might identify the out-of-reference alignments and abort otherwise.
 rule Picard_sortFinal:
      input:
-        bam = rules.postaling_RSS.output.bam if dataType == "ATACseq" else rules.postalign_MAPQ.output.bam 
+        bam = rules.postaling_RSS.output.bam if dataType == "ATACseq" else rules.postalign_MAPQ.output.bam
      output:
         bam = temp(Picard_sortFinal_outputName + '.bam')
      log: expand('{dir}/Picard_sortFinal.{{sample}}.log', dir = LOG_BENCHMARK_dir)
@@ -2002,7 +2009,6 @@ if nSamplesUnique > 1 and doPeakCalling:
             minOverlapValues = ",".join(map(str, rangeOverlap))
         script: script_PCA
 
-# @RIM here? if annotatePeaks
 
 if annotatePeaks and doPeakCalling:
     rule annotatePeaks:
@@ -2014,7 +2020,7 @@ if annotatePeaks and doPeakCalling:
             consensus_stringent    = consensusPeaks_stringent,
             consensus_nonStringent = consensusPeaks_nonStringent
             #pooled = rule.poolPeaksReplicateSamples.output.pooledPeaks,
-            #replicates = rule.poolPeaksReplicateSamples.output.replicatePeaks 
+            #replicates = rule.poolPeaksReplicateSamples.output.replicatePeaks
         output:
             annotation_individual_encode = annotation_individualPeaksEncode,
             annotation_individual_stringent = annotation_individualPeaksStringent,
@@ -2022,14 +2028,14 @@ if annotatePeaks and doPeakCalling:
             annotation_consensus_encode = annotation_consensusPeaks_ENCODE,
             annotation_consensus_stringent = annotation_consensusPeaks_stringent,
             annotation_consensus_nonStringent = annotation_consensusPeaks_nonStringent
-        log: expand('{dir}/annotatePeaks.R.log', dir = LOG_BENCHMARK_dir) 
+        log: expand('{dir}/annotatePeaks.R.log', dir = LOG_BENCHMARK_dir)
         singularity: "shub://chrarnold/Singularity_images:atac_seq_r"
         message: "{ruleDisplayMessage} Annotating peaks for all peak files with the script {script_annoPeaks}..."
         threads: 1
         params:
             tabulateOutput = "true",
             assemblyVersion = config["par_align"]["assemblyVersion"],
-            tssRegion = "-5000,5000" 
+            tssRegion = "-5000,5000"
         script: script_annoPeaks
 
 
@@ -2059,7 +2065,7 @@ if annotatePeaks and doPeakCalling:
 #         other            = "--plot --use-best-multisummit-IDR",
 #         inputFileType    = "narrowPeak" # File type of --samples and --peak-list
 #     conda: "condaEnvironments/idr.yaml"
-#     shell:    
+#     shell:
 #        """{idr_exec:q} \
 #         --samples {input.sample1Peaks} {input.sample2Peaks} \
 #         --input-file-type {params.inputFileType} \
@@ -2129,35 +2135,54 @@ rule fragment_length_distr:
 #     if encodePeaks:
 #         consensusPeaks_ENCODE_frip = expand("{dir1}/{dir2}/consensusPeaks{GCBias}{peakType}.minOverlap{overlap}_frip.pdf", dir1 = REPORTS_dir_frip, dir2 = "Encode", GCBias = GCBiasStr, peakType = [".broadPeak", ".gappedPeak", ".narrowPeak"], overlap = rangeOverlap)
 #         fripOutput.append(consensusPeaks_ENCODE_frip)
-    
-
 
-    
+#
+#
+# expand("{dir}/consensusPeaks{GCBias}.minOverlap{overlap}.bed.gz", dir = PEAK_STRINGENT_dir, GCBias = GCBiasStr, overlap = rangeOverlap)
+#
+#
+#
+# def generateInputFiles (wildcards):
+#     if len(getSampleBasenamesForIndividual(wildcards.individual)) == 0:
+#         raise AssertionError("Cannot determine sample basenames for wildcard " + wildcards.individual + ": " + str(len(getSampleBasenamesForIndividual(wildcards.individual))))
+#     return expand('{dir}/{samples}{GCBiasStr}final.{analysisType}.{peakType}Peak.filtered.bed.gz', dir = wildcards.dir, samples = getSampleBasenamesForIndividual(wildcards.individual), GCBiasStr = wildcards.GCBiasStr, analysisType = wildcards.analysisType, peakType = wildcards.peakType)
+#
+#
+# def getSampleBasenamesForIndividual(individual):
+#     """text"""
+#     sampleBasenames = numpy.asarray(samplesData.loc[samplesData["individual"] == individual, "sampleName"])
+#     return sampleBasenames
+#
+#
+#
+#         lambda wildcards: expand('{dir}/{samples}.final.bam', dir = FINAL_OUTPUT_dir, samples = getSampleBasenamesForIndividual(wildcards.individual))
 
-    
-# rule deeptools_plotEnrichment:
-#     input:
-#         bams = expand('{dir}/{samples}.final.bam', dir = FINAL_OUTPUT_dir, samples = allSamplesUnique),
-#         bed = expand('{dir}/{subdir}/{{basename}}.bed.gz', dir = PEAKCALLING_dir)
-#     output:
-#         fig = expand( '{dir}/{{subdir}}_{{basename}}_frip.pdf', dir = REPORTS_dir_frip)
-#     log : expand('{dir}/deeptools_plotEnrichment_{{subdir}}_{{basename}}.log', dir = LOG_BENCHMARK_dir)
-#     message: "{ruleDisplayMessage} Computing FRiP scores"
-#     threads: 1
-#     conda: "condaEnvironments/deepTools.yaml"
-#     singularity: "shub://chrarnold/Singularity_images:atac_seq_conda"
-#     params:
-#         plotTitle = "FRiP"
-#     shell:
-#         """
-#             plotEnrichment \
-#                 -b {input.bams} \
-#                 --BED {input.bed} \
-#                 -o {output.fig} \
-#                 -p {threads} \
-#                 --smartLabels \
-#                 --plotTitle {params.plotTitle} 2> {log:q}
-#         """
+rule deeptools_plotEnrichment:
+    input:
+        bams = expand('{dir}/{samples}.final.bam', dir = FINAL_OUTPUT_dir, samples = allSamplesUnique),
+        # bed = expand('{dir}/{subdir}/{{basename}}.bed.gz', dir = PEAKCALLING_dir)
+        beds = expand('{dir}/{{basename}}Peak' + '.filtered.bed.gz', dir = PEAKCALLING_dir)
+    output:
+        fig       = expand('{dir}/{{basename}}_frip.pdf', dir = REPORTS_dir_frip),
+        rawCounts = expand('{dir}/{{basename}}_frip.rawCounts', dir = REPORTS_dir_frip),
+    log : expand('{dir}/deeptools_plotEnrichment_{{basename}}.log', dir = LOG_BENCHMARK_dir)
+    message: "{ruleDisplayMessage} Computing FRiP scores for basename {wildcars.basename}"
+    threads: 1
+    conda: "condaEnvironments/deepTools.yaml"
+    singularity: "shub://chrarnold/Singularity_images:atac_seq_conda"
+    params:
+        plotTitle = "FRiP score"
+    shell:
+        """
+            plotEnrichment \
+                -b {input.bams} \
+                --BED {input.beds} \
+                -o {output.fig} \
+                -p {threads} \
+                --outRawCounts {output.rawCounts} \
+                --smartLabels \
+                --plotTitle {params.plotTitle} 2> {log:q}
+        """
 
 # bam = FINAL_OUTPUT_dir + '/{sample}.final.bam',
 # bed = PEAKCALLING_dir + '/{basename}.bed.gz'
@@ -2165,7 +2190,7 @@ rule fragment_length_distr:
 
 rule stats:
     input:
-        rules.deeptools_plotEnrichment.output, # test?
+        # rules.deeptools_plotEnrichment.output, # test?
         expand('{dir}/allSamples_fragmentLengthDistr.pdf',
             dir = REPORTS_dir_summary),
         expand('{dir}/{sample}.s.bam',
@@ -2553,7 +2578,7 @@ rule deepTools_plotPCA_genomeWide:
         else
             touch {output.data} {output.plot}
             echo "input file {input.coverage} is empty" > {log}
-        fi  
+        fi
         """
 
 # rule deepTools_plotPCA_peaksOnly:
@@ -2590,7 +2615,7 @@ rule deepTools_plotPCA_genomeWide:
 #         else
 #             touch {output.data} {output.plot}
 #             echo "input file {input.coverage} is empty" > {log}
-#         fi       
+#         fi
 #        """
 
 use rule deepTools_plotPCA_genomeWide as deepTools_plotPCA_peaksOnly with:
@@ -2604,7 +2629,7 @@ use rule deepTools_plotPCA_genomeWide as deepTools_plotPCA_peaksOnly with:
     params:
         other     = "--transpose",
         titlePlot = "PCA plot for all samples (peaks only)"
-    singularity: 
+    singularity:
         "shub://chrarnold/Singularity_images:atac_seq_conda"
 
 # Enforce that it runs at the very end