sbx_mummer.smk

import pandas as pd
from itertools import product
import os
import glob

configfile: "config.yaml"

# --- Function for finding absolute path to a genome ---
def find_fasta_path(parent_dir, sample_name):
    matches = []
    
    # Check if absolute path is provided and exists
    if os.path.exists(parent_dir) and os.path.isfile(parent_dir):
        return parent_dir
    
    # Standard sbx_sga pipeline output structure
    pattern1 = os.path.join(parent_dir, "**", "shovill", sample_name, f"{sample_name}.fa")
    matches = sorted(glob.glob(pattern1, recursive=True))

    # Non-standard output structure
    if not matches:
        pattern2 = os.path.join(parent_dir, "**", sample_name, f"{sample_name}.fa")
        matches = sorted(glob.glob(pattern2, recursive=True))
    
    # Non-standard path/extension format
    if not matches:
        pattern3 = os.path.join(parent_dir, "**", f"{sample_name}.f*a*")
        matches = sorted(glob.glob(pattern3, recursive=True))
        
    if matches:
        if len(matches) > 1:
            print(f"WARNING: Multiple matches for {sample_name}. Using: {matches[0]}")
        return matches[0]
    else:
        return f"MISSING_FILE_FOR_{sample_name}"

# --- Read in query and subject CSV files ---
queries = pd.read_csv(config['query_list'])
queries_sample_header = config['query_sample_name_header']
queries_path_header = config['query_path_name_header']

subjects = pd.read_csv(config['subject_list'])
subjects_sample_header = config['subject_sample_name_header']
subjects_path_header = config['subject_path_name_header']

query_genomes = queries.drop_duplicates(subset=[queries_sample_header])
subject_genomes = subjects.drop_duplicates(subset=[subjects_sample_header])

# --- Add queries to subjects so that queries can also be compared to other queries ---
queries_as_subjects = query_genomes[[queries_sample_header, queries_path_header]].rename(
    columns={
        queries_sample_header: subjects_sample_header,
        queries_path_header: subjects_path_header
    }
)

extended_subject_genomes = pd.concat([subject_genomes, queries_as_subjects]).drop_duplicates(
    subset=[subjects_sample_header]
)

# --- Get all samples ---
query_samples = query_genomes[queries_sample_header].tolist()
subject_samples = extended_subject_genomes[subjects_sample_header].tolist()

# Ensure we use ref/query ordering consistent with the rules
all_pairs = list(product(subject_samples, query_samples))

# --- Get path to all genomes using function above ---
path_map = {}

for _, row in query_genomes.iterrows():
    s_name = row[queries_sample_header]
    p_dir = row[queries_path_header]
    path_map[s_name] = find_fasta_path(p_dir, s_name)

for _, row in extended_subject_genomes.iterrows():
    s_name = row[subjects_sample_header]
    p_dir = row[subjects_path_header]
    path_map[s_name] = find_fasta_path(p_dir, s_name)

# --- SNAKEMAKE RULES ---
rule all:
    input:
        "all_mummer_summary.csv"

rule run_mummer:
    input:
        ref = lambda w: path_map[w.ref_name],
        query = lambda w: path_map[w.query_name],
    output:
        report = os.path.join(config['output_dir'], "{ref_name}_{query_name}.report")
    conda:
        "envs/mummer.yml"
    params:
        prefix = "{ref_name}_{query_name}",
        out_dir = config['output_dir']
    shell:
        """
        # Run dnadiff with a temporary local prefix
        dnadiff {input.ref} {input.query} -p {params.prefix} && \
        
        # Move the required report to the output directory designated by Snakemake
        mv {params.prefix}.report {output.report} && \
        
        # Clean up any other collateral files produced by dnadiff
        find . -maxdepth 1 -type f -name "{params.prefix}*" -delete
        """
	
rule summarize_mummer:
    input:
        [f"{config['output_dir']}/{ref}_{q}.report" for ref, q in all_pairs],
    output:
        "all_mummer_summary.csv",
    script:
        "scripts/parse_dnadiff.py"