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

# 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.






#O script recebe(obrigatoriamente):
#1) um arquivo multifasta com reads (sequencias) de shotgun 
#2) um arquivo multifasta com sequencias-consulta (ESTs, proteinas ou clusters "full length" de ESTs) 

#O script recebe(opcionalmente):
#3) o nome do diretorio de saida usado pelo genseed 
#4) o nome do diretorio de saida usado pelo pipeline 
#5) o numero maximo de rodadas de montagem (usado pelo genseed) 
#6) a porcentagem de identidade minima exigida nos alinhamentos (usada na filtragem dos mapeamentos do blat)
#7) a cobertura minima exigida das shotguns  (usada na filtragem dos mapeamentos do blat)
#8) a cobertura minima exigida das sequencias-consulta (usada na filtragem dos mapeamentos do blat)      
#9) o numero de nucleos da CPU
#10) o tipo de montador que se deseja utilizar: mira ou phrap (phrap é o padrão)
#11) a tecnologia de sequenciamento utilizada: sanger, 454, iontor, solexa ou solid (sanger é o padrão)

#e produz:

#1) um arquivo do blat no formato PSL, resultante do mapeamento das shotguns nas sequencias-consulta
#2) os mesmos arquivos do phrap (que sao o resultado da montagem das shotguns mapeadas - pelo blat - para a mesma sequencia-consulta)
#3) arquivo multifasta com os contigs resultantes da montagem (usados como sementes para o genseed)
#4) arquivo multifasta com as shotguns que nao foram mapeadas em #1 e os singlets do phrap obtidos em #2 (formando o "database" do genseed) 
#5) os mesmos arquivos dos genseed (resultantes da montagem de todas as shotguns, usando como sementes os contigs obtidos em #2 e o "database" descrito em #4)

my $nome_arquivo_queries = "";
my $fim_nome_arquivo_queries;
my $arquivo_queries_fh;
my $nome_arquivo_queries_selecionadas;

my $nome_arquivo_reads = "";
my $fim_nome_arquivo_reads;
my $arquivo_reads_fh;
my @nomes_arquivos_reads = ();
my $nome_saida_blat;
my $diretorio_saida_genseed = "output_genseed";
my $diretorio_saida_pipeline = "output_pipeline";
my $numero_maximo_rodadas_genseed = "100";
my $montador = "newbler";
my $sequenciamento = "454";

my @sementes;#guarda os nomes dos contigs gerados pelo phrap (serao usados posteriormente pelo genseed)
my @arquivos_banco;#guarda o nome dos arquivos que irao compor o banco (tambem sera usado pelo genseed)

my %query = (); #hash cujas chaves sao os nomes de queries e cujos valores sao as proprias sequencias das queries
my %read = (); #hash cujas chaves sao nomes de sequencias de read e cujos valores sao as proprias sequencias
my %reads_mapeados_para = (); #hash cujas chaves sao nomes de ESTs|proteinas e cujos valores sao os nomes das sequencias de shotgun que foram mapeados para aquele EST|proteina

my @nomes_reads_mapeados = (); #armazena os nomes de quais shotguns foram mapeadas para algum gene
my @sementes_nao_montadas; #vetor com nomes de arquivos que possuam as shotguns mapeadas para o mesmo gene

my $cobertura_minima_read =  0.9;
my $identidade_minima = 0.9;#porcentagem minima de identidade exigida entre as porcoes das sequencias que foram alinhadas
my $cobertura_minima_query = 0.9;

my $numero_cores = 1;
my $modelo_alinhamento;
my $tipo_query; #EST, proteina ou cluster de ESTs
my $pula_genseed = 0;

my $nome_arquivo_sementes = "regioes_genicas.seeds";
my $nome_arquivo_banco = "reads_nao_mapeados.dbs";


my $help = 
"Uso: perl $0 -g <shotgun_reads.mfasta> -p|e|c <queries.mfasta>\n
Parametros obrigatorios:

-g     Arquivo em formato multi-fasta com as sequencias genomicas (reads de shotgun)

-p     Arquivo em formato multi-fasta com as proteinas
-e     Arquivo em formato multi-fasta com as ESTs
-c     Arquivo em formato multi-fasta com os clusters 'full length' de ESTs

Parametros opcionais:

-og    Nome do diretorio de saida usado pelo genseed (padrao: $diretorio_saida_genseed)
-op    Nome do diretorio de saida usado pelo pipeline (padrao: $diretorio_saida_pipeline)   
-r     Numero maximo de rodadas de montagem no genseed (padrao: $numero_maximo_rodadas_genseed)
-cs    Cobertura minima (numero entre 0 e 1) exigida dos reads (padrao: $cobertura_minima_read)
-cq    Cobertura minima (numero entre 0 e 1) exigida das sequencias-consulta (ESTs ou proteinas - padrao: $cobertura_minima_query)
-nn    Numero de nucleos da CPU a serem utilizados (numero inteiro maior que 0 - padrao: $numero_cores)
-pi    Porcentagem de identidade minima exigida nos alinhamentos (numero entre 0 e 1 - padrao: $identidade_minima)
-pg    Pula genseed (nao executa o genseed, que e a ultima fase do pipeline)
-m     O tipo de montador que se deseja utilizar (mira, phrap ou newbler - padrao: $montador) 
-s     A tecnologia de sequenciamento que gerou os reads (sanger, 454, iontor, solexa ou solid - padrao: $sequenciamento)\n\n";

