#!/usr/bin/perl 
use strict;
use warnings;


# Copyright (C) 2013 by Pedro Faria 
#
# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to deal
# in the Software without restriction, including without limitation the rights
# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
# copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
#
# The above copyright notice and this permission notice shall be included in
# all copies or substantial portions of the Software.
#
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
# THE SOFTWARE.



my $nome_arquivo_contigs;
my $fim_nome_arquivo_contigs;

my $nome_arquivo_contigs_comparacao;
my $fim_nome_arquivo_contigs_comparacao;

my $nome_arquivo_queries;
my $fim_nome_arquivo_queries;

my $nome_arquivo_confiaveis;
my $fim_nome_arquivo_confiaveis;


#obs.: nao faz muito sentido  (em relacao ao resultado final) fazer exigencias em relacao a cobertura ou identidade, ja que apenas os melhores alinhamentos (ou seja, os de maior pontuacao) serao selecionados (porem, restringir a identidade pode acelerar a velocidade dos alinhamentos)
my $cobertura_minima_queries = 0;
my $cobertura_minima_targets = 0;
my $identidade_minima = 0.50;

my $nome_tabela_saida;

my %melhor_alinhamento_queries_vs_contigs;
my $parametros_alinhamento_queries_vs_contigs = "-q=prot -t=dnax";
my $hash_key_type_alinhamento_queries_vs_contigs = "target"; #indica que o hash com os melhores alinhamentos tera como chave o nome dos targets (contigs)

my %melhor_alinhamento_contigs_vs_confiaveis;
my $parametros_alinhamento_contigs_vs_confiaveis = "-fastMap";
my $hash_key_type_alinhamento_contigs_vs_confiaveis = "query"; #indica que o hash com os melhores alinhamentos tera como chave o nome das queries (contigs)



my %melhor_alinhamento_queries_vs_contigs_comparacao;
my $parametros_alinhamento_queries_vs_contigs_comparacao = "-q=prot -t=dnax";
my $hash_key_type_alinhamento_queries_vs_contigs_comparacao = "target"; #indica que o hash com os melhores alinhamentos tera como chave o nome dos targets (contigs de comparacao)

my %melhor_alinhamento_contigs_comparacao_vs_confiaveis;
my $parametros_alinhamento_contigs_comparacao_vs_confiaveis = "-fastMap";
my $hash_key_type_alinhamento_contigs_comparacao_vs_confiaveis = "query"; #indica que o hash com os melhores alinhamentos tera como chave o nome das queries (contigs de comparacao)



my %melhor_alinhamento_contigs_vs_contigs_comparacao;
my $parametros_alinhamento_contigs_vs_contigs_comparacao = "-fastMap";
my $hash_key_type_alinhamento_contigs_vs_contigs_comparacao = "query"; #indica que o hash com os melhores alinhamentos tera como chave o nome das queries (contigs)

my $help = 
"Uso: perl $0 -c <contigs_interesse.mfasta> -q <queries.mfasta> -g <contigs_confiaveis.mfasta> -t <contigs_comparacao.mfasta>\n
Parametros obrigatorios:

-c     Arquivo em formato multi-fasta com os 'contigs' que se deseja analisar (i.e., os contigs de interesse)
-q     Arquivo em formato multi-fasta com as 'queries' (proteinas, ESTs, cDNAs) usadas para montar os contigs que se deseja analisar
-g     Arquivo em formato multi-fasta com os contigs confiaveis (usados como 'gabarito')
-t     Arquivo em formato multi-fasta com os contigs de comparacao ('targets'), que serao comparados com os contigs que se deseja analisar

Parametros opcionais:

-cq    Cobertura minima (numero entre 0 e 1) exigida das sequencias utilizadas como 'queries' nos alinhamentos 
-ct    Cobertura minima (numero entre 0 e 1) exigida das sequencias utilizadas como 'targets' nos alinhamentos
-pi    Porcentagem de identidade minima exigida nos alinhamentos (entre 0 e 1)\n\n";


