Final bug fixes before submitting the paper

.gitlab-ci.yml

@@ -0,0 +1,3 @@
+# This file is a template, and might need editing before it works on your project.
+include:
+  template: Jobs/Code-Quality.gitlab-ci.yml

Snakefile

@@ -1,47 +1,105 @@
 
 configfile: "config.yaml"
+from snakemake.utils import Paramspace
+import pandas as pd
+import os
+import random
+df = pd.read_csv(os.path.join(config["input_dir"],config["barcodes"]), sep='\t', dtype="object").set_index('Sample')
+SAMPLES = df.index
+SAMPLE = SAMPLES[0]
+flowCell = df.Flowcell[0]
+lane = df.Lane[0]
+projectName=random.randint(1,1000000) #To ensure non overlapping tmp directories
+
+
+paramspace = Paramspace(pd.read_csv("src/parameter_test/paramTest.tsv", sep="\t"))
+
 
 
 if config["mode"] == "reference":
     include: "src/rules/fastqc-ref.rules"
-    include: "src/rules/with_reference.rules"
+    include: "src/rules/reference.rules"
     include: "src/rules/demultiplex.rules"
-    include: "src/rules/report-ref.rules"
+    include: "src/rules/trimming.rules"
+
 
 if config["mode"] == "denovo":
-    include: "src/rules/denovo_reference.rules"
     include: "src/rules/demultiplex.rules"
+    include: "src/rules/denovo.rules"
+    include: "src/rules/trimming.rules"
     include: "src/rules/fastqc.rules"
-    include: "src/rules/report.rules"
+
+if config["mode"] == "paramTest":
+    include: "src/rules/demultiplex.rules"
+    include: "src/rules/paramTest.rules"
+    include: "src/rules/trimming.rules"
+
 
 if config["mode"]== "reference":
+    rule all:
+        input: expand("{out}/output_demultiplex/clone-stacks/{sample}-Watson.1.fq.gz \
+            {out}/output_demultiplex/clone-stacks/{sample}-Watson.2.fq.gz \
+            {out}/output_demultiplex/clone-stacks/{sample}-Crick.1.fq.gz \
+            {out}/output_demultiplex/clone-stacks/{sample}-Crick.2.fq.gz \
+            {out}/fastqc/ \
+            {out}/multiQC_report.html \
+		    {out}/cutadapt/{sample}_trimmed_filt_merged.1.fq.gz \
+            {out}/cutadapt/{sample}_trimmed_filt_merged.2.fq.gz \
+		    {out}/alignment/{sample}_trimmed_filt_merged.1_bismark_bt2_pe.bam \
+		    {out}/methylation_calling/{sample}_bismark_bt2_pe.CX_report.txt.gz \
+            {out}/snp_calling/snp.vcf.gz".split(),out=config["output_dir"],sample=SAMPLES)
+
+if config["mode"]== "denovo":
     rule all:
         input: expand("{out}/output_demultiplex/Watson_R2.fq.gz \
             {out}/output_demultiplex/Watson_R1.fq.gz \
             {out}/output_demultiplex/Crick_R2.fq.gz \
             {out}/output_demultiplex/Crick_R1.fq.gz \
-            {out}/multiQC_report.html \
             {out}/fastqc/ \
-            {out}/trimmed/Watson_R1_trim_fastqc.html \
-            {out}/trimmed/Crick_R2_trim_fastqc.html \
-            {out}/report.html \
-            {out}/mapping/watson.bam \
-            {out}/mapping/crick.bam \
-            {out}/mapping/methylation.bed \
-            {out}/mapping/snp.vcf.gz".split(), out=config["output_dir"])
+            {out}/multiQC_report.html \
+            {out}/output_denovo/consensus_cluster.renamed.fa \
+		    {out}/cutadapt/{sample}_trimmed_filt_merged.1.fq.gz \
+            {out}/cutadapt/{sample}_trimmed_filt_merged.1.fq.gz \
+		    {out}/alignment/{sample}_trimmed_filt_merged.1_bismark_bt2_pe.bam \
+		    {out}/methylation_calling/{sample}_bismark_bt2_pe.CX_report.txt.gz \
+            {out}/snp_calling/snp.vcf.gz".split(),out=config["output_dir"],sample=SAMPLES)
 
