auto16s

auto16s takes paired-end 16S rRNA amplicon reads from raw FASTQ files to diversity statistics, differential abundance results and a written report.

It is designed for a typical study: samples sequenced in one Illumina run and compared between groups. Each step uses a widely used tool (DADA2, SILVA, phyloseq, vegan, MaAsLin 3) with fixed settings (listed under Parameters), so there is little to configure. What usually has to be worked out by hand (which primers are in the reads, where to truncate them, how deep to rarefy) is determined from the data, and each decision is shown in the report.

The only required input is a samplesheet.

Contents

• Pipeline summary
• Quick start
• Input
• Parameters and defaults
• What the pipeline decides from the data
• Outputs
• Statistics
• Two variables (e.g. disease and sex)
• Figures
• Reading the report
• Scope and limitations
• Testing
• Citations
• License

Pipeline summary

1. Read QC of the raw reads (FastQC)

2. Primer detection: the read starts are matched against a library of common 16S primers (assets/primers_16s.tsv).

3. Primer removal, per read and by sequence. This handles variable-length spacers before the primer and pairs in reverse orientation. Pairs without both primers are dropped. Afterwards every read is checked again for the primer, so the report shows whether removal worked. When the amplicon is shorter than the reads (e.g. V4 sequenced 2×300), reads run past its end into the opposite primer and the sequencing adapter; they are cut there.

4. ASV inference (DADA2):

   • quality filtering with automatic truncLen and maxEE;
   • error models for R1 and R2;
   • denoising;
   • mergePairs;
   • consensus chimera removal.

   The filtered and truncated reads go through FastQC a second time, so the report shows the quality before and after filtering side by side.

5. Taxonomy: DADA2 assignTaxonomy against SILVA 138.1.

6. phyloseq object (phyloseq). Unassigned, eukaryotic, chloroplast and mitochondrial ASVs are removed, then rarefaction at an automatically chosen depth.

7. Composition: phylum, family and genus bar plots per sample and per group, plus a genus heatmap (all reads).

8. Alpha diversity: Observed, Shannon and Simpson, with Wilcoxon or Kruskal-Wallis tests (rarefied).

9. Beta diversity: Bray-Curtis PCoA, PERMANOVA, betadisper and pairwise PERMANOVA (vegan; rarefied).

10. Differential abundance and prevalence at genus level (MaAsLin 3; all reads).

With a second grouping column (--strata, e.g. sex), steps 8–10 run for each comparison: the main variable overall (adjusted, with the interaction test) and within each level of the second variable, and the second variable within each group (details). 11. QC summary (MultiQC): FastQC before and after filtering, primer detection and removal, DADA2 steps, warnings and the main figures. 12. Final report: one HTML file that reads like a short study report: the samples analysed, a methods section written with the values and software versions of this run (ready to adapt for a paper), and the results with all key figures and tables. Every table and figure is also delivered as a separate file.

Quick start

1. Install Nextflow (≥ 25.04, Java 17+) and one of the software stacks below. On Windows, use WSL.

2. Test the installation (8 small samples, a few minutes). Nextflow downloads the pipeline from GitHub:

   nextflow run ozlemsagiroglu/auto16s -profile test,docker      # or test,singularity / test,conda

3. Create the samplesheet from your FASTQ folder, then fill in the group column (details):

   nextflow run ozlemsagiroglu/auto16s --input raw/

4. Run your data:

   nextflow run ozlemsagiroglu/auto16s -profile docker --input samplesheet.csv --outdir results

5. Open results/report/auto16s_report.html for the samples, methods and results. For quality control, open results/multiqc/multiqc_report.html and read the "DADA2 summary and warnings" table first.

Add -resume to rerun after a change; finished steps are reused. Add -r v0.6.0 to run a fixed release.

Software stacks: Docker is not required

All tools come either from containers or from Conda. Choose the profile that matches what is installed:

