Nanopore barcoding analysis

Run name: its-run1

Date: 12 December, 2025

Analyzed with DadaNanoBC, version: v0.2.0

1 The workflow (in short)

1. Search the primers and short sample indexes (located up-/downstream of the primer sequences)
2. Infer the barcode sequences with DADA2 and/or fixed-threshold clustering (for low-coverage variants), align and infer the consensus
3. Auto-assign the taxonomy and compare with the name given by morphological identification (if present) to validate and check for contamination.
4. Compare with already known sequences (if present)
5. Export summary table, which can be further manually curated

See documentation for more details.

2 Configuration

The following settings are used (contents of analysis/its-run1/config.yaml):

cores: 4

demultiplex:
  primer_max_err: 0.2
  idx_max_diffs: 0
  min_read_length: 200
  error_threshold: 2.5

cluster:
  max_sample_depth: 2000

taxonomy:
    # UNITE database (https://doi.org/10.15156/BIO/3301243, all Eukaryotes, QIIME-formatted)
    db_url: https://s3.hpc.ut.ee/plutof-public/original/0ee4ff53-bd04-454b-a937-af1f466bfe8b.tgz
    db_type: unite
    known_contaminants:
      family:
        - Aspergillaceae
        - Capnodiaceae
      genus:
        - Vishniacozyma
    kingdom: Fungi

3 Software

Required R packages are automatically installed.

Other software needs to be installed and accessible either in the system PATH or in the bin subdirectory.

4 Input files

4.1 Sequencing data

Base-called sequences need to be present in analysis/<run>/reads.fastq.gz.

4.2 Samples

The run has 1 amplicon(s) with the following sample numbers:

amplicon	# samples	# with morphological identification	# with known sequence
ITS5_ITS4	1344	1333	54

Note: primers are searched in the order that amplicons appear in the sample sheet. With nested amplicons, the shorter one should be placed at the end.

4.3 Primers

The primers are defined in the primers sheet of the meta.xlsx input file. Here is a summary of the primers and the attached sample indexes.

amplicon	primer	seq	index length	min. index dist.	median dist.	max. dist.
ITS5_ITS4	ITS5	GGAAGTAAAAGTCGTAACAAGG	9	3	6	8
ITS5_ITS4	ITS4	TCCTCCGCTTATTGATATGC	9	3	6	9

The forward and reverse sample indexes have a length of 9 bp and a minimum edit distance of 3 (see also Figure 4.1). The maximum allowed number of sequencing errors in barcodes is 0. This value should always be smaller than the minimum distance between any sample index).

Note: Using a lower idx_max_diffs threshold reduces false assignments of reads to samples, but also leads to less sequences being included in the analysis. On the other hand, an increased error tolerance increases the read count. You can check the unused barcode combinations (unspecific in below figures) for any “background noise”.

Figure 4.1: Distribution of pairwise edit distances between sample indexes

5 De-multiplexing

Steps:

• Search for the primers (up to 20 % error rate)
• Search for the sample indexes next to the primers (max. 0 differences)
• Trim primers, filter by quality (max. 2.5 expected errors) and distribute into different files by barcode (de-multiplex)

5.1 Primer trimming

Assuming that the basecalling was run without adapter trimming (--no-trim Dorado option), there are always small tails of to remove before the primer starts:

amplicon	avg. amplicon length	avg. tail length
ITS5_ITS4	672.2 bp	46.1 / 36 bp

A small part of the reads has tail lengths around the mean fragment length (or multiples of it), which indicates the presence of concatenated artifacts (Figures 5.1, 5.2).

Figure 5.1: Distribution of the primer trimmed tail lengths (note the non-linear scale)

Trimmed reads are therefore checked for the presence of any primer sequence and removed if any was found. The average length of these fragments is usually longer than the mean length (Figure 5.2).

Figure 5.2: Distribution of trimmed fragment lengths (top panel), and fragment lengths for putative concatenated artifacts (bottom) containing primer sequences, which were removed.

5.2 Quality filtering

The quality filtering should remove highly errorneous sequencing reads, while still retaining enough reads for inferring the barcode sequences. The default threshold of 2.5 expected errors per read can be tweaked (error_threshold setting, use Figure 5.3 for guidance).

Figure 5.3: Read counts by category

According to Figure 5.4, reads with many primer or sample index mismatches also tend to have lower quality scores. With the current setting (error_threshold), reads with >2.5 errors and reads shorter than 200 bp (min_read_length setting) are discarded.

Figure 5.4: Mismatch statistics/quality statistics, based on the quality scores from the base-called FASTQ file. (a) Comparison of index/primer mismatches, colored by average read quality; (b/c) Distribution of expected errors per read, colored by primer/index mismatches.

6 Infer barcode sequences

Before the DADA2 denoising, it is necessary to determine how the quality scores in the sequencing reads correlate with the substitution error rates. Usually, the relationship looks fairly log-linear, although there can be a plateau/increase at Q > 40 (Figure 6.1).

Figure 6.1: Q-scores vs. error frequencies learned by DADA2

From each of the de-multiplexed and filtered read files, a maximum of 2000 reads (max_sample_depth setting) are analyzed (Figure 6.2a).

The DADA2 denoising method relies on presence of at least a few error-free reads in the filtered Nanopore data. This should usually be the case with recent (R10.4) sequencing data, and assuming that amplicons are not extremely long. Specifically, the proportion of singletons¹ has to be lower than 100 % (Figure 6.2b), as duplicated sequences are likely to represent true sequences in the sample, while singletons often contain sequencing errors.

