MergeRibo

Script: mergeribo.sh Package: prok Class: MergeRibo.java

Merges ribosomal sequence files to keep ranked members per TaxID. By default, a consensus is generated per TaxID, then the sequence best matching that consensus is used: First, all sequences per TaxID are aligned to a reference consensus. Second, the best-matching sequence is used as a seed, and all sequences for that TaxID are aligned to the seed to generate a new consensus. Third, in 'consensus' mode, that consensus is simply output. In 'best' mode (default), all sequences are aligned again to the new consensus, and up to maxpertaxid members are output (one by default). All copies, including duplicates, contribute to the per-TaxID consensus. Optional exact deduplication occurs after ranking and before selecting outputs.

Basic Usage

mergeribo.sh in=<file,file> out=<file>

MergeRibo processes ribosomal RNA sequences to produce ranked representative sequences per taxonomic ID. Supported modes are 16S (default), 18S, ITS, LSU (23S/25S/26S/28S), 5S, and 5.8S. This is particularly useful for creating non-redundant databases or consensus sequences from multiple rRNA sequences.

Parameters

Parameters are organized by their function in the sequence merging process. The tool follows a multi-step algorithm: alignment to reference consensus, seed selection, individual consensus generation, and best sequence selection.

Standard parameters

in=<file,file>
Comma-delimited list of input FASTA files. Headers may contain tid|N, tid_N, ncbi|N or ncbi_N anywhere in the header. Invalid or missing positive TaxIDs are skipped. Inputs are unpaired.
out=<file>
Output FASTA file. Contains up to maxpertaxid selected sequences per taxID (one by default).
out2=<file>
Not used by MergeRibo, which reads unpaired input and writes only out=.
overwrite=t
(ow) Set to false to force the program to abort rather than overwrite an existing file.
t=auto
Number of worker threads.
reads=-1
Maximum reads from each input file; -1 means unlimited.
verbose=f
Print additional loading, alignment and selection diagnostics.
ordered=f
Use the ordered output writer. Taxon order still depends on worker completion, rather than input or taxonomic order.
showspeed=t
(ss) Set to 'f' to suppress display of processing speed statistics during execution.
ziplevel=2
(zl) Set to 1 (lowest) through 9 (max) to change compression level; lower compression is faster.
fastawrap=70
FASTA line wrap length. 4000 is recommended to minimize filesize.

Taxonomy parameters

tree=f
Load a TaxTree (default f); tree=t, tree=auto, or tree=<path> enables it.
usetree=f
Alias for tree=f. dada2=t or level=<level> requires a tree; explicit f fails.
treefile=auto
Select the TaxTree source without changing the local gate.
TaxIDs (header format)
Not a parameter. Headers may contain tid|N, tid_N, ncbi|N or ncbi_N anywhere in the header. Invalid or missing positive TaxIDs are skipped.

Processing parameters