+if config["mode"]== "paramTest":
+    rule all:
+        input: expand("{out}/output_demultiplex/clone-stacks/{samples}-Watson.1.fq.gz \
+            {out}/output_demultiplex/clone-stacks/{samples}-Watson.2.fq.gz \
+            {out}/output_demultiplex/clone-stacks/{samples}-Crick.1.fq.gz \
+            {out}/output_demultiplex/clone-stacks/{samples}-Crick.2.fq.gz \
+            {out}/output_demultiplex/Watson_R2.fq.gz \
+            {out}/output_demultiplex/Watson_R1.fq.gz \
+            {out}/output_demultiplex/Crick_R2.fq.gz \
+            {out}/output_demultiplex/Crick_R1.fq.gz \
+            {out}/paramTest/{params}/alignment/{sample}_trimmed_filt_merged.1_bismark_bt2_pe.bam \
+            {out}/paramTest/{params}/output_denovo/consensus_cluster.renamed.fa \
+            {out}/paramTest/Assembled.txt \
+            {out}/paramTest/averageDepth.txt \
+            {out}/paramTest/denovoParameter.tsv \
+            {out}/paramTest/denovoParameter.tiff \
+		    {out}/cutadapt/{sample}_trimmed_filt_merged.1.fq.gz \
+            {out}/cutadapt/{sample}_trimmed_filt_merged.2.fq.gz".split(),out=config["output_dir"],params=paramspace.instance_patterns,sample=SAMPLE,samples=SAMPLES)
+
 
-if config["mode"]== "denovo":
+if config["mode"] == "legacy":
+    include: "src/rules/legacy.rules"
+    include: "src/rules/demultiplex.rules"
+    include: "src/rules/report.rules"
+
+
+if config["mode"]== "legacy":
     rule all:
         input: expand("{out}/output_demultiplex/Watson_R2.fq.gz \
             {out}/output_demultiplex/Watson_R1.fq.gz \
             {out}/output_demultiplex/Crick_R2.fq.gz \
             {out}/output_demultiplex/Crick_R1.fq.gz \
-            {out}/fastqc/ \
-            {out}/multiQC_report.html \
-            {out}/report.html \
             {out}/mapping/watson.bam \
             {out}/mapping/crick.bam \
             {out}/output_denovo/consensus_cluster.renamed.fa \
             {out}/mapping/methylation.bed \
             {out}/mapping/snp.vcf.gz".split(), out=config["output_dir"])
+
+

config.yaml

@@ -17,10 +17,7 @@ barcodes: "barcodes.tsv"
 # the pipeline produces some temporary files. Please indicate the tmp location on your server (in most cases /tmp)
 tmpdir        : "/tmp"
 
-# some of the steps can be run in parallel. Please set the number of available computing threads on your system
-threads: "12"
-
-# mode of running pipeline (set denovo or reference)
+# mode of running pipeline (set denovo, reference or legacy. PLEASE NOTE: legacy is not supported)
 mode: "reference"
 
 # genome directory (leaave it blank in denovo mode)
@@ -32,20 +29,10 @@ genome: "consensus_cluster.renamed.fa"
 # advanced users have the possibility to change different parameter, leave them blank or write "default" to run them in default mode
 
 # parameters in the denovo reference creation:
-# identity: percentage of sequence identity in the last clustering step, in decimal number e.g. for 90% identity write 0.90, default 0.95
-# min-depth: minimal cluster depth in the first clustering step to include a cluster, default 0
-# max-depth: maximal cluster depth in the first clustering step to include a cluster, default 0
+# identity: percentage of sequence identity in the last clustering step, in decimal number e.g. for 95% identity write 0.95, default 0.97
+# min-depth: minimal cluster depth in the first clustering step to include a cluster, default 10
+# max-depth: maximal cluster depth in the first clustering step to include a cluster, default 10000
 param_denovo:
   identity: ""
   min-depth: ""
   max-depth: ""