le_parametros(\$nome_arquivo_contigs, \$nome_arquivo_queries, \$nome_arquivo_confiaveis, \$nome_arquivo_contigs_comparacao,
              \$cobertura_minima_queries, \$cobertura_minima_targets, \$identidade_minima);


die $help if not($nome_arquivo_contigs and $nome_arquivo_queries and $nome_arquivo_confiaveis);

($fim_nome_arquivo_queries) = ($nome_arquivo_queries =~  /(([^\/]*\/)*)(.*)$/ ) [2];   
($fim_nome_arquivo_contigs) = ($nome_arquivo_contigs =~  /(([^\/]*\/)*)(.*)$/ ) [2];   
($fim_nome_arquivo_contigs_comparacao) = ($nome_arquivo_contigs_comparacao =~  /(([^\/]*\/)*)(.*)$/ ) [2];   
($fim_nome_arquivo_confiaveis) = ($nome_arquivo_confiaveis =~  /(([^\/]*\/)*)(.*)$/ ) [2];   


$nome_tabela_saida = "contigs_".$fim_nome_arquivo_contigs."_vs_contigs_".$fim_nome_arquivo_confiaveis."_vs_contigs_".$fim_nome_arquivo_contigs_comparacao.".tsv";

print "Montando histograma para $fim_nome_arquivo_contigs...\n";
monta_histograma_tamanho_contigs($nome_arquivo_contigs);
print "Histograma para $fim_nome_arquivo_contigs terminado.\n";


#alinhamentos para os contigs que se deseja analisar

print "Fazendo alinhamentos $fim_nome_arquivo_queries vs $fim_nome_arquivo_contigs  ...\n";
encontra_melhores_alinhamentos(\%melhor_alinhamento_queries_vs_contigs, $nome_arquivo_queries, $nome_arquivo_contigs, $cobertura_minima_queries, $cobertura_minima_targets,$identidade_minima, $parametros_alinhamento_queries_vs_contigs, $hash_key_type_alinhamento_queries_vs_contigs);
print "Alinhamentos $fim_nome_arquivo_queries vs $fim_nome_arquivo_contigs terminados.\n";


print "Fazendo alinhamentos $fim_nome_arquivo_contigs vs $fim_nome_arquivo_confiaveis ...\n";
encontra_melhores_alinhamentos(\%melhor_alinhamento_contigs_vs_confiaveis, $nome_arquivo_contigs, $nome_arquivo_confiaveis, $cobertura_minima_queries, $cobertura_minima_targets, $identidade_minima, $parametros_alinhamento_contigs_vs_confiaveis, $hash_key_type_alinhamento_contigs_vs_confiaveis);
print "Alinhamentos $fim_nome_arquivo_contigs vs $fim_nome_arquivo_confiaveis terminados.\n";



#alinhamentos para os contigs de comparacao


print "Fazendo alinhamentos $fim_nome_arquivo_queries vs $fim_nome_arquivo_contigs_comparacao ...\n";
encontra_melhores_alinhamentos(\%melhor_alinhamento_queries_vs_contigs_comparacao, $nome_arquivo_queries, $nome_arquivo_contigs_comparacao, $cobertura_minima_queries, $cobertura_minima_targets,$identidade_minima, $parametros_alinhamento_queries_vs_contigs_comparacao, $hash_key_type_alinhamento_queries_vs_contigs_comparacao);
print "Alinhamentos $fim_nome_arquivo_queries vs $fim_nome_arquivo_contigs_comparacao terminados.\n";


print "Fazendo alinhamentos $fim_nome_arquivo_contigs_comparacao vs $fim_nome_arquivo_confiaveis ...\n";
encontra_melhores_alinhamentos(\%melhor_alinhamento_contigs_comparacao_vs_confiaveis, $nome_arquivo_contigs_comparacao, $nome_arquivo_confiaveis, $cobertura_minima_queries, $cobertura_minima_targets, $identidade_minima, $parametros_alinhamento_contigs_comparacao_vs_confiaveis, $hash_key_type_alinhamento_contigs_comparacao_vs_confiaveis);
print "Alinhamentos $fim_nome_arquivo_contigs_comparacao vs $fim_nome_arquivo_confiaveis terminados.\n";