alt=<file>
Lower priority data, used only for TaxIDs with no accepted representative in the primary input.
best=t
Output ranked real representatives per TaxID. Rank is identity times min(length,ideal)/max(length,ideal), with longer sequences first on ties. The FIRST reference sets the ideal length, even when another reference matches better.
maxpertaxid=1
Maximum representatives per TaxID; must be at least 1.
dedupe=f
Remove exact duplicate sequences within each TaxID AFTER consensus building and ranking, keeping the highest-ranked copy. Use maxpertaxid=10 dedupe=t for up to 10 distinct copies. With dedupe=f, repeated sequences may occupy multiple outputs.
consensus=f
Output a consensus per taxID instead of the best input sequence. Mutually exclusive with best; requires maxpertaxid=1. Groups with fewer than 3 members output their top-ranked real sequence instead (see Consensus vs Best Mode).
fast=f
Rank against the reference set without building a per-TaxID consensus; still honors maxpertaxid and dedupe. Fast mode and groups with fewer than 3 members use global-seed ranking; larger non-fast groups use their own consensus.
minid=0.62
Ignore sequences with identity lower than this to the global consensus (the reference set; each sequence uses its maximum identity across references).
maxns=-1
Ignore sequences with more than this many N bases, if non-negative. Set to -1 to disable filtering.
minlen=1
Ignore sequences shorter than this length in bases. Mode defaults: 1 for 16S, 18S, and ITS; 80 for 5S; 50 for 5.8S. An explicit minlen overrides the mode default in any argument order.
maxlen=4000
Ignore sequences longer than this length in bases. Mode defaults: 4000 for 16S, 18S, and ITS; 6000 for LSU; 200 for 5S; 300 for 5.8S. An explicit maxlen overrides the mode default in any argument order. Typical full-length 16S is ~1500bp.
16S=t
(process16S) Use the default 16S reference, or ref= if supplied. Modes are mutually exclusive; the last enabled mode wins.
18S=f
(process18S) Use the default 18S reference, or ref= if supplied. Modes are mutually exclusive; the last enabled mode wins.
ITS=f
(processITS) Use the default ITS lineage references (ITS_fungi, ITS_plant, ITS_animal, ITS_other), or ref= if supplied. Modes are mutually exclusive; the last enabled mode wins.
LSU=f
(processLSU; also 23S, 25S, 26S, 28S) Large-subunit mode. Aligns to the prokaryotic 23S consensus, or to ref= if given (required for eukaryotic LSU, which is too divergent from 23S). Sets maxlen=6000 unless overridden. Modes are mutually exclusive; the last enabled mode wins.
5S=f
(process5S) Named 5S mode; use ALL shipped 5S reference records, or ALL records in ref= if supplied. Defaults to minlen=80 maxlen=200. Modes are mutually exclusive; the last enabled mode wins.
5.8S=f
(process5.8S; also 58S) 5.8S mode; requires ref= (no configured MergeRibo default). Sets minlen=50 maxlen=300 unless overridden. Modes are mutually exclusive; the last enabled mode wins.
ref=<file>
Override the reference set in EVERY mode, including 16S, 18S, and ITS. ALL records are used; each input receives its maximum identity across them. The first record still sets the ideal length for ranking. Empty or unreadable reference sets are errors.
level=
(taxlevel) If set to a term like 'species' or 'genus', nodes will be promoted to that level, minimum, before consensus generation. Requires taxonomy tree; the tree loads automatically unless tree=f is explicitly set, which is an error.
dada2=f
Output headers in DADA2 format with complete taxonomic lineages instead of simple sequence IDs. Requires a taxonomy tree; the tree loads automatically unless tree=f is explicitly set, which is an error.

Java Parameters

-Xmx
This will set Java's memory usage, overriding autodetection. -Xmx20g will specify 20 gigs of RAM, and -Xmx200m will specify 200 megs. The max is typically 85% of physical memory.
-eoom
This flag will cause the process to exit if an out-of-memory exception occurs. Requires Java 8u92+.
-da
Disable assertions.

Examples

Basic SSU Merging

mergeribo.sh in=ssu_sequences.fa out=merged_ssu.fa

Merge SSU sequences from a single file, producing one representative 16S sequence per taxonomic ID.

Multiple Input Files

mergeribo.sh in=file1.fa,file2.fa,file3.fa out=consensus.fa minid=0.8

Process multiple input files with higher identity threshold for more stringent filtering.

18S Processing

mergeribo.sh in=euk_ssu.fa out=merged_18s.fa 18S=t 16S=f

Process 18S rRNA sequences from eukaryotes instead of bacterial 16S.

5S Mode with Multiple Ranked Outputs

mergeribo.sh in=5s.fa out=best5s.fa 5S=t maxpertaxid=10 dedupe=t

Select up to 10 distinct ranked 5S sequences per TaxID. Uses all shipped 5S reference records and the 5S default length window (minlen=80, maxlen=200). Duplicates are removed after consensus building and ranking.

5.8S Mode with a Custom Reference

mergeribo.sh in=its_region.fa out=reps58.fa 5.8S=t ref=ref58.fa

5.8S mode requires ref=. All records in ref58.fa are used, and the 5.8S length window (minlen=50, maxlen=300) applies unless overridden.

Consensus Mode with Quality Filtering

mergeribo.sh in=raw_ssu.fa out=consensus_ssu.fa consensus=t best=f minlen=1200 maxlen=1600 maxns=10

Generate consensus sequences instead of best representatives, with length filtering for full-length 16S and N-base filtering. Groups with fewer than 3 members output their top-ranked real sequence instead of a consensus.

Species-Level Grouping with DADA2 Output

mergeribo.sh in=amplicons.fa out=species_reps.fa level=species dada2=t

Group sequences at species level and output with DADA2-compatible taxonomic headers.

Algorithm Details

