delly.py

#!/usr/bin/env python3

import sys
import os
import shutil
import argparse
from pathlib import Path
from tempfile import mkdtemp
from shutil import rmtree, copyfile
from subprocess import run

import pysam


def eprint(*args, **kwargs):
    ''' A simple to sys.stderr printer wrapper '''
    print(*args, file=sys.stderr, **kwargs)


def existing_file_abs_path(f):
    try:
        Path(f).resolve()
    except FileNotFoundError:
        print("File %s is not an existing and valid file" % f)
        exit(1)
    else:
        my_abs_path = os.path.abspath(f)
    return my_abs_path


def file_abs_path(f):
    return os.path.abspath(f)


def get_bamfiles(bamlist):
    with open(bamlist, "r") as fin:
        bamfiles = [existing_file_abs_path(f.strip()) for f in fin]
    return bamfiles


def get_all_chromosomes(bamfile):
    with pysam.AlignmentFile(bamfile, "rb") as samfile:
        chromosomes = [contig['SN'] for contig in samfile.header['SQ']]
    return chromosomes


def get_unwanted_chroms(chroms_all, chromosome):
    if chromosome not in chroms_all:
        eprint("Selected chromosome is not a valid chromosome")
    chroms_all.remove(chromosome)
    return chroms_all


def append_chroms_to_exclude(exclusion_file, unwanted_chromosomes, odir="."):
    new_exclusion = os.path.join(odir, "excluded.tsv")
    copyfile(exclusion_file, new_exclusion)
    with open(new_exclusion,  "a") as fout:
        for chrom in unwanted_chromosomes:
            fout.write(chrom+"\n")
    return new_exclusion


def runDelly(args):
    bamlist = args.bamlist
    genome = args.genome
    chrom = args.chrom
    svtype = args.svtype
    excluded = existing_file_abs_path(args.excluded)
    output = file_abs_path(args.output)
    template = file_abs_path(args.template)

    bamfiles = get_bamfiles(bamlist)

    chromosomes_all = get_all_chromosomes(bamfiles[0])

    # removing the selected chrom from chromosomes_all
    unwanted_chromosomes = get_unwanted_chroms(chromosomes_all, chrom)

    tempdir = mkdtemp(dir=".")
    oldir = os.getcwd()
    os.chdir(tempdir)

    # add unwanted chromosomes to the excluded file
    exclusion_file = append_chroms_to_exclude(excluded, unwanted_chromosomes)

    delly_raw = "delly_raw.bcf"

    # construct delly string
    # TODO warning the genologin delly is not multiprocess
    delly_exe = "/home/faraut/save/Softwares/delly_v7.9_src/delly/src/delly"
    # delly_call_str = "delly call "
    delly_call_str = delly_exe + " call "
    delly_call_str += "-g %s " % genome
    delly_call_str += "-x %s " % exclusion_file
    delly_call_str += "-t %s " % svtype
    delly_call_str += "-o %s " % delly_raw
    delly_call_str += " ".join(bamfiles)

    completed = run(delly_call_str, shell=True)

    if completed:
        if os.path.exists(delly_raw):
            # construct delly string
            delly_filter_str = "delly filter "
            delly_filter_str += "-f germline "
            delly_filter_str += "-m 100 "
            delly_filter_str += "-t %s " % svtype
            delly_filter_str += "-o %s " % output
            delly_filter_str += "%s " % delly_raw
            completed = run(delly_filter_str, shell=True)

        else:
            command = "bcftools convert -O b -o {output}.bcf {template};\n".format(output=os.path.splitext(output)[0],
                                                                                   template=template)
            completed = run(command, shell=True)

    os.chdir(oldir)
    rmtree(tempdir)

    return completed.returncode


def parse_arguments():
    parser = argparse.ArgumentParser(
        prog="delly.py",
        description="Delly wrapper for a specific chromosome  "
                    "  ")
    parser.add_argument("-b", "--bamlist", required=True,
                        help="A file with a list of BAM files")
    parser.add_argument("-c", "--chrom", required=True,
                        help="a chromosome")
    parser.add_argument("-g", "--genome", required=True,
                        help="the genome")
    parser.add_argument("-x", "--excluded", required=True,
                        help="File describing regions to exclude")
    parser.add_argument("-t", "--svtype", required=True,
                        help="the svtype")
    parser.add_argument("-o", "--output", required=True,
                        help="outputfile")
    parser.add_argument("-e", "--template", required=True,
                        help="empty template file (used if no variant was found)")

    args = parser.parse_args()

    return args


def main():
    args = parse_arguments()
    return runDelly(args)


if __name__ == "__main__":
    sys.exit(main())