<commit_before>#!/usr/bin/env python3

"""
Converts the btab output of AAT (nap/gap2) to GFF3 format with option-dependent modeling.

Example input:
jcf7180000787896        Aug 28 2013     16242   /usr/local/packages/aat/nap     /usr/local/projects/mucormycosis/protein_alignments/fungi_jgi/fungi_jgi.faa     jgi|Rhior3|10928|RO3G_01207     1       103     1162    1198    25.000000       41.666667       15      1       1               -1      Plus    1484
jcf7180000787896        Aug 28 2013     16242   /usr/local/packages/aat/nap     /usr/local/projects/mucormycosis/protein_alignments/fungi_jgi/fungi_jgi.faa     jgi|Rhior3|10928|RO3G_01207     870     1729    1198    1484    91.958042       93.356643       1384    1       2               -1      Plus    1484

The output here depends on the mode requested.  If --export_mode=model (default), then it is:

jcf7180000787896        nap     gene    1       1729    .       +       .       ID=gene.000001
jcf7180000787896        nap     mRNA    1       1729    .       +       .       ID=mRNA.000001;Parent=gene.000001
jcf7180000787896        nap     exon    1       103     .       +       .       ID=exon.000001;Parent=mRNA.000001
jcf7180000787896        nap     CDS     1       103     .       +       0       ID=exon.000001;Parent=mRNA.000001
jcf7180000787896        nap     exon    870     1729    .       +       .       ID=exon.000002;Parent=mRNA.000001
jcf7180000787896        nap     CDS     870     1729    .       +       0       ID=exon.000002;Parent=mRNA.000010

On the other hand, if you specify --export_mode=match you'll get:


Note that chain membership is (intentionally) lost in this mode, and all segments appear as
individual, distinct matches.

"""

import argparse
import os
import re

next_ids = {'gene':1, 'mRNA':1, 'exon':1, 'CDS':1, 'match':1 }
MATCH_TERM = 'nucleotide_to_protein_match'

def main():
    parser = argparse.ArgumentParser( description='AAT (nap/gap2) converter to GFF3 format')

    ## output file to be written
    parser.add_argument('-i', '--input_file', type=str, required=True, help='Path to an input file to parse' )
    parser.add_argument('-o', '--output_file', type=str, required=True, help='Path to an output file to be created' )
    parser.add_argument('-p', '--name_prefix', type=str, required=False, help='So that resulting GFF3 files from multiple runs can be merged, you can supply a prefix for the IDs generated here' )
    parser.add_argument('-d', '--perc_identity_cutoff', type=float, required=False, help='Filters on the percent identity of over the length of the alignment' )
    parser.add_argument('-s', '--perc_similarity_cutoff', type=float, required=False, help='Filters on the percent similarity of over the length of the alignment' )
    parser.add_argument('-m', '--max_intron_cutoff', type=int, required=False, help='Excludes alignment chains which propose an intron greater than this value' )
    parser.add_argument('-e', '--export_mode', type=str, required=False, default='model', help='Controls whether chains are modeled as genes or match regions' )
    
    args = parser.parse_args()
    algorithm = None
    total_chains = 0
    chains_exported = 0

    ofh = open(args.output_file, 'w')
    ofh.write("##gff-version 3\n")
    
    current_chain_number = 1
    gene_segments = []

    for line in open(args.input_file, 'r'):
        cols = line.split("\t")
        chain_num = int(cols[13]);
        
        segment = { 'contig_id':cols[0], 'contig_start':int(cols[6]), 'contig_end':int(cols[7]), \
                    'hit_id':cols[5], 'hit_start':int(cols[8]), 'hit_end':int(cols[9]), \
                    'pct_id':float(cols[10]), 'pct_sim':float(cols[11]), \
                    'strand':cols[17], 'product':cols[15], 'score':cols[18] }

        ## aat product fields are a bit of mess, so lets tease a few parts out:
        #   | organism=Neospora_caninum_Liverpool | product=Iron regulatory protein-like protein, related | location=FR823391:1231501-1237765(-) | length=986 | sequence_SO=chromosome | SO=protein_coding
        #   | organism=Toxoplasma_gondii_ME49 | product=aconitate hydratase ACN/IRP | location=TGME49_chrX:1400716-1407887(-) | length=1055 | sequence_SO=chromosome | SO=protein_coding
        m = re.search("product=(.+?) \|", segment['product'])
        if m:
            hit_product = m.group(1)
        else:
            hit_product = None
        
        if algorithm is None:
            algorithm = os.path.basename(cols[3])

        if args.export_mode == 'match':
            segment_min = min( segment['contig_start'], segment['contig_end'] )
            segment_max = max( segment['contig_start'], segment['contig_end'] )
            hit_min = min( segment['hit_start'], segment['hit_end'] )
            hit_max = max( segment['hit_start'], segment['hit_end'] )
            
            if segment['strand'] == 'Minus':
                segment_strand = '-'
            else:
                segment_strand = '+'

            match_id = get_next_id('match', args.name_prefix)
            data_column = "ID={0};Target={1} {2} {3};Name={4}".format(match_id, segment['hit_id'], hit_min, hit_max, hit_product)
            ofh.write("{0}\t{1}\t{2}\t{3}\t{4}\t{5}\t{6}\t{7}\t{8}\n".format( \
                    segment['contig_id'], algorithm, MATCH_TERM, segment_min, segment_max, '.', segment_strand, \
                    '.', data_column ))
        else:
            if chain_num != current_chain_number:
                chain_pct_id = global_pct_id( gene_segments )

                # this is a new chain<commit_msg>Added support for perc_id and perc_sim cutoffs when in match export mode<commit_after>#!/usr/bin/env python3