-
-# parameters in the mpileup step (variant callin)
-# max-depth: At a position, read maximally INT reads per input file, default 10000000, check if mpileup automatically reduced the maximal depth even more. You might have to re-run your analysis with a lower max-depth value to avoid this.
-# min-MQ: Minimum mapping quality for an alignment to be used, default 0
-# min-BQ: Minimum base quality for a base to be considered, default 15
-
-param_SNPcalling:
-  max-depth: ""
-  min-MQ: ""
-  min-BQ: ""

src/demultiplex/barcode_stacks.py

@@ -118,7 +118,7 @@ def run_subprocess(cmd, args, log_message):
 cmd = "process_radtags -1 %s -2 %s " % (read1,
                                         read2)
 cmd += "-b %s -o %s " % (os.path.join(args.outputdir, "barcode_stacks.tsv"), os.path.join(args.outputdir, "clone-stacks"))
-cmd += "-r -D --inline_inline --barcode_dist_2 0 "
+cmd += "-r -D --inline_inline --barcode_dist_1 0 "
 cmd += "--renz_1 %s --renz_2 %s --retain_header --barcode_dist_2 0" % (line2split[bc_dict["ENZ_R1"][0]],
                                                      line2split[bc_dict["ENZ_R2"][0]])
 log = "Run process_radtags from stacks"

src/env/SNP-calling.yaml

@@ -0,0 +1,29 @@
+name: SNP-calling
+channels:
+  - conda-forge
+  - bioconda
+  - defaults
+  - r
+dependencies:
+  - _libgcc_mutex=0.1=conda_forge
+  - _openmp_mutex=4.5=1_gnu
+  - ca-certificates=2020.12.5=ha878542_0
+  - certifi=2020.12.5=py36h5fab9bb_1
+  - libgcc-ng=9.3.0=h2828fa1_18
+  - libgomp=9.3.0=h2828fa1_18
+  - ncurses=5.9=10
+  - openssl=1.0.2u=h516909a_0
+  - pigz=2.3.4=hed695b0_1
+  - pip=20.3.3=pyhd8ed1ab_0
+  - python=3.6.3=4
+  - python_abi=3.6=1_cp36m
+  - readline=7.0=0
+  - samtools=1.3.1=0
+  - setuptools=49.6.0=py36h5fab9bb_3
+  - sqlite=3.20.1=2
+  - tabix=0.2.6=ha92aebf_0
+  - tk=8.6.10=h21135ba_1
+  - wheel=0.36.2=pyhd3deb0d_0
+  - xz=5.2.5=h516909a_1
+  - zlib=1.2.11=h516909a_1010
+prefix: /home/NIOO.INT/fleurg/.conda/envs/SNP-calling

src/env/SNP_calling_epiDiverse.yaml