#leitura da entrada
le_parametros(\$nome_arquivo_reads, \$diretorio_saida_genseed, \$numero_maximo_rodadas_genseed,
              \$nome_arquivo_queries, \$modelo_alinhamento, \$cobertura_minima_read,
              \$cobertura_minima_query, \$numero_cores,\$identidade_minima, \$tipo_query,
              \$diretorio_saida_pipeline, \$nome_arquivo_queries_selecionadas, \$pula_genseed,
              \$fim_nome_arquivo_queries, \$fim_nome_arquivo_reads, \$montador,
              \$sequenciamento);

die $help if not($nome_arquivo_queries and $nome_arquivo_reads);

mkdir($diretorio_saida_pipeline) unless -d $diretorio_saida_pipeline;

print "vou popular hash reads\n";
$arquivo_reads_fh = popula_hash(\%read, $nome_arquivo_reads);
print "hash reads populado\n";

print "vou popular hash queries\n";
$arquivo_queries_fh = popula_hash(\%query, $nome_arquivo_queries);
print "hash queries populado\n";

print "vou dividir arquivo com os reads\n";
divide_arquivo_com_os_reads(\%read, \@nomes_arquivos_reads, $numero_cores,
                            $nome_arquivo_reads, $arquivo_reads_fh, $fim_nome_arquivo_reads);
print "arquivo com os reads dividido\n";

system("date >data_inicio_blat_$fim_nome_arquivo_reads.$montador.$numero_cores");

print "vou rodar o blat\n";
roda_e_processa_saida_blat($modelo_alinhamento, $nome_arquivo_queries, \@nomes_arquivos_reads,
                           $tipo_query, \%reads_mapeados_para, \@nomes_reads_mapeados,
                           $identidade_minima, $cobertura_minima_query, $cobertura_minima_read,
                           \%read, $diretorio_saida_pipeline, $numero_cores,
                           $fim_nome_arquivo_queries, $fim_nome_arquivo_reads);
print "saidas do blat criadas e processadas\n";

system("date >data_fim_blat_$fim_nome_arquivo_reads.$montador.$numero_cores");

print "vou montar arquivo com as queries selecionadas\n";
gera_arquivo_queries_selecionadas(\%reads_mapeados_para, \%query,
                                   $nome_arquivo_queries_selecionadas, $arquivo_queries_fh);
print "arquivo com os queries selecionadas montado\n";


print "vou montar arquivos com os conjuntos de reads selecionados\n";
gera_arquivos_sementes_nao_montadas(\%reads_mapeados_para, \@sementes_nao_montadas,
                                     \%read, $arquivo_reads_fh);
print "arquivos com os conjuntos de reads selecionados montados\n";

print "nomes das sementes nao montadas: @sementes_nao_montadas\n";

print "vou montar arquivo com os reads nao mapeados\n";
gera_arquivo_reads_nao_mapeados(\@nomes_reads_mapeados, \%read,
                                 \@arquivos_banco, $arquivo_reads_fh);
print "arquivo com os reads nao mapeados montado\n";


system("date >data_inicio_montador_$fim_nome_arquivo_reads.$montador.$numero_cores");

roda_montador($montador, $sequenciamento, $numero_cores, \@sementes, \@sementes_nao_montadas);

system("date >data_fim_montador_$fim_nome_arquivo_reads.$montador.$numero_cores");


if (not $pula_genseed) {
    monta_arquivos_genseed(\@sementes, \@arquivos_banco, $diretorio_saida_genseed, $numero_maximo_rodadas_genseed,
                       $nome_arquivo_sementes, $nome_arquivo_banco);
    roda_genseed($nome_arquivo_sementes, $nome_arquivo_banco, $diretorio_saida_genseed, $numero_maximo_rodadas_genseed);
    system("date >data_fim_genseed_$fim_nome_arquivo_reads.$numero_cores");
}

close($arquivo_reads_fh);
close($arquivo_queries_fh);

exit;


#############################
#                           # 
#  sub-rotinas principais   #
#                           #
#############################