"""
Converts the btab output of AAT (nap/gap2) to GFF3 format with option-dependent modeling.

Example input:
jcf7180000787896        Aug 28 2013     16242   /usr/local/packages/aat/nap     /usr/local/projects/mucormycosis/protein_alignments/fungi_jgi/fungi_jgi.faa     jgi|Rhior3|10928|RO3G_01207     1       103     1162    1198    25.000000       41.666667       15      1       1               -1      Plus    1484
jcf7180000787896        Aug 28 2013     16242   /usr/local/packages/aat/nap     /usr/local/projects/mucormycosis/protein_alignments/fungi_jgi/fungi_jgi.faa     jgi|Rhior3|10928|RO3G_01207     870     1729    1198    1484    91.958042       93.356643       1384    1       2               -1      Plus    1484

The output here depends on the mode requested.  If --export_mode=model (default), then it is:

jcf7180000787896        nap     gene    1       1729    .       +       .       ID=gene.000001
jcf7180000787896        nap     mRNA    1       1729    .       +       .       ID=mRNA.000001;Parent=gene.000001
jcf7180000787896        nap     exon    1       103     .       +       .       ID=exon.000001;Parent=mRNA.000001
jcf7180000787896        nap     CDS     1       103     .       +       0       ID=exon.000001;Parent=mRNA.000001
jcf7180000787896        nap     exon    870     1729    .       +       .       ID=exon.000002;Parent=mRNA.000001
jcf7180000787896        nap     CDS     870     1729    .       +       0       ID=exon.000002;Parent=mRNA.000010

On the other hand, if you specify --export_mode=match you'll get:


Note that chain membership is (intentionally) lost in this mode, and all segments appear as
individual, distinct matches.

"""

import argparse
import os
import re

next_ids = {'gene':1, 'mRNA':1, 'exon':1, 'CDS':1, 'match':1 }
MATCH_TERM = 'nucleotide_to_protein_match'

def main():
    parser = argparse.ArgumentParser( description='AAT (nap/gap2) converter to GFF3 format')

    ## output file to be written
    parser.add_argument('-i', '--input_file', type=str, required=True, help='Path to an input file to parse' )
    parser.add_argument('-o', '--output_file', type=str, required=True, help='Path to an output file to be created' )
    parser.add_argument('-p', '--name_prefix', type=str, required=False, help='So that resulting GFF3 files from multiple runs can be merged, you can supply a prefix for the IDs generated here' )
    parser.add_argument('-d', '--perc_identity_cutoff', type=float, required=False, help='Filters on the percent identity of over the length of the alignment' )
    parser.add_argument('-s', '--perc_similarity_cutoff', type=float, required=False, help='Filters on the percent similarity of over the length of the alignment' )
    parser.add_argument('-m', '--max_intron_cutoff', type=int, required=False, help='Excludes alignment chains which propose an intron greater than this value' )
    parser.add_argument('-e', '--export_mode', type=str, required=False, default='model', help='Controls whether chains are modeled as genes or match regions' )
    
    args = parser.parse_args()
    algorithm = None
    total_chains = 0
    chains_exported = 0

    ofh = open(args.output_file, 'w')
    ofh.write("##gff-version 3\n")
    
    current_chain_number = 1
    gene_segments = []

    for line in open(args.input_file, 'r'):
        cols = line.split("\t")
        chain_num = int(cols[13]);
        
        segment = { 'contig_id':cols[0], 'contig_start':int(cols[6]), 'contig_end':int(cols[7]), \
                    'hit_id':cols[5], 'hit_start':int(cols[8]), 'hit_end':int(cols[9]), \
                    'pct_id':float(cols[10]), 'pct_sim':float(cols[11]), \
                    'strand':cols[17], 'product':cols[15], 'score':cols[18] }

        ## aat product fields are a bit of mess, so lets tease a few parts out:
        #   | organism=Neospora_caninum_Liverpool | product=Iron regulatory protein-like protein, related | location=FR823391:1231501-1237765(-) | length=986 | sequence_SO=chromosome | SO=protein_coding
        #   | organism=Toxoplasma_gondii_ME49 | product=aconitate hydratase ACN/IRP | location=TGME49_chrX:1400716-1407887(-) | length=1055 | sequence_SO=chromosome | SO=protein_coding
        m = re.search("product=(.+?) \|", segment['product'])
        if m:
            hit_product = m.group(1)
        else:
            hit_product = None
        
        if algorithm is None:
            algorithm = os.path.basename(cols[3])

        if args.export_mode == 'match':
            if args.perc_identity_cutoff is not None and segment['pct_id'] < args.perc_identity_cutoff:
                continue

            if args.perc_similarity_cutoff is not None and segment['pct_sim'] < args.perc_similarity_cutoff:
                continue
            
            segment_min = min( segment['contig_start'], segment['contig_end'] )
            segment_max = max( segment['contig_start'], segment['contig_end'] )
            hit_min = min( segment['hit_start'], segment['hit_end'] )
            hit_max = max( segment['hit_start'], segment['hit_end'] )
            
            if segment['strand'] == 'Minus':
                segment_strand = '-'
            else:
                segment_strand = '+'

            match_id = get_next_id('match', args.name_prefix)
            data_column = "ID={0};Target={1} {2} {3};Name={4}".format(match_id, segment['hit_id'], hit_min, hit_max, hit_product)
            ofh.write("{0}\t{1}\t{2}\t{3}\t{4}\t{5}\t{6}\t{7}\t{8}\n".format( \
                    segment['contig_id'], algorithm, MATCH_TERM, segment_min, segment_max, '.', segment_strand, \
                    '.', data_column ))
        else:
            if chain_num != current_chain_number:
                chain_pct_id = global_pct_id( gene_segments )

                # this is a new chain
{'fraction_non_alphanumeric': 0.11883252258512857, 'fraction_numerical': 0.07348853370396108, 'mean_word_length': 2.4615153337342153, 'pattern_counts': {'":': 0, '<': 5, '<?xml version=': 0, '>': 3, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 12, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '1186602', 'n_tokens_mistral': 4312, 'n_tokens_neox': 3676, 'n_words': 864}