The Role of BAM Preprocessing in Variant Calling: Marking Duplicates, Base Quality Recalibration, and Indel Realignment

By Dr. Zubair Khalid, DVM, MS, PhD ·

The Role of BAM Preprocessing in Variant Calling: Marking Duplicates, Base Quality Recalibration, and Indel Realignment

Key Takeaways

  • BAM preprocessing is critical for reducing systematic artifacts that mimic true genetic variation, ensuring variant calls reflect biology rather than library preparation or sequencing noise. Key steps include marking duplicates, base quality recalibration (BQSR), and indel realignment.
  • Marking duplicates addresses PCR and optical artifacts by flagging reads originating from the same DNA fragment, preventing inflated read depth and false variant calls, particularly in WGS, WES, and targeted panels that utilize PCR amplification.
  • Base Quality Recalibration (BQSR) corrects systematic errors in base quality scores, which are often miscalibrated in specific sequence contexts (e.g., read ends, CpG motifs, homopolymers), by comparing observed base calls to known variant sites.
  • Indel realignment, largely superseded by modern aligners and variant callers like GATK HaplotypeCaller, corrects local alignment artifacts around insertions and deletions that can manifest as false SNVs.
  • The decision to skip preprocessing steps depends on the sequencing strategy and variant caller; for instance, UMI-based protocols obviate duplicate marking, and DeepVariant does not require BQSR as it models errors internally.
  • Reproducibility in preprocessing is achieved through version-controlled tools, containerization (e.g., nf-core pipelines), and workflow managers (e.g., Snakemake, Nextflow), ensuring consistent results across different computational environments.

BAM preprocessing is the sequence of file-level operations applied to aligned sequencing reads before variant calling. These steps, primarily marking duplicates, base quality recalibration, and indel realignment, exist to reduce systematic artifacts that mimic true genetic variation. For researchers moving from raw sequencing data to a list of candidate variants, preprocessing determines whether the final call set reflects biology or library preparation noise. This article explains what each preprocessing step does, why it matters for germline and somatic calling, when it can be skipped, and how to implement it using GATK tools within reproducible workflows.

The scope here covers DNA sequencing data, specifically whole-genome sequencing (WGS), whole-exome sequencing (WES), and targeted amplicon panels. The reader is assumed to have aligned reads in BAM format and to be deciding which preprocessing operations to run before invoking a variant caller. The practical outcome is a defensible preprocessing strategy that matches your study design, sequencing platform, and downstream analysis goals.

The Purpose of BAM Preprocessing in the Variant Calling Workflow

Variant calling begins with a BAM file that contains read alignments to a reference genome. Between alignment and variant calling, the BAM file can be modified to correct for known sources of error. These errors fall into two broad categories: those introduced during library preparation and sequencing, and those introduced during the alignment step itself.

Library preparation errors include PCR duplicates, where multiple sequencing reads originate from the same original DNA fragment. During PCR amplification, fragments are copied, and the sequencer may read the same fragment multiple times. If these duplicate reads are counted as independent observations, they inflate read depth at specific positions and can cause false variant calls, particularly for variants that were present in the original fragment before amplification.

Sequencing errors include base quality score miscalibration. Every base call from a sequencer carries a quality score that estimates the probability of error. These scores are generated by the instrument's software and can be systematically wrong, especially in contexts such as the ends of reads, motifs like CpG dinucleotides, or runs of identical bases. If the variant caller trusts these scores, systematic errors can be mistaken for true variants.

Alignment errors include reads that map incorrectly around insertions or deletions (indels). When a read spans an indel, the alignment algorithm may place mismatching bases near the indel boundary. These mismatches look like single nucleotide variants (SNVs) to the caller. Realignment corrects these local alignment artifacts.

The preprocessing steps address these three error sources. Marking duplicates flags reads that are likely PCR or optical duplicates. Base quality recalibration adjusts quality scores to match empirical error rates. Indel realignment, now largely superseded by modern aligners and callers, locally realigns reads to minimize mismatches around indels.

The decision to run these steps depends on the variant caller, the sequencing strategy, and the biological question. The Genome Analysis Toolkit (GATK) best practices workflow includes duplicate marking and base quality recalibration for germline short variant discovery. Somatic workflows, such as those used in cancer genomics, may modify these steps because tumor samples often have low purity and high heterogeneity.

At a Glance: Preprocessing Steps and Their Role

Preprocessing StepPrimary Artifact AddressedTypical ApplicationWhen to Consider Skipping
Mark DuplicatesPCR and optical duplicates inflating read depthWGS, WES, and targeted panels with PCR amplificationUnique molecular identifier (UMI) based protocols where duplicates are collapsed during consensus generation
Base Quality Recalibration (BQSR)Systematic base quality score errorsGermline WGS and WES with GATK HaplotypeCallerTargeted panels with known high-quality reference calls, or when using callers that model base quality errors internally
Indel RealignmentMismatches near indel boundaries causing false SNVsLegacy pipelines with older alignersModern aligners and variant callers that perform local realignment internally, such as GATK HaplotypeCaller

The table above summarizes the three core preprocessing operations. Each step targets a distinct artifact class, and the decision to run or skip a step should be based on the specific error profile of your data and the capabilities of your variant caller.

Marking Duplicates: Why Duplicate Reads Distort Variant Calls

Duplicate reads are sequencing reads that originate from the same DNA fragment. During library preparation, PCR amplification creates multiple copies of each fragment. The sequencer then reads these copies, producing reads that align to the same genomic coordinates. Optical duplicates arise when the sequencer's imaging system misidentifies clusters, producing reads that appear identical.

The problem with duplicates is statistical. Variant callers use read depth and allele fraction to determine whether a variant is real. If a library contains a PCR error in one original fragment, amplification creates many copies of that error. The sequencer reads these copies, and the variant caller sees high depth at the error position with a consistent alternate allele. Without duplicate marking, this artifact looks like a true heterozygous or homozygous variant.

Duplicate marking identifies reads that share the same alignment start position and orientation. The marking algorithm assigns a flag to all but one read in each duplicate set. The variant caller then ignores flagged reads or treats them as a single observation.

For whole-genome and whole-exome sequencing, duplicate marking is standard practice. The GATK best practices workflow includes MarkDuplicates as a required step before variant calling. The tool calculates metrics such as the duplicate rate, which indicates library complexity. A high duplicate rate suggests that the library has low complexity, meaning many fragments were lost during preparation and the remaining fragments were over-amplified.

For targeted amplicon panels, duplicate marking requires careful consideration. Amplicon-based protocols generate reads that all start at the same position because the primers define the fragment boundaries. These reads are not true duplicates in the biological sense, but they share alignment coordinates. Marking them as duplicates would discard most of the data. Instead, amplicon pipelines often use read family information or primer coordinates to distinguish true duplicates from independent reads.

UMI-based protocols offer an alternative. Unique molecular identifiers are short random sequences attached to each fragment before amplification. Reads sharing the same UMI and alignment coordinates are true duplicates. The pipeline collapses these reads into a consensus sequence, which corrects PCR errors. In this case, duplicate marking is replaced by consensus generation, and the resulting BAM contains one read per original fragment.

The practical decision for duplicate marking depends on your library preparation method. If you used PCR amplification without UMIs, mark duplicates. If you used a PCR-free protocol, duplicate rates are low, and marking may still be useful for optical duplicates. If you used UMIs, follow the consensus generation workflow specific to your UMI tool.

Base Quality Recalibration: Correcting Systematic Quality Score Errors

Base quality scores are the foundation of variant calling confidence. Each base in a read carries a Phred-scaled quality score, where a score of Q20 means a 1 in 100 chance of error and Q30 means a 1 in 1000 chance. Variant callers use these scores to weigh evidence for or against a variant.

Sequencers estimate quality scores from the intensity of the fluorescence signal and the surrounding sequence context. These estimates are imperfect. Systematic errors occur in specific contexts, such as:

  • The last 10 to 20 bases of a read, where signal quality degrades
  • Motifs such as CpG dinucleotides, where methylation-related chemistry affects base incorporation
  • Homopolymer runs, where the polymerase stutters
  • Specific dinucleotide contexts that vary by sequencer model and chemistry version

Base quality recalibration (BQSR) addresses these systematic errors by comparing observed base calls to a known truth set. The GATK BaseRecalibrator tool takes the BAM file and a set of known variant sites, typically from dbSNP and the 1000 Genomes project. At these known sites, the tool assumes that the reference allele is correct. Any base that disagrees with the reference is an error. The tool builds a model of error rates as a function of covariates, such as read group, machine cycle, and dinucleotide context.

The model produces a recalibration table. The ApplyBQSR tool then adjusts the quality scores in the BAM file. Bases in error-prone contexts receive lower quality scores, and bases in clean contexts may receive higher scores. The variant caller then uses these recalibrated scores, reducing the chance that a systematic error passes the quality threshold.

BQSR requires a known sites resource. For human data, the standard resource is dbSNP plus the 1000 Genomes phase 1 sites. For non-human organisms, a comparable resource may not exist. In that case, BQSR cannot be run in its standard form. Some pipelines skip BQSR for non-model organisms, relying instead on the variant caller's internal error modeling.

The decision to run BQSR also depends on the variant caller. GATK HaplotypeCaller was designed to work with recalibrated base qualities. DeepVariant, a deep learning based caller, does not require BQSR because it learns error patterns from the raw read data. The GermVarX workflow integrates both GATK HaplotypeCaller and DeepVariant, allowing researchers to compare call sets and build consensus calls. In this workflow, BQSR is applied for the GATK path, while DeepVariant uses the un-recalibrated BAM.

For targeted panels with high depth, BQSR has less impact. The high depth provides many independent observations at each position, and the variant caller can overcome individual base quality errors. However, systematic errors in specific contexts can still bias allele fractions. If your panel targets regions with known difficult contexts, such as GC-rich regions or homopolymers, BQSR may still be beneficial.

Indel Realignment: Addressing Local Alignment Artifacts

Indel realignment is the oldest of the three preprocessing steps. It was developed to correct a specific problem: reads that span an insertion or deletion are often misaligned, with bases placed incorrectly around the indel boundary. These misalignments create mismatches that look like SNVs.

The GATK IndelRealigner tool performed local realignment by identifying regions with excess mismatches, then realigning reads in those regions to minimize the total number of mismatches. The tool used known indel sites and the reads themselves to determine the optimal alignment.

Modern aligners, such as BWA-MEM, handle indels more accurately than older aligners. Modern variant callers, particularly GATK HaplotypeCaller, perform local assembly of the reads in each region. HaplotypeCaller reassembles the reads into haplotypes, which naturally handles indels without a separate realignment step. For this reason, the GATK best practices workflow no longer includes IndelRealigner. The tool is deprecated.

The practical implication is that you should not run IndelRealigner unless you are using an older pipeline that requires it. If your aligner is BWA-MEM and your caller is HaplotypeCaller, DeepVariant, or Mutect2, skip indel realignment. The caller handles local alignment internally.

For legacy pipelines that use older callers, such as UnifiedGenotyper, indel realignment may still be necessary. If you are maintaining an older workflow, verify that your caller does not perform local realignment internally before deciding to include IndelRealigner.

The GATK Best Practices Workflow: A Reference Point

The GATK best practices workflow provides a reference implementation for germline short variant discovery. The workflow proceeds through these stages:

  1. Align reads to the reference genome with BWA-MEM
  2. Mark duplicates with MarkDuplicates
  3. Base quality recalibration with BaseRecalibrator and ApplyBQSR
  4. Variant calling with HaplotypeCaller
  5. Variant filtering with Variant Quality Score Recalibration (VQSR) or hard filters

The workflow is designed for human WGS and WES data. It assumes a high-quality reference genome, a known sites resource for BQSR, and a large cohort for VQSR. The Galaxy Training Network provides accessible tutorials that walk through this workflow step by step, making it possible to learn the process without command-line expertise.

The nf-core documentation describes community pipelines that implement the GATK best practices in a reproducible manner. These pipelines use containerization to ensure that the same tool versions run consistently across different computing environments. For researchers who want to avoid building a pipeline from scratch, nf-core pipelines offer a tested starting point.

The Bioconductor project provides R packages for downstream analysis of variant call sets. After preprocessing and variant calling, Bioconductor packages can filter, annotate, and visualize variants. The integration of preprocessing, calling, and downstream analysis is essential for a complete workflow.

Somatic Variant Calling: Different Preprocessing Considerations

Somatic variant calling identifies mutations present in tumor tissue but absent from normal tissue. The goal is to find variants that drive cancer, which may be present at low allele fractions due to tumor heterogeneity and normal cell contamination.

Somatic workflows differ from germline workflows in several ways. The Onkopipe pipeline integrates quality control, read alignment, BAM preprocessing, and variant calling to detect SNVs, copy number variants, and structural variants without matched normal samples. This design is important for clinical settings where normal tissue may not be available.

For somatic calling with matched tumor-normal pairs, the preprocessing steps are similar to germline calling, but the quality control thresholds differ. Tumor samples often have lower quality due to formalin fixation and paraffin embedding (FFPE), which introduces deamination artifacts. These artifacts appear as C to T transitions and can be mistaken for true mutations.

Duplicate marking is essential for somatic calling because tumor samples often have low input DNA, requiring more PCR amplification cycles. The resulting high duplicate rate can mask true variants if duplicates are not marked.

BQSR is also important for somatic calling, but the known sites resource must be chosen carefully. Using germline dbSNP sites is standard, but somatic mutation databases should not be used for BQSR because the goal is to preserve true somatic mutations.

The VariantMedium caller demonstrates that machine learning based callers can achieve high sensitivity in somatic SNV calling, particularly in regions with high sequencing error rates. The caller was trained on experimentally confirmed variants and uses a 3D DenseNet architecture. For researchers using such callers, the preprocessing requirements may differ from GATK-based workflows. The VariantMedium pipeline is available as an end-to-end solution, meaning the preprocessing steps are built into the pipeline.

The metapipeline-DNA workflow analyzes both germline and somatic genomes from targeted and whole-genome sequencing data. It runs preprocessing, feature detection by multiple algorithms, quality control, and data visualization. Each step can be run independently, allowing researchers to customize the preprocessing for their specific needs.

When to Skip Preprocessing Steps

The decision to skip a preprocessing step should be based on evidence, not convenience. The following scenarios justify skipping specific steps.

Skip duplicate marking when using UMI-based protocols. The UMI consensus generation step replaces duplicate marking because it collapses reads with the same UMI into a single consensus read. Running MarkDuplicates after UMI consensus generation would be redundant and could discard useful information.

Skip BQSR when using DeepVariant. DeepVariant learns error patterns from the raw read data and does not require recalibrated base qualities. Running BQSR before DeepVariant adds computational cost without improving accuracy. The GermVarX workflow demonstrates this by running DeepVariant on un-recalibrated BAM files while running GATK HaplotypeCaller on recalibrated BAM files.

Skip BQSR when no known sites resource exists. For non-model organisms without a comprehensive variant database, BQSR cannot be calibrated. Some pipelines use a self-calibration approach, but this is less reliable than using known sites.

Skip indel realignment when using modern aligners and callers. BWA-MEM and HaplotypeCaller handle indels internally. Running IndelRealigner adds computational cost and may introduce errors.

Skip all preprocessing when using a pipeline that integrates preprocessing internally. The AgrOmicSo interface provides a one-step mode that runs quality control, read mapping, variant calling, and annotation automatically. The metapipeline-DNA workflow similarly integrates preprocessing into the pipeline. In these cases, the user does not need to run preprocessing separately.

Implementing Preprocessing with GATK Tools

The GATK tools for preprocessing are MarkDuplicates, BaseRecalibrator, ApplyBQSR, and IndelRealigner. The following steps describe a standard implementation for germline WGS or WES data.

Step 1: Mark Duplicates

Run MarkDuplicates on the aligned BAM file. The tool requires the input BAM and produces an output BAM with duplicates flagged and a metrics file. The metrics file contains the duplicate rate and other library statistics.

The command structure is:

gatk MarkDuplicates \
  -I input.bam \
  -O marked_duplicates.bam \
  -M metrics.txt

Review the metrics file after running. A duplicate rate above 30 percent for WGS suggests low library complexity. For WES, higher duplicate rates are expected due to the capture step. Record the duplicate rate in your laboratory notebook or analysis log.

Step 2: Base Quality Recalibration

Run BaseRecalibrator to generate the recalibration table. The tool requires the BAM file, the reference genome, and a known sites resource.

gatk BaseRecalibrator \
  -I marked_duplicates.bam \
  -R reference.fasta \
  --known-sites known_sites.vcf.gz \
  -O recalibration_table.table

Review the recalibration table to understand the error model. The table shows error rates by covariate, such as read group and cycle. Large discrepancies between reported and empirical quality scores indicate significant miscalibration.

Run ApplyBQSR to create the recalibrated BAM:

gatk ApplyBQSR \
  -I marked_duplicates.bam \
  -R reference.fasta \
  --bqsr-recal-file recalibration_table.table \
  -O recalibrated.bam

Step 3: Variant Calling

Run HaplotypeCaller on the recalibrated BAM:

gatk HaplotypeCaller \
  -I recalibrated.bam \
  -R reference.fasta \
  -O raw_variants.vcf.gz

The output VCF contains raw variant calls that require filtering before downstream analysis.

Step 4: Variant Filtering

Filter the raw variants using VQSR or hard filters. VQSR requires a large cohort and known sites resources. Hard filters use thresholds for quality metrics such as QD, FS, and MQ.

For small cohorts or targeted panels, hard filters are more practical. The Galaxy Training Network provides tutorials on both filtering approaches.

Records and Measurements: What to Track During Preprocessing

Preprocessing generates metrics that document data quality and inform downstream decisions. The following records should be maintained for each sample.

The duplicate rate from MarkDuplicates indicates library complexity. Record the percentage of reads marked as duplicates. Compare this rate across samples in the same batch to identify library preparation failures.

The base quality recalibration table shows the error model. Record the number of bases used for calibration and the overall error rate. A high error rate before recalibration suggests poor sequencing quality.

The alignment metrics from the aligner, such as the percentage of reads mapped and the percentage of reads properly paired, provide context for the preprocessing metrics. These metrics are typically generated by the aligner and recorded in the BAM header or a separate metrics file.

The EMBL-EBI training resources provide guidance on interpreting sequencing quality metrics. The NCBI data resources offer access to reference genomes and variant databases used in preprocessing.

For reproducible analysis, record the exact tool versions and parameters used. The nf-core documentation emphasizes the importance of version control and containerization for reproducibility. The Carpentries lessons teach the shell, Git, and programming skills needed to manage reproducible workflows.

Common Failure Patterns in Preprocessing

Several failure patterns recur in preprocessing. Recognizing these patterns helps diagnose problems early.

High duplicate rates indicate low library complexity. This can result from insufficient input DNA, excessive PCR cycles, or inefficient adapter ligation. If the duplicate rate exceeds 50 percent, consider whether the library preparation protocol needs adjustment.

BQSR failure occurs when the known sites resource is missing or incompatible with the reference genome. The BaseRecalibrator tool will error if the known sites VCF uses a different reference build. Verify that the reference genome and known sites resource match.

Indel realignment failure occurs when the tool cannot find known indel sites. The IndelRealigner requires a known indels resource, typically from Mills and 1000 Genomes. Without this resource, the tool may not realign effectively.

Preprocessing order errors occur when steps are run in the wrong sequence. MarkDuplicates must run before BQSR because duplicate reads should not contribute to the error model. BQSR must run before variant calling because the caller uses the recalibrated qualities.

Sample mix-ups are detected by comparing the reported sex, ancestry, or relatedness to the expected values. Preprocessing does not detect sample mix-ups, but the downstream variant calls can reveal them. The GermVarX workflow includes sample and cohort level quality control to detect such issues.

Limitations of Preprocessing

Preprocessing cannot correct all sequencing artifacts. The following limitations should be understood before interpreting variant calls.

Duplicate marking cannot distinguish true biological duplicates from PCR duplicates in all cases. Reads from different fragments that happen to align to the same coordinates may be incorrectly marked as duplicates. This is more likely in targeted panels where fragments are short and start positions are constrained.

BQSR cannot correct errors that are not captured by the covariates. If the error model does not include a relevant covariate, such as a specific sequence motif, the recalibration will not correct that error. The model is only as good as the covariates and the known sites resource.

Indel realignment cannot fix alignment errors in repetitive regions. Reads in tandem repeats and homopolymers may align ambiguously, and realignment may not find the correct placement. Variant callers may still produce false calls in these regions.

Preprocessing does not address mapping errors. Reads that map to the wrong location in the genome will produce false variants at the mapped location and missed variants at the true location. Preprocessing operates on the aligned reads and cannot correct alignment errors.

The VariantMedium study highlights that certain genomic regions have high sequencing error rates that challenge even advanced callers. The study found that machine learning based callers can achieve higher accuracy in these regions, but the limitation is inherent to the sequencing technology, not the preprocessing.

Quality Control and Validation After Preprocessing

After preprocessing, verify that the BAM file is ready for variant calling. The following checks provide confidence in the preprocessing output.

Check the duplicate rate and confirm it is within the expected range for your library type. For WGS, a duplicate rate below 20 percent is typical for PCR-free libraries. For PCR-amplified libraries, rates of 10 to 30 percent are common.

Check the base quality scores before and after recalibration. The mean base quality should increase after recalibration if the original scores were underestimated. The distribution of quality scores should be smooth, without spikes at specific values.

Check the coverage distribution across the target regions. For WES, the coverage should be relatively uniform across the captured regions. For targeted panels, the coverage should be consistent across all amplicons.

Check for contamination using tools that estimate the contamination rate from the allele fractions at known polymorphic sites. High contamination can cause false variant calls.

The AgrOmicSo interface provides a step-by-step mode that allows users to review each stage of the analysis. This mode is useful for quality control because it lets the user inspect intermediate files before proceeding.

Professional Escalation Criteria

Some situations require escalation to a bioinformatics specialist or laboratory manager. The following criteria indicate that preprocessing results are not reliable.

Escalate when the duplicate rate exceeds 50 percent and the library cannot be re-prepared. High duplicate rates reduce the effective sequencing depth and may compromise variant calling sensitivity.

Escalate when BQSR produces an error model that does not converge. This may indicate a problem with the known sites resource or the BAM file.

Escalate when the variant caller produces an unexpectedly high number of variants. This may indicate a preprocessing failure or a sample quality issue.

Escalate when sample identity checks fail. If the reported sex or ancestry does not match the expected values, the sample may be mislabeled or contaminated.

Escalate when the preprocessing pipeline fails repeatedly with the same error. This may indicate a software bug, a resource issue, or a data format problem.

The metapipeline-DNA workflow includes automated failure recovery and consistent verification of inputs, outputs, and parameters. These features help identify the source of failures and reduce the need for manual intervention.

Reproducibility and Workflow Management

Reproducibility requires that the same input data and parameters produce the same output. Preprocessing steps must be version controlled, and the tool versions must be recorded.

Containerization ensures that the same tool versions run consistently across different computing environments. The nf-core documentation describes how community pipelines use containers to achieve reproducibility. The Onkopipe pipeline is containerized and provides reproducibility, parallelization, and customization features.

Workflow managers such as Snakemake and Nextflow orchestrate the preprocessing steps. These tools handle dependencies, parallelization, and error recovery. The GermVarX workflow uses Nextflow DSL2, while Onkopipe uses Snakemake.

For researchers who prefer a graphical interface, AgrOmicSo provides a client-server architecture that allows users to manage and execute complex pipelines from a local computer. The interface supports GATK, DeepVariant, and VarScan variant calling algorithms.

The Carpentries lessons teach the foundational skills needed to work with command-line tools and workflow managers. These skills include shell scripting, version control with Git, and data management.

The Bioconductor project provides R packages for reproducible genomic analysis. These packages can document the analysis steps and generate reports that capture the preprocessing parameters and results.

A Practical Decision Framework for Preprocessing Choices

Choosing which preprocessing steps to run is not a single decision but a sequence of decisions that depend on your sequencing strategy, variant caller, and biological question. The framework below organizes these decisions into a structured assessment that can be applied before you commit computational resources to preprocessing. This framework is designed to prevent the common error of applying a default workflow without considering whether each step is appropriate for your specific data.

Step 1: Characterize Your Library Preparation Method

The first decision point is the library preparation protocol. This determines whether duplicate marking is necessary and how it should be implemented.

Ask whether your protocol used PCR amplification. If the protocol is PCR-free, such as the Illumina TruSeq DNA PCR-Free kit, the duplicate rate will be low, and marking duplicates primarily addresses optical duplicates. If the protocol used PCR amplification, duplicates will be present and must be marked.

Ask whether your protocol used unique molecular identifiers. If UMIs were attached to fragments before amplification, the preprocessing workflow changes fundamentally. Instead of marking duplicates, you should collapse reads sharing the same UMI and alignment coordinates into consensus sequences. This consensus generation corrects PCR errors and produces one read per original fragment. Running MarkDuplicates after UMI consensus generation is redundant.

Ask whether your protocol is amplicon-based. Amplicon panels generate reads that share the same start positions because primers define fragment boundaries. These reads are not true biological duplicates, and marking them as duplicates would discard most of the data. Amplicon pipelines should use primer coordinates or read family information to distinguish independent reads from true duplicates.

Step 2: Assess Your Variant Caller Requirements

The second decision point is the variant caller. Different callers have different requirements for base quality recalibration and indel realignment.

If you are using GATK HaplotypeCaller, the caller was designed to work with recalibrated base qualities. The GATK best practices workflow includes BQSR as a required step. Running HaplotypeCaller without BQSR may produce more false positives because the caller trusts the original quality scores.

If you are using DeepVariant, the caller learns error patterns from the raw read data and does not require BQSR. The GermVarX workflow demonstrates this by running DeepVariant on un-recalibrated BAM files while running GATK HaplotypeCaller on recalibrated BAM files. Running BQSR before DeepVariant adds computational cost without improving accuracy.

If you are using a somatic caller such as Mutect2, the preprocessing requirements are similar to germline calling, but the known sites resource for BQSR must be chosen carefully. Germline dbSNP sites are appropriate, but somatic mutation databases should not be used because the goal is to preserve true somatic mutations.

If you are using a machine learning based caller such as VariantMedium, the preprocessing steps may be built into the pipeline. The VariantMedium pipeline is available as an end-to-end solution, meaning the preprocessing steps are integrated and do not need to be run separately.

Step 3: Evaluate Known Sites Resource Availability

The third decision point is the availability of a known sites resource for BQSR. This resource is a VCF file containing known variant sites, typically from dbSNP and the 1000 Genomes project for human data.

If a known sites resource exists for your organism and reference build, BQSR can be run in its standard form. Verify that the resource matches the reference genome build. A mismatch between the known sites VCF and the reference genome will cause BaseRecalibrator to error or produce an incorrect error model.

If no known sites resource exists, BQSR cannot be calibrated reliably. For non-model organisms without a comprehensive variant database, some pipelines use a self-calibration approach, but this is less reliable than using known sites. In this case, skip BQSR and rely on the variant caller's internal error modeling.

Step 4: Determine Whether Indel Realignment Is Needed

The fourth decision point is indel realignment. This step is deprecated in the GATK best practices workflow because modern aligners and callers handle indels internally.

If your aligner is BWA-MEM and your caller is HaplotypeCaller, DeepVariant, or Mutect2, skip indel realignment. These callers perform local assembly or learning that naturally handles indels.

If you are maintaining a legacy pipeline that uses an older caller such as UnifiedGenotyper, indel realignment may still be necessary. Verify that your caller does not perform local realignment internally before deciding to include IndelRealigner.

Step 5: Document Your Decision and Rationale

The final step is documentation. Record the preprocessing decisions for each sample or batch, including the rationale for each choice. This documentation is essential for reproducibility and for interpreting downstream results.

Record the library preparation protocol, including whether PCR amplification and UMIs were used. Record the variant caller and version. Record whether BQSR was run and which known sites resource was used. Record whether indel realignment was run and why.

The nf-core documentation emphasizes the importance of version control and containerization for reproducibility. The Carpentries lessons teach the shell, Git, and programming skills needed to manage reproducible workflows.

A Comparison Table for Preprocessing Decisions

The following table summarizes the decision framework for common sequencing scenarios. Use this table as a starting point and adjust based on your specific data and requirements.

ScenarioMark DuplicatesBQSRIndel RealignmentRationale
Germline WGS, PCR-free, HaplotypeCallerYesYesNoPCR-free libraries have low duplicate rates, but marking is still useful for optical duplicates. BQSR is required for HaplotypeCaller. Modern callers handle indels internally.
Germline WES, PCR-amplified, HaplotypeCallerYesYesNoPCR amplification creates duplicates that must be marked. BQSR corrects systematic quality score errors.
Germline WGS, DeepVariantYesNoNoDeepVariant learns error patterns from raw reads and does not require BQSR.
Targeted amplicon panel, HaplotypeCallerConditionalOptionalNoUse primer coordinates or read family information instead of MarkDuplicates. BQSR has less impact at high depth but may help in difficult contexts.
UMI-based protocol, any callerNoOptionalNoUMI consensus generation replaces duplicate marking. BQSR may still be beneficial depending on the caller.
Somatic tumor-normal, Mutect2YesYesNoTumor samples often have high duplicate rates due to low input DNA. BQSR uses germline known sites only.
Non-model organism, no known sitesYesNoNoBQSR cannot be calibrated without a known sites resource. Rely on the caller's internal error modeling.

Implementing the Decision Framework in Practice

To implement this framework, create a preprocessing decision checklist for each sample or batch. The checklist should include the following items:

  1. Library preparation protocol and whether PCR amplification was used
  2. Whether UMIs were used and the consensus generation tool
  3. Whether the protocol is amplicon-based and how duplicates will be handled
  4. The variant caller and version
  5. Whether a known sites resource is available and matches the reference build
  6. Whether indel realignment is needed based on the aligner and caller
  7. The preprocessing steps selected and the rationale for each

This checklist serves as a record that can be reviewed by collaborators or included in a methods section. The EMBL-EBI training resources provide guidance on documenting bioinformatics workflows.

For researchers using a graphical interface, AgrOmicSo provides both a one-step mode for rapid automated batch processing and a step-by-step mode for detailed customized analyses. The step-by-step mode allows users to review each stage of the analysis and apply the decision framework at each point.

For researchers using a workflow manager, the metapipeline-DNA workflow analyzes targeted and whole-genome sequencing data from raw reads through preprocessing, feature detection, quality control, and data visualization. Each step can be run independently, allowing the decision framework to be applied at each stage.

Common Mistakes in Applying the Decision Framework

Several mistakes recur when researchers apply preprocessing decisions. Recognizing these mistakes helps avoid them.

The first mistake is applying a germline workflow to somatic data without adjustment. Somatic workflows require different quality control thresholds and known sites resources. The Onkopipe pipeline demonstrates that somatic pipelines can detect SNVs, copy number variants, and structural variants without matched normal samples, but the preprocessing must be configured for tumor data.

The second mistake is running BQSR without verifying the known sites resource. A mismatch between the known sites VCF and the reference genome build will produce an incorrect error model. Always verify that the resource matches the reference build before running BaseRecalibrator.

The third mistake is running MarkDuplicates on UMI-based data after consensus generation. This discards useful information and reduces the effective depth. If UMIs were used, the consensus generation step replaces duplicate marking.

The fourth mistake is running IndelRealigner on modern data. This adds computational cost and may introduce errors. Modern aligners and callers handle indels internally.

The fifth mistake is skipping preprocessing entirely because a pipeline integrates it internally. The AgrOmicSo interface and metapipeline-DNA workflow integrate preprocessing into the pipeline, but the user should still understand what preprocessing steps are being run and why.

Recording Preprocessing Decisions for Reproducibility

The decision framework produces a record that supports reproducibility. For each sample, record the following information in a structured format:

The library preparation protocol, including the kit name and whether PCR amplification was used. The UMI strategy, if any, and the consensus generation tool. The variant caller and version. The known sites resource and reference build. The preprocessing steps run and the tool versions. The metrics generated by each preprocessing step, such as the duplicate rate and the BQSR error model.

The Bioconductor project provides R packages for reproducible genomic analysis. These packages can document the analysis steps and generate reports that capture the preprocessing parameters and results.

The Galaxy Training Network provides accessible tutorials that walk through the preprocessing workflow step by step. These tutorials are useful for learning the decision framework and for training new team members.

The NCBI data resources offer access to reference genomes and variant databases used in preprocessing. The EMBL-EBI training resources provide guidance on interpreting sequencing quality metrics and making preprocessing decisions.

Escalation Criteria for Preprocessing Decisions

Some situations require escalation to a bioinformatics specialist or laboratory manager. The following criteria indicate that the decision framework cannot be applied without additional expertise.

Escalate when the library preparation protocol is unclear or undocumented. Without knowing whether PCR amplification or UMIs were used, the preprocessing decisions cannot be made reliably.

Escalate when the known sites resource is unavailable for the organism or reference build. A bioinformatics specialist may be able to identify an alternative resource or a self-calibration approach.

Escalate when the variant caller has unusual preprocessing requirements. Some callers have specific requirements that are not covered by the standard decision framework.

Escalate when the preprocessing decisions produce unexpected results, such as a very high duplicate rate or a BQSR error model that does not converge. These results may indicate a problem with the data or the decisions.

The metapipeline-DNA workflow includes automated failure recovery and consistent verification of inputs, outputs, and parameters. These features help identify the source of failures and reduce the need for manual intervention.

Frequently Asked Questions

What is the difference between marking duplicates and removing duplicates?

Marking duplicates assigns a flag to duplicate reads in the BAM file. The reads remain in the file but are ignored by the variant caller. Removing duplicates deletes the reads from the file. Marking is preferred because it preserves the option to include the reads in other analyses, such as coverage estimation. The GATK MarkDuplicates tool marks duplicates by default.

Does base quality recalibration improve variant calling accuracy for all organisms?

Base quality recalibration improves accuracy when a reliable known sites resource is available. For human data, dbSNP and 1000 Genomes provide this resource. For non-model organisms, a comparable resource may not exist, and BQSR cannot be run in its standard form. In such cases, the variant caller's internal error modeling may be sufficient.

Can I run variant calling without any preprocessing?

You can run variant calling without preprocessing, but the call set will contain more artifacts. Duplicate reads inflate depth and can cause false calls. Uncalibrated base qualities can cause systematic errors to pass quality thresholds. The magnitude of the problem depends on the library preparation method and the variant caller. DeepVariant is more tolerant of raw BAM files than GATK HaplotypeCaller.

Why is indel realignment no longer recommended in GATK best practices?

Modern aligners such as BWA-MEM place reads around indels more accurately than older aligners. Modern variant callers, particularly GATK HaplotypeCaller, perform local assembly of reads in each region. This assembly naturally handles indels without a separate realignment step. IndelRealigner is deprecated and should only be used in legacy pipelines with older callers.

How do UMI-based protocols change the preprocessing workflow?

UMI-based protocols attach a unique molecular identifier to each fragment before amplification. Reads sharing the same UMI and alignment coordinates are true duplicates. The pipeline collapses these reads into a consensus sequence, which corrects PCR errors. This consensus generation replaces duplicate marking. The resulting BAM contains one read per original fragment, and the duplicate rate is effectively zero.

What metrics should I record from the preprocessing steps?

Record the duplicate rate from MarkDuplicates, the number of bases used for BQSR calibration, and the overall error rate from the recalibration table. Record the tool versions and parameters for reproducibility. Record the alignment metrics from the aligner, such as the percentage of reads mapped and properly paired. These metrics provide context for interpreting the preprocessing results.

How does preprocessing differ for somatic variant calling?

Somatic variant calling requires careful handling of tumor samples, which often have low quality due to FFPE preservation. FFPE introduces deamination artifacts that appear as C to T transitions. Duplicate marking is essential because tumor samples often require more PCR amplification. BQSR should use germline known sites, not somatic mutation databases, to preserve true somatic mutations.

What should I do if the duplicate rate is very high?

A high duplicate rate indicates low library complexity. If the rate exceeds 50 percent, the effective sequencing depth is reduced, and variant calling sensitivity may be compromised. Consider whether the library preparation protocol needs adjustment, such as increasing input DNA or reducing PCR cycles. If the library cannot be re-prepared, document the limitation and interpret the variant calls with caution.

Related Bioinformatics Guides

Related Clinical & Scientific Guides

References and Further Reading

This article is educational and does not replace validated analysis plans, institutional policy, clinical interpretation, or specialist review.