RestoreBases

Overview

Many short and long-read aligners output secondary and supplementary alignments with SEQ=* and QUAL=*, reducing them to position-only records. This makes them unsuitable for variant calling tools like CallVariants, which require actual base sequences to evaluate evidence.

RestoreBases solves this by matching each read by name, finding its primary (non-supplementary) alignment which carries the true bases and qualities, and copying those fields onto all secondary/supplementary records. The restoration is reference-independent: you get the actual sequenced bases, not bases reconstructed from the MD tag, which would inherit any reference bias in the alignment.

Key advantage: For a true heterozygous variant at allele frequency 0.50, MD-tag reconstruction drops the estimated AF to ~0.17 due to reference bias. RestoreBases preserves the true AF because it uses actual sequenced bases.

Basic Usage

restorebases.sh in=input.bam out=output.bam

The input file must be a SAM or BAM file containing secondary and/or supplementary alignments. RestoreBases will read all records, match them by name, restore bases/quals onto non-primary alignments from their primary counterpart, and write a complete SAM/BAM output with correct SO:unsorted header tag.

How It Works

RestoreBases uses a two-phase approach to bound memory usage:

1. Partition Phase: Stream the input once and partition records into ways temporary headerless sam.gz subfiles, routed by qname.hashCode()%ways. This ensures all alignments of a read land in the same subfile.
2. Restore Phase: For each subfile: load into memory, name-sort, group by (qname, pairnum), then for each group copy the primary's seq/qual onto its secondaries/supplementaries. Write results to the final output.

With ways=1, partitioning is skipped (no temp files)—the entire input is loaded, sorted, and restored at once. This is recommended for bacterial BAMs and other inputs that fit comfortably in RAM.

Handling Hard-Clipped Supplementary Alignments

Some aligners hard-clip supplementary records (H in CIGAR), meaning the record's SEQ does not span the full read. RestoreBases offers two strategies in fix mode:

• convert (default): Convert H→S so the CIGAR spans the full read, then copy the full seq/qual. Larger output, but simpler semantics.
• truncate: Keep the H cigar unchanged and copy only the sub-range of the read that this record's own segment covers. Gives identical variant calls but smaller output (avoids quadratic blowup).

Note: if the primary itself has hard-clips, that group cannot be restored (the full read is not available). Such groups are counted in the statistics.

Parameters

RestoreBases organizes parameters into input/output, processing mode, and memory/performance categories.

Common Use Cases

Preparing minimap2 Output for Variant Calling

minimap2 is a fast long-read aligner that emits supplementary alignments for reads mapping to multiple loci. Its secondary/supplementary records carry SEQ=* to save space:

minimap2 -a -x map-ont ref.fa reads.fastq > reads.sam
samtools view -b reads.sam > reads.bam
samtools sort reads.bam > reads.sorted.bam
restorebases.sh in=reads.sorted.bam out=reads.restored.bam
callvariants.sh in=reads.restored.bam ref=ref.fa out=vars.vcf

After restoration, CallVariants can use all alignments (primary, secondary, supplementary) for evidence, improving variant discovery on repetitive regions and multi-mapping reads.

Restore with Truncate Mode for Space Efficiency

If output file size is a concern and you have supplementary hard-clips, use hardclip=truncate to copy only the relevant segment of each read:

restorebases.sh in=multi.bam out=restored.bam hardclip=truncate

This avoids quadratic seq blowup (where a single read might be repeated N times in different supplementary records) while still providing all bases needed for variant calling.

Testing and Debugging

Use the reads parameter to process a small subset for validation:

restorebases.sh in=large.bam out=test.bam reads=1000

Check the output and statistics before running on the full file.

Large-Scale Human Genome Restoration

For whole-genome BAMs that don't fit in memory, tune ways to control per-subfile memory:

restorebases.sh in=NA12878.bam out=NA12878.restored.bam ways=100 tmpdir=/scratch/restbases_tmp -Xmx8000m

Increase ways to reduce per-subfile memory and prevent OOM. Temporary files are created and cleaned up automatically.

Attaching Tags Instead of Overwriting SEQ/QUAL

In tag mode, the original SEQ=* is preserved and the true bases appear in OS:Z: and OQ:Z: tags:

restorebases.sh in=reads.bam out=reads.tagged.bam mode=tag

Useful if you want to preserve the original BAM structure but have the true bases available for downstream analysis.

Algorithm and Implementation Details