#comparacao entre os contigs que se deseja analisar e os contigs de comparacao


print "Fazendo alinhamentos $fim_nome_arquivo_contigs vs $fim_nome_arquivo_contigs_comparacao ...\n";
encontra_melhores_alinhamentos(\%melhor_alinhamento_contigs_vs_contigs_comparacao, $nome_arquivo_contigs, $nome_arquivo_contigs_comparacao, $cobertura_minima_queries, $cobertura_minima_targets, $identidade_minima, $parametros_alinhamento_contigs_vs_contigs_comparacao, $hash_key_type_alinhamento_contigs_vs_contigs_comparacao);
print "Alinhamentos $fim_nome_arquivo_contigs vs $fim_nome_arquivo_contigs_comparacao terminados.\n";


print "Montando tabela de comparacao...\n";
monta_tabela_comparacao(\%melhor_alinhamento_queries_vs_contigs, \%melhor_alinhamento_contigs_vs_confiaveis, \%melhor_alinhamento_queries_vs_contigs_comparacao, \%melhor_alinhamento_contigs_comparacao_vs_confiaveis, \%melhor_alinhamento_contigs_vs_contigs_comparacao,$nome_tabela_saida);
print "Tabela de comparacao terminada.\n";



sub monta_histograma_tamanho_contigs {
    my ($nome_arquivo) = @_;
    my $fim_nome_arquivo  = ($nome_arquivo =~  /(([^\/]*\/)*)(.*)$/ ) [2];
    my $nome_script_gnuplot = "script_histograma_".$fim_nome_arquivo.".p";
    my $nome_arquivo_tamanhos = "tamanhos_".$fim_nome_arquivo.".dat";
    my $nome_arquivo_histograma = qq("histograma_para_$fim_nome_arquivo.png");
    my $min;
    my $max;
    my $width = 10; #tamanho dos intervalos ("bins")
    my $total_bases_montadas = 0;

    undef ($min);
    undef ($max);

    open (my $fh, "<", $nome_arquivo) or die "Erro na abertura do arquivo $nome_arquivo: $!\n";
    open (my $fh2, ">", $nome_arquivo_tamanhos) or die "Erro na abertura do arquivo $nome_arquivo_tamanhos: $!\n";
    {
        local $/ = ">";
        <$fh>; #descarta o primeiro ">"
        while (my $fasta = <$fh>) {
            chomp($fasta);
            my ($header, $sequence) = split("\n", $fasta, 2);
            my ($name) = ($header =~ /(\S+)/);
            $sequence =~ s/\s+//g; descarta todos os espacos em branco
            die("Sequencia nula: $name\n") if $sequence eq "";             
            my $length = length($sequence);
            #print "Nome: |$name|, header:|$header|\n";
            #print "Seq : |$sequence|\n";
            print {$fh2} ($length, "\n");
            $total_bases_montadas += $length;
            if (not defined $min or $length < $min) {
                $min = $length;
            }
            if (not defined $max or $length > $max) {
                $max = $length;
            }
        }
    }
    close($fh);
    close($fh2);

    $max = $max + $width;

    open (my $fh3, ">", $nome_script_gnuplot) or die "Erro na abertura do arquivo $nome_script_gnuplot: $!\n";
    print {$fh3} "reset\n";
    print {$fh3} "hist(x,width)=width*floor(x/width)+width/2.0\n";
    print {$fh3} "set term png\n";
    print {$fh3} "set output $nome_arquivo_histograma\n";
    print {$fh3} "set xrange [$min:$max]\n";
    print {$fh3} "set yrange [0:]\n"; 

    print {$fh3} "set xtics $min, ($max-$min)/5, $max\n";
    print {$fh3} "set boxwidth $width*0.9\n"; 
    print {$fh3} "set style fill solid 0.5\n";
    print {$fh3} "set tics out nomirror\n";
    print {$fh3} "set xlabel \"Tamanhos dos contigs\"\n"; 
    print {$fh3} "set ylabel \"Frequencia absoluta\"\n";
    print {$fh3} "set title \"Histograma dos tamanhos - numero de bases montadas = $total_bases_montadas\"\n";
    print {$fh3} "plot \"$nome_arquivo_tamanhos\" u (hist(\$1,$width)):(1.0) smooth freq w boxes lc rgb \"green\" notitle\n";
    close($fh3);
    system("gnuplot $nome_script_gnuplot");
}