sub le_parametros {
    my ($nome_arquivo_reads_ref, $diretorio_saida_genseed_ref, $numero_maximo_rodadas_genseed_ref,
        $nome_arquivo_queries_ref, $modelo_alinhamento_ref, $cobertura_minima_read_ref,
        $cobertura_minima_query_ref, $numero_cores_ref, $identidade_minima_ref,
        $tipo_query_ref, $diretorio_saida_pipeline_ref,$nome_arquivo_queries_selecionadas_ref,               
        $pula_genseed_ref, $fim_nome_arquivo_queries_ref, $fim_nome_arquivo_reads_ref,
        $montador_ref, $sequenciamento_ref) = @_;



    for (my $i = 0; $i < scalar @ARGV; $i++){
        if ($ARGV[$i] eq "-g"){ 
            $i++;
            die "Arquivo com os reads de shotgun nao foi fornecido" if $i >= scalar @ARGV;
            ${$nome_arquivo_reads_ref} = $ARGV[$i];  
            (${$fim_nome_arquivo_reads_ref}) = ($ARGV[$i] =~  /(([^\/]*\/)*)(.*)$/ ) [2];   
        }
        elsif ($ARGV[$i] eq "-p") { 
            $i++;
            die "Arquivo com as proteinas nao foi fornecido" if $i >= scalar @ARGV;    
            ${$nome_arquivo_queries_ref} = $ARGV[$i]; 
            (${$fim_nome_arquivo_queries_ref}) = ($ARGV[$i] =~  /(([^\/]*\/)*)(.*)$/ ) [2];   

            ${$nome_arquivo_queries_selecionadas_ref} = ${$fim_nome_arquivo_queries_ref}.".selecionadas";               
            ${$modelo_alinhamento_ref} = "-q=prot -t=dnax";
            ${$tipo_query_ref} = "proteina";
        }
        elsif ($ARGV[$i] eq "-e") { 
            $i++;
            die "Arquivo com as ESTs nao foi fornecido" if $i >= scalar @ARGV;    
            ${$nome_arquivo_queries_ref} = $ARGV[$i];    
            (${$fim_nome_arquivo_queries_ref}) = ($ARGV[$i] =~  /(([^\/]*\/)*)(.*)$/ ) [2];   
            ${$nome_arquivo_queries_selecionadas_ref} = ${$fim_nome_arquivo_queries_ref}.".selecionadas";               
            ${$modelo_alinhamento_ref} = ""; 
            ${$tipo_query_ref} = "est";   
        }
        elsif ($ARGV[$i] eq "-c") { 
            $i++;
            die "Arquivo com os clusters 'full length' de ESTs nao foi fornecido" if $i >= scalar @ARGV;    
            ${$nome_arquivo_queries_ref} = $ARGV[$i];   
            (${$fim_nome_arquivo_queries_ref}) = ($ARGV[$i] =~  /(([^\/]*\/)*)(.*)$/ ) [2];    
            ${$nome_arquivo_queries_selecionadas_ref} = ${$fim_nome_arquivo_queries_ref}.".selecionadas";               
            ${$modelo_alinhamento_ref} = ""; 
            ${$tipo_query_ref} = "full-length-est";   
        }
        elsif ($ARGV[$i] eq "-og"){ 
            $i++;
            die "Diretorio de saida do genseed nao foi fornecido" if $i >= scalar @ARGV;
            ${$diretorio_saida_genseed_ref}  = "$ARGV[$i]";     
        }

        elsif ($ARGV[$i] eq "-op"){ 
            $i++;
            die "Diretorio de saida do pipeline nao foi fornecido" if $i >= scalar @ARGV;
            ${$diretorio_saida_pipeline_ref}  = $ARGV[$i];     
        }

        elsif ($ARGV[$i] eq "-r"){ 
            $i++;
            die "Numero maximo de rodadas do genseed nao foi fornecido" if $i >= scalar @ARGV;
            die "Numero maximo de rodadas do genseed deve ser um valor maior que 0" if $ARGV[$i] <= 0;
            ${$numero_maximo_rodadas_genseed_ref} = "$ARGV[$i]";     
        }
        elsif ($ARGV[$i] eq "-cs"){ 
            $i++;
            die "Cobertura minima da shotgun nao foi fornecida" if $i >= scalar @ARGV;
            die "Cobertura minima da shotgun deve ser um valor entre 0 e 1"  if not (0 <= $ARGV[$i] and $ARGV[$i] <= 1);
            ${$cobertura_minima_read_ref} = $ARGV[$i];     
        }
        elsif ($ARGV[$i] eq "-cq"){ 
            $i++;
            die "Cobertura minima da query nao foi fornecida" if $i >= scalar @ARGV;
            die "Cobertura minima da query deve ser um valor entre 0 e 1"  if not (0 <= $ARGV[$i] and $ARGV[$i] <= 1);
            ${$cobertura_minima_query_ref} = $ARGV[$i];     
        }
        elsif ($ARGV[$i] eq "-nn"){ 
            $i++;
            die "Numero de nucleos da CPU nao foi fornecido" if $i >= scalar @ARGV;
            die "Numero de nucleos da CPU deve ser maior ou igual a 1"  if $ARGV[$i] < 1;
            ${$numero_cores_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];     
        }
        elsif ($ARGV[$i] eq "-pg"){
            ${$pula_genseed_ref} = 1;     
        }
        elsif ($ARGV[$i] eq "-m"){ 
            $i++;
            die "Montador nao foi fornecido" if $i >= scalar @ARGV;
            die "Montador deve ser 'phrap', 'mira' ou 'newbler'"  if not ($ARGV[$i] =~ /^(mira|phrap|newbler)$/);
            ${$montador_ref} = $ARGV[$i];     
        }
        elsif ($ARGV[$i] eq "-s"){ 
            $i++;
            die "Tecnologia de sequeciamento nao foi fornecida" if $i >= scalar @ARGV;
            die "Tecnologia de sequenciamento deve ser 'sanger', '454', 'iontor', 'solexa' ou 'solid'"  if not ($ARGV[$i] =~ /^(sanger|454|iontor|solexa|solid)$/);
            ${$sequenciamento_ref} = $ARGV[$i];     
        }
        else {
            die "Argumento desconhecido: '$ARGV[$i]'\n$help";
        }    
    }
}