Profile	Needs	Notes
docker	Docker	Simplest on a laptop or workstation
singularity / apptainer	Singularity or Apptainer	No Docker needed. Usual choice on HPC clusters; the same images are pulled and converted automatically
conda	Conda / Mamba / Miniforge	No containers at all. Environments are created on the first run (this takes a while once)

The tools of the R analysis steps (phyloseq, vegan, MaAsLin 3, ggplot2) are in the image ghcr.io/ozlemsagiroglu/auto16s-r, built from containers/r-analysis and published automatically. Nothing has to be built by hand. For offline clusters, the image can be built without Docker:

cd containers/r-analysis && apptainer build auto16s-r.sif auto16s-r.def
nextflow run ozlemsagiroglu/auto16s -profile singularity --r_container $PWD/auto16s-r.sif --input samplesheet.csv

Without any container or Conda, the pipeline also runs on a machine where FastQC, MultiQC and R are installed. It needs R ≥ 4.4 with dada2, phyloseq, vegan, maaslin3, ggplot2 and ggrepel; the R packages can be installed with Rscript containers/r-analysis/install_r_packages.R (plus BiocManager::install("dada2")). Then run without -profile.

Input

The pipeline needs one table, the samplesheet. It tells the pipeline which FASTQ files belong to which sample and which group each sample is in. One row per sample.

Let the pipeline write it

Point --input to the folder with your FASTQ files:

nextflow run ozlemsagiroglu/auto16s --input raw/

The pipeline:

• pairs the R1 and R2 files by their names;
• writes samplesheet.csv in the current folder, with the group column left empty;
• stops.

Fill in the group of every sample (for example control or patient), then start the analysis:

nextflow run ozlemsagiroglu/auto16s -profile docker --input samplesheet.csv

File names it recognises: S1_R1.fastq.gz/S1_R2.fastq.gz, S1_1.fq.gz/S1_2.fq.gz, S1.R1.fastq, and Illumina's S1_S1_L001_R1_001.fastq.gz; the last one becomes sample S1. Sample names are cleaned automatically: characters such as ş, ç or spaces are replaced, and a name that starts with a digit gets an S in front. Files without a partner are listed and left out.

Or write it yourself

Example: a project folder

project/
├── samplesheet.csv
└── raw/
    ├── K1_R1.fastq.gz    K1_R2.fastq.gz
    ├── K2_R1.fastq.gz    K2_R2.fastq.gz
    ├── P1_R1.fastq.gz    P1_R2.fastq.gz
    └── P2_R1.fastq.gz    P2_R2.fastq.gz

and its samplesheet.csv:

sample,fastq_1,fastq_2,group,age,sex
K1,raw/K1_R1.fastq.gz,raw/K1_R2.fastq.gz,control,34,F
K2,raw/K2_R1.fastq.gz,raw/K2_R2.fastq.gz,control,41,M
P1,raw/P1_R1.fastq.gz,raw/P1_R2.fastq.gz,patient,38,F
P2,raw/P2_R1.fastq.gz,raw/P2_R2.fastq.gz,patient,45,M

Column	What to write	Rules
sample	Name of the sample; it is used in every table and figure	Must start with a letter; English letters, digits, ., _ and - only (no spaces). Each name once
fastq_1	Path to the sample's R1 (forward) file	.fastq.gz, .fq.gz, .fastq or .fq
fastq_2	Path to the sample's R2 (reverse) file	same
group	Group the sample belongs to (control, patient, treatment, …)	At least 2 groups; 3 or more samples per group recommended. Spelling must match exactly: Control and control are two different groups. Another column can be used with --group_col
other columns	Optional, e.g. age or sex	Not used in the analysis; kept in the phyloseq object for your own analyses in R

The first line must contain the column names exactly as above.

File paths can be written in two ways:

• relative to the samplesheet's folder, as in the example (raw/K1_R1.fastq.gz); the project folder can then be moved anywhere;
• absolute, e.g. /home/user/project/raw/K1_R1.fastq.gz. On Windows with WSL, write /mnt/c/Users/..., not C:\....