sub monta_tabela_comparacao {

    my ($melhor_alinhamento_queries_vs_contigs_hashref, $melhor_alinhamento_contigs_vs_confiaveis_hashref, $melhor_alinhamento_queries_vs_contigs_comparacao_hashref, $melhor_alinhamento_contigs_comparacao_vs_confiaveis, $melhor_alinhamento_contigs_vs_contigs_comparacao_hashref, $nome_tabela_saida) = @_;

    my @linhas_tabela = ();



    foreach (keys %{$melhor_alinhamento_queries_vs_contigs_hashref}) {
        #   0        1      2       3       4        5       6     7         8          9          10         11                12 
        #$qName, $qStart, $qEnd, $qSize, $tName, $tStart, $tEnd, $tSize, $identidade, $_5_size, $_3_size, $cobertura_query, $cobertura_target


        my ($queryName, $queryStart, $queryEnd, $querySize, $contigName, $contigStart, $contigEnd, $contigSize, $id_query_contig, $_5_size, $_3_size, $cobertura_query, $cobertura_contig) = @{${$melhor_alinhamento_queries_vs_contigs_hashref}{$_}};
        my ($contigName_c, $contigStart_c, $contigEnd_c, $contigSize_c, $cName, $cStart, $cEnd, $cSize, $id_contig_c, $_5_size_c, $_3_size_c, $cobertura_contig_c, $cobertura_c) = @{${$melhor_alinhamento_contigs_vs_confiaveis_hashref}{$_}};




        my ($contigName_comp, $contigStart_comp, $contigEnd_comp, $contigSize_comp, $compName, $compStart, $compEnd, $compSize, $id_contig_comp, $_5_size_contig_comp, $_3_size_contig_comp, $cobertura_contig_comp, $cobertura_comp_contig) = @{${$melhor_alinhamento_contigs_vs_contigs_comparacao_hashref}{$_}};






        my ($queryName_comp, $queryStart_comp, $queryEnd_comp, $querySize_comp, $compName_q, $compStart_q, $compEnd_q, $compSize_q, $id_query_comp, $_5_size_comp, $_3_size_comp, $cobertura_query_c, $cobertura_comp_q) = @{${$melhor_alinhamento_queries_vs_contigs_comparacao_hashref}{$compName}};
        my ($compName_c, $compStart_c, $compEnd_c, $compSize_c, $cName_comp, $cStart_comp, $cEnd_comp, $cSize_comp, $id_comp_c, $_5_size_comp_c, $_3_size_comp_c, $cobertura_comp_c, $cobertura_c_comp) = @{${$melhor_alinhamento_contigs_comparacao_vs_confiaveis}{$compName}};





        my $contig_estende_query;#usado para indicar se o contig conseguiu estender a query em algum sentido (5' ou 3')

        if ($_5_size == -1 and $_3_size == -1) {
            $contig_estende_query = "N"; #"NAO" (isto é, o contig não estende a query) 
        }
        else {
            $contig_estende_query = "S"; #"SIM" (isto é, o contig estende a query) 
        }



        my $contig_contido_contig_comparacao;#usado para indicar se o contig esta (quase) totalmente contido no contig de comparacao

        if ($contigStart_comp < 10 and $contigSize_comp - $contigEnd_comp < 10) {
            $contig_contido_contig_comparacao = "S"; #"SIM" (isto é, o contig esta (quase) totalmente contido no contig de comparacao) 
        }
        else {
            $contig_contido_contig_comparacao = "N"; #"NAO" (isto é, o contig nao esta (quase) totalmente contido no contig de comparacao) 
        }




        push(@linhas_tabela, [$contigName, $compName,
                              $queryName, $queryStart, $queryEnd, $querySize, $contigStart, $contigEnd, $contigSize, $_5_size, $_3_size, $cobertura_query, $id_query_contig, $cobertura_contig_c, $id_contig_c,
                              $queryName_comp, $queryStart_comp, $queryEnd_comp, $querySize_comp, $compStart_q, $compEnd_q , $compSize_q  ,$_5_size_comp, $_3_size_comp, $cobertura_query_c, $id_query_comp, $cobertura_comp_c, $id_comp_c,
                              $cobertura_contig_comp, $id_contig_comp]);




    }

    #ordena de forma crescente de acordo com $queryName e depois por $queryStart (isso agrupa os contigs mapeados para a mesma query, permitindo visualizar melhor quais contigs  - os que foram mapeados para regioes intermediarias da query - não teriam chance de pegar as regioes 5' ou 3' da query)                           
    @linhas_tabela = sort { ${$a}[2] cmp ${$b}[2]  || ${$a}[3] <=> ${$b}[3]} @linhas_tabela;


    open (my $fh, ">", $nome_tabela_saida) or die "Erro na abertura do arquivo $nome_tabela_saida: $!\n";

    print {$fh} "contigName\tcompName\tpipe_queryName\tpipe_queryStart\tpipe_queryEnd\tpipe_querySize\tpipe_contigStart\tpipe_contigEnd\tpipe_contigSize\tpipe_5_size\tpipe_3_size\tpipe_cobertura_query\tpipe_id_query\tpipe_cobertura_contig_no_correto\tpipe_id_contig_no_correto\tcomp_queryName\tcomp_queryStart\tcomp_queryEnd\tcomp_querySize\tcomp_contigStart\tcomp_contigEnd\tcomp_contigSize\tcomp_5_size\tcomp_3_size\tcomp_cobertura_query\tcomp_id_query\tcomp_cobertura_comp_no_correto\tcomp_id_comp_no_correto\tpipecomp_cobertura_contig_no_comp\tpipecomp_id_contig_no_comp\n";    

    for (my $i = 0; $i < scalar @linhas_tabela; $i++) {
        my $linha = join ("\t", @{$linhas_tabela[$i]} );
        printf {$fh} "$linha\n";
    }
    close($fh);
}


