Eco-Flow/nanoporemetabarcoding

Introduction

nf-core/nanoporemetabarcoding is a bioinformatics pipeline for processing nanopore metabarcoding data.

Overview:

Steps:

1. Filtering and trimming (NanoFilt)
2. Read QC (NanoPlot)
   • Ran on both raw and filtered and trimmed reads
3. Tags+primer based demultiplexing and trimming (Cutadapt). Divided in two steps:
   1. First, demultiplexing based on the forward tags-primers
   2. Then, demultiplexing based on matching reverse tags-primers (combination of forward and reverse tags-primers)
4. Group amplicons reads into “species” (consensus sequences) (amplicon_sorter)
5. Consensus sequence correction (Medaka)
   • Correction using the consensus sequence from aplicon_sorter as reference and the grouped amplicon reads as the basecalled data
6. Create custom BLAST database (makeblastdb)
7. Consensus sequence annotation based on database (blastn)
   • Annotation is based on the best blast hit per consensus. And best blast hit is based on:
     1. First on the bit score
     2. Second on the e-value
8. Assign taxonomy to blast hits using taxonomizr (taxonomizr). Only works with NCBI accessions (GenBank and RefSeq). If an ASV has multiple hits with the same top bitscore, e-value, and percent identity, the lowest common taxonomic rank across all hits is assigned.

Usage

[!NOTE] If you are new to Nextflow and nf-core, please refer to this page on how to set-up Nextflow. Make sure to test your setup with -profile test before running the workflow on actual data.

To use nanoporemetabarcoding, first clone this repo:

git clone https://github.com/Eco-Flow/nanoporemetabarcoding.git

Before running nanoporemetabarcoding on your data, you can test if the pipeline is suited to your setup by running:

nextflow run main.nf \
   -profile test,<docker/singularity/conda/.../institute> \
   --outdir <OUTDIR>

The running time should be around 8 minutes on a machine with 12 cpus/threads.

For running the pipleine on your data, prepare a samplesheet with your input FASTQs that looks as follows:

id,fastq
fastq_id_1,path/to/reads1.fastq.gz
fastq_id_2,path/to/reads2.fastq.gz
...

The first column represents the FASTQ id (e.g. plate number if the FASTQ represents a single plate), and the second the location to the FASTQ file.

Now, you can run the pipeline using:

nextflow run main.nf \
   -profile <docker/singularity/.../institute> \
   --input samplesheet.csv \
   --outdir <OUTDIR>

[!WARNING] Please provide pipeline parameters via the CLI or Nextflow -params-file option. Custom config files including those provided by the -c Nextflow option can be used to provide any configuration except for parameters; see docs.

Parameters

Check the config file nextflow.config in the params block and fill in the parameters with the proper values:

Input options:

You can optionally provide a metadata file in CSV format with the primer-tag combination in the first column and the corresponding sample name in the second column:

    // Metadata with primer combinations
    metadata                   = path/to/metadata.csv

The strucutre of the CSV should be as follows:

id,primer_comb,sample
fastq_id_1,<forward-primer-tag-id>_<forward-primer-tag-id>,sample1
fastq_id_2,<forward-primer-tag-id>_<forward-primer-tag-id>,sample2
...

[!NOTE] If and/or in the metadata samplesheet don't match the ids in the primer-tag FASTA files, demultiplexing won't work. Make sure they have the exact same name.

The id of the metadata should match the id of the samplesheet (see usage), and values for the primer_comb and sample fields should be unique within id groups. See ./test_data/metadata.csv for an example. If set to null (value by default), the pipeline will use the primer-tag combination as the sample ID in the final ASV table.

Demultiplex options:

1. List of forward primer-tag combinations in FASTA format:

    // Cutadapt options
    tags_f                      = 'path/to/forward/primer-tag.fasta' // List of forward primer-tag combinations in fasta format

1. List of reverse primer-tag combinations in FASTA format:

    tags_r                      = 'path/to/reverse/primer-tag.fasta' // List of forward primer-tag combinations in fasta format

1. Set the error rate for adapter removal. This can be a value between 0 and 1 (1 not included) if a maximum error rate wants to be applied to all adapters, or it can be equal or greater than 1, in which case it will be converted to a maximum error rate depending on the adapter length. Check cutadatpt documentation for more information:

    error_rate                  = 2 // Error rate for adapter removal

1. Primer-tag overlap. To calculate mismatches with the --error-rate option, cutadapt considers only the aligned region (overlap) between primer and read. The minimum length of this overlap is set by the --overlap option. For example, with an error rate of 0.2 (20%) and a 20 bp primer, if the required overlap is 10 bp and there are 3 mismatches in that region, the match is rejected because the mismatch rate is 3/10 = 30%, which is above the 20% threshold. But if the overlap is 20 bp and there are 3 mismatches, the match is accepted because the mismatch rate is 3/20 = 15%, which is below the 20% threshold. If your primers are very similar, it is often best to set --overlap to at least the length of the longest primer to avoid spurious matches.

Nanofilt options:

Nanofilt is run before demultiplexing and tag+primer trimming:

    // Nanofilt options
    nano_quality                = null // Minimun read quality (phred score)
    nano_read_length            = 250 // Minimum reads length

1. Filter FASTQs with a mumber of reads equal or lower than (changing this value is not recommended):

    filt_fastq                  = 0

Check ./test_data/primers_f.fasta and ./test_data/primers_r.fasta for examples of tags_f and tags_r files, respectively.

