Version Compatibility

Reference examples tested with: DADA2 1.30+, cutadapt 4.4+

Before using code patterns, verify installed versions match. If versions differ:

• R: packageVersion('<pkg>') then ?function_name to verify parameters
• CLI: <tool> --version then <tool> --help to confirm flags

If code throws ImportError, AttributeError, or TypeError, introspect the installed package and adapt the example to match the actual API rather than retrying.

eDNA Metabarcoding

"Process my eDNA samples to identify species" → Transform raw amplicon reads (COI, 12S, rbcL, ITS) into species occurrence tables through primer removal, denoising (DADA2 ASVs), taxonomy assignment against reference databases (BOLD, MIDORI2), and contamination filtering with decontam.

• CLI: cutadapt for primer removal, obi (OBITools3) for paired-end assembly
• R: dada2::filterAndTrim() → dada() → assignTaxonomy() for ASV pipeline

Processes environmental DNA amplicon reads into species occurrence tables with taxonomy assignment, contamination filtering, and occupancy modeling.

Primer Removal with Cutadapt

Goal: Remove amplicon primers from paired-end eDNA reads while discarding untrimmed read pairs.

Approach: Use cutadapt linked adapter trimming with marker-specific primer sequences and minimum overlap enforcement.

Linked adapter trimming removes primer pairs while discarding reads lacking primers:

# COI primers (mlCOIintF / jgHCO2198)
cutadapt -g 'GGWACWGGWTGAACWGTWTAYCCYCC;min_overlap=20' \
         -G 'TAIACYTCIGGRTGICCRAARAAYCA;min_overlap=20' \
         --discard-untrimmed --pair-filter=any \
         -o trimmed_R1.fastq.gz -p trimmed_R2.fastq.gz \
         raw_R1.fastq.gz raw_R2.fastq.gz

# 12S MiFish-U primers
cutadapt -g 'GTCGGTAAAACTCGTGCCAGC;min_overlap=18' \
         -G 'CATAGTGGGGTATCTAATCCCAGTTTG;min_overlap=18' \
         --discard-untrimmed --pair-filter=any \
         -o trimmed_R1.fastq.gz -p trimmed_R2.fastq.gz \
         raw_R1.fastq.gz raw_R2.fastq.gz

# ITS primers (ITS1F / ITS2)
cutadapt -g 'CTTGGTCATTTAGAGGAAGTAA;min_overlap=18' \
         -G 'GCTGCGTTCTTCATCGATGC;min_overlap=18' \
         --discard-untrimmed --pair-filter=any \
         -o trimmed_R1.fastq.gz -p trimmed_R2.fastq.gz \
         raw_R1.fastq.gz raw_R2.fastq.gz

OBITools3 Pipeline

Goal: Process paired-end eDNA reads through the full OBITools3 workflow to produce a taxonomy-assigned species table.

Approach: Import reads, align pairs, filter by score/length, demultiplex, dereplicate, denoise with obi clean, assign taxonomy with ecotag, and export.

Full pipeline from paired-end reads to taxonomy table:

# Import paired FASTQ into OBITools3 DMS
obi import --fastq-input raw_R1.fastq.gz EDNA/reads1
obi import --fastq-input raw_R2.fastq.gz EDNA/reads2

# Paired-end alignment
obi alignpairedend -R EDNA/reads2 EDNA/reads1 EDNA/aligned

# Filter by alignment score and length
# score >= 50: removes poorly overlapping pairs
# length 100-500: typical COI amplicon range
obi grep -p 'sequence["score"] >= 50' EDNA/aligned EDNA/filtered
obi grep -p 'len(sequence) >= 100 and len(sequence) <= 500' \
    EDNA/filtered EDNA/length_filtered

# Demultiplex if needed (with ngsfilter file)
obi ngsfilter -t ngsfilter.txt -u EDNA/unassigned \
    EDNA/length_filtered EDNA/demux

# Dereplicate (obi uniq creates merged_sample attribute automatically)
obi uniq EDNA/demux EDNA/derep

# Remove singletons (count >= 2 removes PCR/sequencing errors)
obi grep -p 'sequence["count"] >= 2' EDNA/derep EDNA/no_singletons

# Denoise (remove PCR/sequencing errors)
# obi clean: denoises by merging low-abundance variants into parent sequences
# -s merged_sample: per-sample denoising (critical for multi-sample datasets)
# -r 0.05: ratio threshold; -H: keep only head sequences (discard variants)
obi clean -s merged_sample -r 0.05 -H EDNA/no_singletons EDNA/denoised

# Taxonomy assignment against reference database
obi ecotag -R EDNA/refdb --taxonomy EDNA/taxonomy EDNA/denoised EDNA/assigned

# Export to tab-separated format
obi export --tab-output EDNA/assigned > species_table.tsv

DADA2 Pipeline for eDNA

Goal: Generate amplicon sequence variants (ASVs) from eDNA reads with error correction and chimera removal.

Approach: Run the DADA2 workflow: filter/trim, learn error rates, denoise, merge pairs, remove chimeras, and assign taxonomy against a reference database.

DADA2 provides an alternative ASV-based approach:

library(dada2)

# After cutadapt primer removal, get file paths
filt_fwd <- file.path('filtered', basename(fwd_reads))
filt_rev <- file.path('filtered', basename(rev_reads))

