Gatekeeper Pipeline

This pipeline should be used to perform check quality and trimming in sequence data with FASTP and filtering reads based on Kraken2 classification.

Requirements

• Nextflow (tested on version 23.04)
• dependencies in env.yml
• Kraken2 Standard dataset (tested with Standard-8 from 2023/06/05).
• For sylph: pre-built indexes and associated metadata.

Running Pipeline

Parameters:

• input_dir = A directory containing input fastq.gz files
• seq_platform = Sequencing platform illumina or ont
• kraken2_db_path = Path to the kraken2 dataset
• workflow = speciate, flu, or myco
• mem_map (Optional) = Use Kraken 2 memory mapping (requires preloading kraken2 db cp -r <path_to_kraken2_db> /dev/shm) true or false
• publish_dir (Optional) If set will output files to directory provided.
• sylph_db = Path to sylph database. See here
• sylph_db_metadata = Path to sylph database metadata file for taxonomic report

See parameters with:

nextflow run main.nf --help

nextflow run main.nf --workflow myco --input_dir ${INPUT_PATH} --kraken2_db_path ${KRAKEN2_DB_PATH} --seq_platform <illumina or ont>

A simple test for running locally can be done with

make run_local

The input directory can contain multiple samples. By default it will look for files in the input directory based on the following params:

• params.input_paired_suffix = "*_{1,2}.fastq.gz"
• params.input_single_suffix = "*.fastq.gz"

but these can be overriden. e.g.

nextflow run ... --input_paired_suffix "tb_sample*_{1,2}.fna.gz"

Memory mapping

You can put a kraken db into Ram disk (/dev/shm) and then use that with kraken using the --memory-mapping option to save on loading each time. To do this you may need to dedicate more ram space to /dev/shm

sudo mount -o remount,size=150G /dev/shm

Then

cp -r kraken_index /dev/shm/
kraken2 --db /dev/shm/kraken_index --memory-mapping ...

Sylph

Sylph is a k-mer based classification system that processes all the reads in one go, rather than classifying per read. It can accept multiple databases, so to incorporate fungi or viruses an extra database can be provided rather than needing to rebuild a single db.

Parameters

There are two parameters used in the Sylph process:

• sylph_query_ani_threshold: sets the minimum ani threshold from the query step. Sylph authors recommend 90.
• sylph_profile_ani_threshold: sets the minimum ani threshold from the profile step. Sylph authors recommend 95.

Databases for use

Preference is to go with NCBI over GTDB as it is more familiar to clinical users. Likewise with Virus the focus should be on human pathogenic viruses rather than phages etc.

See sylph_references for details of creating databases from RefSeq (NCBI).

To be used in competitive mapping downstream also require that there is a folder of the reference genomes to make a manifest from. As such that pre-built Sylph dbs are less useful except for GTDB which matches certain releases.

Data for tests

Sylph tests expect a folder sylph_dbs containing the sylph pre-built databases for testing. These are included in the knowledge bucket in the sylph directory. Specifically uses:

• gtdb-r220-c1000-dbv1.syldb and gtdb_r220_metadata.tsv.gz
• virus_imgvr_c200_v0.3.0.syldb and virus_imgvr_4.1_metadata.tsv.gz

The c1000 means that the k-mers are sampled every 1000bp making it smaller than c200, and so useful for testing.

Development/Install

The python side of this repo forms a python project. Due to a seqtk dependency it needs conda to install

conda env create -f env.yml
conda activate gatekeeper
pip install -e .[dev]
pytest

The nextflow processes require a container to run. By default will use the currently released version but can be overrided with a local version as follows:

# build local version
docker build -t test_container_gatekeeper .

# run using this container
nextflow run main.nf -profile local_docker ...
# or
nextflow run main.nf --test_container_gatekeeper test_container_gatekeeper ...

