argparser.py

import argparse
import os
import subprocess
from .getSegmentsCoordinates import seg_coord
from pathlib import Path
import sys
from tqdm import tqdm

def dir_path(path):
    if os.path.isdir(path):
        return path
    else:
        raise argparse.ArgumentTypeError(f"readable_dir:{path} is not a valid path")
    
def range_type(astr, min=0, max=100):
    value = int(astr)
    if min<= value <= max:
        return value
    else:
        raise argparse.ArgumentTypeError('value not in range %s-%s'%(min,max))
        
def arg():

    first_parser=argparse.ArgumentParser(add_help=False)
    first_parser.add_argument('-pav','--pav_matrix',help="output a presence-absence variation matrix to recapitulate the annotation transfers; requires the -ann/--annotation option",action='store_true')
    # annotation transfer options
    first_parser.add_argument('-id', '--identity', metavar="[0-100]", type=range_type, choices=range(0,101), help="minimum identity percentage required for the feature to be transfered, between 0 and 100; default is 80")
    first_parser.add_argument('-cov', '--coverage', metavar="[0-100]", type=range_type, choices=range(0,101), help="minimum coverage percentage required for the feature to be transfered, between 0 and 100; default is 80")
    first_parser.add_argument('-sh','--source_haplotype',type=str,default='.',help="haplotype of the source genome that the annotation file reffers to; by default the first haplotype will be selected; requires the -ht/--haplotype option")
    first_args, _ = first_parser.parse_known_args()


    parser = argparse.ArgumentParser(prog='grannot',description='Annotation transfer between genome and pangenome graph.',parents=[first_parser],usage='%(prog)s  graph_file  source_annotation_file  source_genome_name  [options]')

    # input files
    parser.add_argument('graph', metavar='graph.gfa', type=argparse.FileType('r'), help='pangenome graph file in GFA format containing the genome that the annotation refers to and the target genomes. the GFA file must have W-lines to describe the genomes walks in the graph, and not P-lines')
    parser.add_argument('-ht','--haplotype',action="store_true",required=(first_args.source_haplotype!='.'),help="use if there are several haplotypes in the graph; default is false")
    parser.add_argument('gff', metavar='annotation.gff', type=argparse.FileType('r') ,help='annotation file in GFF format containing the annotations to be transfered')

    # optionnal input file
    parser.add_argument('-coord', '--segment_coordinates_path', metavar="path/to/files", type=dir_path, default=".", help="path to the files given by the preprocessing of the graph; if not given the program will do the preprocessing; recommended if the program will run several times")

    # genome names
    parser.add_argument('source_genome', metavar="source_genome", type=str, help="name of the annotated genome's path in the GFA (2nd field in the W line or first part (separated by '#') of the 1st field in a P line), that the annotation file refers to")
    parser.add_argument('-t', '--target', metavar="target_genome", nargs='*', type=str, default=[], help="name of the target genomes' paths in the GFA (2nd field in the W line or first part (separated by '#') of the 1st field in a P line) for the annotation; if the option is not used the programm will transfer the annotation on all the genomes")


    # graph annotation options
    parser.add_argument('-gff','--graph_gff',help="output the annotation of the graph in GFF format (not recommended)",action='store_true')
    parser.add_argument('-gaf','--graph_gaf',help="output the annotation of the graph in GAF format",action='store_true')

    # genome annotation options
    parser.add_argument('-ann','--annotation',help="output the annotation transfer in GFF format",action='store_true',required=(first_args.pav_matrix))
    parser.add_argument('-var','--variation',help="output the detail of the variations for the transfered features",action='store_true')
    parser.add_argument('-aln','--alignment',help="output the alignment of the transfered features",action='store_true')

    # general options
    parser.add_argument('-o','--outdir',metavar="outdir",help="directory so store the output files; if not given, the output will go in the current directory",type=str,default='.')
    parser.add_argument('-v','--verbose',action='store_true',help="write more information about the program's operations as it runs",default=False)


    # version
    parser.add_argument('--version', action='version', version='GrAnnoT 1.0.0')

    return parser.parse_args()