#popula um hash a partir de um arquivo FASTA
sub popula_hash {
    my ($hash_ref, $nome_arquivo) = @_;

    open (my $fh, "<", $nome_arquivo) or die "Erro na abertura do arquivo $nome_arquivo: $!\n";
    {
        local $/ = ">";
        <$fh>; #descarta o primeiro ">"
        my $filepos = tell($fh);
        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
            #print "Nome: |$name|, header:|$header|\n";
            #print "Seq : |$sequence|\n";
            if ($sequence ne "") { #se a sequencia nao e nula, entao a posicao do arquivo onde ela esta e armazenada no hash           
                ${$hash_ref}{$name} = $filepos;
            }
            $filepos = tell($fh);
            #die("Sequencia nula: $name\n") if $sequence eq "";             
        }
    }
    return $fh;
}


sub divide_arquivo_com_os_reads {
    my ($read_hashref, $nomes_arquivos_reads_vecref, $numero_cores,
        $nome_arquivo_reads, $arquivo_reads_fh, $fim_nome_arquivo_reads) = @_;
    my $numero_reads_por_arquivo = round((scalar keys %{$read_hashref}) / $numero_cores);

    my $numero_arquivo = 1;
    my $numero_reads_arquivo_atual = 0;
    my $numero_atual_reads = 0;
    my $numero_total_reads = scalar keys %{$read_hashref};
    my $nome_arquivo = $fim_nome_arquivo_reads.".part".$numero_arquivo;
    

    push (@{$nomes_arquivos_reads_vecref}, $nome_arquivo);
    open (FASTA, ">", $nome_arquivo) or die "Erro na abertura do arquivo $nome_arquivo: $!\n";
    foreach my $nome_read (keys %{$read_hashref}) {
        print FASTA ">$nome_read\n";
        print FASTA (get_sequence_by_name($nome_read, $read_hashref, $arquivo_reads_fh));
        $numero_reads_arquivo_atual++;
        $numero_atual_reads++;
        if ($numero_reads_arquivo_atual == $numero_reads_por_arquivo) {
            $numero_reads_arquivo_atual = 0;
            $numero_arquivo++;
            close(FASTA);
            if ($numero_atual_reads < $numero_total_reads) {
                $nome_arquivo = $fim_nome_arquivo_reads.".part".$numero_arquivo;
                push (@{$nomes_arquivos_reads_vecref}, $nome_arquivo);
                open (FASTA, ">", $nome_arquivo) or die "Erro na abertura do arquivo $nome_arquivo: $!\n";
            }
        }
    } 
}

sub roda_e_processa_saida_blat {
    my ($modelo_alinhamento, $nome_arquivo_queries, $nomes_arquivos_reads_vecref, $tipo_query,
        $hashref_reads_mapeados_para, $vecref_nomes_reads_mapeados,
        $identidade_minima, $cobertura_minima_query, $cobertura_minima_read,
        $hashref_read, $diretorio_saida_pipeline, $numero_cores,
        $fim_nome_arquivo_queries, $fim_nome_arquivo_reads) = @_;

    my %query_mapeada_para = (); #hash cujas chaves sao nomes de sequencias de shotgun (targets) e cujo valor e um vetor com nome da query para o qual ela foi mapeada e a pontuacao do alinhamento
    my @nomes_arquivos_blat = ();


    #roda o blat de forma paralela
    my $numero_filhos_ativos = 0;
    foreach my $nome_arquivo_reads (@{$nomes_arquivos_reads_vecref}) {
        my $nome_arquivo_blat = "alinhamentos_".$fim_nome_arquivo_queries."_em_".$nome_arquivo_reads.".psl";
        push(@nomes_arquivos_blat, $nome_arquivo_blat);

        if ($numero_filhos_ativos >= $numero_cores) { #espera um dos filhos terminar
            wait();
            $numero_filhos_ativos--;
        }        

        #cria mais um processo filho
        my $pid = fork();
        if ($pid == 0) { #processo filho
            exec("blat $nome_arquivo_reads $nome_arquivo_queries $nome_arquivo_blat $modelo_alinhamento") or print STDERR "Erro na execucao do blat: $!";
        }
        elsif(not defined $pid) {
            print STDERR "Erro: fork() nao funcionou\n";    
            exit(1);
        }
        else { #processo pai
            $numero_filhos_ativos++;
        }   
    }
    
    #espera todos os filhos terminarem 
    while(wait() != -1) {}


    #move os arquivos de saida do blat (gerados por cada filho) para o diretorio de saida do pipeline
    foreach my $nome_arquivo_blat (@nomes_arquivos_blat) {
        system("mv $nome_arquivo_blat $diretorio_saida_pipeline");   
    }
    chdir($diretorio_saida_pipeline);   

 
    my $queryEhProteina = ($tipo_query eq "proteina") ? 1 : 0; 

    #popula o hash %reads_mapeados_para, analisando todos os alinhamentos obtidos pelos processos filhos
    foreach my $nome_arquivo_blat (@nomes_arquivos_blat) {
        my $pula_linhas = 0;

        open (ARQUIVO_BLAT, $nome_arquivo_blat) or die "Erro na abertura do arquivo $nome_arquivo_blat: $!\n";
        while  (my $linha = <ARQUIVO_BLAT>) {
            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);

            ($qName) = ($qName =~ /([^\s\|]+)/);

            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, $queryEhProteina)/100;
            my $pontuacao = pslScore($matches, $repMatches, $misMatches,
                                     $qNumInsert, $tNumInsert, $queryEhProteina);
            my $cobertura_query = ($qEnd - $qStart) / $qSize;
            my $cobertura_read = ($tEnd - $tStart) / $tSize;

            if ($identidade >= $identidade_minima) {

                if ( $cobertura_query >= $cobertura_minima_query
                     or $cobertura_read >= $cobertura_minima_read
                     or possui_splice_site_bordas($tStart, $tEnd, $tSize, $qStart, $qEnd,$qSize,  
                           \@strand, get_sequence_by_name($tName,$hashref_read, $arquivo_reads_fh),
                           $tipo_query)
                    ) {  

                    #print("query: $qName, target: $tName, score: $pontuacao\n");              

                    if(exists $query_mapeada_para{$tName}) {
                        my ($qName_antiga, $pontuacao_antiga) = @{$query_mapeada_para{$tName}};
                        if ($pontuacao > $pontuacao_antiga) {
                            @{$query_mapeada_para{$tName}} = ($qName, $pontuacao);
                            if ($qName_antiga ne $qName) {
                                remove_string($tName,
                                              ${$hashref_reads_mapeados_para}{$qName_antiga});

                                if (not exists ${$hashref_reads_mapeados_para}{$qName}) {
                                    ${$hashref_reads_mapeados_para}{$qName} = [];
                                }
                                adiciona_string($tName, $reads_mapeados_para{$qName});
                                   
                            }
                        }

                    }
                    else { 
                        $query_mapeada_para{$tName} = [$qName, $pontuacao]; 
                        if (not exists ${$hashref_reads_mapeados_para}{$qName}) {
                            ${$hashref_reads_mapeados_para}{$qName} = [];
                        }                   
                        adiciona_string($tName, ${$hashref_reads_mapeados_para}{$qName});   
                        adiciona_string($tName, $vecref_nomes_reads_mapeados);                    
                    }

                }

            }
        }
        close (ARQUIVO_BLAT);
    }    

}


sub gera_arquivo_queries_selecionadas {
    my ($hashref_reads_mapeados_para, $hashref_query,
        $nome_arquivo_queries_selecionadas, $arquivo_queries_fh) = @_;

    open(FASTA, ">", $nome_arquivo_queries_selecionadas) or die "Erro na abertura do arquivo $nome_arquivo_queries_selecionadas : $!";

    foreach my $nome_query (keys %{$hashref_reads_mapeados_para}) {
        if (scalar @{${$hashref_reads_mapeados_para}{$nome_query}} == 0) {
            next;
        }
        print FASTA ">$nome_query\n";
        print FASTA (get_sequence_by_name($nome_query, $hashref_query, $arquivo_queries_fh));
    }

    close (FASTA);
}

#monta os arquivos multifasta com as sequencias de shotgun que foram mapeadas para a mesma sequencia-consulta
sub gera_arquivos_sementes_nao_montadas {
    my ($hashref_reads_mapeados_para, $arrayref_sementes_nao_montadas,
        $hashref_read, $arquivo_reads_fh) = @_;

    foreach my $nome_query (keys %{$hashref_reads_mapeados_para}) {
        if (scalar @{${$hashref_reads_mapeados_para}{$nome_query}} == 0) {
            next;
        }

        # print "$nome_query @{${$hashref_reads_mapeados_para}{$nome_query}}\n";

        #nome do arquivo que contem todos os reads que foram mapeados para a sequencia $nome_query
        $nome_query =~ s/_//g;
        my $nome_arquivo_query = "dataset:$nome_query";
        open(FASTA, ">",  $nome_arquivo_query) or die "Erro na abertura do arquivo $nome_arquivo_query : $!";
        foreach my $nome_read (@{${$hashref_reads_mapeados_para}{$nome_query}}) {
            print FASTA ">$nome_read|dataset:$nome_query\n";
            print FASTA (get_sequence_by_name($nome_read, $hashref_read, $arquivo_reads_fh));
        }
        close (FASTA);
        push (@{$arrayref_sementes_nao_montadas}, $nome_arquivo_query);
    }
}

#monta os arquivos multifasta com as sequencias de shotgun que nao foram mapeadas para qualquer sequencia-consulta
sub gera_arquivo_reads_nao_mapeados {
    my ($arrayref_nomes_reads_mapeados, $hashref_read,
        $arrayref_arquivos_banco, $arquivo_reads_fh) = @_;
    my $nome_arquivo_reads_nao_mapeados = "reads_nao_mapeados.mfasta";

    open(FASTA, ">", $nome_arquivo_reads_nao_mapeados) or die "Erro na abertura do arquivo $nome_arquivo_reads_nao_mapeados : $!";
    foreach my $nome_read (keys %{$hashref_read}) {
        if (posicao_string($nome_read, $arrayref_nomes_reads_mapeados) == -1) { #read nao foi mapeado
            print FASTA ">$nome_read\n";
            print FASTA (get_sequence_by_name($nome_read, $hashref_read, $arquivo_reads_fh));
        }
    }
    close (FASTA);
    push(@{$arrayref_arquivos_banco}, $nome_arquivo_reads_nao_mapeados); #usar qq("$nome_arquivo_reads_nao_mapeados") se for para usar o cat

}

#usa o phrap para montar cada uma das provaveis regioes genicas separadamente
sub roda_montador {
    my ($montador, $sequenciamento, $numero_cores,
        $sementes_ref, $sementes_nao_montadas_ref) = @_;
    my $numero_filhos_ativos = 0;


    #usa o phrap para montar cada uma das provaveis regioes genicas separadamente
    foreach my $nome_arquivo_reads_genicos (@{$sementes_nao_montadas_ref}) {
        my $saida_montador = $nome_arquivo_reads_genicos.".$montador"."_output";  
        $saida_montador = qq("$saida_montador");
        my $nome_arquivo_semente_1;
        my $nome_arquivo_semente_2;
        my $diretorio_saida_mira;
        my $arquivo_entrada_mira;
        my $diretorio_saida_newbler;
        if ($montador eq "phrap") {
            $nome_arquivo_semente_1 = $nome_arquivo_reads_genicos.".contigs";
            $nome_arquivo_semente_2 = $nome_arquivo_reads_genicos.".singlets"; 

            push(@{$sementes_ref}, $nome_arquivo_semente_2);#usar qq("$nome_arquivo_semente_2") se for para usar o cat depois 
        }
        elsif ($montador eq "mira") { 
            $diretorio_saida_mira = $nome_arquivo_reads_genicos."_mira_assembly";
            $arquivo_entrada_mira = $nome_arquivo_reads_genicos."_in.$sequenciamento.fasta";
            $nome_arquivo_semente_1 = "$diretorio_saida_mira/$nome_arquivo_reads_genicos"."_assembly/$nome_arquivo_reads_genicos"."_d_results/$nome_arquivo_reads_genicos"."_out.unpadded.fasta";
        }
        else { #$montador eq "newbler"
            $diretorio_saida_newbler = $nome_arquivo_reads_genicos."_newbler_assembly";
            $nome_arquivo_semente_1 = "$diretorio_saida_newbler/454AllContigs.fna";
        }    
        $nome_arquivo_reads_genicos = qq("$nome_arquivo_reads_genicos");  

        push(@{$sementes_ref}, $nome_arquivo_semente_1);#usar qq("$nome_arquivo_semente_1") se for para usar o cat depois


        #print "|$nome_arquivo_reads_genicos|";
        #print "\n";    
        #print qq("$nome_arquivo_semente_1");
        #print "\n";    
        #print qq("$nome_arquivo_semente_2");
        #print "\n";    
        if ($numero_filhos_ativos >= $numero_cores) { #espera um dos filhos terminar
            wait();
            $numero_filhos_ativos--;
        }        

        #cria mais um filho
        my $pid = fork();
        if ($pid == 0) { #processo filho
            if ($montador eq "phrap") {
                exec("phrap $nome_arquivo_reads_genicos > $saida_montador") or
                      print STDERR "Erro na execucao do phrap: $!";
            }
            elsif ($montador eq "mira") {
                mkdir($diretorio_saida_mira) unless -d $diretorio_saida_mira;
                system ("mv ".$nome_arquivo_reads_genicos." ".qq("$diretorio_saida_mira/$arquivo_entrada_mira")); 
                chdir($diretorio_saida_mira);
                exec("mira --project=".$nome_arquivo_reads_genicos." --job=denovo,genome,accurate,$sequenciamento --fasta --noqualities ".uc($sequenciamento)."_SETTINGS -LR:mxti=no -AS:bdq=30 -OUT:sssip=yes > $saida_montador") or
                      print STDERR "Erro na execucao do mira: $!";
            }
            else { #$montador eq "newbler"
                $diretorio_saida_newbler =~ s/\|/\\|/g;
                $diretorio_saida_newbler = qq("$diretorio_saida_newbler"); 

                print("vou executar |runAssembly -cpu 1 -o $diretorio_saida_newbler -force $nome_arquivo_reads_genicos > $saida_montador|\n");
                exec("runAssembly -cpu 1 -o $diretorio_saida_newbler -force $nome_arquivo_reads_genicos > $saida_montador") or
                      print STDERR "Erro na execucao do newbler: $!";
            }
        }
        elsif(not defined $pid) {
            print STDERR "Erro: fork() nao funcionou\n";    
            exit(1);
        }
        else { #processo pai
            $numero_filhos_ativos++;
        }
    }

    #espera todos os filhos terminarem
    while(wait() != -1) {}

}

sub monta_arquivos_genseed {
    my($sementes_ref, $arquivos_banco_ref, $diretorio_saida_genseed, $numero_maximo_rodadas_genseed) = @_;

    {
        local $/;#ativa "slurp mode" (le um arquivo inteiro a cada <>)

        #gera arquivo com as sementes do genseed (concatenando os arquivos em @{$sementes_ref})
        open (SEMENTES, ">", $nome_arquivo_sementes) or die "Erro na abertura do arquivo $nome_arquivo_sementes: $!\n";
        foreach my $semente (@{$sementes_ref}) {
            open(SEMENTE, "<", $semente) or die "Erro na abertura do arquivo $semente: $!\n";;
            my $conteudo = <SEMENTE>;
            print SEMENTES $conteudo;
            close(SEMENTE);
        }
        close(SEMENTES);

        #gera arquivo com o banco do genseed (concatenando os arquivos em @{$arquivos_banco_ref})

        open (BANCO, ">", $nome_arquivo_banco) or die "Erro na abertura do arquivo $nome_arquivo_banco: $!\n";
        foreach my $item_banco (@{$arquivos_banco_ref}) {
            open(ITEM, "<", $item_banco) or die "Erro na abertura do arquivo $item_banco: $!\n";;
            my $conteudo = <ITEM>;
            print BANCO $conteudo;
            close(ITEM);
        }
        close(BANCO);
    }
}

sub roda_genseed {
    my ($nome_arquivo_sementes, $nome_arquivo_banco, $diretorio_saida_genseed, $numero_maximo_rodadas_genseed) = @_;
    system ("genseed -s $nome_arquivo_sementes -d $nome_arquivo_banco -o $diretorio_saida_genseed -r $numero_maximo_rodadas_genseed");
    #system ("/home/pedro/genseedHMM/genseed-hmm_31_07.pl -seed $nome_arquivo_sementes -db $nome_arquivo_banco -assembler_name Newbler -output $diretorio_saida_genseed -max_number_rounds $numero_maximo_rodadas_genseed -next_seed_size 10");
}

#############################
#                           # 
#  sub-rotinas auxiliares   #
#                           #
#############################

#devolve uma sequencia de um arquivo fasta indexado por um hash que possui as posicoes das sequencias no arquivo
sub get_sequence_by_name{
    my ($seqid, $hash_ref, $arquivo_fh) = @_;
    my $fasta;
    my ($header, $sequence);
    {
        local $/ = ">";
        seek($arquivo_fh, ${$hash_ref}{$seqid},0);
        $fasta = <$arquivo_fh>;
        chomp($fasta);
        ($header, $sequence) = split("\n", $fasta, 2);
    }
    die "Erro: sequencia $seqid esta vazia\n" if $sequence eq "";
    return $sequence;
}

#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;
}

#funcao auxiliar no calculo da porcentagem de identidade, retirada de http://genome.ucsc.edu/FAQ/FAQblat.html#blat4 (Replicating web-based Blat percent identity and score calculations)
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 pslScore {
    my ($match, $repMatch, $misMatch, $qNumInsert, $tNumInsert, $isQueryProtein) = @_;
    my $sizeMul = $isQueryProtein ? 3 : 1;

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

#devolve TRUE (1) se um alinhamento possuir DNA genomico nao mapeado que seja compativel com
#regioes intronicas, e FALSE (0) caso contrario
sub possui_splice_site_bordas {

    my ($tstart, $tend, $tsize, $qstart,$qend,$qsize, $vecref_strand, $target, $tipo_query) = @_;

    my $possui_ss_borda_esquerda = 0;
    my $possui_ss_borda_direita = 0;
    my $indicio_ss_borda_esquerda = 0;
    my $indicio_ss_borda_direita = 0;

    #valores abaixo retirados (ou deduzidos) do script "blat2hints.pl" (http://augustus.gobics.de/binaries/scripts/blat2hints.pl)
    my $maxDNAoverhang = 10000;
    my $minQUERYoverhang = ($tipo_query eq "full-length-est")? 40 : 40/3;


    my $ss_5linha = "GT";
    my $ss_3linha = "AG";

    if ($tipo_query eq "est") {
        return 0; #nao faz sentido exigir sitios de splicing na borda do alinhamento se a query for uma EST
    }

    if(scalar @{$vecref_strand} == 2 and ${$vecref_strand}[1] eq "-") { #target possui orientacao negativa 
        $ss_5linha = "CT";
        $ss_3linha = "AC";
    }

    # indicios de splice sites quando o "cDNA e' maior que o gDNA (DNA genomico)":
    # - existe cDNA nao alinhado 
    # - o DNA genomico termina assim que o alinhamento acaba

    # ponta genomica esquerda
    if ($tstart < $maxDNAoverhang && $tstart > 5 &&
        ((${$vecref_strand}[0] eq '+' && $qstart > $minQUERYoverhang) ||
         (${$vecref_strand}[0] eq '-' && $qsize-$qend > $minQUERYoverhang))){
        $indicio_ss_borda_esquerda  = 1;
        # pequena janela [-4,+4], que permite ate 4 "matches" aleatorios ou 4 "mismatches"
        my ($inicio, $fim) = ($tstart-4, $tstart+4);   
        $possui_ss_borda_esquerda = (substr($target, $inicio, $fim - $inicio)=~ /$ss_5linha/i);
    }

	# ponta genomica direita
	if ($tsize-$tend < $maxDNAoverhang && $tsize-$tend > 5 &&
	    ((${$vecref_strand}[0] eq '+' && $qsize-$qend > $minQUERYoverhang) ||
         (${$vecref_strand}[0] eq '-' && $qstart > $minQUERYoverhang))){
        $indicio_ss_borda_direita  = 1;
        # pequena janela [-4,+4], que permite ate 4 "matches" aleatorios ou 4 "mismatches"
        my ($inicio, $fim) = ($tend-4, $tend+4);  
        $possui_ss_borda_direita = (substr($target, $inicio, $fim - $inicio)=~ /$ss_3linha/i);
    }
    
     
    if ($indicio_ss_borda_esquerda and $indicio_ss_borda_direita) {
        return $possui_ss_borda_esquerda and $possui_ss_borda_direita;
    }
    elsif ($indicio_ss_borda_esquerda and not $indicio_ss_borda_direita) {
        return $possui_ss_borda_esquerda; 
    }
    elsif (not $indicio_ss_borda_esquerda and $indicio_ss_borda_direita) {
        return $possui_ss_borda_direita; 
    }
    else { # not $indicio_ss_borda_esquerda and not $indicio_ss_borda_direita
        return 0;
    }

}

#devolve o complemento reverso de uma sequencia de DNA
sub revcom {
    my($dna) = @_;

    # First reverse the sequence
    my $revcom = reverse $dna;

    # Next, complement the sequence, dealing with upper and lower case
    # A->T, T->A, C->G, G->C
    $revcom =~ tr/ACGTacgt/TGCAtgca/;

    return $revcom;
}

#remove a string "string" do vetor
sub remove_string {
    my ($string, $array_ref) = @_;
    splice (@{$array_ref}, posicao_string($string, $array_ref), 1);
}


#adiciona uma string a um vetor de forma a manter a ordenacao
sub adiciona_string {
    my ($string, $array_ref) = @_;
    splice (@{$array_ref}, posicao_insercao_string($string, $array_ref) , 0, $string);
}


#procura a posicao de insercao de uma string num vetor ordendado 
sub posicao_insercao_string {
    my ($string, $vecref) = @_;   

    my $inicio = 0;
    my $fim = $#{$vecref};
    my $meio;
    while($inicio <= $fim){
        $meio = int(($inicio + $fim)/2);
        if($string eq ${$vecref}[$meio] ){
            return $meio;
        }
        elsif($string lt ${$vecref}[$meio]){
            $fim = $meio - 1;
        }
        else{
            $inicio = $meio + 1;
        }
    }
    return $inicio;
}

#procura uma string num vetor ordendado, devolvendo sua posicao caso ela esteja no vetor ou -1
#caso contrario
sub posicao_string {
    my ($string, $vecref) = @_;      

    my $inicio = 0;
    my $fim = $#{$vecref};
    my $meio;
    while($inicio <= $fim){
        $meio = int(($inicio + $fim)/2);
        if($string eq ${$vecref}[$meio]){
            return $meio;
        }
        elsif($string lt ${$vecref}[$meio]){
            $fim = $meio - 1;
        }
        else{
            $inicio = $meio + 1;
        }
    }
    return -1;
}
