created new ANALYSIS_RULES and workflow for incorporating the cd-hit and quast analyses

2019_GDB/envs/bbmap.yaml

+name: bbmap
+channels:
+  - imp
+  - bioconda
+  - conda-forge
+  - defaults
+dependencies:
+  - _libgcc_mutex=0.1
+  - _openmp_mutex=4.5
+  - alsa-lib=1.1.5
+  - bbmap=38.79
+  - fontconfig=2.13.1
+  - freetype=2.10.2
+  - giflib=5.2.1
+  - icu=67.1
+  - jpeg=9c
+  - lcms2=2.9
+  - libgcc-ng=9.2.0
+  - libgomp=9.2.0
+  - libiconv=1.15
+  - libpng=1.6.37
+  - libstdcxx-ng=9.2.0
+  - libtiff=4.1.0
+  - libuuid=2.32.1
+  - libwebp-base=1.1.0
+  - libxcb=1.13
+  - libxml2=2.9.10
+  - lz4-c=1.9.2
+  - openjdk=11.0.1
+  - pthread-stubs=0.4
+  - xorg-fixesproto=5.0
+  - xorg-inputproto=2.3.2
+  - xorg-kbproto=1.0.7
+  - xorg-libx11=1.6.9
+  - xorg-libxau=1.0.9
+  - xorg-libxdmcp=1.1.3
+  - xorg-libxext=1.3.4
+  - xorg-libxfixes=5.0.3
+  - xorg-libxi=1.7.10
+  - xorg-libxrender=0.9.10
+  - xorg-libxtst=1.2.3
+  - xorg-recordproto=1.14.2
+  - xorg-renderproto=0.11.1
+  - xorg-xextproto=7.3.0
+  - xorg-xproto=7.0.31
+  - xz=5.2.5
+  - zlib=1.2.11
+  - zstd=1.4.4
+prefix: /home/users/sbusi/apps/miniconda3/envs/bbmap

2019_GDB/envs/cd-hit.yaml

@@ -7,24 +7,24 @@ channels:
   - etetoolkit
   - au-eoed
 dependencies:
-  - _libgcc_mutex=0.1=conda_forge
-  - _openmp_mutex=4.5=0_gnu
-  - bzip2=1.0.8=h516909a_2
-  - ca-certificates=2020.4.5.1=hecc5488_0
-  - cd-hit=4.8.1=h8b12597_3
-  - gawk=5.1.0=h516909a_0
-  - gettext=0.19.8.1=hc5be6a0_1002
-  - libffi=3.2.1=he1b5a44_1007
-  - libgcc-ng=9.2.0=h24d8f2e_2
-  - libgomp=9.2.0=h24d8f2e_2
-  - libidn2=2.3.0=h516909a_0
-  - libstdcxx-ng=9.2.0=hdf63c60_2
-  - libunistring=0.9.10=h14c3975_0
-  - llvm-openmp=8.0.1=hc9558a2_0
-  - mmseqs2=11.e1a1c=h2d02072_0
-  - openmp=8.0.1=0
-  - openssl=1.1.1f=h516909a_0
-  - wget=1.20.1=h22169c7_0
-  - zlib=1.2.11=h516909a_1006
+  - _libgcc_mutex=0.1
+  - _openmp_mutex=4.5
+  - bzip2=1.0.8
+  - ca-certificates=2020.4.5.1
+  - cd-hit=4.8.1
+  - gawk=5.1.0
+  - gettext=0.19.8.1
+  - libffi=3.2.1
+  - libgcc-ng=9.2.0
+  - libgomp=9.2.0
+  - libidn2=2.3.0
+  - libstdcxx-ng=9.2.0
+  - libunistring=0.9.10
+  - llvm-openmp=8.0.1
+  - mmseqs2=11.e1a1c
+  - openmp=8.0.1
+  - openssl=1.1.1g
+  - wget=1.20.1
+  - zlib=1.2.11
 prefix: /home/users/sbusi/apps/miniconda3/envs/cd-hit
 

2019_GDB/envs/diamond.yaml

@@ -2,23 +2,3 @@ channels:
  - bioconda
 dependencies:
  - diamond=0.9.25
-
-# if running on esb-compute-01 use the below configuration
-#name: diamond
-#channels:
-#  - bioconda
-#  - conda-forge
-#  - defaults
-#dependencies:
-#  - _libgcc_mutex=0.1
-#  - _openmp_mutex=4.5
-#  - boost-cpp=1.70.0
-#  - bzip2=1.0.8
-#  - diamond=0.9.32
-#  - icu=64.2
-#  - libgcc-ng=9.2.0
-#  - libgomp=9.2.0
-#  - libstdcxx-ng=9.2.0
-#  - xz=5.2.5
-#  - zlib=1.2.11
-#prefix: /home/susheel.busi/miniconda3/envs/diamond
\ No newline at end of file

2019_GDB/envs/flye_v2_7.yaml

+name: flye
+channels:
+  - conda-forge
+  - bioconda
+  - defaults
+  - r
+  - etetoolkit
+  - au-eoed
+dependencies:
+  - _libgcc_mutex=0.1
+  - _openmp_mutex=4.5
+  - ca-certificates=2020.4.5.1
+  - certifi=2020.4.5.1
+  - flye=2.7.1
+  - ld_impl_linux-64=2.34
+  - libffi=3.2.1
+  - libgcc-ng=9.2.0
+  - libgomp=9.2.0
+  - libstdcxx-ng=9.2.0
+  - ncurses=6.1
+  - openssl=1.1.1g
+  - pip=20.1
+  - python=3.7.6
+  - python_abi=3.7
+  - readline=8.0
+  - setuptools=46.1.3
+  - sqlite=3.30.1
+  - tk=8.6.10
+  - wheel=0.34.2
+  - xz=5.2.5
+  - zlib=1.2.11
+prefix: /home/users/sbusi/apps/miniconda3/envs/flye
+

2019_GDB/envs/quast.yaml

+name: quast
+channels:
+  - bioconda
+  - conda-forge
+  - defaults
+dependencies:
+  - blast=2.6.0
+  - ca-certificates=2020.4.5.1
+  - certifi=2020.4.5.1
+  - circos=0.69.8
+  - cycler=0.10.0
+  - expat=2.2.9
+  - fontconfig=2.13.1
+  - freetype=2.10.2
+  - giflib=5.2.1
+  - glimmerhmm=3.0.4
+  - icu=64.2
+  - joblib=0.14.1
+  - jpeg=9c
+  - kiwisolver=1.2.0
+  - libblas=3.8.0
+  - libcblas=3.8.0
+  - libcxx=10.0.0
+  - libffi=3.2.1
+  - libgd=2.2.5
+  - libgfortran=4.0.0
+  - libiconv=1.15
+  - liblapack=3.8.0
+  - libopenblas=0.3.9
+  - libpng=1.6.37
+  - libtiff=4.1.0
+  - libwebp=1.0.2
+  - libxml2=2.9.10
+  - llvm-openmp=10.0.0
+  - lz4-c=1.9.2
+  - matplotlib=3.2.1
+  - matplotlib-base=3.2.1
+  - ncurses=6.1
+  - numpy=1.18.4
+  - openjdk=11.0.1
+  - openssl=1.1.1g
+  - perl=5.26.2
+  - perl-autoloader=5.74
+  - perl-carp=1.38
+  - perl-clone=0.42
+  - perl-config-general=2.63
+  - perl-digest-perl-md5=1.9
+  - perl-dynaloader=1.25
+  - perl-exporter=5.72
+  - perl-exporter-tiny=1.002001
+  - perl-extutils-makemaker=7.36
+  - perl-font-ttf=1.06
+  - perl-gd=2.71
+  - perl-io-string=1.08
+  - perl-list-moreutils=0.428
+  - perl-list-moreutils-xs=0.428
+  - perl-math-bezier=0.01
+  - perl-math-round=0.07
+  - perl-math-vecstat=0.08
+  - perl-module-implementation=0.09
+  - perl-module-runtime=0.016
+  - perl-number-format=1.75
+  - perl-params-validate=1.29
+  - perl-pathtools=3.75
+  - perl-readonly=2.05
+  - perl-regexp-common=2017060201
+  - perl-scalar-list-utils=1.52
+  - perl-set-intspan=1.19
+  - perl-statistics-basic=1.6611
+  - perl-svg=2.84
+  - perl-text-format=0.59
+  - perl-threaded=5.26.0
+  - perl-time-hires=1.9760
+  - perl-try-tiny=0.30
+  - perl-xml-parser=2.44
+  - perl-xsloader=0.24
+  - pip=20.1
+  - pyparsing=2.4.7
+  - python=3.6.10
+  - python-dateutil=2.8.1
+  - python_abi=3.6
+  - quast=5.0.2
+  - readline=8.0
+  - setuptools=46.3.0
+  - simplejson=3.8.1
+  - six=1.14.0
+  - sqlite=3.30.1
+  - tk=8.6.10
+  - tornado=6.0.4
+  - wheel=0.34.2
+  - xz=5.2.5
+  - zlib=1.2.11
+  - zstd=1.4.4
+prefix: /home/users/sbusi/apps/miniconda3/envs/quast

2019_GDB/rules/ANALYSIS_RULES

+# For assessing the levels of partial and unique gene calls after prodigal analyses on assemblies
+
+import os
+from tempfile import TemporaryDirectory
+
+configfile: "config/CONFIG.yaml"
+DATA_DIR = config["data_dir"]
+RESULTS_DIR = config["results_dir"]
+DB_DIR=config["db_dir"]
+BARCODES=config["barcodes"]
+# ASSEMBLERS=config["assemblers"]
+MAPPERS=["bwa", "mmi"]
+# SAMPLES=config["samples"]
+SAMPLES=["ONT3_MG_xx_Rashi_S11"]
+BINNING_SAMPLES=config["binning_samples"]
+HYBRID_ASSEMBLER=config["hybrid_assembler"]
+
+
+#############################
+######## CD-HIT  ############
+#############################
+rule prodigal_partial_gene_comparison:
+    input:
+        int1=expand(os.path.join(RESULTS_DIR, "annotation/proteins/{assembler}/{sample}/final.contigs.faa"), assembler="metaspades", sample=SAMPLES),
+        int2=expand(os.path.join(RESULTS_DIR, "annotation/proteins/{assembler}/lr_{barcode}-sr_{sample}/contigs.faa"), assembler="metaspades_hybrid", barcode=BARCODES, sample=SAMPLES)
+    output:
+        dat1=expand(os.path.join(RESULTS_DIR, "analysis/cdhit/{assembler}_partial_counts.txt"), assembler="metaspades"),
+        dat2=expand(os.path.join(RESULTS_DIR, "analysis/cdhit/{assembler}_partial_counts.txt"), assembler="metaspades_hybrid"),
+        dat3=os.path.join(RESULTS_DIR, "analysis/cdhit/partial_gene_counts.txt")
+    log: os.path.join(RESULTS_DIR, "analysis/cdhit/compare.log")
+    shell:
+        """
+        (date &&\
+        grep -E -o "partial.{0,3}" {input.int1} | grep '10\|01\|11' | wc -l > {output.dat1} &&\
+        grep -E -o "partial.{0,3}" {input.int2} | grep '10\|01\|11' | wc -l > {output.dat2} &&\
+        cat {output.dat1} {output.dat2} | paste - - > {output.dat3} &&\
+        date) &> >(tee {log})
+        """
+
+rule bbmap_rename:
+    input:
+        mm1=rules.prodigal_partial_gene_comparison.input.int1,
+        mm2=rules.prodigal_partial_gene_comparison.input.int2
+    output:
+        mo1=os.path.join(RESULTS_DIR, "analysis/cdhit/spades.faa"),
+        mo2=os.path.join(RESULTS_DIR, "analysis/cdhit/hybrid.faa"),
+    log: os.path.join(RESULTS_DIR, "analysis/cdhit/bbmap.rename.log")
+    threads: config["mmseq2"]["threads"]
+    conda: "../envs/bbmap.yaml"
+    shell:
+        """
+        (date &&\
+        rename.sh in={input.mm1} out={output.mo1} prefix=spades ignorejunk=t &&\
+        rename.sh in={input.mm2} out={output.mo2} prefix=hybrid ignorejunk=t &&\
+        date) &> >(tee {log})
+        """
+
+rule cdhit_compare:
+    input:
+        co1=rules.bbmap_rename.output.mo1,
+        co2=rules.bbmap_rename.output.mo2
+    output:
+        form1=os.path.join(RESULTS_DIR, "analysis/cdhit/spades_hybrid.fa"),
+        form2=os.path.join(RESULTS_DIR, "analysis/cdhit/hybrid_spades.fa"),
+        form3=os.path.join(RESULTS_DIR, "analysis/cdhit/cdhit_unique_counts.txt"),
+        form4=os.path.join(RESULTS_DIR, "analysis/cdhit/cdhit_compare.done")
+    log: os.path.join(RESULTS_DIR, "analysis/cdhit/cdhit.unique.log")
+    conda: "../envs/cd-hit.yaml"
+    shell:
+        """
+        (date &&\
+        cd-hit-2d -i {input.co1} -i2 {input.co2} -o {output.form1} -c 0.9 -n 5 -d 0 -M 16000 -T 8 &&\
+        cd-hit-2d -i {input.co2} -i2 {input.co1} -o {output.form2} -c 0.9 -n 5 -d 0 -M 16000 -T 8 &&\
+        echo {output.form1} > tmp && grep -c '>' {output.form1} >> tmp &&\
+        echo {output.form2} >> tmp && grep -c '>' {output.form2} >> tmp &&\
+        cat tmp | paste - - > {output.form3} &&\
+        touch {output.form4}
+        date) &> >(tee {log})
+        """
+
+rule quast:
+    input:
+        expand(os.path.join(RESULTS_DIR, "assembly/{assembler}.fa"), assembler=["flye", "megahit", "metaspades", "metaspades_hybrid"])
+    output: 
+        os.path.join(RESULTS_DIR, "analysis/quast/report.txt")
+    log: os.path.join(RESULTS_DIR, "analysis/quast/quast.assemblies.log")
+    threads: config["mmseq2"]["threads"]
+    conda: "../envs/quast.yaml"
+    shell:
+        """
+        (date &&\
+        metaquast.py --max-ref-number 0 --threads {threads} {input} -o $(dirname {output}) &&\
+        date) &> >(tee {log})
+        """