@@ -0,0 +1,126 @@
+name: snps
+channels:
+  - conda-forge
+  - bioconda
+  - defaults
+dependencies:
+  - _libgcc_mutex=0.1=conda_forge
+  - _openmp_mutex=4.5=1_gnu
+  - _r-mutex=1.0.1=anacondar_1
+  - bcftools=1.10.2=h4f4756c_3
+  - binutils_impl_linux-64=2.34=h2122c62_9
+  - binutils_linux-64=2.34=h47ac705_27
+  - bwidget=1.9.14=0
+  - bz2file=0.98=py_0
+  - bzip2=1.0.8=h516909a_3
+  - c-ares=1.16.1=h516909a_3
+  - ca-certificates=2020.6.20=hecda079_0
+  - cairo=1.16.0=h3fc0475_1005
+  - certifi=2020.6.20=py38h32f6830_0
+  - curl=7.71.1=he644dc0_5
+  - cycler=0.10.0=py_2
+  - expat=2.2.9=he1b5a44_2
+  - fontconfig=2.13.1=h1056068_1002
+  - freebayes=1.3.2=py38h40864fe_2
+  - freetype=2.10.2=he06d7ca_0
+  - fribidi=1.0.10=h516909a_0
+  - gcc_impl_linux-64=7.5.0=hda68d29_13
+  - gcc_linux-64=7.5.0=hf34d7eb_27
+  - gettext=0.19.8.1=hc5be6a0_1002
+  - gfortran_impl_linux-64=7.5.0=h1104b78_16
+  - gfortran_linux-64=7.5.0=ha781d05_27
+  - glib=2.66.0=h0dae87d_0
+  - graphite2=1.3.13=he1b5a44_1001
+  - graphviz=2.42.3=h0511662_0
+  - gsl=2.6=h294904e_0
+  - gxx_impl_linux-64=7.5.0=h64c220c_13
+  - gxx_linux-64=7.5.0=ha781d05_27
+  - harfbuzz=2.7.2=hee91db6_0
+  - htslib=1.10.2=hd3b49d5_1
+  - icu=67.1=he1b5a44_0
+  - jpeg=9d=h516909a_0
+  - kernel-headers_linux-64=2.6.32=h77966d4_13
+  - khmer=3.0.0a3=py38hfcc51a0_1
+  - kiwisolver=1.2.0=py38hbf85e49_0
+  - krb5=1.17.1=hfafb76e_3
+  - kwip=0.2.0=h8b12597_3
+  - lcms2=2.11=hbd6801e_0
+  - ld_impl_linux-64=2.34=hc38a660_9
+  - libblas=3.8.0=17_openblas
+  - libcblas=3.8.0=17_openblas
+  - libcurl=7.71.1=hcdd3856_5
+  - libdeflate=1.6=h516909a_0
+  - libedit=3.1.20191231=he28a2e2_2
+  - libev=4.33=h516909a_1
+  - libffi=3.2.1=he1b5a44_1007
+  - libgcc-ng=9.3.0=h24d8f2e_16
+  - libgfortran-ng=7.5.0=hdf63c60_16
+  - libgomp=9.3.0=h24d8f2e_16
+  - libiconv=1.16=h516909a_0
+  - liblapack=3.8.0=17_openblas
+  - libnghttp2=1.41.0=h8cfc5f6_2
+  - libopenblas=0.3.10=pthreads_hb3c22a3_4
+  - libpng=1.6.37=hed695b0_2
+  - libssh2=1.9.0=hab1572f_5
+  - libstdcxx-ng=9.3.0=hdf63c60_16
+  - libtiff=4.1.0=hc7e4089_6
+  - libtool=2.4.6=h516909a_1005
+  - libuuid=2.32.1=h14c3975_1000
+  - libwebp-base=1.1.0=h516909a_3
+  - libxcb=1.13=h14c3975_1002
+  - libxml2=2.9.10=h68273f3_2
+  - lz4-c=1.9.2=he1b5a44_3
+  - make=4.3=h516909a_0
+  - matplotlib=3.3.1=1
+  - matplotlib-base=3.3.1=py38h91b0d89_1
+  - ncurses=6.2=he1b5a44_1
+  - numpy=1.19.1=py38hbc27379_2
+  - olefile=0.46=py_0
+  - openssl=1.1.1g=h516909a_1
+  - pango=1.42.4=h7062337_4
+  - parallel=20200722=0
+  - pcre=8.44=he1b5a44_0
+  - pcre2=10.35=h2f06484_0
+  - perl=5.30.3=h516909a_1
+  - pillow=7.2.0=py38h9776b28_1
+  - pip=20.2.3=py_0
+  - pixman=0.38.0=h516909a_1003
+  - pthread-stubs=0.4=h14c3975_1001
+  - pyparsing=2.4.7=pyh9f0ad1d_0
+  - pysam=0.16.0.1=py38hbdc2ae9_1
+  - python=3.8.5=h1103e12_7_cpython
+  - python-dateutil=2.8.1=py_0
+  - python_abi=3.8=1_cp38
+  - r-base=4.0.2=he766273_1
+  - r-dotcall64=1.0_0=r40h9bbef5b_1005
+  - r-fields=11.4=r40h9bbef5b_0
+  - r-maps=3.3.0=r40hcdcec82_1004
+  - r-spam=2.5_1=r40hbf399a0_2
+  - readline=8.0=he28a2e2_2
+  - samtools=1.10=h2e538c0_3
+  - screed=1.0.4=pyh9f0ad1d_0
+  - sed=4.8=hbfbb72e_0
+  - setuptools=49.6.0=py38h32f6830_0
+  - six=1.15.0=pyh9f0ad1d_0
+  - sqlite=3.33.0=h4cf870e_0
+  - sysroot_linux-64=2.12=h77966d4_13
+  - tk=8.6.10=hed695b0_0
+  - tktable=2.10=h555a92e_3
+  - tornado=6.0.4=py38h1e0a361_1
+  - wheel=0.35.1=pyh9f0ad1d_0
+  - xorg-kbproto=1.0.7=h14c3975_1002
+  - xorg-libice=1.0.10=h516909a_0
+  - xorg-libsm=1.2.3=h84519dc_1000
+  - xorg-libx11=1.6.12=h516909a_0
+  - xorg-libxau=1.0.9=h14c3975_0
+  - xorg-libxdmcp=1.1.3=h516909a_0
+  - xorg-libxext=1.3.4=h516909a_0
+  - xorg-libxpm=3.5.13=h516909a_0
+  - xorg-libxrender=0.9.10=h516909a_1002
+  - xorg-libxt=1.1.5=h516909a_1003
+  - xorg-renderproto=0.11.1=h14c3975_1002
+  - xorg-xextproto=7.3.0=h14c3975_1002
+  - xorg-xproto=7.0.31=h14c3975_1007
+  - xz=5.2.5=h516909a_1
+  - zlib=1.2.11=h516909a_1009
+  - zstd=1.4.5=h6597ccf_2

