pandora_lib_promethion/src/callers/savana.rs

use crate::{
    annotation::{Annotation, Annotations, Caller, CallerCat},
    collection::{vcf::Vcf, HasOutputs, Initialize, Version},
    commands::{
        bcftools::{bcftools_keep_pass, BcftoolsConfig},
        longphase::{LongphaseConfig, LongphaseHap, LongphasePhase},
    },
    config::Config,
    helpers::force_or_not,
    io::vcf::read_vcf,
    runners::{run_wait, CommandRun, Run},
    variant::{
        variant::{RunnerVariants, Variants},
        variant_collection::VariantCollection,
    },
};
use anyhow::Context;
use log::info;
use rayon::prelude::*;
use std::{fs, path::Path};

#[derive(Debug)]
pub struct Savana {
    pub id: String,
    pub passed_vcf: String,
    pub phased_germline_vcf: String,
    pub normal_hp_bam: String,
    pub tumoral_hp_bam: String,
    pub config: Config,
    pub log_dir: String,
}

impl Initialize for Savana {
    fn initialize(id: &str, config: Config) -> anyhow::Result<Self> {
        let log_dir = format!("{}/{}/log/savana", config.result_dir, id);
        if !Path::new(&log_dir).exists() {
            fs::create_dir_all(&log_dir)
                .context(format!("Failed  to create {log_dir} directory"))?;
        }

        // Check for phased germline vcf
        let phased_germline_vcf = config.constit_phased_vcf(id);
        if !Path::new(&phased_germline_vcf).exists() {
            LongphasePhase::initialize(id, config.clone())?.run()?;
        }

        let normal_hp_bam = config.normal_haplotagged_bam(id);
        let tumoral_hp_bam = config.tumoral_haplotagged_bam(id);

        let passed_vcf = config.savana_passed_vcf(id);
        fs::create_dir_all(config.savana_output_dir(id))?;

        Ok(Self {
            id: id.to_string(),
            passed_vcf,
            config,
            log_dir,
            phased_germline_vcf,
            normal_hp_bam,
            tumoral_hp_bam,
        })
    }
}

impl Run for Savana {
    fn run(&mut self) -> anyhow::Result<()> {
        force_or_not(
            &self.config.savana_output_vcf(&self.id),
            self.config.savana_force,
        )?;
        info!("Running Savana v{}", Savana::version(&self.config)?);

        let id = &self.id;
        let output_vcf = &self.config.savana_output_vcf(id);
        if !Path::new(&output_vcf).exists() {
            // Check for haplotagged bam
            if !Path::new(&self.normal_hp_bam).exists() {
                LongphaseHap::new(
                    id,
                    &self.config.normal_bam(id),
                    &self.phased_germline_vcf,
                    LongphaseConfig::default(),
                )
                .run()?;
            }

            if !Path::new(&self.tumoral_hp_bam).exists() {
                LongphaseHap::new(
                    id,
                    &self.config.tumoral_bam(id),
                    &self.phased_germline_vcf,
                    LongphaseConfig::default(),
                )
                .run()?;
            }

            let savana_args = [
                // "run",
                "--tumour",
                &self.tumoral_hp_bam,
                "--normal",
                &self.normal_hp_bam,
                "--outdir",
                &self.config.savana_output_dir(id),
                "--ref",
                &self.config.reference,
                "--phased_vcf",
                &self.phased_germline_vcf,
                "--no_blacklist",
                "--threads",
                &self.config.savana_threads.to_string(),
            ];
            let args = [
                "-c",
                &format!(
                    "source {} && conda activate savana && {} {}",
                    self.config.conda_sh,
                    self.config.savana_bin,
                    savana_args.join(" ")
                ),
            ];
            let mut cmd_run = CommandRun::new("bash", &args);
            let report = run_wait(&mut cmd_run).context(format!(
                "Error while running `severus.py {}`",
                args.join(" ")
            ))?;

            let log_file = format!("{}/savana_", self.log_dir);
            report
                .save_to_file(&log_file)
                .context(format!("Error while writing logs into {log_file}"))?;
        }

        // Keep PASS
        if !Path::new(&self.passed_vcf).exists() && Path::new(output_vcf).exists() {
            let report =
                bcftools_keep_pass(output_vcf, &self.passed_vcf, BcftoolsConfig::default())
                    .context(format!(
                        "Error while running bcftools keep PASS for {output_vcf}"
                    ))?;
            let log_file = format!("{}/bcftools_pass", self.log_dir);
            report
                .save_to_file(&log_file)
                .context(format!("Error while writing logs into {log_file}"))?;
        }
        Ok(())
    }
}

impl HasOutputs for Savana {
    fn has_output(&self, id: &str) -> (bool, bool) {
        let output = &self.config.savana_passed_vcf(id);
        let passed = &self.config.savana_passed_vcf(id);
        (Path::new(output).exists(), Path::new(passed).exists())
    }
}

impl Version for Savana {
    fn version(config: &Config) -> anyhow::Result<String> {
        let savana_args = ["--version"];
        let args = [
            "-c",
            &format!(
                "source {} && conda activate savana && {} {}",
                config.conda_sh,
                config.savana_bin,
                savana_args.join(" ")
            ),
        ];
        let mut cmd_run = CommandRun::new("bash", &args);
        let report = run_wait(&mut cmd_run)
            .context(format!("Error while running `savana {}`", args.join(" ")))?;
        let log = report.log;
        let start = log
            .find("Version ")
            .context("Failed to find 'Version ' in the log")?;
        let start_index = start + "Version ".len();
        let end = log[start_index..]
            .find('\n')
            .context("Failed to find newline after 'Version '")?;
        Ok(log[start_index..start_index + end]
            .to_string()
            .trim()
            .to_string())
    }
}

impl CallerCat for Savana {
    fn caller_cat(&self) -> (Caller, Annotation) {
        (Caller::Savana, Annotation::Somatic)
    }
}

impl Variants for Savana {
    fn variants(&self, annotations: &Annotations) -> anyhow::Result<VariantCollection> {
        let (caller, category) = self.caller_cat();
        let add = vec![Annotation::Callers(caller.clone()), category.clone()];
        info!(
            "Loading variants from Savana {} with annotations: {:?}",
            self.id, add
        );
        let variants = read_vcf(&self.passed_vcf)?;

        variants.par_iter().for_each(|v| {
            annotations.insert_update(v.hash(), &add);
        });
        Ok(VariantCollection {
            variants,
            vcf: Vcf::new(self.passed_vcf.clone().into())?,
            caller,
            category,
        })
    }
}

impl RunnerVariants for Savana {}