2019_GDB/rules/BINNING_RULES

@@ -9,7 +9,7 @@ RESULTS_DIR = config["results_dir"]
 DB_DIR=config["db_dir"]
 BARCODES=config["barcodes"]
 ASSEMBLERS=config["assemblers"]
-MAPPERS=["bwa", "mm"]
+MAPPERS=["bwa", "mmi"]
 #SAMPLES=config["samples"]
 SAMPLES=["flye", "megahit", "metaspades_hybrid"]
 BINNING_SAMPLES=config["binning_samples"]

2019_GDB/rules/TAXONOMY_RULES

@@ -9,7 +9,7 @@ RESULTS_DIR = config["results_dir"]
 DB_DIR=config["db_dir"]
 BARCODES=config["barcodes"]
 ASSEMBLERS=config["assemblers"]
-MAPPERS=["bwa", "mm"]
+MAPPERS=["bwa", "mmi"]
 #SAMPLES=config["samples"]
 SAMPLES=["flye", "megahit", "metaspades_hybrid"]
 BINNING_SAMPLES=config["binning_samples"]

2019_GDB/workflows/analysis.smk

+# For running the ANALYSIS "workflow"
+
+import os
+from tempfile import TemporaryDirectory
+
+configfile: "config/CONFIG.yaml"
+DATA_DIR = config["data_dir"]
+RESULTS_DIR = config["results_dir"]
+DB_DIR=config["db_dir"]
+BARCODES=config["barcodes"]
+#ASSEMBLERS=config["assemblers"]
+MAPPERS=["bwa", "mmi"]
+# SAMPLES=config["samples"]
+#SAMPLES=["flye", "megahit", "metaspades_hybrid"]
+BINNING_SAMPLES=config["binning_samples"]
+HYBRID_ASSEMBLER=config["hybrid_assembler"]
+
+
+# specify which rules to run
+include:
+    '../rules/ANALYSIS_RULES'
+
+# Rule all for running the analysing partial gene counts in prodigal and unique gene counts with cd-hit on the proteins
+rule CDHIT:    
+    input:
+        expand(os.path.join(RESULTS_DIR, "analysis/cdhit/cdhit_compare.done")), 
+        expand(os.path.join(RESULTS_DIR, "analysis/quast/report.txt"))
+    output:
+        touch('analysis.done')

envs/bbmap.yaml