Avoid spaces in folder and file names.

Excel works too. Columns may be separated by commas, semicolons (what Excel writes with Turkish and most European language settings) or tabs; the separator is recognised automatically. "Save as CSV" from Excel is fine.

Checks before anything runs. The pipeline stops with a clear message on:

• missing columns, duplicate or invalid sample names;
• files that do not exist or are not FASTQ;
• empty groups, or fewer than 2 groups.

It warns about groups with fewer than 3 samples, and prints the number of samples per group at the start.

Parameters and defaults

Required: --input: the samplesheet, or a folder of FASTQ files to create one (see Input).

Optional, only if needed:

Parameter	Default	Description
--outdir	results	Output folder
--group_col	group	Samplesheet column with the groups
--strata	none	A second samplesheet column, e.g. sex; see Two variables
--maaslin_reference	alphabetically first group	Reference group for MaAsLin 3 coefficients
--silva_db	SILVA 138.1 training set, downloaded from Zenodo	Local copy, to avoid downloading in every new run
--fw_primer, --rv_primer	detected	Primer sequences (5'→3', IUPAC), only if your primers are not in the library
--trunc_len_f, --trunc_len_r	automatic	Fixed truncation lengths
--rarefy_depth	automatic	Fixed rarefaction depth
--max_cpus, --max_memory, --max_time	all CPUs and all memory of the machine, 24.h	Upper limits per task. Tasks never ask for more than these, so the pipeline also runs on a laptop

Settings used in every run. The pipeline runs with the following values (stored in nextflow.config):

Step	Setting	Value	What it does	Why this value
Primer removal	Primer position	within the first 12 bases	Where in the read the primer may start	Leaves room for the variable-length spacers some library preparations add before the primer
	Mismatches	1 (primer < 20 nt), 2 (≥ 20 nt); IUPAC codes honoured	Differences tolerated when matching the primer	Tolerates sequencing errors in the primer without matching unrelated sequences
Adapter removal	Read-through at the 3' end	cut before the reverse complement of the opposite primer, or before an Illumina adapter (TruSeq AGATCGGAAGAGC, Nextera CTGTCTCTTATACACATCT); partial matches of ≥ 10 bases at the read end count	Removes primer and adapter bases that follow the end of the amplicon	Happens when the amplicon is shorter than the reads; the opposite primer marks the amplicon end for any library kit, the adapter sequences catch reads with a damaged primer
DADA2 filterAndTrim	truncLen	chosen from the data	Cuts every read at a fixed length	DADA2 needs reads of equal length; the position keeps the most read pairs while R1 and R2 still overlap (see below)
	maxEE	2 (R1 and R2)	Removes reads with more than 2 expected errors	The value used in the DADA2 tutorial; removes error-rich reads and keeps most of the rest
	truncQ	2	Cuts a read at the first base with quality ≤ 2	Q2 marks unusable bases on Illumina; such reads become too short and are removed
	maxN	0	Removes reads that contain an N	DADA2's denoising does not accept N bases
	rm.phix	TRUE	Removes PhiX reads	PhiX is a control added to Illumina runs, not part of the sample
DADA2 mergePairs	minOverlap	12	Minimum overlap needed to merge R1 and R2	DADA2 default; mismatches in the overlap are not allowed
	trimOverhang	TRUE	Removes bases that extend past the start of the partner read	Prevents leftover primer or adapter bases from ending up in the ASVs
Chimera removal	removeBimeraDenovo	consensus	Checks every sample for chimeras and removes ASVs flagged in most of the samples that contain them	DADA2 default
Taxonomy	assignTaxonomy minBoot	50	Minimum bootstrap confidence to assign a rank	DADA2 default; 50 % is recommended for reads shorter than 250 bp (Wang et al. 2007)
Rarefaction	Depth	chosen from the data	Subsamples every sample to the same number of reads	Alpha and beta diversity need equal depth; the rule keeps one failed sample from lowering all others (see below)
	Seed	711	Starting value of the random subsampling	Makes the results reproducible