Blast options:

1. Path to an already built blast database:

    blast_db                   = 'path/to/blast.db' // Path to already built database

1. Or build a local BLAST database from a collection of sequences:

    custom_db                  = 'path/to/local/database.fasta' // Path to database to be built

1. The E value describes the number of one can “expect” to see by chance when searching a database of a particular size. The lower the value, the more significant a match is:

    evalue                     = 0.001 // evalue cutoff

1. Maximum number of hits per query (consensus sequence):

    max_target_seqs            = 5 // Maximum number of hits per query

1. Maximum number of high-scoring segment pairs (HSPs) per subject (sequence in the database):

    max_hsps                   = null // Maximum number of HSPs per subject (subjects in the database)

1. Minimum query (consensus sequence) coverage per HSP:

    qcov_hsp_perc              = 90 // Minimum query coverage per HSP

For more details, check the appendix section of the blast documentation.

Assign taxonomy options:

1. Assign taxonomy at different levels based on an identity threshold:

    // Assign taxonomy options
    spident                     = null // Identity threshold (in %) for taxonomy assignment at species level. If not set, the value by default will be 99. Set to > 100 if not wanting to assign species
    gpident                     = null // Identity threshold (in %) for taxonomy assignment at genus level. If not set, the value by default will be 90
    fpident                     = null // Identity threshold (in %) for taxonomy assignment at family level. If not set, the value by default will be 80
    opident                     = null // Identity threshold (in %) for taxonomy assignment at order level. If not set, the value by default will be 70

1. Path to taxonomizr SQL database:

    sql_db                      = 'path/to/taxonomizr/database.sqlite' // Path to the taxonomizr SQLite database

Command line arguments:

Alternatively, instead of modifying the nextflow.config file, you can provide/modify pipeline options by specifying them as command line arguments. For example:

nextflow run main.nf \
   -profile <docker/singularity/.../institute> \
   --input samplesheet.csv \
   --outdir <OUTDIR>
   --nano_quality 10
   --evalue 1e-10
   --tags_f path/to/forward/tags+primers.fasta
   --tags_r path/to/reverse/tags+primers.fasta
   ...

Pipeline output

The files listed below will be created in the results directory (set by --outdir) after the pipeline has finished. All paths are relative to the top-level results directory.

Quality reports

Nanoplot

NanoPlot is a tool that outputs QC reports for raw, filtered -after nanofilt-, and demultiplexed reads. Nanoplot reports are store inside the nanoplot folder, and come in different formats.

Multiqc

MultiQC is a visualization tool that generates a single HTML report summarising all of the NanoPlot QC results for raw, filtered, demultiplexed reads. Most of the pipeline QC results are visualised in the report and further statistics are available in the report data directory.

Taxonomy assignment

Each consensus sequence (ASV) is assigned a taxon using a two-step approach:

1. Rank assignment by percent identity: each BLAST hit is assigned to a taxonomic rank based on percent identity thresholds (species, genus, family, order) (see options Assign taxonomy options: in the Parameters section).
2. Last Common Rank (LCR) consensus: because each ASV can have multiple BLAST hits, the final taxon is resolved by finding the most specific rank at which all hits agree. If hits disagree at the assigned rank, the pipeline falls back to progressively coarser ranks (genus → family → order → class → phylum) until a consensus is reached. If no consensus can be found, the ASV is labelled Unassigned.

The final output is ASV_taxa_final.csv, one row per ASV, with the following columns:

• ASV: ASV identifier.
• sample_name: Sample name, derived from the combination of FASTQ id and sample id.
• read_count: Number of reads assigned to the ASV.
• pident: Mean percent identity of the best BLAST alignment.
• length: Mean alignment length.
• mismatch: Mean number of mismatched positions in the alignment (excluding gaps).
• evalue: Expect value — the lower, the more significant the match. Depends on database size.
• bitscore: Bit score — the higher, the better the alignment. Independent of database size.
• taxaId: NCBI taxonomy ID of the matched sequence.
• phylum, class, order, family, genus, species: Full taxonomic lineage of the matched sequence.
• Resolved.taxon: Final consensus taxon assigned to the ASV. A trailing * indicates that one or more BLAST hits for the same ASV had no taxonomy in the database and were excluded from the LCR resolution. This ASV can be classified as Unasigned instead for a more conservative approach.

ASV_table_pre-assigned.csv contains the per-hit taxonomy assignments before LCR consensus resolution — one row per BLAST hit rather than one row per ASV.

Credits

nf-core/nanoporemetabarcoding was originally written by Fernando Duarte, Chris Wyatt.

We thank the following people for their extensive assistance in the development of this pipeline:

Contributions and Support

If you would like to contribute to this pipeline, please see the contributing guidelines.

For further information or help, don’t hesitate to get in touch on the Slack #nanoporemetabarcoding channel (you can join with this invite).

Citations

An extensive list of references for the tools used by the pipeline can be found in the CITATIONS.md file.

You can cite the nf-core publication as follows:

The nf-core framework for community-curated bioinformatics pipelines.

Philip Ewels, Alexander Peltzer, Sven Fillinger, Harshil Patel, Johannes Alneberg, Andreas Wilm, Maxime Ulysse Garcia, Paolo Di Tommaso & Sven Nahnsen.

Nat Biotechnol. 2020 Feb 13. doi: 10.1038/s41587-020-0439-x.