Natrix is an open-source bioinformatics pipeline for the preprocessing of raw sequencing data. The need for a scalable, reproducible workflow for the processing of environmental amplicon data led to the development of Natrix. It is divided into quality assessment, read assembly, dereplication, chimera detection, split-sample merging, ASV or OTU-generation and taxonomic assessment. The pipeline is written in Snakemake (Köster and Rahmann 2018), a workflow management engine for the development of data analysis workflows. Snakemake ensures reproducibility of a workflow by automatically deploying dependencies of workflow steps (rules) and scales seamlessly to different computing environments like servers, computer clusters or cloud services. While Natrix was only tested with 16S and 18S amplicon data, it should also work for other kinds of sequencing data. The pipeline contains seperate rules for each step of the pipeline and each rule that has additional dependencies has a seperate conda environment that will be automatically created when starting the pipeline for the first time. The encapsulation of rules and their dependencies allows for hassle-free sharing of rules between workflows.

DAG of an example workflow: each node represents a rule instance to be executed. The direction of each edge represents the order in which the rules are executed, which dashed lines showing rules that are exclusive to the OTU version and dotted lines rules exclusive to the ASV variant of the workflow. Disjoint paths in the DAG can be executed in parallel. Below is a schematic representation of the main steps of the pipeline, the color coding represents which rules belong to which main step.

If you use Natrix, please cite: Welzel, M., Lange, A., Heider, D. et al. Natrix: a Snakemake-based workflow for processing, clustering, and taxonomically assigning amplicon sequencing reads. BMC Bioinformatics 21, 526 (2020). https://doi.org/10.1186/s12859-020-03852-4

Dependencies

• Conda

• GNU screen (optional)

Conda can be downloaded as part of the Anaconda or the Miniconda plattforms (Python 3.7). We recommend to install miniconda3. Using Linux you can get it with:

$ wget https://repo.anaconda.com/miniconda/Miniconda3-latest-Linux-x86_64.sh
$ bash Miniconda3-latest-Linux-x86_64.sh

GNU screen can be found in the repositories of most Linux distributions:

• Debian / Ubuntu based: apt-get install screen

• RHEL based: yum install screen

• Arch based: pacman -S screen

All other dependencies will be automatically installed using conda environments and can be found in the corresponding environment.yaml files in the envs folder and the natrix.yaml file in the root directory of the pipeline.

Getting Started

To install Natrix, you'll need the open-source package management system conda and, if you want to try Natrix using the accompanying pipeline.sh script you'll need GNU screen . After cloning this repository to a folder of your choice, it is recommended to create a general natrix conda environment with the accompanying natrix.yaml . In the main folder of the cloned repository, execute the following command:

$ conda env create -f natrix.yaml

This will create a conda environment containing all dependencies for Snakemake itself.

With Natrix comes an example primertable example_data.csv , configfile example_data.yaml and an example amplicon dataset in the folder example_data .

To try out Natrix using the example data, type in the following command:

$ ./pipeline.sh
Please enter the name of the project
$ example_data

The pipeline will then start a screen session using the project name (here, example_data ) as session name and will beginn downloading dependencies for the rules. To detach from the screen session, press Ctrl+a, d ( first press Ctrl+a and then d ). To reattach to a running screen session, type in:

$ screen -r

When the workflow has finished, you can press Ctrl+a, k ( first press Ctrl+a and then k ). This will end the screen session and any processes that are still running.

Tutorial

Prerequisites: dataset, primertable and configuration file

The FASTQ files need to follow a specific naming convention:

samplename_unit_direction.fastq.gz

with:

• samplename as the name of the sample, without special characters.

• unit , identifier for split-samples ( A , B ). If the split-sample approach is not used, the unit identifier is simply A .

• direction , identifier for forward ( R1 ) and reverse ( R2 ) reads of the same sample. If the reads are single-end, the direction identifier is R1 .

A dataset should look like this (two samples, paired-end, no split-sample approach):

S2016RU_A_R1.fastq.gz
S2016RU_A_R2.fastq.gz
S2016BY_A_R1.fastq.gz
S2016BY_A_R2.fastq.gz

Besides the FASTQ data from the sequencing process Natrix needs a primertable containing the sample names and, if they exists in the data, the length of the poly-N tails, the sequence of the primers and the barcodes used for each sample and direction. Besides the sample names all other information can be omitted if the data was already preprocessed or did not contain the corresponding subsequence. Natrix also needs a configuration file in YAML format, specifying parameter values for tools used in the pipeline.

The primertable, configfile and the folder containing the FASTQ files all have to be in the root directory of the pipeline and have the same name (with their corresponding file extensions, so project .yaml, project .csv and the project folder containing the FASTQ files). The first configfile entry ( filename ) also needs to be the name of the project.

Running Natrix with the pipeline.sh script

IF everything is configured correctly, you can start the pipeline by typing in the following commands into your terminal emulator:

$ ./pipeline.sh
Please enter the name of the project
$ example_data

The pipeline will then start a screen session using the project name as session name and will beginn downloading dependencies for the rules. To detach from the screen session, press Ctrl+a, d ( first press Ctrl+a and then d ). To reattach to a running screen session, type in:

$ screen -r

When the workflow has finished, you can press Ctrl+a, k ( first press Ctrl+a and then k ). This will end the screen session and any processes that are still running.

Running Natrix with Docker or docker-compose

Pulling the image from Dockerhub

Natrix can be run inside a Docker-container. Therefore, Docker has to be installed. Please have a look at the Docker website to find out how to install Docker and set up an environment if you have not used it before.

The easiest way to run the docker container is to download the pre-build container from dockerhub .

$ docker pull mw55/natrix

The docker container has all environments pre-installed, eliminating the need for downloading the environments during first-time initialization of the workflow. To connect to the shell inside the docker container, input the following command:

docker run -it --label natrix_container -v </host/database>:/app/database -v </host/results>:/app/results -v </host/input_folder>:/app/input -v </host/demultiplexed>:/app/demultiplexed mw55/natrix bash

/host/database is the full path to a local folder, in which you wish to install the database (SILVA or NCBI). This part is optional and only needed if you want to use BLAST for taxonomic assignment.

/host/results is the full path to a local folder in which the results of the workflow should be stored for the container to use.

/host/input_folder is the full path to a local folder in which the input (the project folder, the project .yaml and the project .csv) should be saved.

/host/demultiplexed is the full path to a local folder in which the demultiplexed data, or, if demultiplexing is turned off, the input data will be saved.

After you connected to the container shell, you can follow the running Natrix manually tutorial.

Directly starting the workflow using docker-compose

Alternatively, you can start the workflow using the the docker-compose command in the root directory of the workflow (it will pull the latest natrix image from DockerHub):

$ PROJECT_NAME="<project>" docker-compose up (-d)

with project being the name of your project. e.g.:

$ PROJECT_NAME="example_data" docker-compose up # sudo might be needed

all output folders will be available at /srv/docker/natrix/ make sure to copy your project folder, project .yaml and project .csv files to /srv/docker/natrix/input/ or create a new volume-mapping using the docker-compose.yml file. By default the container will wait until the input files exist. At first launch the container will download the required databases to /srv/docker/natrix/databases/, this process might take a while.

Building the container yourself

If you prefer to build the docker container yourself from the repository (for example, if you modified the source code of Natrix) the container can be build and started directly: ( host folders have to be changed!)

docker build . --tag natrix
docker run -it --label natrix_container -v </host/database>:/app/database -v </host/results>:/app/results -v </host/input_folder>:/app/input -v </host/demultiplexed>:/app/demultiplexed natrix bash # -v /host/database:/app/database is optional

You will then be at the command prompt inside the docker container, from there you can follow the tutorial for running Natrix manually .

Running Natrix manually

If you prefer to run the preperation script and snakemake manually, you have to start by activating the snakemake environment:

$ conda activate natrix

Followed by running the preperation script, with project being the name of your project:

$ python3 create_dataframe.py <project>.yaml

This command will create the units.tsv file, containing the file information in a way that Natrix can use it.

To start the main pipeline, type in:

snakemake --use-conda --configfile <project>.yaml --cores <cores>

with project being the name of your project and cores being the amount of cores you want to allocate for Natrix to use.

Should the pipeline prematurely terminate (either because of an error or by deliberatly stopping it) running the command above again will start the pipeline from the point it was terminated.

Cluster execution

Natrix can be run easily on a cluster system using either conda or the docker container. Adding --cluster to the start command of Natrix, together with a command to submit jobs (e. g. qsub) is enough for most cluster computing environments. An example command would be:

$ snakemake -s <path/to/Snakefile> --use-conda --configfile <path/to/configfile.yaml> --cluster "qsub -N <project name> -S /bin/bash/" --jobs 100

Further qsub arguments including brief explanations can be found under qsub arguments . For additional commands that should be executed for each job the argument --jobscript path/to/jobscript.sh can be used. A simple jobscript that sources before the execution of each job the bashrc and activates the snakemake environment looks like this:

#!/usr/bin/env bash
source ~/.bashrc
conda activate natrix
{exec_job}

Instead of directly passing the cluster submission arguments to the snakemake command it is also possible to write a profile that contains cluster commands and arguments. The use of profiles allows the assignment of rule-specific hardware requirements. For example, the BLAST rule benefits from a large amount of CPU cores, while other rules, like AmpliconDuo, do not. With profiles, rules could be assigned a low amount of CPU cores as default, with rules like BLAST being assigned a larger amount of cores. This allows optimal usage of cluster resources and shorter waiting times. The creation of profiles is largely dependant on the software and hardware available on the cluster. With a profile Natrix can simply be run with

$ snakemake -s <path/to/Snakefile> --profile myprofile 

The Snakemake documentation contains a tutorial for profile creation and the Snakemake profiles GitHub page contains example profiles for different cluster systems.

Output

After the workflow is finished, the original data can be found under Natrix-Pipeline/demultiplexed/ , while files created during the workflow can be found under Natrix-Pipeline/results/ .

Output file hierarchy, blue nodes represent folders, orange nodes represent files that are created in both variants of the workflow, green nodes are files exclusive to the OTU variant and purple nodes are files exclusive to the ASV variant of the workflow.

Steps of the Pipeline

Initial demultiplexing

The sorting of reads according to their barcode is known as demultiplexing.

Quality control

For quality control the pipeline uses the programs FastQC (Andrews 2010), MultiQC (Ewels et al. 2016) and PRINSEQ (Schmieder and Edwards 2011).

FastQC

FastQC generates a quality report for each FASTQ file, containing information such as the per base and average sequence quality (using the Phred quality score), overrepresented sequences, GC content, adapter and the k-mer content of the FASTQ file.

MultiQC

MultiQC aggregates the FastQC reports for a given set of FASTQ files into a single report, allowing reviews of all FASTQ files at once.

PRINSEQ

PRINSEQ is used to filter out sequences with an average quality score below the threshold that can be defined in the configuration file of the pipeline.

Read assembly

Define primer

The define_primer rule specifies the subsequences to be removed by the assembly rule, specified by entries of the configuration file and a primer table that contains information about the primer and barcode sequences used and the length of the poly-N subsequences. Besides removing the subsequences based on their nucleotide sequence, it is also possible to remove them based solely on their length using an offset. Using an offset can be useful if the sequence has many uncalled bases in the primer region, which could otherwise hinder matches between the target sequence defined in the primer table and the sequence read.

Assembly & removal of undesired subsequences (OTU-variant)