PERMANOVA, betadisper	Permutations	999	Number of random permutations for the p-value	vegan default; the smallest possible p-value is 0.001
MaAsLin 3	Models	abundance (linear) and prevalence (logistic)	Tests separately whether a genus is more abundant and whether it is present more often	MaAsLin 3 default; absence and low abundance are different effects
	Normalisation, transform	TSS, LOG	Relative abundance, then log	MaAsLin 3 defaults
	min_prevalence	0.1	Tests only genera found in at least 10 % of samples	Rarer genera cannot be tested reliably and only add to the multiple-testing correction
	Significance	q < 0.05 (Benjamini-Hochberg)	Threshold for reporting a genus	A common threshold; stricter than the MaAsLin 3 default of 0.1
Bar plots	Taxa shown	top 12 by mean relative abundance, rest as "Other"	Number of taxa with their own colour	More colours can no longer be told apart

What the pipeline decides from the data

Primers

The first reads of every sample are matched against assets/primers_16s.tsv. The list has 20 widely used 16S primers, among them 27F, 341F, 515F, 515F-Y, 799F, 805R, 806R, 806RB, 926R and 1193R. The forward/reverse pair found in most read pairs is selected and removed from every read by sequence, which also handles:

• spacers: variable-length bases before the primer ("heterogeneity spacers", phased primers);
• mixed orientation: pairs whose R1 starts with the reverse primer, as in some ligation-based libraries, are swapped.

Pairs that do not carry both primers are removed; they are usually off-target products. The detected primers, the share of read pairs that carry them and a per-sample summary are in primers/ and in the report.

If no known primer pair is found, the pipeline assumes the primers were already removed by the sequencing provider and uses the reads as they are. The report then shows the most common read start, so an unlisted primer can be recognised. Give such a primer with --fw_primer/--rv_primer, or add one line to the TSV.

Truncation length (truncLen)

The truncation length is computed from the reads themselves. FastQC is only for viewing. On the first reads of every sample, after primer removal, the pipeline:

1. estimates the amplicon length by aligning R1 to the reverse complement of R2 (the longest common length is used, so variable-length regions such as V3-V4 are covered);
2. for every pair of truncation positions, counts how many read pairs would pass the quality filter (maxEE = 2, no N, no base with Q ≤ 2);
3. keeps only the pairs that leave R1 and R2 overlapping by at least 20 bp;
4. picks the pair that keeps the most read pairs. Among pairs within 1 percentage point of the best, the longest reads are preferred.

Truncating too late loses reads, because errors at the 3' end push them over the expected-error limit. Truncating too early prevents merging, the most common reason why read pairs silently fail in DADA2. This procedure finds the balance from the data, following the idea of FIGARO (Weinstein et al. 2019). On the test data it keeps 83 % of the read pairs, against 68 % with a rule based on the median quality. The quality profile, the chosen positions, the expected overlap and the expected share of read pairs passing are shown in dada2/dada2_quality_truncation.png and in the report.

Rarefaction depth

Rarefaction is used for alpha and beta diversity only. The depth is the smallest depth among samples with at least 1,000 reads and at least 10 % of the median depth. One failed sample therefore cannot force all others down to a few hundred reads. Samples below the depth are left out of alpha and beta diversity and listed in the report. Composition and MaAsLin 3 use all samples and all reads.

Outputs

All results are in --outdir. Besides the figures, every number behind a figure is delivered as a table, and the full data as phyloseq objects, so you can make your own figures or analyses.