+name: bbmap
+channels:
+  - imp
+  - bioconda
+  - conda-forge
+  - defaults
+dependencies:
+  - _libgcc_mutex=0.1
+  - _openmp_mutex=4.5
+  - alsa-lib=1.1.5
+  - bbmap=38.79
+  - fontconfig=2.13.1
+  - freetype=2.10.2
+  - giflib=5.2.1
+  - icu=67.1
+  - jpeg=9c
+  - lcms2=2.9
+  - libgcc-ng=9.2.0
+  - libgomp=9.2.0
+  - libiconv=1.15
+  - libpng=1.6.37
+  - libstdcxx-ng=9.2.0
+  - libtiff=4.1.0
+  - libuuid=2.32.1
+  - libwebp-base=1.1.0
+  - libxcb=1.13
+  - libxml2=2.9.10
+  - lz4-c=1.9.2
+  - openjdk=11.0.1
+  - pthread-stubs=0.4
+  - xorg-fixesproto=5.0
+  - xorg-inputproto=2.3.2
+  - xorg-kbproto=1.0.7
+  - xorg-libx11=1.6.9
+  - xorg-libxau=1.0.9
+  - xorg-libxdmcp=1.1.3
+  - xorg-libxext=1.3.4
+  - xorg-libxfixes=5.0.3
+  - xorg-libxi=1.7.10
+  - xorg-libxrender=0.9.10
+  - xorg-libxtst=1.2.3
+  - xorg-recordproto=1.14.2
+  - xorg-renderproto=0.11.1
+  - xorg-xextproto=7.3.0
+  - xorg-xproto=7.0.31
+  - xz=5.2.5
+  - zlib=1.2.11
+  - zstd=1.4.4
+prefix: /home/users/sbusi/apps/miniconda3/envs/bbmap

envs/quast.yaml

+name: quast
+channels:
+  - bioconda
+  - conda-forge
+  - defaults
+dependencies:
+  - blast=2.6.0
+  - ca-certificates=2020.4.5.1
+  - certifi=2020.4.5.1
+  - circos=0.69.8
+  - cycler=0.10.0
+  - expat=2.2.9
+  - fontconfig=2.13.1
+  - freetype=2.10.2
+  - giflib=5.2.1
+  - glimmerhmm=3.0.4
+  - icu=64.2
+  - joblib=0.14.1
+  - jpeg=9c
+  - kiwisolver=1.2.0
+  - libblas=3.8.0
+  - libcblas=3.8.0
+  - libcxx=10.0.0
+  - libffi=3.2.1
+  - libgd=2.2.5
+  - libgfortran=4.0.0
+  - libiconv=1.15
+  - liblapack=3.8.0
+  - libopenblas=0.3.9
+  - libpng=1.6.37
+  - libtiff=4.1.0
+  - libwebp=1.0.2
+  - libxml2=2.9.10
+  - llvm-openmp=10.0.0
+  - lz4-c=1.9.2
+  - matplotlib=3.2.1
+  - matplotlib-base=3.2.1
+  - ncurses=6.1
+  - numpy=1.18.4
+  - openjdk=11.0.1
+  - openssl=1.1.1g
+  - perl=5.26.2
+  - perl-autoloader=5.74
+  - perl-carp=1.38
+  - perl-clone=0.42
+  - perl-config-general=2.63
+  - perl-digest-perl-md5=1.9
+  - perl-dynaloader=1.25
+  - perl-exporter=5.72
+  - perl-exporter-tiny=1.002001
+  - perl-extutils-makemaker=7.36
+  - perl-font-ttf=1.06
+  - perl-gd=2.71
+  - perl-io-string=1.08
+  - perl-list-moreutils=0.428
+  - perl-list-moreutils-xs=0.428
+  - perl-math-bezier=0.01
+  - perl-math-round=0.07
+  - perl-math-vecstat=0.08
+  - perl-module-implementation=0.09
+  - perl-module-runtime=0.016
+  - perl-number-format=1.75
+  - perl-params-validate=1.29
+  - perl-pathtools=3.75
+  - perl-readonly=2.05
+  - perl-regexp-common=2017060201
+  - perl-scalar-list-utils=1.52
+  - perl-set-intspan=1.19
+  - perl-statistics-basic=1.6611
+  - perl-svg=2.84
+  - perl-text-format=0.59
+  - perl-threaded=5.26.0
+  - perl-time-hires=1.9760
+  - perl-try-tiny=0.30
+  - perl-xml-parser=2.44
+  - perl-xsloader=0.24
+  - pip=20.1
+  - pyparsing=2.4.7
+  - python=3.6.10
+  - python-dateutil=2.8.1
+  - python_abi=3.6
+  - quast=5.0.2
+  - readline=8.0
+  - setuptools=46.3.0
+  - simplejson=3.8.1
+  - six=1.14.0
+  - sqlite=3.30.1
+  - tk=8.6.10
+  - tornado=6.0.4
+  - wheel=0.34.2
+  - xz=5.2.5
+  - zlib=1.2.11
+  - zstd=1.4.4
+prefix: /home/users/sbusi/apps/miniconda3/envs/quast