src/env/bcftools.yaml

@@ -0,0 +1,37 @@
+name: bcftools
+channels:
+  - conda-forge
+  - bioconda
+  - r
+  - defaults
+dependencies:
+  - _libgcc_mutex=0.1=conda_forge
+  - _openmp_mutex=4.5=1_gnu
+  - bcftools=1.11=h7c999a4_0
+  - bzip2=1.0.8=h7f98852_4
+  - c-ares=1.11.0=h470a237_1
+  - ca-certificates=2020.12.5=ha878542_0
+  - gsl=2.6=he838d99_2
+  - htslib=1.11=hd3b49d5_2
+  - krb5=1.17.2=h926e7f8_0
+  - libblas=3.9.0=7_openblas
+  - libcblas=3.9.0=7_openblas
+  - libcurl=7.71.1=hcdd3856_8
+  - libdeflate=1.7=h7f98852_5
+  - libedit=3.1.20191231=he28a2e2_2
+  - libev=4.33=h516909a_1
+  - libgcc-ng=9.3.0=h2828fa1_18
+  - libgfortran-ng=9.3.0=hff62375_18
+  - libgfortran5=9.3.0=hff62375_18
+  - libgomp=9.3.0=h2828fa1_18
+  - libnghttp2=1.41.0=hab1572f_1
+  - libopenblas=0.3.12=pthreads_h4812303_1
+  - libssh2=1.9.0=hab1572f_5
+  - libstdcxx-ng=9.3.0=h6de172a_18
+  - ncurses=6.2=h58526e2_4
+  - openssl=1.1.1i=h7f98852_0
+  - perl=5.32.0=h36c2ea0_0
+  - tk=8.6.10=h21135ba_1
+  - xz=5.2.5=h516909a_1
+  - zlib=1.2.11=h516909a_1010
+prefix: /home/NIOO.INT/maartenp/.conda/envs/bcftools