Read Grouping and Matching

RestoreBases groups alignments by (qname, pairnum)—i.e., read name and pair number (0 for single-end or R1, 1 for R2). Within each group, the first record after name-sorting is assumed to be the primary alignment (non-supplementary, non-secondary, with valid bases). All subsequent records in that group are treated as secondary or supplementary.

If no usable primary is found (e.g., all records are supplementary, or the primary has SEQ=*), the entire group is skipped and counted in "Groups with no primary" statistics.

Multi-Way Partitioning Strategy

The partition phase uses qname.hashCode()%ways to route each record to a subfile, ensuring all alignments of a read go to the same subfile. This allows each subfile to be processed independently with bounded memory. Records are buffered in chunks (200 per buffer by default) before writing to reduce I/O overhead.

Hard-Clip Handling in Truncate Mode

In truncate mode, the code extracts the leading hard-clip count from the CIGAR, calculates the query length (soft+hard-inclusive bases), and determines the trailing hard-clip as the remainder. For minus-strand records, the CIGAR is reverse-complemented relative to the read in the forward strand, so the segment start is calculated accordingly. Only the relevant portion of the primary's seq/qual is copied to avoid duplication.

Reference Bias Avoidance

Unlike MD-tag reconstruction (which reconstructs bases by comparing the SEQ to the MD tag and reference), RestoreBases copies the actual sequenced bases from the primary record. This preserves the true allele frequencies: a 50/50 het remains 50/50, not degraded to ~17% as would happen with MD-tag restoration on a misaligned supplementary.

Header Patching

The output header is modified to include or replace the SO: (sort order) field with SO:unsorted, indicating that records are not in any standard sort order. If no @HD line exists, one is created with VN:1.6 and SO:unsorted.

Statistics and Output

RestoreBases prints a summary of work performed:

• Mode: fix or tag
• Hardclip: convert or truncate
• Ways: Number of subfiles (partitioning factor)
• Alignments in: Total input records processed
• Restore attempts: Groups where restoration was attempted
• Records restored: Actual successful restorations
• Hardclips converted: Number of H→S conversions (fix/convert mode)
• Hardclips truncated: Number of truncations applied (fix/truncate mode)
• Skipped (hardclipped primary): Groups where the primary itself was hard-clipped (cannot restore)
• Skipped (length mismatch): Groups where supplementary CIGAR length doesn't match primary seq length
• Skipped (no cigar/unmapped): Supplementaries without a valid CIGAR or marked unmapped
• Groups with no primary: Read groups with no usable primary alignment

Performance Considerations

Memory Management

The ways parameter directly controls peak memory usage. Each subfile is held in memory sequentially, so peak memory ≈ (total input size) / ways. For a 100 GB BAM with ways=31, expect ~3.2 GB per subfile.

Speed

RestoreBases performs I/O-bound streaming and sorting, with computational overhead from SAM parsing and name-sorting. For large inputs, temporary file I/O dominates. Using an SSD or high-speed storage for tmpdir significantly improves performance.

Disk Space

During multi-way processing, temporary files consume space roughly equal to the input size. Ensure tmpdir has sufficient free space. Temporary files are cleaned up immediately after each subfile is processed, so total temp space never exceeds input size.

Compressed vs. Uncompressed Output

BAM (binary, compressed) output is smaller and faster to read than SAM (text). Use out=output.bam unless you specifically need text format.

Troubleshooting

Out-of-Memory Errors

If you encounter "java.lang.OutOfMemoryError", either increase the heap size with -Xmx or increase ways to reduce per-subfile memory. For example:

restorebases.sh in=input.bam out=output.bam ways=100 -Xmx8000m

Unexpected Output Size

In fix/convert mode, hard-clips are converted to soft-clips and the full read is copied, which can increase output file size. If this is undesirable, use hardclip=truncate instead.

Partial Restoration

If the statistics show a high number of "Skipped (length mismatch)" records, the input BAM may be corrupted or the supplementary alignments were generated by a non-standard aligner. Check that the CIGAR strings are valid and that the primary and supplementary records truly match.

Temporary Files Not Cleaned Up

If the process crashes, temporary files in tmpdir may remain. You can safely delete them manually; they are headerless SAM.gz files prefixed part_*.sam.gz.

Support

For questions and support:

• Email: bbushnell@lbl.gov
• Documentation: bbmap.org
• GitHub: BBTools Repository