sub le_parametros {

    my ($nome_arquivo_contigs_ref, $nome_arquivo_queries_ref, $nome_arquivo_confiaveis_ref, $nome_arquivo_contigs_comparacao_ref,
        $cobertura_minima_queries_ref, $cobertura_minima_targets_ref, $identidade_minima_ref) = @_;

    for (my $i = 0; $i < scalar @ARGV; $i++){
        if ($ARGV[$i] eq "-c"){ 
            $i++;
            die "Arquivo com os 'contigs' que se deseja analisar nao foi fornecido" if $i >= scalar @ARGV;
            ${$nome_arquivo_contigs_ref} = $ARGV[$i];     
        }
        elsif ($ARGV[$i] eq "-q"){ 
            $i++;
            die "Arquivo com as 'queries' usadas para montar os contigs que se deseja analisar nao foi fornecido" if $i >= scalar @ARGV;
            ${$nome_arquivo_queries_ref} = $ARGV[$i];     
        }
        elsif ($ARGV[$i] eq "-g"){ 
            $i++;
            die "Arquivo com os contigs confiaveis (usados como 'gabarito') nao foi fornecido" if $i >= scalar @ARGV;
            ${$nome_arquivo_confiaveis_ref} = $ARGV[$i];     
        }

        elsif ($ARGV[$i] eq "-t"){ 
            $i++;
            die "Arquivo com os contigs de comparacao (usados como 'targets' no alinhamento contra os contigs que se deseja analisar) nao foi fornecido" if $i >= scalar @ARGV;
            ${$nome_arquivo_contigs_comparacao_ref} = $ARGV[$i];     
        }

        elsif ($ARGV[$i] eq "-cq"){ 
            $i++;
            die "Cobertura minima das 'queries' nao foi fornecida" if $i >= scalar @ARGV;
            die "Cobertura minima das 'queries' deve ser um valor entre 0 e 1"  if not (0 <= $ARGV[$i] and $ARGV[$i] <= 1);
            ${$cobertura_minima_queries_ref} = $ARGV[$i];     
        }
        elsif ($ARGV[$i] eq "-ct"){ 
            $i++;
            die "Cobertura minima dos 'targets' nao foi fornecida" if $i >= scalar @ARGV;
            die "Cobertura minima dos 'targets' deve ser um valor entre 0 e 1"  if not (0 <= $ARGV[$i] and $ARGV[$i] <= 1);
            ${$cobertura_minima_targets_ref} = $ARGV[$i];     
        }
        elsif ($ARGV[$i] eq "-pi"){ 
            $i++;
            die "Identidade minima nao foi fornecida" if $i >= scalar @ARGV;
            die "Identidade minima deve ser um valor entre 0 e 1"  if not (0 <= $ARGV[$i] and $ARGV[$i] <= 1);
            ${$identidade_minima_ref} = $ARGV[$i];     
        }

        else {
            die "Argumento desconhecido: '$ARGV[$i]'\n$help";
        }    
    }
}


