A Snakemake-wrapper for evaluating de novo bacterial isolate genome assemblies, e.g. from Oxford Nanopore (ONT) or Illumina sequencing, using multiple programs. The results are summarized in a HTML report.

The workflow is published in Snakemake workflows for long-read bacterial genome assembly and evaluation in GigaByte.

Following programs are included in the workflow:

• pomoxis assess_assembly and assess_homopolymers
• dnadiff from the mummer package
• NucDiff
• QUAST
• BUSCO
• ideel, which uses prodigal and diamond
• bakta

Clone repository, for example:

git clone https://github.com/pmenzel/score-assemblies.git /opt/software/score-assemblies

Create a new conda environment containing all necessary programs:

conda env create -n score-assemblies --file /opt/software/score-assemblies/env/environment.yaml

and activate the environment:

conda activate score-assemblies

First, prepare a data folder, which must contain subfolders assemblies/ containing the assemblies.
Additionally, the sub-foldersreferences/ and references-protein/ can contain reference genomes and reference proteins with which the assemblies and predicted proteins will be compared.
For example:

.
├── assemblies
│   ├── example-mtb_flyehq4.fa
│   ├── example-mtb_flyehq4+medaka.fa
│   ├── example-mtb_flyehq.fa
│   ├── example-mtb_flyehq+racon4.fa
│   ├── example-mtb_flyehq+racon4+medaka.fa
│   ├── example-mtb_raven4.fa
│   ├── example-mtb_raven4+medaka.fa
│   ├── example-mtb_raven4+medaka+pilon.fa
│   └── example-mtb_unicycler.fa
├── references
│   └── AL123456.3.fa
└── references-protein
    └── AL123456.3.faa

NB: The assembly and reference FASTA files need to have the .fa extension and protein reference FASTA files need to have the extension .faa.

This is the same folder structure used by ont-assembly-snake, i.e. score-assemblies can be run directly in the same folder.

To run the workflow, e.g. with 20 threads, use this command:

snakemake -s /opt/software/score-assemblies/Snakefile --cores 20 --use-conda

Output files of each program will be written to various folders in score-assemblies-data/.

If no references are supplied, then only ideel and BUSCO are done, otherwise score-assemblies will run these programs on each assembly:

Each assembly will be compared against each reference genome using the assess_assembly and assess_homopolymers scripts from pomoxis. Additionally to the tables and plots generated by these programs, summary plots for each reference genome will be plotted in score-assemblies-data/pomoxis/<reference>_assess_assembly_all_meanQ.pdf.

Set the lineage via the snakemake call:

snakemake -s /opt/software/score-assemblies/Snakefile --cores 20 --config busco_lineage=bacillales

If not set, the default lineage bacteria will be used. Available datasets can be listed with busco --list-datasets

The number of complete, fragmented and missing BUSCOs per assembly is tabulated in the file score-assemblies-data/busco/all_stats.tsv and also drawn as dotplot in score-assemblies-data/busco/busco_stats.pdf.

Each assembly is compared with each reference and the output files will be located in score-assemblies-data/dnadiff/<reference>/<assembly>-dnadiff.report. The values for AvgIdentity (from 1-to-1 alignments) and TotalIndels are extracted from these files and are plotted for each reference in score-assemblies-data/dnadiff/<reference>_dnadiff_stats.pdf.

Each assembly is compared with each reference and the output files will be located in the folder score-assemblies-data/nucdiff/<reference>/<assembly>-nucdiff/. The values for Insertions, Deletions, and Substitutions are extracted from the file results/nucdiff_stat.out and are drawn for each reference in score-assemblies-data/nucdiff/<reference>_nucdiff_stats.pdf.

One QUAST report is generated for each reference genome, containing the results for all assemblies. The report files are located in score-assemblies-data/quast/<reference>/report.html. The main report file score-assemblies-report.html also links the these individual reports.

Open reading frames are predicted from each assembly via Prodigal and are search in the Uniprot sprot database with diamond, retaining the best alignment for each ORF. For each assembly, the distribution of the ratios between length of the ORF and the matching database sequence are plotted to ideel/ideel_uniprot_histograms.pdf and ideel/ideel_uniprot_boxplots.pdf.

Additionally, diamond alignments are done between the predicted ORFs and the supplied reference proteins and ratios are plotted to score-assemblies-data/ideel/<reference>_ideel_histograms.pdf and score-assemblies-data/ideel/<reference>_ideel_boxplots.pdf.

bakta is only run when specified as extra config argument in the snakemake call:

snakemake -s /opt/software/score-assemblies/Snakefile --cores 20 --use-conda --config bakta=1

The bakta outfiles files are written to in the folder score-assemblies-data/bakta/<assembly>/.

NB: It takes a long time to download the bakta database and run bakta on all assemblies.

All measurements are summarized in a HTML page in score-assemblies-report.html.