nf-core/ampliseq

Path to sample sheet, either tab-separated (.tsv), comma-separated (.csv), or in YAML format (.yml/.yaml), that points to compressed fastq files.

The sample sheet must have two to four tab-separated columns/entries with the following headers:

• sampleID (required): Unique sample IDs, must start with a letter, and can only contain letters, numbers or underscores
• forwardReads (required): Paths to (forward) reads zipped FastQ files
• reverseReads (optional): Paths to reverse reads zipped FastQ files, required if the data is paired-end
• run (optional): If the data was produced by multiple sequencing runs, any string

Related parameters are:

• --pacbio and --iontorrent if the sequencing data is PacBio data or IonTorrent data (default expected: paired-end Illumina data)
• --single_end if the sequencing data is single-ended Illumina data (default expected: paired-end Illumina data)
• Choose an appropriate reference taxonomy for the type of amplicon (16S/18S/ITS/CO1) (default: DADA2 assignTaxonomy and 16S rRNA sequence database)

Path to fasta format file with sequences that will be taxonomically classified. The fasta file input option can be used to taxonomically classify previously produced ASV/OTU sequences.

The fasta sequence header line may contain a description, that will be kept as part of the sequence name. However, tabs will be changed into spaces.

Related parameters are:

• Choose an appropriate reference taxonomy for the type of amplicon (16S/18S/ITS/CO1) (default: DADA2 assignTaxonomy and 16S rRNA sequence database)

Path to folder containing compressed fastq files.

Example for input data organization from one sequencing run with two samples, paired-end data:

data
  ├─sample1_1_L001_R1_001.fastq.gz
  ├─sample1_1_L001_R2_001.fastq.gz
  ├─sample2_1_L001_R1_001.fastq.gz
  └─sample2_1_L001_R2_001.fastq.gz

Please note the following requirements:

1. The path must be enclosed in quotes
2. The folder must contain gzip compressed demultiplexed fastq files. If the file names do not follow the default ("/*_R{1,2}_001.fastq.gz"), please check --extension.
3. Sample identifiers are extracted from file names, i.e. the string before the first underscore _, these must be unique
4. If your data is scattered, produce a sample sheet
5. All sequencing data should originate from one sequencing run, because processing relies on run-specific error models that are unreliable when data from several sequencing runs are mixed. Sequencing data originating from multiple sequencing runs requires additionally the parameter --multiple_sequencing_runs and a specific folder structure.

Related parameters are:

• --pacbio and --iontorrent if the sequencing data is PacBio data or IonTorrent data (default expected: paired-end Illumina data)
• --single_end if the sequencing data is single-ended Illumina data (default expected: paired-end Illumina data)
• --multiple_sequencing_runs if the sequencing data originates from multiple sequencing runs
• --extension if the sequencing file names do not follow the default ("/*_R{1,2}_001.fastq.gz")
• Choose an appropriate reference taxonomy for the type of amplicon (16S/18S/ITS/CO1) (default: DADA2 assignTaxonomy and 16S rRNA sequence database)

In amplicon sequencing methods, PCR with specific primers produces the amplicon of interest. These primer sequences need to be trimmed from the reads before further processing and are also required for producing an appropriate classifier. Do not use here any technical sequence such as adapter sequences but only the primer sequence that matches the biological amplicon.

For example:

--FW_primer "GTGYCAGCMGCCGCGGTAA" --RV_primer "GGACTACNVGGGTWTCTAAT"

In amplicon sequencing methods, PCR with specific primers produces the amplicon of interest. These primer sequences need to be trimmed from the reads before further processing and are also required for producing an appropriate classifier. Do not use here any technical sequence such as adapter sequences but only the primer sequence that matches the biological amplicon.

For example:

--FW_primer GTGYCAGCMGCCGCGGTAA --RV_primer GGACTACNVGGGTWTCTAAT