Folder	Content
report/	auto16s_report.html: final report (samples, methods, results); open it first
multiqc/	multiqc_report.html: QC summary (FastQC, primers, DADA2 steps, warnings)
phyloseq/	phyloseq_object.rds (all reads), phyloseq_object_rarefied.rds, ASV tables (raw and rarefied), taxonomy, ASV sequences (FASTA), sample metadata, depths, rarefaction curves
composition/	Counts and relative abundances per phylum, family and genus; bar plots; genus heatmap
alpha_diversity/	Values per sample, tests, figure
comparisons/	Only with --strata: one folder per comparison, each with alpha_diversity/, beta_diversity/ and maaslin3/ (these then replace the top-level folders)
beta_diversity/	Bray-Curtis distance matrix, PCoA coordinates, PERMANOVA, betadisper, interpretation, figures
maaslin3/	All results and significant genera (CSV), volcano, coefficient plot, abundance and prevalence plots, raw MaAsLin 3 output
primers/	Detected primers, detection plot, per-sample primer statistics
dada2/	Summary and warnings, read tracking, truncation, error models, ASV lengths, sequence table
taxonomy/	SILVA assignment per ASV sequence
fastqc/	FastQC reports: raw/ before and filtered/ after DADA2 quality filtering
pipeline_info/	Nextflow execution report, timeline, trace

Each file is described in docs/output.md. Figures are saved as PNG (300 dpi) and PDF.

Continuing in R

library(phyloseq)
ps      <- readRDS("results/phyloseq/phyloseq_object.rds")            # all reads
ps_rare <- readRDS("results/phyloseq/phyloseq_object_rarefied.rds")   # rarefied

sample_data(ps)                       # your samplesheet columns
plot_bar(tax_glom(ps, "Phylum"), fill = "Phylum")

genus <- read.delim("results/composition/genus_relative_abundance.tsv", row.names = 1, check.names = FALSE)

Statistics

Analysis	Data	Method
Alpha diversity	rarefied	2 groups: Wilcoxon rank-sum. More than 2: Kruskal-Wallis and pairwise Wilcoxon, Benjamini-Hochberg
Beta diversity	rarefied	Bray-Curtis; PERMANOVA (adonis2) and betadisper. More than 2 groups: pairwise PERMANOVA, Benjamini-Hochberg
Differential abundance and prevalence	all reads, genus level	MaAsLin 3 (TSS, LOG; abundance and prevalence models), Benjamini-Hochberg q-values
With --strata (overall comparison)	as above	Alpha: linear model of the group adjusted for the second variable. Beta: PERMANOVA with marginal tests (by = "margin"). MaAsLin 3: second variable as covariate. Interaction: linear model, PERMANOVA with sequential terms, MaAsLin 3 with an interaction term

PERMANOVA cannot tell a shift in composition from a difference in within-group spread. The pipeline therefore reads it together with betadisper and writes the interpretation into the report. Rarefying for alpha and beta diversity, but not for differential abundance, follows current recommendations (Schloss 2024).

Why MaAsLin 3 reports two kinds of results. A genus can differ between groups in two ways: it is present in fewer samples of one group (prevalence), or it is present everywhere but less abundant (abundance). Older methods, including MaAsLin2, mix the two, because absent genera enter the model as zeros. MaAsLin 3 tests them separately: the abundance model uses only the samples in which the genus was found, the prevalence model only whether it was found. Every significant result is therefore labelled abundance or prevalence. Following the MaAsLin 3 defaults, abundance coefficients are tested against the median coefficient of all genera rather than against zero, which corrects for the compositional nature of relative abundances.

Two variables (e.g. disease and sex)

With a second column in the samplesheet, one run tests the main variable overall and within each level of the second, and the second within each group of the main variable:

nextflow run ozlemsagiroglu/auto16s -profile docker --input samplesheet.csv --group_col disease --strata sex

Comparison	Samples	Model
disease	all	disease adjusted for sex; plus the disease × sex interaction
disease_within_sex-female, disease_within_sex-male	one sex	disease
sex_within_disease-healthy, sex_within_disease-patient	one disease group	sex

Each comparison gets its own alpha diversity, beta diversity and MaAsLin 3 results in comparisons/<name>/, and its own section in the final report. Composition plots show the combined groups (e.g. "patient / female").

