adding Clair3 and bcftools mpileup calling

Author: pmorrell

.github/workflows/long-read-mapping-sanity.yml

@@ -0,0 +1,14 @@
+name: Long Read Mapping Sanity
+
+on:
+  push:
+    branches: [dev]
+  pull_request:
+
+jobs:
+  sanity:
+    runs-on: ubuntu-latest
+    steps:
+      - uses: actions/checkout@v4
+      - name: Run long-read sanity checks
+        run: bash ci/sanity_long_read_mapping.sh

Contributing.md

@@ -2,7 +2,7 @@
 
 # `sequence_handling` Design Principles
 
-This document records the guiding principles for developing and extending the [`sequence_handling`](https://www.google.com/search?q=%5Bhttps://github.com/MorrellLAB/sequence_handling%5D\(https://github.com/MorrellLAB/sequence_handling\)) pipeline. It is intended as an enduring reference for contributors, particularly during the ongoing modernization of the pipeline to support current GATK versions, long-read sequencing technologies (ONT and PacBio HiFi), and updated tooling.
+This document records the guiding principles for developing and extending the [`sequence_handling`](https://github.com/MorrellLAB/sequence_handling) pipeline. It is intended as an enduring reference for contributors, particularly during the ongoing modernization of the pipeline to support current GATK versions, long-read sequencing technologies (ONT and PacBio HiFi), and updated tooling.
 
 These principles apply to all new handlers, to revisions of existing handlers, and to supporting infrastructure such as config files and Slurm job scripts.
 
@@ -43,7 +43,7 @@ Consistent terminology prevents ambiguity across documentation, code, and issues
 - **Config files** (e.g., `Config`) define all parameters for a run. A handler should source a single config and rely entirely on the variables defined there.
 - **HelperScripts** (in `HelperScripts`) are scripts for handling multiple samples.
 - **Sequence_Accessories** (e.g., `PanDepthCoverage.sh`) are tools that may be used occasionally or that supplement `sequence_handling`.
-- This separation of concerns -- job logic in Handlers, scheduling logic in SlurmJobScripts, parameters in Config -- should be preserved as the pipeline grows.
+- This separation of concerns - job logic in Handlers, scheduling logic in SlurmJobScripts, parameters in Config - should be preserved as the pipeline grows.
 
 When adding support for new sequencing platforms or tools, follow this multi-layer structure. Do not mix scheduling directives into handler logic.
 
@@ -54,7 +54,7 @@ When adding support for new sequencing platforms or tools, follow this multi-lay
 All shell code must be safe, readable, and linter-compliant.
 
 - Use strict shell options at the top of every handler: `set -euo pipefail`. This ensures the script exits on errors, treats unset variables as errors, and catches failures in pipelines.
-- All handlers should pass [`shellcheck`](https://www.google.com/search?q=%5Bhttps://www.shellcheck.net/%5D\(https://www.shellcheck.net/\)) without warnings. `shellcheck` compliance is the baseline standard.
+- All handlers should pass [`shellcheck`](https://www.shellcheck.net/) without warnings. `shellcheck` compliance is the baseline standard.
 - Where applicable, code should also satisfy [DeepSource](https://deepsource.com/) static analysis alerts.
 - Variables that come from the config should be validated before use (e.g., check that a file path is non-empty and the file exists before passing it to a tool).
 - Avoid `eval`, unquoted variable expansions in paths, and other patterns that are fragile or unsafe in HPC environments.
@@ -102,7 +102,7 @@ Code that runs silently and produces incorrect results is worse than code that f
 Handlers should produce output that is useful for both humans and AI-assisted troubleshooting.
 
 - Log the tool version, key parameters, input files, and output files at the start and end of each handler run. This information should appear in both stdout and the Slurm log.
-- Capture metrics that are meaningful for QC: read counts, mapping rates, coverage depth, duplication rates, variant counts, and so forth -- whichever are appropriate for the handler's task.
+- Capture metrics that are meaningful for QC: read counts, mapping rates, coverage depth, duplication rates, variant counts, and so forth - whichever are appropriate for the handler's task.
 - Use structured, grep-friendly log lines where possible (e.g., `[Fastplong] Sample: ${SAMPLE} | Reads before: ${N_BEFORE} | Reads after: ${N_AFTER}`).
 - Do not suppress stderr from tools unless you have explicitly handled the error conditions. Suppressed errors are a common source of silent failures in pipelines.
 - Metadata captured at runtime (tool versions, parameters, timestamps) should be written to a per-sample or per-run log file that persists after the job completes.
@@ -128,7 +128,7 @@ All active development occurs on the `dev` branch.
 
 - New handlers, refactored handlers, and bug fixes should be developed on `dev` (or on feature branches that merge into `dev`). The `main` branch reflects the last stable release.
 - `dev` is the integration target: it should remain functional and pass basic testing at all times, even if individual features are incomplete.
-- When a meaningful set of changes has accumulated -- for example, the addition of long-read support, or a GATK version bump -- prepare a new versioned release. Update the changelog, tag the release, and archive via Zenodo so the version is citable.
+- When a meaningful set of changes has accumulated - for example, the addition of long-read support, or a GATK version bump - prepare a new versioned release. Update the changelog, tag the release, and archive via Zenodo so the version is citable.
 - Pull requests into `dev` should include: updated config stanzas (if new parameters are introduced), handler-level comments explaining non-obvious logic, and a brief description in the PR of what changed and why.
 - Breaking changes to the config format or handler interface should be clearly flagged in the changelog and in a migration note for existing users.
 
@@ -139,7 +139,7 @@ The `dev` branch is where the next version of `sequence_handling` is being built
 ## 8. Documentation
 
 - Documentation of the workflow is provided in the `README.md` file. More specifics are available in a [wiki](https://github.com/MorrellLAB/sequence_handling/wiki).
-- [`Dependencies`](https://www.google.com/search?q=%5Bhttps://github.com/MorrellLAB/sequence_handling/wiki/Dependencies%5D\(https://github.com/MorrellLAB/sequence_handling/wiki/Dependencies\)) are listed on a dedicated page, which should be updated so users can easily identify the tools needed for successful execution.
+- [`Dependencies`](https://github.com/MorrellLAB/sequence_handling/wiki/Dependencies) are listed on a dedicated page, which should be updated so users can easily identify the tools needed for successful execution.
 - Where necessary, clarifying comments should be included in code, particularly for more complex operations.
 * * *
 
@@ -157,21 +157,32 @@ _Last updated: March 2026. To be revised as the pipeline evolves._
 
 * * *
 
-## Next steps
+## Recently implemented
 
 - Update from GATK v4.1 to GATK v4.6.
 
 - Add an accessory script to generate an AllSites VCF for [pixy v2.0](https://github.com/ksamuk/pixy). The specific examples for using bcftools mpileup and GATK are [here](https://pixy.readthedocs.io/en/latest/generating_invar/generating_invar.html).
 
-- Document all changes to the new version in a Release file similar to that for [v3.0.0](https://github.com/MorrellLAB/sequence_handling/releases/tag/v3.0.0). Recent updates to the dev branch include replacing' fastp' and' fastplong' with `fastp` and `fastplong` for quality assessment and adapter trimming. This replaces the full front end of the workflow. We have also added `minimap2` for long-read mapping and updated the `config` to include read-mapping presets for PacBio (HiFi reads) and ONT (Q20 reads).
+- Document all changes to the new version in a Release file similar to that for [v3.0.0](https://github.com/MorrellLAB/sequence_handling/releases/tag/v3.0.0). Recent updates to the `dev` branch include replacing `fastp` and `fastplong` for quality assessment and adapter trimming. This replaces the full front end of the workflow. We have also added `minimap2` for long-read mapping and updated the `config` to include read-mapping presets for PacBio (HiFi reads) and ONT (Q20 reads).
 
 - Make sure that the concatenation of gzipped fastq files uses `zcat` rather than `cat`. They don't produce the same results.
 
 - Determine if new GATK indel or SV callers require any changes in our workflow.
 
+
     - Assessment (March 2026): no immediate mandatory workflow changes for germline SNP/indel calling. The current `Haplotype_Caller -> Genomics_DB_Import -> Genotype_GVCFs` path remains valid under GATK 4.6.
     - Indel-specific note: GATK 4.6 includes HaplotypeCaller fixes (including long-deletion edge cases), but does not introduce a replacement germline indel caller that requires handler redesign.
     - SV-specific note: GATK 4.6 includes SV tooling improvements (for annotation/concordance), but full production SV calling is still typically handled by dedicated SV workflows/tools (e.g., GATK-SV WDL stack, pbsv, Sniffles2). Integrating a new SV-calling branch in `sequence_handling` is optional future work, not a blocker for this release.
     - Recommended follow-up: add a dedicated design note before implementing any optional SV branch (inputs, caller choice, output normalization, and filtering strategy).
 
+## Next steps
+
+- Implement a long-read variant caller, probably [`Clair3`](https://github.com/HKU-BAL/Clair3), that will pick appropriate error models given our various read types already specified in the updated `config` file. The read types we know we need to handle include ONT R9 reads, ONT Q20 reads, and PacBio HiFi. For `Clair3`, it may also be important to track which basecaller was used for ONT reads; that would need to be added to the `config` file. This long-read path should bypass most of the GATK germline workflow and land on a VCF that enters a dedicated downstream filtering step. The recommended joint calling path for Clair3 outputs is [`GLnexus`](https://github.com/dnanexus-rnd/GLnexus) rather than GATK's Genomics_DB_Import/Genotype_GVCFs (since GLnexus is designed for consolidating VCFs from non-GATK callers). If needed, we should also be able to create an "AllSites VCF" similar to that for pixy.
+
+- There is also a need to integrate code for some steps in handling Ultima Genomics UG100 resequencing data. For now, this could probably be exclusively in [sequence_accessories](https://github.com/MorrellLAB/sequence_accessories/tree/master). This involves joint variant calling with [GLnexus](https://github.com/dnanexus-rnd/GLnexus) using [GLnexus.sh](https://github.com/pmorrell/Utilities/blob/030effbd0599dd0a0d823cfa19c3bf90bd5e150c/variant_calling/GLnexus.sh#L15) and filtering of those variants using [UG100_filter.sh](https://github.com/MorrellLAB/sequence_accessories/blob/master/Accessories/UG100_filter.sh).
+
+- Integrate existing code for `bcftools mpileup` calling as an alternative for "SNP-only" calling from short-read sequencing. The [bcftools_mpileup.sh](https://github.com/pmorrell/Utilities/blob/030effbd0599dd0a0d823cfa19c3bf90bd5e150c/bcftools_mpileup.sh) script could be added as an additional branch in the workflow after BAM files are sorted and indexed.
+
+- In `sequence_handling_fastp`, we should probably trim the opening number selection to eliminate " 13 | GBS_Demultiplex (in progress)" and all the Nanopore Workflow options. We are integrating long read protocols into the main workflow as side channels. We can move handlers not being deployed to a "Deprecated" directory.
+
 * * *

Handlers/Bcftools_Mpileup_Variant_Calling.sh

@@ -0,0 +1,103 @@
+#!/bin/bash
+
+#   This script performs SNP-only variant calling using bcftools mpileup
+#   as an alternative lightweight path after BAM files are sorted and indexed.
+#   This handler generates only SNP variants (no indels) for short-read data.
+
+set -e
+set -o pipefail
+
+#   What are the dependencies for Bcftools_Mpileup_Variant_Calling?
+declare -a Bcftools_Mpileup_Variant_Calling_Dependencies=(bcftools)
+
+function Bcftools_Mpileup_Variant_Calling() {
+    local sample_list="$1" # What is our sample list (BAM files)?
+    local out_dir="$2" # Where are we storing our results?
+    local ref="$3" # Where is the reference sequence?
+    local threads="$4" # How many threads can we use?
+    local tmp="${5:-}" # Temporary directory (optional)
+    
+    declare -a sample_array=($(grep -E ".bam" "${sample_list}")) # Turn the list into an array
+    
+    # Create output directories
+    mkdir -p "${out_dir}/Bcftools_Mpileup_Variant_Calling"
+    mkdir -p "${out_dir}/Bcftools_Mpileup_Variant_Calling/Intermediates"
+    
+    # Check if temp variable is provided
+    if [[ -z "${tmp}" ]]; then
+        echo "Not using temp directory, proceeding..."
+    else
+        echo "Making temp directory: ${tmp}"
+        mkdir -p "${tmp}"
+    fi
+    
+    # Check if reference genome is indexed for bcftools
+    if ! [[ -f "${ref}.fai" ]]; then
+        echo "Indexing reference genome for bcftools..."
+        bcftools faidx "${ref}"
+    fi
+    
+    # Verify all BAM files are indexed
+    for bam in "${sample_array[@]}"; do
+        if ! [[ -f "${bam}.bai" ]]; then
+            echo "Indexing BAM file: ${bam}"
+            samtools index "${bam}"
+        fi
+    done
+    
+    # Perform mpileup variant calling
+    local mpileup_vcf="${out_dir}/Bcftools_Mpileup_Variant_Calling/Intermediates/mpileup_raw.vcf.gz"
+    local filtered_vcf="${out_dir}/Bcftools_Mpileup_Variant_Calling/mpileup_snps_only.vcf.gz"
+    local sample_log="${out_dir}/Bcftools_Mpileup_Variant_Calling/variant_calling.log"
+    
+    echo "[Bcftools_Mpileup] Starting SNP-only variant calling with bcftools mpileup" | tee -a "${sample_log}"
+    echo "[Bcftools_Mpileup] Input BAM files:" | tee -a "${sample_log}"
+    for bam in "${sample_array[@]}"; do
+        echo "[Bcftools_Mpileup]   ${bam}" | tee -a "${sample_log}"
+    done
+    echo "[Bcftools_Mpileup] Reference genome: ${ref}" | tee -a "${sample_log}"
+    
+    # Run bcftools mpileup
+    if bcftools mpileup \
+        --fasta-ref "${ref}" \
+        --output "${mpileup_vcf}" \
+        --output-type z \
+        --threads "${threads}" \
+        --annotate FORMAT/AD,FORMAT/DP \
+        "${sample_array[@]}" >> "${sample_log}" 2>&1; then
+        
+        echo "[Bcftools_Mpileup] Successfully completed mpileup" | tee -a "${sample_log}"
+        
+        # Filter to SNPs only (exclude indels)
+        echo "[Bcftools_Mpileup] Filtering to SNPs only..." | tee -a "${sample_log}"
+        
+        if bcftools view \
+            "${mpileup_vcf}" \
+            --types snps \
+            --output "${filtered_vcf}" \
+            --output-type z >> "${sample_log}" 2>&1; then
+            
+            echo "[Bcftools_Mpileup] Successfully filtered to SNPs only" | tee -a "${sample_log}"
+            echo "[Bcftools_Mpileup] Output VCF: ${filtered_vcf}" | tee -a "${sample_log}"
+            
+            # Index the final VCF
+            if bcftools index "${filtered_vcf}" >> "${sample_log}" 2>&1; then
+                echo "[Bcftools_Mpileup] Successfully indexed VCF" | tee -a "${sample_log}"
+            else
+                echo "[WARN] Failed to index VCF file" >&2
+            fi
+        else
+            echo "[ERROR] Failed to filter VCF to SNPs only" >&2
+            echo "[ERROR] Check log file: ${sample_log}" >&2
+            return 1
+        fi
+    else
+        echo "[ERROR] Bcftools mpileup failed" >&2
+        echo "[ERROR] Check log file: ${sample_log}" >&2
+        return 1
+    fi
+    
+    echo "[Bcftools_Mpileup] SNP-only variant calling complete"
+}
+
+export -f Bcftools_Mpileup_Variant_Calling