def read_args(args):
    # check bedtools installation
    bedtools_test=subprocess.run("bedtools --version",shell=True,check=False,stdout=subprocess.DEVNULL,stderr=sys.stderr)
    if int(bedtools_test.returncode)>0:
        raise SystemExit('bedtools not installed')

    # create outdir
    args.outdir=Path(args.outdir).resolve()
    args.outdir.mkdir(exist_ok=True)

    # set default values
    if args.identity==None:
        args.identity=80
    if args.coverage==None:
        args.coverage=80

    # convert input files in Path objects
    args.outdir=Path(args.outdir).resolve() # resolve() finds the absolute path, and adds the '/' if needed
    seg_coord_path=Path(args.segment_coordinates_path).resolve()
    args.graph=Path(args.graph.name).resolve()
    args.gff=Path(args.gff.name).resolve()

    # run getSegCoord on the target genomes + the source genome if needed
    if not(args.annotation or args.variation or args.alignment or args.pav_matrix): # if we only need the graph annotation, no need to get SegCoord for all the genomes
        path_list=[]
        transfer=False
    else:
        transfer=True
        path_list=args.target[:]
    path_list.append(args.source_genome)
    if args.segment_coordinates_path==".":
        print("Computing the segments coordinates on the genomes")
        seg_coord_path=args.outdir.joinpath("seg_coord")
        seg_coord_path.mkdir(exist_ok=True) # create dir to store bed files
        seg_coord(args.graph.as_posix(),path_list,seg_coord_path,args.verbose,transfer,args.haplotype)
    args.segment_coordinates_path=seg_coord_path

    # intersect between gff and all the bed files from the source genome given by seg_coord
    run_intersect(args,seg_coord_path)

# get the genome name that corresponds to the file_name (file_name if genome_name+walk_name). not universal..?
def get_genome_name(target_genomes,file_name):
    for genome in target_genomes:
        if genome in file_name:
            return genome
    return ""

def run_intersect(args,seg_coord_path):
    print("Computing the intersection between the annotations and the graph segments")

    if args.haplotype:
        if args.source_haplotype=='.':
            files=list(seg_coord_path.glob(f"*{args.source_genome}*.bed")) # get all bed files from the source genome in seg_coord
            # select haplotype
            source_haplotypes=[]
            for file in files:
                file_haplotype=file.name.split("#")[1]
                source_haplotypes.append(file_haplotype)
            first_source_haplo=min(source_haplotypes)
            args.source_haplotype=first_source_haplo
            haplo=f"#{first_source_haplo}#"
            files=list(seg_coord_path.glob(f"*{args.source_genome}*{haplo}*.bed")) # get all bed files from the source genome haplotype in seg_coord
            message=f'No bed files corresponding to the source genome haplotype {first_source_haplo} found in {seg_coord_path}'
        else:
            haplo = f"#{args.source_haplotype}#"
            files=list(seg_coord_path.glob(f"*{args.source_genome}*{haplo}*.bed")) # get all bed files from the source genome in seg_coord
            message=f'No bed files corresponding to the source genome haplotype found in {seg_coord_path}'
        if len(files)==0:
            sys.exit(message)
        if args.verbose:
            print(f'     Found {len(files)} paths corresponding to the source genome haplotype {args.source_haplotype}')
    else:
        files=list(seg_coord_path.glob(f"*{args.source_genome}*.bed")) # get all bed files from the source genome in seg_coord
        message=f'No bed files corresponding to the source genome found in {seg_coord_path}'
        if len(files)==0:
            sys.exit(message)
        if args.verbose:
            print(f'     Found {len(files)} paths corresponding to the source genome')

    intersect_path=args.outdir.joinpath("intersect")
    command=f"echo -n '' > {intersect_path}" # empty the file "intersect"
    subprocess.run(command,shell=True,timeout=None)
    for file in files:
        if args.verbose:
            print(f'          Building the intersection for the path {file.stem}')
        command=f"bedtools intersect -wo -a {file.as_posix()} -b {args.gff.as_posix()}>>{intersect_path}"
        subprocess.run(command,shell=True,timeout=None)