# can then test nextflow locally.
# NOTE: this will look for container with name "test_container_gatekeeper"
nf-test test tests/nextflow/*.nf.test --profile local_docker

Alternatively if all requirements are installed into your current environment can run using -profile no_docker

The kraken nextflow tests require the kraken index. Ensure you have a folder kraken2_db at the project level. The tests have been written to pass with Standard_8Gb index.

To run all tests can simply use:

make container
make test_local

FASTP default Parameters

Note on Adapter trimming

Fastp has bug with adapter detection on some ONT samples running out of memory. Have written a work around which waiting until the fastp issue is solved: OpenGene/fastp#607

Threshold values were chosen based on expert opinion.

--length_required Reads shorter than length_required will be discarded, default is 15. Illumina uses 50, ONT 100.

--low_complexity_filter enable low complexity filter. The complexity is defined as the percentage of base that is different from its next base (base[i] != base[i+1]). Requires at least 30% complexity.

--qualified_quality_phred & --unqualified_percent_limit A base is unqualified if below qualified_quality_phred, and a read can have at post unqualified_percent_limit proportion of bases unqualified. The defaults are Q15 and 40%. In ONT we use Q8.

Adapter Trimming is performed for Illumina but not ONT. guppy/dorado should perform adapter trimming already. There seems to be an issue with some ONT fastqs that crash during the adapter detection step.

per read cutting by quality (Illumina only)

Threshold values were chosen based on expert opinion. By default cut_mean_quality is Q20. And window size used is 4.

--cut_right move a sliding window from front to tail, if meet one window with mean quality < cut_mean_quality, drop the bases in the window and the right part, and then stop.

--cut_tail move a sliding window from tail (3') to front, drop the bases in the window if its mean quality is below cut_mean_quality, stop otherwise.

Data Dictionary

Example output

Parameter Name	Description	Units
unclassified	Number of reads which could not be classified by kraken2	Total single reads
classified	Number of reads which were classified by kraken2	Total single reads
Bacteria	Number of reads which were classified as Bacteria by kraken2	Total single reads
Mycobacteriaceae	Number of reads which were classified as Mycobacteriaceae by kraken2	Total single reads
Homo sapiens	Number of reads which were classified as Homo sapiens by kraken2	Total single reads

Taxa summary

A json summary file of main taxa, from family to species, is also created based on this schema. It has a threshold parameter (currently 0.5%). It will display all species with at least this proportion of reads according to kraken2 (or bracken if available). Can also be set to show the best species from any family with at least this proportion of reads.

A summary csv is also produced which should be easier to aggregate and analyse. It will included all species/genera found at the specified threshold.

Bracken

Bracken is used to give extra results to the taxa summary. Note that this is experimental and has some serious issues that should be resolved if decided that bracken is needed.

• Bracken is designed for short read data, and has a parameter for entering the read length. current 150 is used. It will still work if read lengths are different, but unclear what the effect is
• Bracken doesn't handle a mix of Species and Subspecies re-estimation. So currently set to only re-calculate to species level.
• Bracken will remove unassigned reads from percentage calculations, so can't be directly compared with kraken2 percentages without extra calculation.

Tags, Releases, and Committing

Use conventional commits. This is enforced with commitizen validate action and pre-commit hooks:

pre-commit install

This repo uses a standard gitflow approach, but with some changes to deal with docker containers in nextflow:

• There is a pyproject version which should be updated to match semantic version releases (from main/release branches)

• There is an active_version controlled by version_bumper which allows develop to have commit hash based versions.

• Every push to develop will cause an action to run bumper bump <commit-hash> --no-tag --active. This:

  • bumps the active_version in pyproject.toml
  • bumps the container tag used by nextflow processes
  • leaves the pyproject version as is

  Another workflow then builds and pushes the new container.
  Important: To deploy this you will need to use the hash of the bump commit, not the hash of the merge commit/container.

• In a release branch you can create a release candidate with bumper bump a.b.c-rcX. This also updates the pyproject version. Pushing the changes and new tag (automatically created) will trigger a build action.

• When release branch is ready for main run bumper bump a.b.c --no-tag. Push these changes to main and make a release there to build the container.