Multi-Stage Consensus Generation

MergeRibo implements a multi-stage algorithm for selecting representative sequences per taxonomic ID:

  1. Reference Identity Phase: Each input sequence is filtered by length (minlen, maxlen) and N count (maxns), then aligned against every reference in the selected reference set using SingleStateAlignerFlat2. Each sequence keeps its maximum identity across references. Sequences below the minid threshold are filtered out. Sequences with a valid TaxID are grouped by TaxID (ProcessThread.processRead() and align(), MergeRibo.java lines 706-747).
  2. Seed Selection: Each TaxID group is sorted by its rank product (identity times the length multiplier), with longer sequences first on ties (Ribo.compareTo(), lines 789-795). Groups with fewer than 3 members, and all groups in fast mode, return the top-ranked member directly with no consensus (pickBestInner(), line 670).
  3. Consensus Generation: For non-fast groups with 3 or more members, the top-ranked sequence seeds a BaseGraph, padded toward the full gene length (900 for ITS, 1600 otherwise). Every member of the group, including the seed and identical copies, is aligned to the graph and added as a vote. The graph is then traversed to produce the consensus. BaseGraph settings are set in the constructor (lines 89-94): MAF_sub=0.251, MAF_del=0, MAF_ins=0, MAF_noref=0, trimDepthFraction=0.3, trimNs=true.
  4. Output Selection: In best mode, every member is re-aligned to the consensus, re-scored, and re-sorted; the top maxpertaxid real members are output. In consensus mode, the consensus sequence itself is output as one sequence (consensus=t requires maxpertaxid=1).

Scoring Algorithm

The tool uses a length-weighted identity score implemented in the score(int len, float identity) method: score = lengthMult(length) × identity. The lengthMult(int len) method calculates a penalty ratio as min(len, idealLength) / max(len, idealLength, 1), where idealLength is the length of the first loaded reference sequence (referenceSequences[0]). This penalizes sequences that deviate from the expected gene length; for example, a sequence half the ideal length has a multiplier of 0.5.

Memory Management

The implementation uses several memory optimization strategies:

Consensus vs Best Mode

In best mode (default), the output sequences are real input sequences, ranked by their identity to the generated consensus; their sequence and header are not replaced (the DADA2 header conversion excepted). In consensus mode, the BaseGraph consensus sequence is output for groups with 3 or more members in non-fast mode. Groups with fewer than 3 members, and all groups in fast mode, output their top-ranked real input sequence instead, because no consensus is built for them (pickBestInner(), line 670).

Fast Mode

When fast=t, members are ranked by identity to the reference set and the top maxpertaxid are output, honoring dedupe, without building a per-TaxID consensus (pickBestInner(), line 670; pickBestN(), lines 603-630). Fast mode may produce less accurate representatives for highly divergent sequences within a TaxID, since no per-TaxID consensus is used to rank them.

Taxonomic Level Processing

When the level parameter is specified, the processRead() method calls tree.getIdAtLevelExtended(key, taxLevelE) to promote sequences to the specified taxonomic level before consensus generation. This uses the TaxTree implementation to traverse taxonomic hierarchies. The taxonomy tree is loaded by TaxTree.loadTaxTree(). dada2=t and level= enable the tree automatically; an explicit tree=f or usetree=f with either option causes an error (validateTreeRequirement(), lines 246-254).

Mode Selection

Exactly one mode is active. Enabling a mode disables all others, so the last enabled mode wins. Mode aliases: process16S, process18S, processITS, processLSU (also 23S, 25S, 26S, 28S), process5S, and process5.8S (also 58S). Default length bounds are applied after all arguments are parsed (applyModeLengthDefaults(), lines 233-243), so explicit minlen and maxlen values are never overridden.

Reference Sequence Loading

Without ref=, references are loaded by ProkObject.loadConsensusSequenceType(): "16S", "18S", and "23S" (LSU) load only the first record; "5S" loads all shipped records. ITS loads the first record of each of ITS_fungi, ITS_plant, ITS_animal, and ITS_other that is present. 5.8S has no default and requires ref=. With ref=, every record in the file is loaded in every mode (loadAllSequences()). References are stored as byte[][] referenceSequences and used by SingleStateAlignerFlat2 for identity filtering and scoring. An empty or unreadable reference set is an error (loadReferences(), lines 924-954).

Support

For questions and support: