main.py

# created by Nina Marthe 2023 - nina.marthe@ird.fr
# licensed under MIT

import os
import subprocess
from .Graph_gff import *
from .Functions_output import *
from .argparser import *
#from inference import *
#from GrAnnoT import *
from pathlib import Path


def main():

    args=arg()
    read_args(args) 

    # intersect is in the current directory
    intersect=Path("intersect").resolve()
    gfa=args.graph
    load_intersect(intersect.as_posix())

    segments=args.segment_coordinates_path.joinpath("segments.txt")

    # outputs the gff and gaf of the graph
    if args.graph_gff:
        out_graph_gff=gfa.stem+".gff" # what if there is several suffixes and the last one isn't '.gfa' ?
        graph_gff(out_graph_gff)
    if args.graph_gaf:
        out_graph_gaf=gfa.stem+".gaf"
        graph_gaf(out_graph_gaf,segments)


    # if a transfer on a target genome is asked # what about pav ??
    if args.annotation or args.variation or args.alignment:

        # get list of files in seg_coord
        segment_coord_files=list(args.segment_coordinates_path.glob("*.bed"))

        # get the target genomes if there is none specified
        if len(args.target)==0:
            with open(args.segment_coordinates_path.joinpath("walks.txt"),'r') as walks:
                for line in walks:
                    genome_name=line.split()[1]
                    if (args.source_genome not in genome_name) and (genome_name not in args.target) and ("MINIGRAPH" not in genome_name):
                        args.target.append(genome_name)

        pav_dict={}
        # if a pav matrix is asked, create dictionnary to store the features, and for each feature the value for all the target genomes.
        if args.pav_matrix:
            pav_line=len(args.target)*[1]
            pav_dict["gene_id"]=list(args.target)
            for feat in Features.keys():
                if Features[feat].type=='gene':
                    pav_dict[feat]=pav_line.copy()

        # build a dictionnary with the segment sizes to compute the coverage and id
        seg_size=get_segments_length(segments)

        genome_index=0
        for target_genome in args.target:
            print(f'\n{target_genome} transfer :')
            # create directory to store output files
            command="mkdir "+target_genome # path exist_ok=True, methode mkdir
            subprocess.run(command,shell=True,timeout=None)

            # create dictionnaries with paths and segments positions.
            for file in segment_coord_files:
                genome_name=get_genome_name(args.target,file.name)
                if genome_name==target_genome :
                    print(f'    loading the information from the file {file}')
                    #file_path=args.segment_coordinates_path.joinpath(file)
                    get_segments_positions_on_genome(file.as_posix())
            print(f'    loading the walks for the genome {target_genome}')
            walks_path=args.segment_coordinates_path.joinpath("walks.txt").as_posix()
            target_genome_paths=get_paths(walks_path,target_genome)

            # create output files names
            out_target_gff=target_genome+"/"+target_genome+".gff"
            out_target_var=target_genome+"/"+target_genome+"_var.txt"
            out_aln=target_genome+"/"+target_genome+"_aln.txt"

            list_feat_absent=[]
            # do the annotation transfer (or var/aln)
            transfer_on_target(segments,out_target_gff,out_target_var,out_aln,target_genome,target_genome_paths,list_feat_absent,seg_size,args)
            
            # if pav matrix is asked, add the information of this transfer on the matrix
            if args.pav_matrix:
                for feat in list_feat_absent:
                    pav_dict[feat][genome_index]=0
                genome_index+=1
        
        # print the pav matrix
        if args.pav_matrix:
            print('\ngeneration of the presence-absence matrix for the transfered genes')
            pav_output=''
            for line in pav_dict:
                pav_output+=line
                for field in pav_dict[line]:
                    pav_output+="\t"+str(field)
                pav_output+="\n"

            out_pav="PAV_matrix.txt"
            with open(out_pav,'w') as file_out_pav:
                file_out_pav.write(pav_output)