Reading the results. A difference that is significant in women but not in men does not show that the effect differs between the sexes; the smaller subgroup may simply lack power. Whether the disease effect depends on sex is answered by the interaction test of the first comparison (linear model for alpha diversity, PERMANOVA for composition, MaAsLin 3 interaction term per genus). The interaction test needs more samples than the main effect. Each subgroup comparison uses only part of the samples; the run log and the report show the number of samples per combination, and groups with fewer than 3 samples get a warning. Within-group dispersion (betadisper) is not tested with fewer than 3 samples per group.

Figures

All figures are made with ggplot2. Colours are fixed, so a group or a taxon looks the same in every figure:

Element	Colours
Groups	Okabe-Ito palette (colour-blind safe), in alphabetical order of group names
Taxa in bar plots	20 distinct colours by abundance rank; "Other" in grey
Heatmap	viridis, log10 relative abundance
MaAsLin 3	orange = higher, blue = lower than the reference group

Reading the report

The report starts with DADA2 summary and warnings. Possible warnings:

Warning	Most likely cause
No known 16S primer pair found	Primers were already removed (fine), or a primer that is not in the library. Check the reported read start
Only x % of read pairs carry the primers	Off-target products, or a primer variant not in the library
< 80 % of reads survived chimera removal	Primers still in the reads
Median merge rate < 70 %	Reads too short to overlap after truncation
Primer still found in > 1 % of reads after removal	A primer variant that is not in the library, or very long spacers

Did filtering work? The report has two FastQC sections, raw reads and after DADA2 filtering. After filtering, all reads of a sample have the same length, the low-quality 3' ends are gone and the per-base quality stays high to the end of the read. The section Primer removal check lists, per sample, the share of reads that still contain the primer after removal; it should be close to 0 %.

Merged reads cannot be checked with FastQC, because DADA2 merges the denoised sequences (ASVs), not the reads. The merge is checked by the read tracking table (reads lost at merging) and the ASV length plot below.

Sanity check: ASV length. dada2/dada2_asv_length.png should show one main peak at the length of your region without primers (e.g. V4 ≈ 253 bp, V3-V4 ≈ 400–430 bp). Small peaks far from it are usually off-target products.

Troubleshooting

Message	Solution
Process requirement exceeds available memory	Only with older versions: update with nextflow pull ozlemsagiroglu/auto16s. The pipeline now limits every task to the memory of the machine. To keep memory free for other work, set a lower limit, e.g. --max_memory 6.GB
DADA2 is killed (exit status 137 or 143, Out of memory in journalctl)	DADA2 already uses at most one thread per 2 GB of free memory. If it still runs out, lower the threads with --max_cpus 2, close other programs or, under WSL, raise the memory limit in %UserProfile%\.wslconfig ([wsl2] / memory=12GB) and run wsl --shutdown. Then add -resume
A run stops halfway	Fix the cause and run the same command again with -resume; finished steps are not repeated

Scope and limitations

The workflow is deliberately fixed and assumes:

• One Illumina paired-end run. DADA2's error model is run-specific; samples from different runs should not be combined in one analysis.
• Primers from the built-in list, or given with --fw_primer/--rv_primer.
• No negative controls. If extraction or PCR blanks were sequenced, contaminant removal (e.g. decontam) is standard and is not included.
• No phylogenetic tree, so UniFrac and Faith's PD are not computed.
• A comparison between groups, optionally split by one second variable (--strata). Continuous covariates (e.g. age) and paired or repeated-measures designs are not modelled; use the phyloseq objects in R for those.
• One differential abundance method. Differential abundance methods can give different results on the same data (Nearing et al. 2022). Report MaAsLin 3 results as such.
• Small groups. Permutation tests cannot produce small p-values with very few samples. With 3 vs 2 samples there are only 10 distinct permutations, so p ≥ 0.1.

Testing

