cmp_breakends.py

import argparse
from pathlib import Path
import subprocess

from pyfaidx import Fasta

from conflict_positions import Region

def create_breakend_plot(args):
    reference = Fasta(str(args.reference), as_raw=True, key_function=lambda x: x.split()[0])

    region1 = Region(args.region1)
    region2 = Region(args.region2)

    seq1, seq2 = region1.get_reference_seq(reference), region2.get_reference_seq(reference)

    fname1 = Path(f"{region1:_}.txt")
    fname2 = Path(f"{region2:_}.txt")
    with open(fname1, "w") as tmp1, open(fname2, "w") as tmp2:
        tmp1.write(seq1)
        tmp2.write(seq2) # if tmp1 == tmp2, this just overwrites the other contents

        print(f"Sequences written to {tmp1.name} and {tmp2.name}")

    print("Running dotter on the sequence(s)...")
    if args.save_plot is None:
        cmd = ["dotter", fname1, fname2]
    else:
        cmd = ["dotter", "-b", args.save_plot, fname1, fname2]

    subprocess.run(cmd)

    if args.keep_seqs:
        fname1.unlink()

        if fname1 != fname2:
            fname2.unlink()


if __name__ == "__main__":
    parser = argparse.ArgumentParser()
    parser.add_argument("reference", type=Path)
    parser.add_argument("region1", type=str)
    parser.add_argument("--region2", type=str, default=None)
    parser.add_argument("--save_plot", "-s", type=Path, default=None)
    parser.add_argument("--keep_seqs", "-k", action='store_true', default=None)

    args = parser.parse_args()

    if args.region2 is None:
        args.region2 = args.region1

    create_breakend_plot(args)