Genome Evaluation Pipeline (GEP)

• User-friendly and all-in-one quality control and evaluation pipeline for genome assemblies

• Run multiple genome evaluations in one go (as many as you want!)

• Seamlessly scaled to server, cluster, grid and cloud environments

• Required software stack automatically deployed to any execution environment using snakemake and conda

Getting Started

Step 1. Downloading the workflow

To clone the repository, use the following command:

git clone https://git.imp.fu-berlin.de/cmazzoni/GEP.git

Step 2. Conda management

• Conda (v4.10.3) but may work on older versions

If you already have conda installed on your system, please skip to step 3

Download the linux Miniconda3 installer from the following URL: https://docs.conda.io/en/latest/miniconda.html

Run the miniconda3 installation and check if it worked:

bash /<your_path_to>/Miniconda3-latest-Linux-x86_64.sh
##Follow miniconda3 installation instructions##

source ~/.bashrc

conda update conda

If conda command not found please close and re-open your terminal for conda installation to take effect, and then update.

Step 3. Creating our Snakemake conda environment

The pipeline requires the following software to run:

• snakemake (6.6.1+)
• python (3.9.1+)
• tabulate (0.8.7+)
• beautifulsoup4 (4.9+)
• mamba (0.15.2)

The easiest method to install this software stack is to create a GEP conda environment with the provided installGEP.yaml *Note

conda env create -f /<your_path_to>/GEP/installGEP.yaml

conda activate GEP

##check snakemake installed correctly

snakemake --version

Note If you already have a snakemake installation and would like to avoid installing it again, ensure that the above software are all in your PATH and you have conda installed/activated.

Step 4. Running the pipeline

BELOW NOT COMPLETE

GEP can be run in two modes:

1. Create meryl database

   • Input: Sample sheet outlining either Illumina PE or PacBio HiFi reads
   • Output: (.meryl) k-mer database

2. Run evaluation

   • Input: Sample sheet outlining k-mer database (.meryl) and corresponding assembly (.fa/fasta/fna)

Secondly, we will modify the config.yaml which will look like this:

##Configuration files pathways

Results: "/<PathTo>/desiredDestination/Results_AssemblyEval"

samplesTSV: "/<PathTo>/GEP/configuration/samples.tsv"

busco5Lineage: "vertebrata"

buscoMode: "genome"

1. Path to your Results This will be where all your results are stored for the run. It does not have to be inside the project folder (can be any location accessible by the user). You also do not need to create this folder yourself, it will be handled by snakemake.

2. Path to your samplesTSV. This is the path to the aforemention samples.tsv that was created/modified just above. For now, please keep this file inside the configuration folder, together with this config.yaml

3. busco5Lineage Busco needs a database to be able to run. Here you have a couple of different options.

• Manually download and unpack your desired database from https://busco-data.ezlab.org/v5/data/lineages/ . In this case (or if you already have the database downloaded to a specific location), you can provide the full path:

busco5Lineage: "/<PathTo>/manuallyDownloaded/vertebrata_odb10"

• Alternatively, you can just provide the taxonomy name that you wish to use. In this case, the latest database matching the name provided will be automatically downloaded prior to execution, if it doesn't already exist inside the buscoLineage directory. If it already exists in this buscoLineage either from manual download or from previously automatic download (from previously executed runs), then the pipeline will skip re-downloading.

busco5Lineage: "vertebrata"

4. You can change the busco mode, but considering the scope of this evaluation in it's current state, this option is rather redundant and will be removed/hidden.

Step 5. Running the workflow

If everything is set up correctly, we can run the pipeline very simply.

For now (though it should be a simple fix!), you must run the pipeline while being inside the project folder. In otherwords, you must be inside the folder where the file Snakefile is directly accessible.

Make sure your GEP environment is activated.

First you should run GEP in drymode:

snakemake -n

Which will check to see if some of your parameters/paths have been modified incorrectly. Further, it will install all the necessary environments to be utilised by the workflow, as well as download the busco5 database if it doesn't already exist. Unfortunaly when downloading the busco5 database, there will be lots of output in the terminal - a product of the limitations of the wget command used for downloading.

After the dry-run and downloading has complete, you can simply run the full pipeline with:

snakemake --cores # --use-conda && snakemake --report

Where --cores # is the maximum number of cores (synonomous with threads in this case) you want to utilise.

For example if you run snakemake with the command:

snakemake --cores 32 --use-conda && snakemake --report

Your pipeline will be executed with at most 32 threads per job/process. All jobs defined in the pipeline will use a percentage of these 32 threads allowing for automatic scalability of execution. An example being, if two jobs are ready to be executed and they both are defined as requiring 50% of the total cores avaliable, they can both be run in parralel.

Never use more threads than feasibly available by the architecture used to run the pipeline. Soon I will incorporate further modifications to provide even more scalability and portability (cluster/cloud execution, job scheduler compatability, and more!)

If you want to see some of the options you can use for execution, please see the snakemake help with snakemake -h

For example, you can assign a total amount of memory (e.g. 100GB) allowed by the pipeline by modifying the execution command to:

snakemake --cores 32 --use-conda --resources mem_mb=100000 && snakemake --report

INTRO

This pipeline allows users the ability to produce a wide range of commonly used evaluation metrics for genome assemblies, no matter your level of command-line experience/.