#devolve numero entre 0 e 100 que indica qual a porcentagem de identidade de um alinhamento
sub porcentagem_identidade { 
    my ($match, $repMatch, $misMatch, $tNumInsert, $qNumInsert, $qStart, $qEnd, $blockCount, $tStart, $tEnd, $tSize, $strand_ref, $tStarts_ref, $blockSizes_ref, $isQueryProtein) = @_;
    return 100.0 - pslCalcMilliBad($match, $repMatch, $misMatch, $tNumInsert, $qNumInsert, $qStart,
                                   $qEnd, $blockCount, $tStart, $tEnd, $tSize, $strand_ref,
                                   $tStarts_ref, $blockSizes_ref, $isQueryProtein, 1) * 0.1;
}

sub pslCalcMilliBad {
    my ($match, $repMatch, $misMatch, $tNumInsert, $qNumInsert, $qStart, $qEnd, $blockCount,
        $tStart, $tEnd, $tSize, $strand_ref, $tStarts_ref, $blockSizes_ref, $isQueryProtein,
        $isMrna) = @_;
    my $sizeMul = $isQueryProtein ? 3 : 1;
    my ($qAliSize, $tAliSize, $aliSize);
    my $milliBad = 0;
    my $sizeDif;
    my $insertFactor;
    my $total;

    $qAliSize = $sizeMul * ($qEnd - $qStart);
    $tAliSize = $tEnd - $tStart;
    $aliSize = ($qAliSize < $tAliSize ) ? $qAliSize : $tAliSize;

    if ($aliSize <= 0) {
        return 0;
    }

    $sizeDif = $qAliSize - $tAliSize;

    if ($sizeDif < 0){
        if ($isMrna) {
            $sizeDif = 0;
        }
        else {
            $sizeDif = -$sizeDif;
        }
    }

    $insertFactor = $qNumInsert;
    if (not $isMrna) {
        $insertFactor += $tNumInsert;
    }

    $total = ($sizeMul * ($match + $repMatch + $misMatch));
    if ($total != 0) {
        $milliBad = (1000 * ($misMatch*$sizeMul + $insertFactor +  round(3*log(1+$sizeDif)) ) ) / $total;
    } 
    return $milliBad;
}

#devolve o numero inteiro mais proximo do numero fornecido
sub round {
    my($number) = shift;
    return int($number + 0.5 * ($number <=> 0));
}

sub ceil {
    my($number) = shift;
    return int($number + 0.99* ($number <=> 0));
}

sub floor {
    my($number) = shift;
    return int($number);
}

sub pslScore {
    my ($match, $repMatch, $misMatch, $qNumInsert, $tNumInsert, $isQueryProtein) = @_;
    my $sizeMul = $isQueryProtein ? 3 : 1;

    return $sizeMul * ($match + ($repMatch>>1)) - $sizeMul * $misMatch - $qNumInsert - $tNumInsert;
}