# Filter and trim
# maxEE=2: expected errors threshold balancing sensitivity/specificity
# truncLen: set based on quality profile inspection
out <- filterAndTrim(fwd_reads, filt_fwd, rev_reads, filt_rev,
                     maxN = 0, maxEE = c(2, 2), truncQ = 2,
                     minLen = 100, rm.phix = TRUE, multithread = TRUE)

# Learn error rates
err_fwd <- learnErrors(filt_fwd, multithread = TRUE)
err_rev <- learnErrors(filt_rev, multithread = TRUE)

# Denoise
dada_fwd <- dada(filt_fwd, err = err_fwd, multithread = TRUE)
dada_rev <- dada(filt_rev, err = err_rev, multithread = TRUE)

# Merge paired ends
# minOverlap=12: default; increase for short amplicons
merged <- mergePairs(dada_fwd, filt_fwd, dada_rev, filt_rev, minOverlap = 12)

# Build sequence table and remove chimeras
seqtab <- makeSequenceTable(merged)
# Chimera rate <20% is typical; >30% suggests library prep issues
seqtab_nochim <- removeBimeraDenovo(seqtab, method = 'consensus', multithread = TRUE)

# Taxonomy assignment against MIDORI2 or custom database
# minBoot=80: bootstrap confidence for genus-level; 50 for family
taxa <- assignTaxonomy(seqtab_nochim, 'MIDORI2_DADA2_COI.fasta.gz', minBoot = 80)

Reference Databases

Database	Markers	Format	Source
MIDORI2	COI, 12S, 16S, 18S	DADA2, BLAST, OBITools	midori2.info
MitoFish	12S (fish)	DADA2, BLAST	mitofish.aori.u-tokyo.ac.jp
BOLD	COI	FASTA, API	boldsystems.org
SILVA	18S	DADA2, mothur	silva-project.net
UNITE	ITS (fungi)	DADA2, BLAST	unite.ut.ee

Formatting a custom FASTA for DADA2 taxonomy assignment:

# DADA2 format: >Kingdom;Phylum;Class;Order;Family;Genus;Species
# Download MIDORI2 DADA2-formatted reference:
wget https://reference-midori.info/download/Databases/MIDORI2_DADA2/COI/MIDORI2_LONGEST_NUC_GB259_CO1_DADA2.fasta.gz

Contamination Filtering

Goal: Identify and remove contaminant ASVs introduced during DNA extraction, PCR, or sequencing.

Approach: Apply decontam frequency/prevalence methods using negative controls and DNA concentrations, then remove tag-jumping artifacts with microDecon.

decontam (frequency/prevalence method)

library(decontam)

# Frequency method: uses DNA concentration from extraction blanks
contam_freq <- isContaminant(seqtab_nochim, conc = dna_concentration, method = 'frequency')

# Prevalence method: compares true samples vs negative controls
contam_prev <- isContaminant(seqtab_nochim, neg = is_negative_control, method = 'prevalence',
                             threshold = 0.5)
# threshold=0.5: ASV must be more prevalent in negatives than samples to flag

# Combined method
contam_both <- isContaminant(seqtab_nochim, conc = dna_concentration,
                             neg = is_negative_control, method = 'both')

seqtab_clean <- seqtab_nochim[, !contam_both$contaminant]

microDecon for tag-jumping artifacts

library(microDecon)

# Remove tag-jumping contamination between samples
# numb.blanks: number of blank/negative control columns at start of table
# numb.ind: vector of sample counts per group
decon_result <- decon(otu_table, numb.blanks = 3, numb.ind = c(10, 10, 10))
cleaned_table <- decon_result$decon.table

Occupancy Modeling with occumb

Goal: Correct species occurrence estimates for imperfect detection across replicated eDNA samples.

Approach: Fit a multi-species occupancy model via MCMC (JAGS) on a 3D array of replicated read counts to estimate true occupancy probabilities.

Corrects for imperfect detection in metabarcoding replicates:

library(occumb)

# Requires replicated eDNA samples (biological or PCR replicates)
# y: 3D array [species x sites x replicates] of read counts
data <- occumbData(y = count_array, spec_cov = species_covariates,
                   site_cov = site_covariates)

# Fit multi-species occupancy model (requires JAGS)
# n.chains=4, n.iter=10000: standard MCMC settings
fit <- occumb(data = data, n.chains = 4, n.iter = 10000, n.thin = 5, n.burn = 2500)

# Extract corrected occupancy estimates
summary(fit)

Key Thresholds

Parameter	Value	Rationale
Min read depth per sample	1,000	Below this, rare species detection drops significantly
Singleton removal	count >= 2	Singletons often represent sequencing errors
Chimera rate	<20%	Higher rates suggest library prep problems
COI species identity	>= 97%	Standard COI barcode gap threshold
COI genus identity	>= 95%	Conservative genus-level assignment
12S species identity	>= 98%	12S has less divergence than COI
decontam threshold	0.5	Balanced false positive/negative rate
DADA2 minBoot	80	Genus-level confidence; use 50 for family

Related Skills

• biodiversity-metrics - Diversity analysis from species occurrence tables
• community-ecology - Environmental gradient analysis of communities
• microbiome/amplicon-processing - 16S clinical microbiome alternative
• read-qc/quality-reports - Upstream read quality assessment