The barcode clustering procedure automatically checks the read duplication, and switches from DADA2 denoising to fixed-threshold clustering if there are < 2 identical sequences or < 4 sequences without substitution errors (Figure 6.2b and c; cluster_fixed method). The thresholds can be changed in config.yaml (dada_min_identical and dada_min_n0). In some cases, DADA2-denoised ASVs can be further split into two more or less equally abundant haplotypes (dada_split method in Figure 6.2).

Figure 6.2: Details on the clustering, colored by clustering method. (a) Analyzed read depth of samples (up to 2000 reads); (b) Proportion of singleton reads in relation to sequencing depth; (c) Number of identical sequences supporting the dominant sequence of the top taxon, compared to its abundance. (d) Number of sequences without substitutions (but possibly some InDels) (n0 reported by DADA2, not known for other methods). The vertical line indicates a sequencing depth threshold of 20 reads (low_abund_threshold in config.yaml).

Ambiguous bases in the reported sequence mostly appear in low-depth samples, where the fixed-threshold clustering is applied (cluster_fixed method) and no further separation of haplotypes is attempted (Figure 6.3a/b). These ambiguities may be further inspected manually. They may result from either an unresolved mix of different haplotypes or other sequence polymorphisms, or possibly from sequencing errors (usually homopolymers, long stretches of the same base).

In rare cases, the alignment consensus does not match the inferred DADA2 ASV sequence (Figure 6.3c/d). These barcodes will have a consensus-diffs issue flag in the report, and the consensus sequence is reported (may be inspected manually, but usually the consensus should make sense).

Figure 6.3: (a/b) Distribution of the ambiguous base count in the top consensus sequence (< 65 % in alignment having the same base). (c/d) Number of mismatches between the DADA2 ASV and the consensus sequence. The vertical line indicates a sequencing depth threshold of 20 reads (low_abund_threshold in config.yaml).

7 Haplotypes / polymorphism

The clustering procedure usually yields 1-2 abundant haplotype sequences from the target organism, and few rare/unspecific barcode sequences (Figure 7.1a/b). In addition, there may be additional copies of the target barcode in the genomes. A fixed-threshold clustering at a similarity threshold of 97 % is done to distinguish the target organism from other organisms, which may also be present in the mix (e.g. contamination). If DNA was extracted from the tissue of multiple individuals, there might also be more than the two parental haplotypes. Polymorphisms with < 10 % frequency in an organism group are removed.

The proportion of reads belonging to the final (probably specific) clusters is usually 50-100% at read depths > 1000, but decreases with lower depths (Figure 7.1c). Reasons can be sequencing errors (also in the sample indexes), as well as contamination.

Figure 7.1: (a/b) Distribution of putative haplotype numbers for the top organism in a sample; (c) Proportion of reads from the top organism group relative to the total analyzed read depth

8 Taxonomic assignments

Taxonomic names are assigned with the the SINTAX algorithm using a bootstrap confidence threshold of 80 %.

Caution: these automatic sequence-based assignments are not guaranteed to be correct, it is recommended to further validate them, e.g. using BLAST.

Any blacklisted taxa (known_contaminants in config.yaml) are directly flagged as contaminants, and other barcodes in the same sample are preferred.

Second, morphological identifications (morpho_taxon column in metadata) are compared to the automatic taxonomic assignments. From this, a taxonomic overlap is calculated, which is the proportion of matching ranks in the taxonomic lineages (as far as names are defined at these levels). The lineages are retrieved from the GBIF backbone taxonomy.

A barcode sequence (more precisely, a 97 % group) is automatically classified as contaminant if there exists another sequence (group), which has a better taxonomic overlap with the morphological identification. The overlap is assumed to be better if there exists a less abundant taxon, which has at least three more consistent (matching) taxonomic ranks in the GBIF lineage compared to the top taxon (undefined ranks excluded from the comparison). Figure 8.1 visualizes the final taxonomic assignments. The taxonomic overlap (Figure 8.1b) should mostly be near 100% (bars on right side).

note: the ranking may sometimes be incorrect, or if contaminants are too related to the target taxon, they may not automatically be flagged; a manual investigation of the report is highly recommended.

The following barcodes were found to contain strong contamination:

sample	morpho_taxon	correct taxon	contaminants
sample10	Tricholoma spec.	Tricholoma sp.	Unknown Saccharomycetales
sample99	Cuphophyllus spec.	Cuphophyllus sp.	Cystofilobasidium sp.
sample208	Ramularia sambucina	Ramularia sp.	Sambucus sp., Vishniacozyma sp., Unknown Didymellaceae, Neosetophoma sp.

Figure 8.1: Overview of taxonomic assignments, summarized per plate. (a) Taxonomic composition (based on barcode sequences); (b) Taxonomic overlap between the morphological and sequence-based taxonomic assignments

For samples with putative contaminants, the clustering is re-run without sample depth limit in order to ensure the effective detection of even rare target barcodes.

9 Comparison to known sequences

In case of re-sequencing a specimen, there is the additional possibility of validating the result by calculating the number of differences (mismatches/gaps). Inconsistent barcode sequences might be investigated further. It is still possible, that e.g. the provided Sanger sequences have errors as well, it is best to also investigate the chromatograms in case of inconsistencies.

Figure 9.1: Mismatches to already known sequences (in case of re-sequencing). The vertical red line indicates a critical sequencing depth threshold of 20 reads.

10 Report

The final results are written to an Excel document (output/report.xlsx) containing different sections besides the metadata already present in meta.xlsx. The meaning of different columns is explained in comments attached to the headers.

---

1. Singletons are sequences only present once in a sample (non-duplicated)↩︎