sub encontra_melhores_alinhamentos {
    my ($melhor_alinhamento_queries_vs_targets_hashref, $nome_arquivo_queries, $nome_arquivo_targets,
        $cobertura_minima_queries, $cobertura_minima_targets, $identidade_minima, $parametros_alinhamento,
        $hash_key_type) = @_;
    my $nome_arquivo_blat = "alinhamentos.psl";
    my $porcentagem_identidade_minima = $identidade_minima * 100;

    my $pula_linhas = 0;
    my $queryIsProtein = 0;
    my %melhor_pontuacao;

    system("blat $nome_arquivo_targets $nome_arquivo_queries $nome_arquivo_blat $parametros_alinhamento -minIdentity=$porcentagem_identidade_minima");


    open (my $fh, "<", $nome_arquivo_blat) or die "Erro na abertura do arquivo $nome_arquivo_blat: $!\n";

    while (my $linha = <$fh>) {    
        if ($linha =~ /psLayout/){
	        $pula_linhas = 5;
        }
        if ($pula_linhas > 0) {
            $pula_linhas--;
            next;
        }

        my ($matches, $misMatches, $repMatches, $nCount, $qNumInsert, $qBaseInsert, $tNumInsert,
            $tBaseInsert, $strand, $qName, $qSize, $qStart, $qEnd, $tName, $tSize, $tStart, $tEnd,
            $blockCount, $blockSizes, $qStarts, $tStarts) = split(/\s+/, $linha);

        my @strand = split(//, $strand);
        my @blockSizes = split(/,/, $blockSizes);
        my @qStarts = split(/,/, $qStarts);
        my @tStarts = split(/,/, $tStarts);

        my $identidade = porcentagem_identidade($matches, $repMatches, $misMatches,
                                                $tNumInsert, $qNumInsert, $qStart,
                                                $qEnd, $blockCount, $tStart, $tEnd,
                                                $tSize, \@strand, \@tStarts,
                                                \@blockSizes, $queryIsProtein)/100;
        my $pontuacao = pslScore($matches, $repMatches, $misMatches,
                                 $qNumInsert, $tNumInsert, $queryIsProtein);
        my $cobertura_query = ($qEnd - $qStart) / $qSize;
        my $cobertura_target = ($tEnd - $tStart) / $tSize;

        #analisa a orientacao do target --> nao e' necessario, pois o blat da' as coordenadas da fita  + (mesmo quando o "hit" e' na fita -)    
        #if (scalar @strand == 2 and $strand[1] eq "-") {
        #    ($tStart, $tEnd) = ($tSize - $tEnd, $tSize - $tStart);
        #}


        my $_5_size; 
        my $_3_size;
        my $hash_key;

        if ($hash_key_type eq "query") {
            $hash_key = $qName;
            #analisa o quanto a query conseguiu estender a sequencia target

            if ($tStart <= 10) { #tolerancia de 10 bases
                $_5_size = $qStart;
            }
            else {
                $_5_size = -1;
            }

            if ($tSize - $tEnd <= 10) { #tolerancia de 10 bases
                $_3_size = $qSize - $qEnd;
            }
            else {
                $_3_size = -1;
            }

        }
        else {#$hash_key_type eq "target"
            $hash_key = $tName;
            #analisa o quanto o target conseguiu estender a sequencia query

            if ($qStart <= 10) { #tolerancia de 10 bases
                $_5_size = $tStart;
            }
            else {
                $_5_size = -1;
            }

            if ($qSize - $qEnd <= 10) { #tolerancia de 10 bases
                $_3_size = $tSize - $tEnd;
            }
            else {
                $_3_size = -1;
            }


        }

        if($identidade >= $identidade_minima and
           $cobertura_query >= $cobertura_minima_queries and
           $cobertura_target >= $cobertura_minima_targets) {
            if (not exists ${$melhor_alinhamento_queries_vs_targets_hashref}{$hash_key} or
                $melhor_pontuacao{$hash_key} < $pontuacao) {                
                $melhor_pontuacao{$hash_key} = $pontuacao;
                ${$melhor_alinhamento_queries_vs_targets_hashref}{$hash_key} = [$qName, $qStart, $qEnd, $qSize, $tName, $tStart, $tEnd, $tSize, $identidade, $_5_size, $_3_size, $cobertura_query, $cobertura_target];  
            }
        }

    }
    close ($fh);
}