By harnessing the capabilities of snakemake, we present a workflow that incorporates a number of command-line tools and can be run on multiple independent genome assemblies in parallel. A streamlined user-experience is paramount to the devlopment process of this pipeline, as we strive for three key user-oriented components:

Snakemake will use conda to both install and manage our software packages and required tools. This helps to avoid software dependency conflicts, which will prevent the analysis from being simple to use and easily applied to different hardware. It also means that you, the user, do not have to be concerned with this at all - it is done for you!

#######################################################################

The software/tools used as part of our genome evaluation are as follows:

Pre-Processing (least biased short-read dataset available):

• Trimmomatic (Bolger, A. M., Lohse, M., & Usadel, B. (2014). http://www.usadellab.org/cms/?page=trimmomatic)
• Trim_galore (Felix Krueger bioinformatics.babraham.ac.uk)
• Fastqc (Simon Andrews https://github.com/s-andrews/FastQC
• Multiqc (Ewels, P., Magnusson, M., Lundin, S., Käller, M. (2016). https://doi.org/10.1093/bioinformatics/btw354)

Reference-free Genome Profiling

• GenomeScope2 (Ranallo-Benavidez, T.R., Jaron, K.S. & Schatz, M.C. (2020) https://github.com/tbenavi1/genomescope2.0)

K-mer distribution (copy-number spectra) analysis

• meryl (Rhie, A., Walenz, B.P., Koren, S. et al. (2020). https://doi.org/10.1186/s13059-020-02134-9)
• merqury (Rhie, A., Walenz, B.P., Koren, S. et al. (2020). https://doi.org/10.1186/s13059-020-02134-9)

Assessing quality and annotation completeness with Benchmarking Universal Single-Copy Orthologs (BUSCOs)

• BUSCOv4 (Seppey M., Manni M., Zdobnov E.M. (2019) https://busco.ezlab.org/ )

Scaffold/contig statistics: N# and L# stats, scaffold metrics, sequence counts, GC content, Estimated genome size

• Python scripts (Mike Trizna. assembly_stats 0.1.4 (Version 0.1.4). Zenodo. (2020). http://doi.org/10.5281/zenodo.3968775 )

#######################################################################

How to choose your illumina libraries

Variations in sequencing methods/protocols can lead to an increase in bias in the corresponding raw sequencing libraries. Sequencing a biological sample may often consist of both mate-pair/long-insert (e.g. insert sizes of 5k, 10k, 20k bp, etc.) and short-insert (e.g. insert-sizes 180, 250, 500, 800bp) paired-end libraries, respectively. Usually you can deduce the insert sizes and library types from the metadata found within NCBI or and SRA archive. In order to maintain a little bias as possible whilst maintaining decent coverage, you should ideally use only short-insert paired-end libraries for this evaluation pipeline.

If your library/s was sequenced using 10x barcodes (10X Genomics), you should remove the first 25-30bp of the forward read (R1) only. This will remove all barcode content.

Use trimmomatic

Will be incorporated automatically shortly

Reporting

Instead of or as well as retrieving the result files directly from the locations specified in the Results section (Step 6), the && snakemake --report argument used when running will create an interactive html report upon completion. This .html document will consist of all the relevant key files among other things such as the Directed Acyclic Graph (DAG) that snakemake uses to drive the order of execution, run-times of each individual step, and more (work in progress)

The report will be created in the main project directory, the same location as the Snakefile, where you executed the pipeline from.

Step 6. Results

ALL results can be found at the results directory defined by you in the config.yaml. Within this results folder, you will have a directory for each of the assemblies (assemblyName) you defined in the samples.tsv. The pipeline will produce a large amount of files at each step or for each tool, respectively. Taking this into consideration, the results can be considered in three tiers relative to their importance or ease of viewability.

Tier 3

The full results from all tools, including every file created during the execution of each tool, respectively. These results can be navigated at will, and are separated by their respective intended purposes (i.e. QVstats_merylAndMerqury, assemblystats, etc.)

Tier 2

The key result files (key plots, statistics, ) are aggregated and copied into a separate folder within the assembly folders. e.g./path/to/Results/SpeciesX/keyResults/.

The key result files are:

Multiqc report

• assemblyName_multiqc_report.html

QV Score

• assemblyName_merq.qv

Copy-Number Spectra plots

• assemblyName_merq.spectra-cn.fl.png
• assemblyName_merq.spectra-cn.st.png

Genomescope2 profile plots and summary

• assemblyName_gScope_log_plot.png
• assemblyName_gScope_linear_plot.png
• assemblyName_gScope_summary.txt

Assembly statistics (N#, NG#, L#, LG#, sequence counts, total bps, etc.)

• assemblyName_scaffold_stats.tsv
• assemblyName_contig_stats.tsv

BUSCOv4 results

• short_summary.specific.vertebrata_odb10.assemblyName.txt

Key values from all the above files pulled into aggregated table. Useful for quick glance

• assemblyName_aggregatedResults.tsv

Tier 1 There is a separately created folder within the main results directory (i.e. /path/to/Results/allAssemblies_keyResults )

Within this folder you will find a combined aggregate file (/path/to/Results/allAssemblies_keyResults/key_results.tsv, a tsv that combines the aforementioned key values from each assembly evaluated, respectively, into one single file. This is useful for plotting the key values across multiple assemblies.