Skip to content
Snippets Groups Projects
argparser.py 5.21 KiB
Newer Older
import argparse
import os
import subprocess
from getSegmentsCoordinates import seg_coord

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))
    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; requires the -ann/--annotation option")
    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; requires the -ann/--annotation option")
    first_parser.add_argument('-diff', '--max_difference', metavar="N", type=int, default=0, help="maximum difference accepted for the feature to be transfered (transfer rejected if the feature size on the target genome is N times bigger or smaller than on the source genome; requires the -ann/--annotation 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]')

    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('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 (field ? in the W 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 (field ? in the W 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",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 or first_args.identity!=None or first_args.coverage!=None))
    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')




    # version
    parser.add_argument('--version', action='version', version='%(prog)s 0.1')

    return parser.parse_args()

def read_args(args):
    if args.identity==None:
        args.identity=80
    if args.coverage==None:
        args.coverage=80
    path_list=args.target[:]
    path_list.append(args.source_genome)
    if args.segment_coordinates_path==".":
        print("run getSegCoord.py")
        seg_coord(args.graph.name,path_list)
        args.segment_coordinates_path="seg_coord/"
    elif args.segment_coordinates_path[-1]!="/":
        args.segment_coordinates_path+="/"
    print("run intersect") # intersect between gff and several bed

    list_files=os.scandir(args.segment_coordinates_path)
    command="echo -n "" > intersect"
    subprocess.run(command,shell=True,timeout=None)
    for file in list_files:
        if args.source_genome in file.name:
            command="bedtools intersect -wo -a "+file.path+" -b "+args.gff.name+">>intersect"
            subprocess.run(command,shell=True,timeout=None)

    #command="for f in "+args.segment_coordinates_path+"*"+args.source_genome+"*; do bedtools intersect -wo -a $f -b "+args.gff.name+" >> "+args.segment_coordinates_path+"intersect;done"

def get_genome_name(target_genomes,file_name):
    for genome in target_genomes:
        if genome in file_name:
            return genome
    return ""