To assemble paired-end reads and remove the subsequences described in the pre- vious section PANDAseq (Masella et al. 2012) is used, which uses a probabilistic error correction to assemble overlapping forward- and reverse-reads. After assembly and sequence trimming, it will remove sequences that do not meet a minimal or maximal length threshold, have an assembly quality score below a threshold that can be defined in the configuration file and sequences whose forward- and reverse-read do not have an sufficiently long overlap. The thresholds for each of these procedures can be adjusted in the configuration file. If the reads are single-end, the subsequences (poly-N, barcode and the primer) that are defined in the define_primer rule are removed, followed by the removal of sequences that do not meet a minimal or maximal length threshold as defined in the configuration file.

Removal of undesired subsequences (ASV-variant)

In the ASV variant of the workflow Cutadapt (Martin 2011) is used to remove the undesired subsequences defined in the primer table.

ASV denoising (ASV-Variant)

After the removal of undesired subsequences ASVs are generated using the DADA2 (Callahan et al. 2016) algorithm. It dereplicates the dataset and uses a denoising algorithm. This algorithm infers if a sequence was produced by a different sequence based on the composition, quality, and abundance of a sequence and an Illumina error model. After ASV generation, exactly overlapping forward and reverse reads are assembled. The assembled ASVs are saved as FASTA files for downstream analysis.

Similarity clustering

Conversion of FASTQ files to FASTA files (OTU-variant)

The rule copy_to_fasta converts the FASTQ files to FASTA files to reduce the disc space occupied by the files and to allow the usage of CD-HIT, which requires FASTA formated sequencing files. Clustering of similar sequences The CD-HIT-EST algorithm (Fu et al. 2012) clusters sequences together if they are either identical or if a sequence is a subsequence of another. This clustering approach is known as dereplication. Beginning with the longest sequence of the dataset as the first representative sequence, it iterates through the dataset in order of decreasing sequence length, comparing at each iteration the current query sequence to all representative sequences. If the sequence identity threshold defined in the configuration file is met for a representative sequence, the counter of the representative sequence is increased by one. If the threshold could not be met for any of the existing representative sequences, the query sequence is added to the pool of representative sequences. Cluster sorting The cluster_sorting rule uses the output of the cdhit rule to count the amount of sequences represented by each cluster, followed by sorting the representative sequences in descending order according to the cluster size and adds a specific header to each sequence as required by the UCHIME chimera detection algorithm.

Chimera detection

VSEARCH

VSEARCH is a open-source alternative to the USEARCH toolkit, which aims to functionally replicate the algorithms used by USEARCH for which the source code is not openly available and which are often only rudimentarily described (Rognes et al. 2016). Natrix uses as an alternative to UCHIME the VSEARCH uchime3_denovo algorithm to detect chimeric sequences (further referred to as VSEARCH3). The VSEARCH3 algorithm is a replication of the UCHIME2 algorithm with optimized standard parameters. The UCHIME2 algorithm is described by R. Edgar 2016 as follows:

"Given a query sequence Q , UCHIME2 uses the UCHIME algorithm to construct a model ( M ), then makes a multiple alignment of Q with the model and top hit ( T , the most similar reference sequence). The following metrics are calculated from the alignment: number of differences dQT between Q and T and dQM between Q and M , the alignment score ( H ) using eq. 2 in R. C. Edgar et al. 2011. The fractional divergence with respect to the top hit is calculated as divT = (dQT − dQM)/|Q|. If divT is large, the model is a much better match than the top hit and the query is more likely to be chimeric, and conversely if divT is small, the model is more likely to be a fake."

The difference between the UCHIME2 and UCHIME3 algorithm is that to be selected as a potential parent, a sequence needs to have at least 16 times the abundance of the query sequence in the UCHIME3 algorithm, while it only needs double the abundance of the query sequence to be selected as a potential parent in the UCHIME2 algorithm.

Table creation and filtering

Merging of all FASTA files into a single table