-profile test runs real reads: the V4 example reads (mouse gut, 2×250 bp) and the small example taxonomy training set distributed with the DADA2 R package, with 515F-Y/806RB primers added and their degenerate bases resolved per read. The full SILVA database is not part of the repository; it is downloaded from Zenodo when the pipeline runs (CC BY 4.0). It has 8 samples in 2 groups with a built-in difference in composition. Three other scenarios can be generated with the same script:

Rscript tests/make_test_data.R tests/data_spacer spacer     # 0-7 nt spacers, half of the pairs reverse-oriented
Rscript tests/make_test_data.R tests/data_noprimer none     # primers already removed
Rscript tests/make_test_data.R tests/data_readthrough readthrough   # 2x300 bp: reads run into primer and adapter
nextflow run main.nf -profile docker --input tests/data_spacer/samplesheet.csv \
    --silva_db tests/data/example_train_set.fa.gz

In all of them, the primers are found (or correctly reported as absent), the amplicon length is estimated at 253 bp (V4), and no reads are lost to chimeras. In the read-through scenario every read is cut to exactly 253 bp and no adapter or primer remains in the ASVs.

Citations

If you use auto16s, please cite it (GitHub's "Cite this repository" button gives the reference from CITATION.cff):

Sağıroğlu Ö. (2026). auto16s: automated analysis of paired-end 16S rRNA amplicon data (version 0.6.0). https://github.com/ozlemsagiroglu/auto16s

and the tools it relies on:

• Nextflow: Di Tommaso P. et al. (2017). Nextflow enables reproducible computational workflows. Nat Biotechnol 35:316–319.
• FastQC: Andrews S. (2010). FastQC: a quality control tool for high throughput sequence data.
• DADA2: Callahan B.J. et al. (2016). DADA2: high-resolution sample inference from Illumina amplicon data. Nat Methods 13:581–583.
• SILVA: Quast C. et al. (2013). The SILVA ribosomal RNA gene database project. Nucleic Acids Res 41:D590–D596.
• SILVA for DADA2: McLaren M.R., Callahan B.J. (2021). SILVA 138.1 prokaryotic SSU taxonomic training data formatted for DADA2. Zenodo. doi:10.5281/zenodo.4587955
• phyloseq: McMurdie P.J., Holmes S. (2013). phyloseq. PLoS ONE 8:e61217.
• vegan: Oksanen J. et al. vegan: Community Ecology Package.
• betadisper: Anderson M.J. (2006). Distance-based tests for homogeneity of multivariate dispersions. Biometrics 62:245–253.
• MaAsLin 3: Nickols W.A. et al. (2026). MaAsLin 3: refining and extending generalized multivariate linear models for meta-omic association discovery. Nat Methods.
• MultiQC: Ewels P. et al. (2016). MultiQC. Bioinformatics 32:3047–3048.
• ggplot2: Wickham H. (2016). ggplot2: Elegant Graphics for Data Analysis. Springer.

Background for design choices:

• Wang Q. et al. (2007). Naive Bayesian classifier for rapid assignment of rRNA sequences into the new bacterial taxonomy. Appl Environ Microbiol 73:5261-5267.
• Weinstein M.M. et al. (2019). FIGARO: an efficient and objective tool for optimizing microbiome rRNA gene trimming parameters. bioRxiv 610394.
• Klindworth A. et al. (2013). Evaluation of general 16S ribosomal RNA gene PCR primers. Nucleic Acids Res 41:e1.
• Schloss P.D. (2024). Rarefaction is currently the best approach to control for uneven sequencing effort in amplicon sequence analyses. mSphere.
• Nearing J.T. et al. (2022). Microbiome differential abundance methods produce different results across 38 datasets. Nat Commun 13:342.

License

The auto16s code is released under the MIT License. The tools it runs are separate programs under their own licenses (for example DADA2: LGPL, phyloseq: AGPL-3, vegan: GPL-2, FastQC and MultiQC: GPL-3, MaAsLin 3: MIT). The SILVA reference database is downloaded at run time and is licensed under CC BY 4.0. The test data in tests/data are derived from the example reads distributed with the DADA2 package.