rules/ANALYSIS_RULES

+# For assessing the levels of partial and unique gene calls after prodigal analyses on assemblies
+
+import os
+from tempfile import TemporaryDirectory
+
+configfile: "config/CONFIG.yaml"
+DATA_DIR = config["data_dir"]
+RESULTS_DIR = config["results_dir"]
+DB_DIR=config["db_dir"]
+BARCODES=config["barcodes"]
+# ASSEMBLERS=config["assemblers"]
+MAPPERS=["bwa", "mmi"]
+# SAMPLES=config["samples"]
+SAMPLES=["ONT3_MG_xx_Rashi_S11"]
+BINNING_SAMPLES=config["binning_samples"]
+HYBRID_ASSEMBLER=config["hybrid_assembler"]
+
+
+#############################
+######## CD-HIT  ############
+#############################
+rule prodigal_partial_gene_comparison:
+    input:
+        int1=expand(os.path.join(RESULTS_DIR, "annotation/proteins/{assembler}/{sample}/final.contigs.faa"), assembler="metaspades", sample=SAMPLES),
+        int2=expand(os.path.join(RESULTS_DIR, "annotation/proteins/{assembler}/lr_{barcode}-sr_{sample}/contigs.faa"), assembler="metaspades_hybrid", barcode=BARCODES, sample=SAMPLES)
+    output:
+        dat1=expand(os.path.join(RESULTS_DIR, "analysis/cdhit/{assembler}_partial_counts.txt"), assembler="metaspades"),
+        dat2=expand(os.path.join(RESULTS_DIR, "analysis/cdhit/{assembler}_partial_counts.txt"), assembler="metaspades_hybrid"),
+        dat3=os.path.join(RESULTS_DIR, "analysis/cdhit/partial_gene_counts.txt")
+    log: os.path.join(RESULTS_DIR, "analysis/cdhit/compare.log")
+    shell:
+        """
+        (date &&\
+        grep -E -o "partial.{0,3}" {input.int1} | grep '10\|01\|11' | wc -l > {output.dat1} &&\
+        grep -E -o "partial.{0,3}" {input.int2} | grep '10\|01\|11' | wc -l > {output.dat2} &&\
+        cat {output.dat1} {output.dat2} | paste - - > {output.dat3} &&\
+        date) &> >(tee {log})
+        """
+
+rule bbmap_rename:
+    input:
+        mm1=rules.prodigal_partial_gene_comparison.input.int1,
+        mm2=rules.prodigal_partial_gene_comparison.input.int2
+    output:
+        mo1=os.path.join(RESULTS_DIR, "analysis/cdhit/spades.faa"),
+        mo2=os.path.join(RESULTS_DIR, "analysis/cdhit/hybrid.faa"),
+    log: os.path.join(RESULTS_DIR, "analysis/cdhit/bbmap.rename.log")
+    threads: config["mmseq2"]["threads"]
+    conda: "../envs/bbmap.yaml"
+    shell:
+        """
+        (date &&\
+        rename.sh in={input.mm1} out={output.mo1} prefix=spades ignorejunk=t &&\
+        rename.sh in={input.mm2} out={output.mo2} prefix=hybrid ignorejunk=t &&\
+        date) &> >(tee {log})
+        """
+
+rule cdhit_compare:
+    input:
+        co1=rules.bbmap_rename.output.mo1,
+        co2=rules.bbmap_rename.output.mo2
+    output:
+        form1=os.path.join(RESULTS_DIR, "analysis/cdhit/spades_hybrid.fa"),
+        form2=os.path.join(RESULTS_DIR, "analysis/cdhit/hybrid_spades.fa"),
+        form3=os.path.join(RESULTS_DIR, "analysis/cdhit/cdhit_unique_counts.txt"),
+        form4=os.path.join(RESULTS_DIR, "analysis/cdhit/cdhit_compare.done")
+    log: os.path.join(RESULTS_DIR, "analysis/cdhit/cdhit.unique.log")
+    conda: "../envs/cd-hit.yaml"
+    shell:
+        """
+        (date &&\
+        cd-hit-2d -i {input.co1} -i2 {input.co2} -o {output.form1} -c 0.9 -n 5 -d 0 -M 16000 -T 8 &&\
+        cd-hit-2d -i {input.co2} -i2 {input.co1} -o {output.form2} -c 0.9 -n 5 -d 0 -M 16000 -T 8 &&\
+        echo {output.form1} > tmp && grep -c '>' {output.form1} >> tmp &&\
+        echo {output.form2} >> tmp && grep -c '>' {output.form2} >> tmp &&\
+        cat tmp | paste - - > {output.form3} &&\
+        touch {output.form4}
+        date) &> >(tee {log})
+        """
+
+rule quast:
+    input:
+        expand(os.path.join(RESULTS_DIR, "assembly/{assembler}.fa"), assembler=["flye", "megahit", "metaspades", "metaspades_hybrid"])
+    output: 
+        os.path.join(RESULTS_DIR, "analysis/quast/report.txt")
+    log: os.path.join(RESULTS_DIR, "analysis/quast/quast.assemblies.log")
+    threads: config["mmseq2"]["threads"]
+    conda: "../envs/quast.yaml"
+    shell:
+        """
+        (date &&\
+        metaquast.py --max-ref-number 0 --threads {threads} {input} -o $(dirname {output}) &&\
+        date) &> >(tee {log})
+        """