For further processing the rule unfiltered_table merges all FASTA files into a single, nested dictionary, containing each sequence as the key with another dictionary as value, whose keys are all (split -) samples in which the sequence occurred in and as values the abundance of the sequence in the particular (split - ) sample. For further pipeline processing the dictionary is temporarily saved in JSON format. To ease statistical analysis of the data the dictionary is also exported as a comma separated table. Filtering In the filtering rule of the pipeline all sequences that do not occur in both splitsamples of at least one sample are filtered out. For single-sample data, the filtering rule uses an abundance cutoff value that can be specified in the configuration file to filter out all sequences which have abundances less or equal the specified cutoff value. The filtered data and the filtered out data is subsequently exported as comma separated tables.

Table conversion to FASTA files

As the swarm rule needs FASTA files as input, the resulting table of the filtering is converted to a FASTA file by the rule write_fasta.

AmpliconDuo / Split-sample approach

The pipeline supports both single-sample and split-sample FASTQ amplicon data. The split-sample protocol (Lange et al. 2015) aims to reduce the amount of sequences that are the result of PCR or sequencing errors without the usage of stringent abundance cutoffs, which often lead to the loss of rare but naturally occurring sequences. To achieve this, extracted DNA from a single sample is divided into two split-samples, which are then separately amplified and sequenced. All sequences that do not occur in both split-samples are seen as erroneous sequences and filtered out. The method is therefore based on the idea that a sequence that was created by PCR or sequencing errors does not occur in both samples. A schematic representation of the split-sample method is shown below:

Schematic representation of the split-sample approach: Extracted DNA from a single environmental sample is splitted and separately amplified and sequenced. The filtering rule compares the resulting read sets between two split-samples, filtering out all sequences that do not occur in both split-samples. Image adapted from Lange et al. 2015.

The initial proposal for the split-sample approach by Dr. Lange was accompanied by the release of the R package AmpliconDuo for the statistical analysis of amplicon data produced by the aforementioned split-sample approach. It uses Fishers exact test to detect significantly deviating read numbers between two experimental branches A and B from the sample S. To measure the discordance between two branches of a sample the read-weighted discordance ∆rSθ , which is weigthed by the average read number of sequence i in both branches and the unweighted discordance ∆uSθ for each sequence i are calculated. If ∆uSθ = 0 each branch of sample S contains the same set of sequences, while if ∆rSθ = 0 the read numbers for each sequence in sample S are within the error margin set by the chosen false discovery rate. The results of the discordance calculations are then plotted for visualization purposes and written to an R data file to allow the filtering of significantly deviating sequences.

OTU generation (OTU-variant)

OTUs are generated using the Swarm clustering algorithm (Mahé et al. 2015) in the identically named rule. Swarm clusters sequences into OTUs using an iterative approach with a local threshold: It creates a pool of amplicons from the input file and a empty OTU. Subsequently, it will remove the first amplicon from the pool, which will become the OTU seed. All amplicons left in the pool that differ in their nucleotide composition from the initial seed by a user given threshold (the default threshold used is 1 nucleotide) are removed from the pool and added to the OTU as subseeds. In the next iteration, each amplicon having at most a difference as high as the threshold to any of the subseeds is then removed from the pool and added to the OTU. This iterative process will continue until there are no amplicons left in the pool with a nucleotide difference of at most the threshold to any of the subseeds added in the previous iteration to the OTU, leading to the closure of the OTU and the opening of a new one. This approach to OTU generation circumvents two sources of OTU variability that are inherent to greedy clustering algorithms: the input order dependency, in which the first amplicon in a FASTA file will become the centroid of an OTU and the use of a global threshold, recruiting all amplicons that have less differences to the centroid than a user defined threshold. The iterative approach of Swarm produces a star-shaped minimum spanning tree, with an (usually highly abundant) amplicon as the center, regardless of the chosen first amplicon of the OTU, as visualized in Figure 12. The sequence of the amplicon at the center of each OTU tree is used in subsequent analysis steps as the representative sequence of the corresponding OTU.

Schematic representation of the greedy clustering approach and the iterative Swarm approach. The greedy approach (a) that uses a global clustering threshold t and input order dependent centroid selection can lead to the placement of closely related amplicons into different OTUs. The iterative approach of Swarm (b), with a local threshold d, leads to OTUs containing only closely related amplicons with a natually forming centroid during the iterative growth process of the OTU. Image from Mahé et al. 2015.

