add consensus generation stage (bcftools consensus)

Author: amkram

README.md

@@ -5,12 +5,12 @@ Pangenome-based sequence placement, alignment, and genotyping.
 [Documentation](https://amkram.github.io/panmap/) | [Preprint](https://www.biorxiv.org/content/10.64898/2026.03.29.711974v1)
 
 ## Install
-
+<!-- 
 ```bash
 conda install -c bioconda panmap
 ```
 
-Or with Docker:
+Or with Docker: -->
 
 ```bash
 docker pull alanalohaucsc/panmap:latest
@@ -20,17 +20,17 @@ docker pull alanalohaucsc/panmap:latest
 
 ```bash
 # Place and genotype paired-end reads
-panmap ref.panman reads_R1.fq reads_R2.fq --stop genotype -t 8 -o sample
+panmap ref.panman reads_R1.fq reads_R2.fq -t 8 -o sample
 
 # Metagenomic abundance estimation
-panmap ref.panman reads.fq --meta --index ref.idx -t 8 -o sample
+panmap ref.panman reads.fq --meta -t 8 -o sample
 ```
 
 ## Pipeline
 
 ```
-index  -->  place  -->  align  -->  genotype
- .idx    .placement.tsv  .bam       .vcf
+index  -->  place  -->  align  -->  genotype  -->  consensus
+ .idx    .placement.tsv  .bam       .vcf        .consensus.fa
 ```
 
 By default, panmap stops after placement. Use `--stop` to run further stages.

docs/cli-reference.md

@@ -16,7 +16,7 @@ Use `--help` for common options or `--help-all` for the full list.
 | `-o, --output` | Output file prefix | derived from reads filename |
 | `-t, --threads` | Number of threads | `1` |
 | `-a, --aligner` | `minimap2` or `bwa` | `minimap2` |
-| `--stop` | Stop after: `index`, `place`, `align`, `genotype` | `place` |
+| `--stop` | Stop after: `index`, `place`, `align`, `genotype`, `consensus` | `place` |
 | `--meta` | Enable metagenomic mode | off |
 | `-v, --verbose` | Verbose output | off |
 | `-q, --quiet` | Errors only | off |

docs/quickstart.md

@@ -8,23 +8,23 @@ panmap <panman> [reads1.fq] [reads2.fq] [options]
 
 ## Pipeline
 
-panmap runs four stages in sequence. By default it stops after **placement**. Use `--stop` to run further.
+panmap runs five stages in sequence. By default it stops after **placement**. Use `--stop` to run further.
 
 ```
-index  -->  place  -->  align  -->  genotype
- .idx    .placement.tsv  .bam       .vcf
+index  -->  place  -->  align  -->  genotype  -->  consensus
+ .idx    .placement.tsv  .bam       .vcf        .consensus.fa
 ```
 
 ## Single-sample genotyping
 
-Place reads onto the pangenome, align to the closest reference, and call variants:
+Place reads onto the pangenome, align to the closest reference, call variants, and generate a consensus:
 
 ```bash
 panmap ref.panman reads_R1.fq reads_R2.fq \
-  --stop genotype -t 8 -o sample
+  --stop consensus -t 8 -o sample
 ```
 
-This produces `sample.bam` and `sample.vcf`.
+This produces `sample.bam`, `sample.vcf`, and `sample.consensus.fa`.
 
 ## Metagenomic abundance estimation
 

docs/single-sample.md

@@ -9,7 +9,7 @@ The default mode. Places reads from one sample onto the pangenome tree, aligns t
 panmap ref.panman reads_R1.fq reads_R2.fq -t 8 -o sample
 
 # Run through genotyping
-panmap ref.panman reads_R1.fq reads_R2.fq --stop genotype -t 8 -o sample
+panmap ref.panman reads_R1.fq reads_R2.fq --stop consensus -t 8 -o sample
 ```
 
 ## Pipeline stages
@@ -22,6 +22,7 @@ By default, panmap stops after **placement**. Use `--stop` to control how far th
 | Place | `place` (default) | `<prefix>.placement.tsv` | Places reads onto the pangenome tree |
 | Align | `align` | `<prefix>.bam` | Aligns reads to closest reference |
 | Genotype | `genotype` | `<prefix>.vcf` | Calls variants from alignments |
+| Consensus | `consensus` | `<prefix>.consensus.fa` | Generates consensus FASTA from variants |
 
 !!! note
     When `--stop genotype` is used, `--force-leaf` is enabled automatically.
@@ -54,6 +55,7 @@ panmap examples/data/sars_20000_twilight_dipper.panman \
 | `my_sample.placement.tsv` | Placement on the pangenome tree |
 | `my_sample.bam` | Reads aligned to the closest reference |
 | `my_sample.vcf` | Called variants |
+| `my_sample.consensus.fa` | Consensus sequence |
 
 ### 4. Reuse the index
 
@@ -62,7 +64,7 @@ Once built, the index can be reused across samples:
 ```bash
 panmap examples/data/sars_20000_twilight_dipper.panman \
   sample2_R1.fq.gz sample2_R2.fq.gz \
-  --index my_sample.idx --stop genotype -t 8 -o sample2
+  --index my_sample.idx --stop consensus -t 8 -o sample2
 ```
 
 ## Common options

src/3rdparty/bcftools/bcftools.h

@@ -187,6 +187,7 @@ extern "C" {
 
 int main_mpileup(int argc, char *argv[]);
 int main_vcfcall(int argc, char *argv[]);
+int main_consensus(int argc, char *argv[]);
 
 #ifdef __cplusplus
 }

src/conversion.cpp

@@ -115,6 +115,14 @@ void createVcfWithMutationMatrices(std::string& prefix,
     std::remove(rawVcfFile.c_str());
 }
 
+void createConsensus(const std::string& vcfFileName,
+                     const std::string& refFileName,
+                     const std::string& consensusFileName) {
+    const char* args[] = {
+        "consensus", "-f", refFileName.c_str(), "-o", consensusFileName.c_str(), vcfFileName.c_str()};
+    run_bcftools_in_fork(main_consensus, 6, const_cast<char**>(args));
+}
+
 static uint16_t compute_sam_flags(
     bool is_paired, bool is_read1, uint8_t rev, uint8_t mate_rev, uint8_t proper_frag, bool mate_unmapped) {
     uint16_t flag = 0;

src/conversion.hpp

@@ -34,6 +34,10 @@ void createVcfWithMutationMatrices(std::string& prefix,
                                    std::string& vcfFileName,
                                    const std::vector<std::vector<double>>& substMatrixPhred);
 
+void createConsensus(const std::string& vcfFileName,
+                     const std::string& refFileName,
+                     const std::string& consensusFileName);
+
 // Direct alignment-to-BAM pipeline: parallel minimap2 alignment with direct
 // bam1_t construction (no SAM text intermediate). Writes sorted BAM file.
 void alignAndWriteBam(std::vector<std::string>& readSequences,

src/main.cpp

@@ -100,11 +100,12 @@ inline const char* cyan() {
 }  // namespace color
 
 enum class PipelineStage {
-    Index,     // Build index only
-    Place,     // Placement only
-    Align,     // Placement + Alignment
-    Genotype,  // Placement + Alignment + Genotyping
-    Full       // Full pipeline including assembly
+    Index,      // Build index only
+    Place,      // Placement only
+    Align,      // Placement + Alignment
+    Genotype,   // Placement + Alignment + Genotyping
+    Consensus,  // + Consensus FASTA generation
+    Full        // Full pipeline
 };
 
 struct Config {
@@ -1194,6 +1195,7 @@ int runAlignment(const Config& cfg,
                  const placement::PlacementResult& placement,
                  panmanUtils::Tree* preloadedTree = nullptr);
 int runGenotyping(const Config& cfg);
+int runConsensus(const Config& cfg);
 
 int runBatchPlacement(const Config& cfg) {
     // Parse batch file
@@ -1669,11 +1671,27 @@ int runGenotyping(const Config& cfg) {
     return 0;
 }
 
+int runConsensus(const Config& cfg) {
+    std::string refFileName = cfg.output + ".ref.fa";
+    std::string vcfFileName = cfg.output + ".vcf";
+    std::string consensusFileName = cfg.output + ".consensus.fa";
+
+    auto start = std::chrono::high_resolution_clock::now();
+
+    createConsensus(vcfFileName, refFileName, consensusFileName);
+
+    auto elapsed =
+        std::chrono::duration_cast<std::chrono::milliseconds>(std::chrono::high_resolution_clock::now() - start);
+
+    logging::msg("Consensus complete in {}ms -> {}", elapsed.count(), consensusFileName);
+    return 0;
+}
+
 void printUsage() {
     std::cout << color::bold() << "panmap" << color::reset() << " v" << VERSION << "\n";
     std::cout << "Pangenome-based sequence placement, alignment, and genotyping\n\n";
     std::cout << color::bold() << "Usage:" << color::reset() << "  panmap [options] <panman> [reads.fq] [reads2.fq]\n";
-    std::cout << color::bold() << "Output:" << color::reset() << " <prefix>.vcf, <prefix>.bam, <prefix>.placement.tsv\n"
+    std::cout << color::bold() << "Output:" << color::reset() << " <prefix>.vcf, <prefix>.bam, <prefix>.consensus.fa, ...\n"
               << "        (prefix defaults to reads filename, or use -o)\n\n";
 }
 
@@ -1690,7 +1708,7 @@ int main(int argc, char** argv) {
         ("version,V", "Show version")
         ("output,o", po::value<std::string>(&cfg.output), "Output prefix")
         ("threads,t", po::value<int>(&cfg.threads)->default_value(1), "Threads")
-        ("stop", po::value<std::string>()->default_value("place"), "Stop after: index|place|align|genotype")
+        ("stop", po::value<std::string>()->default_value("place"), "Stop after: index|place|align|genotype|consensus")
         ("meta", po::bool_switch(&cfg.metagenomic), "Metagenomic mode (for more options, see --help-all)")
         ("aligner,a", po::value<std::string>(&cfg.aligner)->default_value("minimap2"), "Aligner: minimap2|bwa")
         ("verbose,v", po::bool_switch(&cfg.verbose), "Verbose output")
@@ -1860,6 +1878,8 @@ int main(int argc, char** argv) {
         cfg.stopAfter = PipelineStage::Align;
     else if (stopStr == "genotype")
         cfg.stopAfter = PipelineStage::Genotype;
+    else if (stopStr == "consensus")
+        cfg.stopAfter = PipelineStage::Consensus;
     else {
         output::error("Invalid stage '{}'", stopStr);
         return 1;
@@ -2064,13 +2084,29 @@ int main(int argc, char** argv) {
             return 1;
         }
         if (signals::check_interrupted()) return 130;
+        if (cfg.stopAfter == PipelineStage::Genotype) {
+            output::done("Variants: " + cfg.output + ".vcf");
+            output::info("");
+            output::info("{}Next steps:{}", color::dim(), color::reset());
+            output::info("  {} View variants: bcftools view {}.vcf | head", color::dim(), cfg.output);
+            output::info("  {} Consensus: panmap {} {} --stop consensus", color::dim(), cfg.panman, cfg.reads1);
+            return 0;
+        }
+
+        // Stage 5: Consensus
+        if (runConsensus(cfg) != 0) {
+            output::error("Consensus generation failed");
+            return 1;
+        }
+        if (signals::check_interrupted()) return 130;
 
         output::info("");
         output::info("{}Pipeline complete.{}", color::green(), color::reset());
         output::info("  Placement:  {}.placement.tsv", cfg.output);
         output::info("  Reference:  {}.ref.fa", cfg.output);
         output::info("  Alignment:  {}.bam", cfg.output);
         output::info("  Variants:   {}.vcf", cfg.output);
+        output::info("  Consensus:  {}.consensus.fa", cfg.output);
         output::info("");
         output::info("{}Next steps:{}", color::dim(), color::reset());
         output::info("  {} View variants: bcftools view {}.vcf | head", color::dim(), cfg.output);