rules/BINNING_RULES

@@ -9,7 +9,7 @@ RESULTS_DIR = config["results_dir"]
 DB_DIR=config["db_dir"]
 BARCODES=config["barcodes"]
 ASSEMBLERS=config["assemblers"]
-MAPPERS=["bwa", "mm"]
+MAPPERS=["bwa", "mmi"]
 #SAMPLES=config["samples"]
 SAMPLES=["flye", "megahit", "metaspades_hybrid"]
 BINNING_SAMPLES=config["binning_samples"]

rules/TAXONOMY_RULES

@@ -9,7 +9,7 @@ RESULTS_DIR = config["results_dir"]
 DB_DIR=config["db_dir"]
 BARCODES=config["barcodes"]
 ASSEMBLERS=config["assemblers"]
-MAPPERS=["bwa", "mm"]
+MAPPERS=["bwa", "mmi"]
 #SAMPLES=config["samples"]
 SAMPLES=["flye", "megahit", "metaspades_hybrid"]
 BINNING_SAMPLES=config["binning_samples"]

workflows/analysis.smk

+# For running the ANALYSIS "workflow"
+
+import os
+from tempfile import TemporaryDirectory
+
+configfile: "config/CONFIG.yaml"
+DATA_DIR = config["data_dir"]
+RESULTS_DIR = config["results_dir"]
+DB_DIR=config["db_dir"]
+BARCODES=config["barcodes"]
+#ASSEMBLERS=config["assemblers"]
+MAPPERS=["bwa", "mmi"]
+# SAMPLES=config["samples"]
+#SAMPLES=["flye", "megahit", "metaspades_hybrid"]
+BINNING_SAMPLES=config["binning_samples"]
+HYBRID_ASSEMBLER=config["hybrid_assembler"]
+
+
+# specify which rules to run
+include:
+    '../rules/ANALYSIS_RULES'
+
+# Rule all for running the analysing partial gene counts in prodigal and unique gene counts with cd-hit on the proteins
+rule CDHIT:    
+    input:
+        expand(os.path.join(RESULTS_DIR, "analysis/cdhit/cdhit_compare.done")), 
+        expand(os.path.join(RESULTS_DIR, "analysis/quast/report.txt"))
+    output:
+        touch('analysis.done')