Sequence comparison

The assignment of taxonomic information to an OTU/ASV is an important part of the processing of environmental amplicon data, as the characteristics of known groups or species can be used to assess the environmental conditions of the samples origin. To find sequences that are similar to the representative sequence of each OTU/ASV the BLAST (basic local alignment search tool) algorithm (Altschul et al. 1990) is used to search for similar sequences in the SILVA database (Pruesse et al. 2007). The SILVA database contains aligned rRNA sequencing data that are curated in a multi-step process. While it has an extensive collection of highquality prokaryotic rRNA sequencing data, it only contains a limited amount of microbial eukaryotic sequencing data. If the database is not locally available, the required files will automatically be downloaded and the database will be build in the rule make_silva_db. The BLAST algorithm itself will be executed in the blast rule. As the aim is to find similar nucleotide sequences in the database for each query sequence the nucleotide-nucleotide BLAST (BLASTn) variation of the BLAST algorithm is used. The tab separated output file of the blast rule contains the following information for each representative sequence if the output of the BLASTn search meets the criteria defined in the configuration file:

Merging of the results

The output data from the write_fasta, swarm and blast rules are merged to a single comma separated table in the merge_results rule, containing for each representative sequence the sequence identification number, the nucleotide sequence, the abundance of the sequence in each sample, the sum of abundances and, if the blast rule found a similar sequence, all information found in the table shown in the previous section.

Example primertable

The primertable should be a .csv file ( project .csv) in the following format:

Configfile

Below are the explainations for the configfile ( project .yaml) entries:

References

• Köster, Johannes and Sven Rahmann (2018). “Snakemake—a scalable bioinformatics workflow engine”. In: Bioinformatics.

• Andrews, S. (2010). FASTQC. A quality control tool for high throughput sequence data.

• Ewels, P et al. (2016). “MultiQC: summarize analysis results for multiple tools and samples in a single report”. In: Bioinformatics 32.19, pp. 3047–3048.

• Callahan, B. J. et al. (2016). “DADA2: High-resolution sample inference from illumina amplicon data.“ In: Nature Methods ,13(7):581–583.

• Schmieder, Robert and Robert A. Edwards (2011). “Quality control and preprocessing of metagenomic datasets.” In: Bioinformatics 27.6, pp. 863–864.

• Masella, Andre P et al. (2012). “PANDAseq: paired-end assembler for illumina sequences”. In: BMC Bioinformatics 13.1 , p. 31.

• Fu, Limin et al. (2012). “CD-HIT: accelerated for clustering the next-generation sequencing data”. In: Bioinformatics 28.23 , pp. 3150–3152.

• Martin, M. (2011). “Cutadapt removes adapter sequences from high-throughput sequencing reads.“ In: EMB- net.journal, 17(1):10.

• Rognes, Torbjørn et al. (2016). “VSEARCH: a versatile open source tool for metagenomics”. In: doi: 10.7287/peerj.preprints.2409v1.

• Edgar, Robert (2016). “UCHIME2: improved chimera prediction for amplicon sequencing”. In: bioRxiv.

• Mahé, Frédéric et al. (2015). “Swarm v2: highly-scalable and high-resolution amplicon clustering”. In: PeerJ 3.

• Lange, Anja et al. (2015). “AmpliconDuo: A Split-Sample Filtering Protocol for High-Throughput Amplicon Sequencing of Microbial Communities”. In: Plos One 10.11.

• Altschul, Stephen F. et al. (1990). “Basic local alignment search tool”. In: Journal of Molecular Biology 215.3 , pp. 403–410.

• Pruesse, E. et al. (2007). “SILVA: a comprehensive online resource for quality checked and aligned ribosomal RNA sequence data compatible with ARB”. In: Nucleic Acids Research 35.21 , pp. 7188–7196.