This is optional, but for performing downstream analysis such as barplots, diversity indices or differential abundance testing, a metadata file is essential.

Related parameter:

• --metadata_category (optional) to choose columns that are used for testing significance

For example:

--metadata "path/to/metadata.tsv"

Please note the following requirements:

1. The path must be enclosed in quotes
2. The metadata file has to follow the QIIME2 specifications (https://docs.qiime2.org/2021.2/tutorials/metadata/)

The first column in the tab-separated metadata file is the sample identifier column (required header: ID) and defines the sample or feature IDs associated with your study. In addition to the sample identifier column, the metadata file is required to have at least one column with multiple different non-numeric values but not all unique. NB: without additional columns there might be no groupings for the downstream analyses.

Sample identifiers should be 36 characters long or less, and also contain only ASCII alphanumeric characters (i.e. in the range of [a-z], [A-Z], or [0-9]), or the dash (-) character. For downstream analysis, by default all numeric columns, blanks or NA are removed, and only columns with multiple different values but not all unique are selected.

The columns which are to be assessed can be specified by --metadata_category. If --metadata_category isn't specified than all columns that fit the specification are automatically chosen.

Path to file with information about sequenced regions, either tab-separated (.tsv), comma-separated (.csv), or in YAML format (.yml/.yaml). This initiates scaffolding multiple regions along a reference.

The file must have four headers:

• region: Unique region identifier
• region_length: Minimal length of region
• FW_primer: Forward primer sequence
• RV_primer: Reverse primer sequence

For more details check the usage documentation.

Related parameters are:

• --sidle_ref_taxonomy to select the reference taxonomic database
• --sidle_ref_tax_custom, --sidle_ref_seq_custom and related parameters for custom reference taxonomic database files

Set this parameter to your e-mail address to get a summary e-mail with details of the run sent to you when the workflow exits. If set in your user config file (~/.nextflow/config) then you don't need to specify this on the command line for every run.

When using a sample sheet with --input containing forward and reverse reads, specifying --single_end will only extract forward reads and treat the data as single ended instead of extracting forward and reverse reads.

This will cause the pipeline to

• not truncate input reads if not --trunclenf and --trunclenr are overwriting defaults
• remove reverse complement primers from the end of reads in case the read length exceeds the amplicon length

From R package 'ShortRead' function 'readFastq': Representation to be used for quality scores, must be one of Auto (infer automatically), FastqQuality (Phred-like base 33 encoding), SFastqQuality (Illumina base 64 encoding).

Expects one sub-folder per sequencing run in the folder specified by --input_folder containing sequencing data of the specific run. Sample identifiers are taken from sequencing files, specifically the string before the first underscore will be the sample ID. Sample IDs across all sequencing runs (all sequencing files) have to be unique. If this is not the case, please use a sample sheet as input instead.

Example for input data organization:

data
  |-run1
  |  |-sample1_1_L001_R1_001.fastq.gz
  |  |-sample1_1_L001_R2_001.fastq.gz
  |  |-sample2_1_L001_R1_001.fastq.gz
  |  |-sample2_1_L001_R2_001.fastq.gz
  |
  |-run2
     |-sample3_1_L001_R1_001.fastq.gz
     |-sample3_1_L001_R2_001.fastq.gz
     |-sample4_1_L001_R1_001.fastq.gz
     |-sample4_1_L001_R2_001.fastq.gz

Example command to analyze this data in one pipeline run:

nextflow run nf-core/ampliseq \
    -profile singularity \
    --input_folder "data" \
    --FW_primer "GTGYCAGCMGCCGCGGTAA" \
    --RV_primer "GGACTACNVGGGTWTCTAAT" \
    --metadata "data/Metadata.tsv" \
    --multiple_sequencing_runs

Indicates the naming of sequencing files (default: "/*_R{1,2}_001.fastq.gz").

Please note:

1. The prepended slash (/) is required
2. The star (*) is the required wildcard for sample names
3. The curly brackets ({}) enclose the orientation for paired end reads, separated by a comma (,).
4. The pattern must be enclosed in quotes

For example for one sample (name: 1) with forward (file: 1_a.fastq.gz) and reverse (file: 1_b.fastq.gz) reads in folder data:

--input_folder "data" --extension "/*_{a,b}.fastq.gz"

Samples with less reads than this threshold at input or after trimming stop the pipeline. Using --ignore_empty_input_files or --ignore_failed_trimming ignores samples with read numbers below the threshold and lets the pipeline continue with less samples.

Ignore input files with less reads than specified by --min_read_counts and continue the pipeline without those samples.

When read sequences are trimmed, untrimmed read pairs are discarded routinely. Use this option to retain untrimmed read pairs. This is usually not recommended and is only of advantage for specific protocols that prevent sequencing PCR primers.

Cutdapt will be run twice, first to remove reads without primers (default), then a second time to remove reads that erroneously contain a second set of primers, not to be used with --retain_untrimmed.

Ignore files with less reads than specified by --min_read_counts after trimming and continue the pipeline without those samples.

Read denoising by DADA2 creates an error profile specific to a sequencing run and uses this to correct sequencing errors. This method prefers when all reads to have the same length and as high quality as possible while maintaining at least 20 bp overlap for merging. One cutoff for the forward read --trunclenf and one for the reverse read --trunclenr truncate all longer reads at that position and drop all shorter reads. If not set, these cutoffs will be determined automatically for the position before the mean quality score drops below --trunc_qmin.

For example:

--trunclenf 180 --trunclenr 120

Please note:

1. Overly aggressive truncation might lead to insufficient overlap for read merging
2. Too little truncation might reduce denoised reads
3. The code choosing these values automatically cannot take the points above into account, therefore checking read numbers is essential

Read denoising by DADA2 creates an error profile specific to a sequencing run and uses this to correct sequencing errors. This method prefers when all reads to have the same length and as high quality as possible while maintaining at least 20 bp overlap for merging. One cutoff for the forward read --trunclenf and one for the reverse read --trunclenr truncate all longer reads at that position and drop all shorter reads. If not set, these cutoffs will be determined automatically for the position before the mean quality score drops below --trunc_qmin.

For example:

--trunclenf 180 --trunclenr 120

Please note:

1. Overly aggressive truncation might lead to insufficient overlap for read merging
2. Too little truncation might reduce denoised reads
3. The code choosing these values automatically cannot take the points above into account, therefore checking read numbers is essential

Automatically determine --trunclenf and --trunclenr before the median quality score drops below --trunc_qmin. The fraction of reads retained is defined by --trunc_rmin, which might override the quality cutoff.

For example:

--trunc_qmin 35

Please note:

1. The code choosing --trunclenf and --trunclenr using --trunc_qmin automatically cannot take amplicon length or overlap requirements for merging into account, therefore use with caution.
2. A minimum value of 25 is recommended. However, high quality data with a large paired sequence overlap might justify a higher value (e.g. 35). Also, very low quality data might require a lower value.
3. If the quality cutoff is too low to include a certain fraction of reads that is specified by --trunc_rmin (e.g. 0.75 means at least 75% percent of reads are retained), a lower cutoff according to --trunc_rmin superseeds the quality cutoff.

Value can range from 0 to 1. 0 means no reads need to be retained and 1 means all reads need to be retained. The minimum lengths of --trunc_qmin and --trunc_rmin are chosen as DADA2 cutoffs.

After truncation, reads with higher than max_ee "expected errors" will be discarded. In case of very long reads, you might want to increase this value. We recommend (to start with) a value corresponding to approximately 1 expected error per 100-200 bp (default: 2)

Remove reads with length less than min_len after trimming and truncation.

Remove reads with length greater than max_len after trimming and truncation. Must be a positive integer.

Ignore files with fewer reads than specified by --min_read_counts after trimming and continue the pipeline without those samples. Please review all quality trimming and filtering options before using this parameter. For example, one sample with shorter sequences than other samples might loose all sequences due to minimum length requirements by read truncation (see --trunclenf).

If samples are treated independent (lowest sensitivity and lowest resources), pooled (highest sensitivity and resources) or pseudo-pooled (balance between required resources and sensitivity).

This parameters specifies how paired-end reads are merged after denoising. By default, read pairs will be merged by overlap. Concatenating read pairs (separated by 10 N's) is an alternative to only analyzing the forward or reverse read in case of non-overlapping paired-end sequencing data, this is not recommended, only if all other options fail. The consensus strategy is merging read pairs by overlap if possible and concatenates non-overlapping read pairs, based on --mergepairs_consensus_* parameters.

This parameter specifies the numerical value added to the alignment score for every pair of bases that match between the forward and reverse reads. A higher value increases the preference for alignments with more matching bases.

This parameter defines the numerical penalty subtracted from the alignment score for each base pair mismatch between the forward and reverse reads. A higher penalty reduces the likelihood of accepting alignments with mismatches.

This parameter sets the numerical penalty subtracted from the alignment score for each gap (insertion or deletion) introduced to align the forward and reverse reads. A higher penalty discourages alignments that require gaps.

This parameter specifies the smallest number of consecutive base pairs that must overlap between the forward and reverse reads for them to be merged. Ensuring sufficient overlap is crucial for accurate merging.

This parameter defines the highest number of mismatched base pairs permitted in the overlap region between forward and reverse reads for a successful merge. Setting this value helps control the stringency of read merging, balancing between sensitivity and accuracy.

ASVs will be clustered with VSEARCH using the id value found in --vsearch_cluster_id.

Lowering or increasing this value can change the number ASVs left over after clustering.

Setting to true adds 'ulimit -s unlimited' to the beginning of the filt_clusters.py command.

Remove ASV that are below the minimum length threshold (default: filter is disabled, otherwise 1). Increasing the threshold might reduce false positive ASVs (e.g. PCR off-targets).

Remove ASV that are above the maximum length threshold (default: filter is disabled, otherwise 1000000). Lowering the threshold might reduce false positive ASVs (e.g. PCR off-targets).

ASVs will be filtered to contain no stop codon in their coding sequence and that their length is a multiple of 3.

By default, when --filter_codons is set, the codons start from the first position of the ASV sequences. The start of the codons can be changed to any position.

By default, when --filter_codons is set, the codons are checked until the end of the ASV sequences. If you would like to change this setting, you can specify until which position of the ASV sequences the codon triplets are checked.

Please note that the length of the ASV from the beginning or from the --orf_start until this position must be a multiple of 3.

By default, when --filter_codons is set, the codons TAA,TAG are set as stop codons. Here you can specify any comma-separated list of codons to be used as stop codons, e.g. --stop_codons "TAA,TAG,TGA"

Choose any of the supported databases, and optionally also specify the version. Database and version are separated by an equal sign (=, e.g. silva=138) . This will download the desired database, format it to produce a file that is compatible with DADA2's assignTaxonomy and another file that is compatible with DADA2's addSpecies.

The following databases are supported:

• GTDB - Genome Taxonomy Database - 16S rRNA
• SBDI-GTDB, a Sativa-vetted version of the GTDB 16S rRNA
• PR2 - Protist Reference Ribosomal Database - 18S rRNA
• RDP - Ribosomal Database Project - 16S rRNA
• SILVA ribosomal RNA gene database project - 16S rRNA
• UNITE - eukaryotic nuclear ribosomal ITS region - ITS
• COIDB - eukaryotic Cytochrome Oxidase I (COI) from The Barcode of Life Data System (BOLD) - COI

Generally, using gtdb, pr2, rdp, sbdi-gtdb, silva, coidb, unite-fungi, or unite-alleuk will select the most recent supported version.

Please note that commercial/non-academic entities require licensing for SILVA v132 database (non-default) but not from v138 on (default).

Overwrites --dada_ref_taxonomy. Either --skip_dada_addspecies (no species annotation) or --dada_ref_tax_custom_sp (species annotation) is additionally required. Consider also setting --dada_assign_taxlevels.

Must be compatible to DADA2's assignTaxonomy function: 'Can be compressed. This reference fasta file should be formatted so that the id lines correspond to the taxonomy (or classification) of the associated sequence, and each taxonomic level is separated by a semicolon.' See also https://rdrr.io/bioc/dada2/man/assignTaxonomy.html

Requires --dada_ref_tax_custom. Must be compatible to DADA2's addSpecies function: 'Can be compressed. This reference fasta file should be formatted so that the id lines correspond to the genus-species binomial of the associated sequence.' See also https://rdrr.io/bioc/dada2/man/addSpecies.html

Typically useful when providing a custom DADA2 reference taxonomy database with --dada_ref_tax_custom. If DADA2's addSpecies is used (default), the last element(s) of the comma separated string must be 'Genus' or 'Genus,Species'.

Expected amplified sequences are extracted from the DADA2 reference taxonomy using the primer sequences, that might improve classification. This is not applied to species classification (assignSpecies) but only for lower taxonomic levels (assignTaxonomy).

Defines the behavior when multiple exact matches against different species are returned. By default only unambiguous identifications are returned. If TRUE, a concatenated string of all exactly matched species is returned.

Reverse-complement of each sequences will be used for classification if it is a better match to the reference sequences than the forward sequence.

Chunks for DADA2's assignTaxonomy and addSpecies can speed up the process and lower required memory.

Headerless, tab-separated, first column with tree leaves, second column with taxonomy ranks separated by semicolon ;. The results take precedence over DADA2 and QIIME2 classifications.

Choose any of the supported databases, and optionally also specify the version. Database and version are separated by an equal sign (=, e.g. silva=138) . This will download the desired database and initiate taxonomic classification with QIIME2 and the chosen database.

If both, --dada_ref_taxonomy and --qiime_ref_taxonomy are used, DADA2 classification will be used for downstream analysis.

The following databases are supported:

• SILVA ribosomal RNA gene database project - 16S rRNA
• UNITE - eukaryotic nuclear ribosomal ITS region - ITS
• Greengenes (only testing!)

Generally, using silva, unite-fungi, or unite-alleuk will select the most recent supported version. For testing purposes, the tiny database greengenes85 (dereplicated at 85% sequence similarity) is available. For details on what values are valid, please either use an invalid value such as x (causing the pipeline to send an error message with all valid values) or see conf/ref_databases.config.

Overwrites --qiime_ref_taxonomy. Either path to tarball (*.tar.gz or *.tgz) that contains sequence (*.fna) and taxonomy (*.tax) data, or alternatively a comma separated pair of filepaths to sequence (*.fna) and taxonomy (*.tax) data (possibly gzipped *.gz).

If you have trained a compatible classifier before, from sources such as SILVA (https://www.arb-silva.de/), Greengenes (http://greengenes.secondgenome.com/downloads) or RDP (https://rdp.cme.msu.edu/).

For example:

--classifier "FW_primer-RV_primer-classifier.qza"

Please note the following requirements:

1. The path must be enclosed in quotes
2. The classifier is a Naive Bayes classifier produced by qiime feature-classifier fit-classifier-naive-bayes (e.g. by this pipeline)
3. The primer pair for the amplicon PCR and the computing of the classifier are exactly the same (or full-length, potentially lower performance)
4. The classifier has to be trained by the same version of scikit-learn as this version of the pipeline uses

Choose any of the supported databases, and optionally also specify the version. Database and version are separated by an equal sign (=, e.g. silve=138) . This will download the desired database and initiate taxonomic classification with Kraken2 and the chosen database.

Consider using --kraken2_confidence to set a confidence score threshold.

The following databases are supported:

• RDP - Ribosomal Database Project - 16S rRNA
• SILVA ribosomal RNA gene database project - 16S rRNA
• Greengenes - 16S rRNA
• Standard Kraken2 database (RefSeq archaea, bacteria, viral, plasmid, human, UniVec_Core) - any amplicon

Generally, using rdp, silva, greengenes, standard will select the most recent supported version.

Please note that commercial/non-academic entities require licensing for SILVA v132 database (non-default) but not from v138 on.

Overwrites --kraken2_ref_taxonomy. Consider also setting --kraken2_assign_taxlevels. Can be compressed tar archive (.tar.gz|.tgz) or folder containing the database. See also https://benlangmead.github.io/aws-indexes/k2.

Typically useful when providing a custom Kraken2 reference taxonomy database with --kraken2_ref_tax_custom. In case a database is given with --kraken2_ref_taxonomy, the default taxonomic levels will be overwritten with --kraken2_assign_taxlevels.

Increasing the threshold will require more k-mers to match at a taxonomic levels and reduce the taxonomic levels shown until the threshold is met.

Choose any of the supported databases, and optionally also specify the version. Database and version are separated by an equal sign (=, e.g. coidb=221216) . This will download the desired database and initiate taxonomic classification with VSEARCH sintax and the chosen database, which if needed is formatted to produce a file that is compatible with VSEARCH sintax.

The following databases are supported:

• COIDB - eukaryotic Cytochrome Oxidase I (COI) from The Barcode of Life Data System (BOLD) - COI
• UNITE - eukaryotic nuclear ribosomal ITS region - ITS

Generally, using coidb, unite-fungi, or unite-alleuk will select the most recent supported version.

If data is long read ITS sequences, that need to be cut to ITS region (full ITS, only ITS1, or only ITS2) for taxonomy assignment.

If using cut_its, this option allows partial ITS sequences, longer than the specified cutoff.

Use with --sidle_ref_seq_custom. Consider also setting --sidle_ref_aln_custom and --sidle_ref_tree_custom. The taxonomy file must be headerless and follow the format of qiime tools import's HeaderlessTSVTaxonomyFormat. Example usage: --sidle_ref_tax_custom 'taxonomy_99_taxonomy.txt'

Use with --sidle_ref_tax_custom. Example usage: --sidle_ref_seq_custom 'rep_set_99.fasta'

May be used with --sidle_ref_tax_custom. Allows sequence reconstruction. Example usage: --sidle_ref_aln_custom 'rep_set_aligned_99.fasta'

Use with --sidle_ref_aln_custom. Recommended with --sidle_ref_tax_custom. Allows phylogenetic tree reconstruction and therefore diversity analysis. Overwrites tree chosen by --sidle_ref_taxonomy

Only effective with --sidle_ref_tax_custom. Recommended with --sidle_ref_tax_custom, default 5 was recommended with SILVA 128. Sets --p-num-degenerates for qiime rescript cull-seqs.

Recommended with --sidle_ref_tax_custom, default '--p-database silva' was recommended with SILVA 128. Therefore, ad-hoc database cleaning will be performed automatically, specifically with regard to the define-missing and ambiguity-handling parameters. Overwrites recommended settings for database chosen by --sidle_ref_taxonomy

Depending on the primers used, PCR might amplify unwanted or off-target DNA. By default sequences originating from mitochondria or chloroplasts are removed. The taxa specified are excluded from further analysis. For example to exclude any taxa that contain mitochondria, chloroplast, or archaea:

--exclude_taxa "mitochondria,chloroplast,archaea"

If you prefer not filtering the data, specify:

--exclude_taxa "none"

Please note the following requirements:

1. Comma separated list enclosed in quotes
2. May not contain whitespace characters
3. Features that contain one or several of these terms in their taxonomical classification are excluded from further analysis
4. The taxonomy level is not taken into consideration
5. Taxa names should be as in Taxonomic database (Default: Silva138), example: 'Bacteria', 'Armatimonadia', 'unidentified', 'p__'
6. Taxon names are case-insensitive and partial match is possible.

Remove entries from the feature table below an absolute abundance threshold (default: 1, meaning filter is disabled). Singletons are often regarded as artifacts, choosing a value of 2 removes sequences with less than 2 total counts from the feature table.

For example to remove singletons choose:

--min_frequency 2

Filtering low prevalent features from the feature table, e.g. keeping only features that are present in at least two samples can be achived by choosing a value of 2 (default: 1, meaning filter is disabled). Typically only used when having replicates for all samples.

For example to retain features that are present in at least two sample:

--min_samples 2

Please note this is independent of abundance.

Here columns in the metadata sheet can be chosen with groupings that are used for diversity indices and differential abundance analysis. By default, all suitable columns in the metadata sheet will be used if this option is not specified. Suitable are columns which are categorical (not numerical) and have multiple different values which are not all unique. For example:

--metadata_category "treatment1,treatment2"

Please note the following requirements:

1. Comma separated list enclosed in quotes
2. May not contain whitespace characters
3. Each comma separated term has to match exactly one column name in the metadata sheet

Here columns in the metadata sheet can be chosen with groupings that are used for average relative abundance barplots. Samples that have empty fields for that column are discarded. For example:

--metadata_category_barplot "treatment1,treatment2"

Please note the following requirements:

1. Comma separated list enclosed in quotes
2. May not contain whitespace characters
3. Each comma separated term has to match exactly one column name in the metadata sheet

Comma separated list of model formula(s), e.g. "treatment1,treatment2". Model formula should contain only independent terms in the sample metadata. These can be continuous variables or factors, and they can have interactions as in a typical R formula. Essentially, columns in the metadata sheet can be chosen that have no empty values, not only unique values, or not only identical values. For example, "treatment1+treatment2" tests whether the data partitions based on "treatment1" and "treatment2" sample metadata. "treatment1*treatment2" test both of those effects as well as their interaction. More examples can be found in the R documentation, https://cran.r-project.org/doc/manuals/r-release/R-intro.html#Formulae-for-statistical-models

Depends on the reference taxonomy database used.

Depends on the reference taxonomy database used. Most default databases have genus level at 6.

ANCOMBC will be performed on all suitable columns in the metadata sheet. Empty values will be removed, therefore it is possible to perform tests on subsets. The reference level will default to highest alphanumeric group (e.g. in alphabetical or numeric order, as applicable) within each metadata column. Formula for specific tests can be supplied with --ancombc_formula.

Comma separated list of model formula(s), e.g. "treatment1,treatment2". The reference level will default to highest alphanumeric group (e.g. in alphabetical or numeric order, as applicable) within each formula term. The reference level can be overwritten by --ancombc_formula_reflvl. Model formula should contain only independent terms in the sample metadata. These can be continuous variables or factors, and they can have interactions as in a typical R formula. Essentially, columns in the metadata sheet can be chosen that have no empty values, not only unique values, or not only identical values. For example, "treatment1+treatment2" tests whether the data partitions based on "treatment1" and "treatment2" sample metadata. "treatment1*treatment2" test both of those effects as well as their interaction. More examples can be found in the R documentation, https://cran.r-project.org/doc/manuals/r-release/R-intro.html#Formulae-for-statistical-models

This will only affect ANCOM-BC started by --ancombc_formula, but for all provided model formula, therefore it might be best to restrict --ancombc_formula to one formula. The syntax is as follows: 'column_name::column_value' or for multiple 'column_name1::column_value1 column_name2::column_value2'