<commit_before>package Tephra::Genome::SoloLTRSearch;

use 5.014;
use Moose;
use Cwd;
use File::Spec;
use File::Find;
use File::Basename;
use File::Copy;
use File::Path          qw(make_path remove_tree);
use IPC::System::Simple qw(capture system);
use Time::HiRes         qw(gettimeofday);
use List::UtilsBy       qw(nsort_by);
use Sort::Naturally;
use Path::Class::File;
use Bio::DB::HTS::Kseq;
use Bio::AlignIO;
use Bio::SearchIO;
use Sort::Naturally;
use Set::IntervalTree;
use Carp 'croak';
use Tephra::Config::Exe;
use namespace::autoclean;
#use Data::Dump::Color;

with 'Tephra::Role::Util';

=head1 NAME

Tephra::Genome::SoloLTRSearch - Find solo-LTRs in a refence genome

=head1 VERSION

Version 0.04.6

=cut

our $VERSION = '0.04.6';
$VERSION = eval $VERSION;

has dir => (
      is       => 'ro',
      isa      => 'Path::Class::Dir',
      required => 1,
      coerce   => 1,
);

has genome => (
    is       => 'ro',
    isa      => 'Path::Class::File',
    required => 1,
    coerce   => 1,
);

has report => (
    is       => 'ro',
    isa      => 'Path::Class::File',
    required => 0,
    coerce   => 1,
);

has outfile => (
    is       => 'ro',
    isa      => 'Path::Class::File',
    required => 1,
    coerce   => 1,
);

has seqfile => (
    is       => 'ro',
    isa      => 'Path::Class::File',
    required => 0,
    coerce   => 1,
);

has percentident => (
    is        => 'ro',
    isa       => 'Num',
    predicate => 'has_percentident',
    lazy      => 1,
    default   => 0.39,
);

has percentcov => (
    is        => 'ro',
    isa       => 'Num',
    predicate => 'has_percentcov',
    lazy      => 1,
    default   => 0.80,
);

has matchlen => (
    is        => 'ro',
    isa       => 'Num',
    predicate => 'has_matchlen',
    lazy      => 1,
    default   => 80,
);

has clean => (
    is       => 'ro',
    isa      => 'Bool',
    required => 0,
    default  => 1,
);

has fullanalysis => (
    is       => 'ro',
    isa      => 'Bool',
    required => 0,
    default  => 0,
);

has numfamilies => (
    is        => 'ro',
    isa       => 'Int',
    predicate => 'has_numfamilies',
    lazy      => 1,
    default   => 20,
);

has threads => (
    is        => 'ro',
    isa       => 'Int',
    predicate => 'has_threads',
    lazy      => 1,
    default   => 1,
);

sub find_soloLTRs {
    my $self = shift;
    my $anno_dir = $self->dir;
    my $genome = $self->genome;

    my ($name, $path, $suffix) = fileparse($genome, qr/\.[^.]*/);                                                                
    my $hmmsearch_summary = $self->report // File::Spec->catfile($path, $name.'_tephra_soloLTRs.tsv'); 

    my @sfs;
    find( sub { push @sfs, $File::Find::name if /_copia\z|_gypsy\z/ }, $anno_dir);
    croak "\nERROR: Could not find the expected sub-directories ending in 'copia' and 'gypsy' please ".
	"check input. Exiting.\n" unless @sfs == 2;

    my $pm = Parallel::ForkManager->new(2);
    local $SIG{INT} = sub {
        warn "Caught SIGINT; Waiting for child processes to finish.";
        $pm->wait_all_children;
        exit 1;
    };

    $pm->run_on_finish( sub { my ($pid, $exit_code, $ident, $exit_signal, $core_dump, $data_ref) = @_;
			      my ($model_dir, $reports) = @{$data_ref}{qw(model_dir reports)};
			      $self->_collate($reports, $hmmsearch_summary);
			      if ($self->clean) {
				  remove_tree( $model_dir, { safe => 1} );
			      }
			} );

    my $dirct = 0;
    for my $dir (nsort @sfs) {
	my $sf = (split /_/, $dir)[-1];
	$pm->start($sf) and next;
	$SIG{INT} = sub { $pm->finish };
	my $data_ref = $self->run_sf_search($sf, $dir, $dirct, $hmmsearch_summary);
	$dirct++;
	$pm->finish(0, $data_ref);
    }
    $pm->wait_all_children;

    if (-s $hmmsearch_summary) {
	$self->write_sololtr_gff($hmmsearch_summary);
    }

    return;
}

sub run_sf_search {
    my $self = shift;
    my ($sf, $dir, $dirct, $hmmsearch_summary) = @_;
    my $genome = $self->genome;
    
    my ($hmmbuild, $hmmsearch) = $self->_find_hmmer;

    #print STDERR "Getting LTR alignments for ",ucfirst($sf),"....";
    my $ltr_aln_files = $self->_get_ltr_alns($dir);
    #say STDERR "done with alignments.";

    ## need masked genome here
    if (!defined $ltr_aln_files || @$ltr_aln_files < 1 || ! -e $genome) {
	croak "\nERROR: No genome was found or the expected alignments files were not found for ",
	    ucfirst($sf),". Exiting.";
    }
        
    #print STDERR "Getting alignment statistics for ",ucfirst($sf),"...";
    my $aln_stats = $self->_get_aln_len($ltr_aln_files); # return a hash-ref
    #say STDERR "done with alignment statistics.";
        
    # make one directory
    my $model_dir = File::Spec->catdir($dir, 'Tephra_LTR_exemplar_models');
    unless ( -e $model_dir ) {
	make_path( $model_dir, {verbose => 0, mode => 0771,} );
    }
        
    #print STDERR "Building LTR exemplar models for ",ucfirst($sf),"...";
    for my $ltr_aln (nsort_by { m/family(\d+)/ and $1 } @$ltr_aln_files) {
	$self->build_model($ltr_aln, $model_dir, $hmmbuild);        
	unlink $ltr_aln;    
    }
    #say STDERR "done with models.";
        
    my @ltr_hmm_files;
    find( sub { push @ltr_hmm_files, $File::Find::name if -f and /\.hmm$/ }, $model_dir);
        
    #print STDERR "Search genome with models for ",ucfirst($sf),"...";
    $self->do_parallel_search($hmmsearch, $genome, \@ltr_hmm_files, $model_dir, $aln_stats);
    #say STDERR "done searching with LTR models.";

    #print STDERR "Collating results and writing GFF...";
    my @reports;
    find( sub { push @reports, $File::Find::name if -f and /\.txt$/ }, $model_dir );

    return { reports => \@reports, model_dir => $model_dir };
}
    

sub do_parallel_search {
    my $self = shift;
    my ($hmmsearch, $genome, $ltr_hmm_files, $model_dir, $aln_stats) = @_;
    my $threads = $self->threads;

    my %reports;
    my $t0 = gettimeofday();
    my $model_num = @$ltr_hmm_files;

    my ($gname, $gpath, $gsuffix) = fileparse($genome, qr/\.[^.]*/);
    my $logfile = File::Spec->catfile($model_dir, 'all_solo-ltr_searches.log');
    open my $log, '>>', $logfile or die "\nERROR: Could not open file: $logfile\n";

    my $thr;
    if ($threads % 2 == 0) {
        $thr = sprintf("%.0f",$threads/2);
    }
    elsif ($threads-1 % 2 == 0) {
        $thr = sprintf("%.0f",$threads-1/2);
    }
    else {
        $thr = 1;
    }

    my $pm = Parallel::ForkManager->new($threads);
    local $SIG{INT} = sub {
        warn "Caught SIGINT; Waiting for child processes to finish.";
        $pm->wait_all_children;
        exit 1;
    };

    $pm->run_on_finish( sub { my ($pid, $exit_code, $ident, $exit_signal, $core_dump, $data_ref) = @_;
			      for my $hmmrep (sort keys %$data_ref) {
				  unlink $hmmrep;
			      }
			      my $t1 = gettimeofday();
			      my $elapsed = $t1 - $t0;
			      my $time = sprintf("%.2f",$elapsed/60);
			      say $log basename($ident),
			        " just finished with PID $pid and exit code: $exit_code in $time minutes";
			} );

    for my $hmm (nsort @$ltr_hmm_files) {
	$pm->start($hmm) and next;
	$SIG{INT} = sub { $pm->finish };
	
	my $hmmsearch_out = $self->search_with_models($gname, $hmm, $hmmsearch, $genome, $aln_stats);
	$reports{$hmmsearch_out} = 1;

	$pm->finish(0, \%reports);
    }

    $pm->wait_all_children;

    my $t2 = gettimeofday();
    my $total_elapsed = $t2 - $t0;
    my $final_time = sprintf("%.2f",$total_elapsed/60);

    say $log "\n========> Finished hmmsearch on $model_num solo-LTR models in $final_time minutes";
    close $log;
}

sub write_sololtr_gff {
    my $self = shift;
    my ($hmmsearch_summary) = @_;
    my $genome  = $self->genome;
    my $outfile = $self->outfile;

    my $seqlen = $self->_get_seq_len($genome);
    open my $in, '<', $hmmsearch_summary or die "\nERROR: Could not open file: $hmmsearch_summary\n";
    open my $out, '>>', $outfile or die "\nERROR: Could not open file: $outfile\n";

    if (! -s $outfile) {
	say $out '##gff-version 3';
	for my $id (nsort keys %$seqlen) {
	    say $out join q{ }, '##sequence-region', $id, '1', $seqlen->{$id};
	}
    }

    my $ct = 0;
    my (%features, %intervals, %ltr_lengths);
    while (my $line = <$in>) {
	chomp $line;
	next if $line =~ /^#/;
	my ($query, $query_length, $number_of_hits, $hit_name, $perc_coverage, $hsp_length, 
	    $hsp_perc_ident, $hsp_query_start, $hsp_query_end, $hsp_hit_start, $hsp_hit_end,
	    $search_type) = split /\t/, $line;
	$query =~ s/_ltrs_muscle-out//;

	if (exists $seqlen->{$hit_name}) {
	    $ct++;
	    say $out join "\t", $hit_name, 'Tephra', 'solo_LTR', $hsp_hit_start, $hsp_hit_end, 
	        '.', '?', '.', "ID=solo_LTR$ct;Parent=$query;Name=solo_LTR;Ontology_term=SO:0001003";
	}
	else {
	    die "\nERROR: $hit_name not found in $genome. This should not happen. Exiting.\n";
	}
    }
    close $in;
    close $out; 

    return;
}

sub build_model {
    my $self = shift;
    my ($aln, $model_dir, $hmmbuild) = @_;

    my $hmmname = basename($aln);
    $hmmname =~ s/\.aln$/\.hmm/;
    my $model_path = File::Spec->catfile($model_dir, $hmmname);
    # hmmer3 has --dna and --cpu options and not -g (global)
    my $hmm_cmd = "$hmmbuild -f --nucleic $model_path $aln";

    my $status = $self->capture_cmd($hmm_cmd);
    unlink $aln && return if $status =~ /failed/i;
}

sub search_with_models {
    my $self = shift;
    my ($gname, $hmm, $hmmsearch, $fasta, $aln_stats) = @_;

    my $hmmsearch_out = $hmm;
    $hmmsearch_out =~ s/\.hmm.*$//;
    $hmmsearch_out .= '_'.$gname.'.hmmer';

    ## no -o option in hmmer2
    my $hmmsearch_cmd = "$hmmsearch --cpu 1 $hmm $fasta > $hmmsearch_out";
    $self->run_cmd($hmmsearch_cmd);

    if (-e $hmmsearch_out) {   
	$self->write_hmmsearch_report($aln_stats, $hmmsearch_out);
    }

    unlink $hmm;
    return $hmmsearch_out;
}

sub write_hmmsearch_report {
    my $self = shift;
    my ($aln_stats, $search_report) = @_;
    my $genome     = $self->genome;
    my $match_pid  = $self->percentident;
    my $match_len  = $self->matchlen;
    my $match_pcov = $self->percentcov;

    my $parsed = $search_report;
    $parsed =~ s/\.hmmer$/\_hmmer_parsed.txt/;
    my $element = $search_report;

    my $model_type = 'local';

    my ($gname, $gpath, $gsuffix) = fileparse($genome, qr/\.[^.]*/);
    my ($name, $path, $suffix)    = fileparse($search_report, qr/\.[^.]*/);
    $element =~ s/_$gname*//;
    open my $out, ">>$parsed" or die "\nERROR: Could not open file: $parsed\n";

    my ($seq, $seqfile);
    if ($self->seqfile) {
	$seqfile = $self->seqfile;
	open $seq, ">>$seqfile" or die "\nERROR: Could not open file: $seqfile";
    }

    my $hmmerin = Bio::SearchIO->new(-file => $search_report, -format => 'hmmer');

    my ($positions, $matches) = (0, 0);
    while ( my $result = $hmmerin->next_result() ) {
	my $query    = $result->query_name();
	my $num_hits = $result->num_hits();
	my $qlen     = $aln_stats->{$query};
	die "\nERROR: Could not determine query length for: $query" 
	    unless defined $qlen;
	while ( my $hit = $result->next_hit() ) {
	    my $hitid = $hit->name();
	    while ( my $hsp = $hit->next_hsp() ) {
		my $percent_q_coverage = sprintf("%.2f", $hsp->length('query')/$qlen * 100);
		my @ident_pos = $hsp->seq_inds('query','identical');
		my $hspgaps   = $hsp->gaps;
		my $hsplen    = $hsp->length('total');
		my $hstart    = $hsp->start('hit');
		my $hstop     = $hsp->end('hit');
		my $qstart    = $hsp->start('query');
		my $qstop     = $hsp->end('query');
		my $qstring   = $hsp->query_string;
		my $qhsplen   = $hsp->length('query');
		
		if (exists $aln_stats->{$query}) {
		    my $percent_coverage = sprintf("%.2f",$hsplen/$aln_stats->{$query});
		    $positions++ for @ident_pos;
		    my $pid = sprintf("%.2f", $positions/$aln_stats->{$query} * 100);
		    if ( $hsplen >= $match_len && $hsplen >= $aln_stats->{$query} * $match_pcov ) {
			my $qid = $query =~ s/_ltrseqs_muscle-out//r; # non-destructive substitution in v5.14+
			if ($pid >= $match_pid) {
			    $matches++;
			    say $out join "\t", $qid, $aln_stats->{$query}, $num_hits, $hitid, 
			        $percent_q_coverage, $hsplen, $pid, $qstart, $qstop, $hstart, $hstop, $model_type;

			    if ($self->seqfile) {
				my $seqid = '>'.$qid.'_'.$hitid.'_'.$hstart.'_'.$hstop; 
				## It makes more sense to show the location of the hit
				## Also, this would pave the way for creating a gff of solo-LTRs
				## my $seqid = ">".$query."|".$hitid."_".$hstart."-".$hstop
				say $seq join "\n", $seqid, $qstring;
			    }
			}
		    }
		    $positions = 0;
		}
	    }
	}
    }
    close $out;
    close $seq;<commit_msg>use join instead of concat for id<commit_after>package Tephra::Genome::SoloLTRSearch;

use 5.014;
use Moose;
use Cwd;
use File::Spec;
use File::Find;
use File::Basename;
use File::Copy;
use File::Path          qw(make_path remove_tree);
use IPC::System::Simple qw(capture system);
use Time::HiRes         qw(gettimeofday);
use List::UtilsBy       qw(nsort_by);
use Sort::Naturally;
use Path::Class::File;
use Bio::DB::HTS::Kseq;
use Bio::AlignIO;
use Bio::SearchIO;
use Sort::Naturally;
use Set::IntervalTree;
use Carp 'croak';
use Tephra::Config::Exe;
use namespace::autoclean;
#use Data::Dump::Color;

with 'Tephra::Role::Util';

=head1 NAME

Tephra::Genome::SoloLTRSearch - Find solo-LTRs in a refence genome

=head1 VERSION

Version 0.04.6

=cut

our $VERSION = '0.04.6';
$VERSION = eval $VERSION;

has dir => (
      is       => 'ro',
      isa      => 'Path::Class::Dir',
      required => 1,
      coerce   => 1,
);

has genome => (
    is       => 'ro',
    isa      => 'Path::Class::File',
    required => 1,
    coerce   => 1,
);

has report => (
    is       => 'ro',
    isa      => 'Path::Class::File',
    required => 0,
    coerce   => 1,
);

has outfile => (
    is       => 'ro',
    isa      => 'Path::Class::File',
    required => 1,
    coerce   => 1,
);

has seqfile => (
    is       => 'ro',
    isa      => 'Path::Class::File',
    required => 0,
    coerce   => 1,
);

has percentident => (
    is        => 'ro',
    isa       => 'Num',
    predicate => 'has_percentident',
    lazy      => 1,
    default   => 0.39,
);

has percentcov => (
    is        => 'ro',
    isa       => 'Num',
    predicate => 'has_percentcov',
    lazy      => 1,
    default   => 0.80,
);

has matchlen => (
    is        => 'ro',
    isa       => 'Num',
    predicate => 'has_matchlen',
    lazy      => 1,
    default   => 80,
);

has clean => (
    is       => 'ro',
    isa      => 'Bool',
    required => 0,
    default  => 1,
);

has fullanalysis => (
    is       => 'ro',
    isa      => 'Bool',
    required => 0,
    default  => 0,
);

has numfamilies => (
    is        => 'ro',
    isa       => 'Int',
    predicate => 'has_numfamilies',
    lazy      => 1,
    default   => 20,
);

has threads => (
    is        => 'ro',
    isa       => 'Int',
    predicate => 'has_threads',
    lazy      => 1,
    default   => 1,
);

sub find_soloLTRs {
    my $self = shift;
    my $anno_dir = $self->dir;
    my $genome = $self->genome;

    my ($name, $path, $suffix) = fileparse($genome, qr/\.[^.]*/);                                                                
    my $hmmsearch_summary = $self->report // File::Spec->catfile($path, $name.'_tephra_soloLTRs.tsv'); 

    my @sfs;
    find( sub { push @sfs, $File::Find::name if /_copia\z|_gypsy\z/ }, $anno_dir);
    croak "\nERROR: Could not find the expected sub-directories ending in 'copia' and 'gypsy' please ".
	"check input. Exiting.\n" unless @sfs == 2;

    my $pm = Parallel::ForkManager->new(2);
    local $SIG{INT} = sub {
        warn "Caught SIGINT; Waiting for child processes to finish.";
        $pm->wait_all_children;
        exit 1;
    };

    $pm->run_on_finish( sub { my ($pid, $exit_code, $ident, $exit_signal, $core_dump, $data_ref) = @_;
			      my ($model_dir, $reports) = @{$data_ref}{qw(model_dir reports)};
			      $self->_collate($reports, $hmmsearch_summary);
			      if ($self->clean) {
				  remove_tree( $model_dir, { safe => 1} );
			      }
			} );

    my $dirct = 0;
    for my $dir (nsort @sfs) {
	my $sf = (split /_/, $dir)[-1];
	$pm->start($sf) and next;
	$SIG{INT} = sub { $pm->finish };
	my $data_ref = $self->run_sf_search($sf, $dir, $dirct, $hmmsearch_summary);
	$dirct++;
	$pm->finish(0, $data_ref);
    }
    $pm->wait_all_children;

    if (-s $hmmsearch_summary) {
	$self->write_sololtr_gff($hmmsearch_summary);
    }

    return;
}

sub run_sf_search {
    my $self = shift;
    my ($sf, $dir, $dirct, $hmmsearch_summary) = @_;
    my $genome = $self->genome;
    
    my ($hmmbuild, $hmmsearch) = $self->_find_hmmer;

    #print STDERR "Getting LTR alignments for ",ucfirst($sf),"....";
    my $ltr_aln_files = $self->_get_ltr_alns($dir);
    #say STDERR "done with alignments.";

    ## need masked genome here
    if (!defined $ltr_aln_files || @$ltr_aln_files < 1 || ! -e $genome) {
	croak "\nERROR: No genome was found or the expected alignments files were not found for ",
	    ucfirst($sf),". Exiting.";
    }
        
    #print STDERR "Getting alignment statistics for ",ucfirst($sf),"...";
    my $aln_stats = $self->_get_aln_len($ltr_aln_files); # return a hash-ref
    #say STDERR "done with alignment statistics.";
        
    # make one directory
    my $model_dir = File::Spec->catdir($dir, 'Tephra_LTR_exemplar_models');
    unless ( -e $model_dir ) {
	make_path( $model_dir, {verbose => 0, mode => 0771,} );
    }
        
    #print STDERR "Building LTR exemplar models for ",ucfirst($sf),"...";
    for my $ltr_aln (nsort_by { m/family(\d+)/ and $1 } @$ltr_aln_files) {
	$self->build_model($ltr_aln, $model_dir, $hmmbuild);        
	unlink $ltr_aln;    
    }
    #say STDERR "done with models.";
        
    my @ltr_hmm_files;
    find( sub { push @ltr_hmm_files, $File::Find::name if -f and /\.hmm$/ }, $model_dir);
        
    #print STDERR "Search genome with models for ",ucfirst($sf),"...";
    $self->do_parallel_search($hmmsearch, $genome, \@ltr_hmm_files, $model_dir, $aln_stats);
    #say STDERR "done searching with LTR models.";

    #print STDERR "Collating results and writing GFF...";
    my @reports;
    find( sub { push @reports, $File::Find::name if -f and /\.txt$/ }, $model_dir );

    return { reports => \@reports, model_dir => $model_dir };
}
    

sub do_parallel_search {
    my $self = shift;
    my ($hmmsearch, $genome, $ltr_hmm_files, $model_dir, $aln_stats) = @_;
    my $threads = $self->threads;

    my %reports;
    my $t0 = gettimeofday();
    my $model_num = @$ltr_hmm_files;

    my ($gname, $gpath, $gsuffix) = fileparse($genome, qr/\.[^.]*/);
    my $logfile = File::Spec->catfile($model_dir, 'all_solo-ltr_searches.log');
    open my $log, '>>', $logfile or die "\nERROR: Could not open file: $logfile\n";

    my $thr;
    if ($threads % 2 == 0) {
        $thr = sprintf("%.0f",$threads/2);
    }
    elsif ($threads-1 % 2 == 0) {
        $thr = sprintf("%.0f",$threads-1/2);
    }
    else {
        $thr = 1;
    }

    my $pm = Parallel::ForkManager->new($threads);
    local $SIG{INT} = sub {
        warn "Caught SIGINT; Waiting for child processes to finish.";
        $pm->wait_all_children;
        exit 1;
    };

    $pm->run_on_finish( sub { my ($pid, $exit_code, $ident, $exit_signal, $core_dump, $data_ref) = @_;
			      for my $hmmrep (sort keys %$data_ref) {
				  unlink $hmmrep;
			      }
			      my $t1 = gettimeofday();
			      my $elapsed = $t1 - $t0;
			      my $time = sprintf("%.2f",$elapsed/60);
			      say $log basename($ident),
			        " just finished with PID $pid and exit code: $exit_code in $time minutes";
			} );

    for my $hmm (nsort @$ltr_hmm_files) {
	$pm->start($hmm) and next;
	$SIG{INT} = sub { $pm->finish };
	
	my $hmmsearch_out = $self->search_with_models($gname, $hmm, $hmmsearch, $genome, $aln_stats);
	$reports{$hmmsearch_out} = 1;

	$pm->finish(0, \%reports);
    }

    $pm->wait_all_children;

    my $t2 = gettimeofday();
    my $total_elapsed = $t2 - $t0;
    my $final_time = sprintf("%.2f",$total_elapsed/60);

    say $log "\n========> Finished hmmsearch on $model_num solo-LTR models in $final_time minutes";
    close $log;
}

sub write_sololtr_gff {
    my $self = shift;
    my ($hmmsearch_summary) = @_;
    my $genome  = $self->genome;
    my $outfile = $self->outfile;

    my $seqlen = $self->_get_seq_len($genome);
    open my $in, '<', $hmmsearch_summary or die "\nERROR: Could not open file: $hmmsearch_summary\n";
    open my $out, '>>', $outfile or die "\nERROR: Could not open file: $outfile\n";

    if (! -s $outfile) {
	say $out '##gff-version 3';
	for my $id (nsort keys %$seqlen) {
	    say $out join q{ }, '##sequence-region', $id, '1', $seqlen->{$id};
	}
    }

    my $ct = 0;
    my (%features, %intervals, %ltr_lengths);
    while (my $line = <$in>) {
	chomp $line;
	next if $line =~ /^#/;
	my ($query, $query_length, $number_of_hits, $hit_name, $perc_coverage, $hsp_length, 
	    $hsp_perc_ident, $hsp_query_start, $hsp_query_end, $hsp_hit_start, $hsp_hit_end,
	    $search_type) = split /\t/, $line;
	$query =~ s/_ltrs_muscle-out//;

	if (exists $seqlen->{$hit_name}) {
	    $ct++;
	    say $out join "\t", $hit_name, 'Tephra', 'solo_LTR', $hsp_hit_start, $hsp_hit_end, 
	        '.', '?', '.', "ID=solo_LTR$ct;Parent=$query;Name=solo_LTR;Ontology_term=SO:0001003";
	}
	else {
	    die "\nERROR: $hit_name not found in $genome. This should not happen. Exiting.\n";
	}
    }
    close $in;
    close $out; 

    return;
}

sub build_model {
    my $self = shift;
    my ($aln, $model_dir, $hmmbuild) = @_;

    my $hmmname = basename($aln);
    $hmmname =~ s/\.aln$/\.hmm/;
    my $model_path = File::Spec->catfile($model_dir, $hmmname);
    # hmmer3 has --dna and --cpu options and not -g (global)
    my $hmm_cmd = "$hmmbuild -f --nucleic $model_path $aln";

    my $status = $self->capture_cmd($hmm_cmd);
    unlink $aln && return if $status =~ /failed/i;
}

sub search_with_models {
    my $self = shift;
    my ($gname, $hmm, $hmmsearch, $fasta, $aln_stats) = @_;

    my $hmmsearch_out = $hmm;
    $hmmsearch_out =~ s/\.hmm.*$//;
    $hmmsearch_out .= '_'.$gname.'.hmmer';

    ## no -o option in hmmer2
    my $hmmsearch_cmd = "$hmmsearch --cpu 1 $hmm $fasta > $hmmsearch_out";
    $self->run_cmd($hmmsearch_cmd);

    if (-e $hmmsearch_out) {   
	$self->write_hmmsearch_report($aln_stats, $hmmsearch_out);
    }

    unlink $hmm;
    return $hmmsearch_out;
}

sub write_hmmsearch_report {
    my $self = shift;
    my ($aln_stats, $search_report) = @_;
    my $genome     = $self->genome;
    my $match_pid  = $self->percentident;
    my $match_len  = $self->matchlen;
    my $match_pcov = $self->percentcov;

    my $parsed = $search_report;
    $parsed =~ s/\.hmmer$/\_hmmer_parsed.txt/;
    my $element = $search_report;

    my $model_type = 'local';

    my ($gname, $gpath, $gsuffix) = fileparse($genome, qr/\.[^.]*/);
    my ($name, $path, $suffix)    = fileparse($search_report, qr/\.[^.]*/);
    $element =~ s/_$gname*//;
    open my $out, ">>$parsed" or die "\nERROR: Could not open file: $parsed\n";

    my ($seq, $seqfile);
    if ($self->seqfile) {
	$seqfile = $self->seqfile;
	open $seq, ">>$seqfile" or die "\nERROR: Could not open file: $seqfile";
    }

    my $hmmerin = Bio::SearchIO->new(-file => $search_report, -format => 'hmmer');

    my ($positions, $matches) = (0, 0);
    while ( my $result = $hmmerin->next_result() ) {
	my $query    = $result->query_name();
	my $num_hits = $result->num_hits();
	my $qlen     = $aln_stats->{$query};
	die "\nERROR: Could not determine query length for: $query" 
	    unless defined $qlen;
	while ( my $hit = $result->next_hit() ) {
	    my $hitid = $hit->name();
	    while ( my $hsp = $hit->next_hsp() ) {
		my $percent_q_coverage = sprintf("%.2f", $hsp->length('query')/$qlen * 100);
		my @ident_pos = $hsp->seq_inds('query','identical');
		my $hspgaps   = $hsp->gaps;
		my $hsplen    = $hsp->length('total');
		my $hstart    = $hsp->start('hit');
		my $hstop     = $hsp->end('hit');
		my $qstart    = $hsp->start('query');
		my $qstop     = $hsp->end('query');
		my $qstring   = $hsp->query_string;
		my $qhsplen   = $hsp->length('query');
		
		if (exists $aln_stats->{$query}) {
		    my $percent_coverage = sprintf("%.2f",$hsplen/$aln_stats->{$query});
		    $positions++ for @ident_pos;
		    my $pid = sprintf("%.2f", $positions/$aln_stats->{$query} * 100);
		    if ( $hsplen >= $match_len && $hsplen >= $aln_stats->{$query} * $match_pcov ) {
			my $qid = $query =~ s/_ltrseqs_muscle-out//r; # non-destructive substitution in v5.14+
			if ($pid >= $match_pid) {
			    $matches++;
			    say $out join "\t", $qid, $aln_stats->{$query}, $num_hits, $hitid, 
			        $percent_q_coverage, $hsplen, $pid, $qstart, $qstop, $hstart, $hstop, $model_type;

			    if ($self->seqfile) {
				my $seqid = join '_', '>'.$qid, $hitid, $hstart, $hstop; 
				## It makes more sense to show the location of the hit
				## Also, this would pave the way for creating a gff of solo-LTRs
				## my $seqid = ">".$query."|".$hitid."_".$hstart."-".$hstop
				say $seq join "\n", $seqid, $qstring;
			    }
			}
		    }
		    $positions = 0;
		}
	    }
	}
    }
    close $out;
    close $seq;
{'fraction_non_alphanumeric': 0.18451413392538557, 'fraction_numerical': 0.008960076121000674, 'mean_word_length': 2.683411214953271, 'pattern_counts': {'":': 0, '<': 9, '<?xml version=': 0, '>': 321, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 4, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '650443', 'n_tokens_mistral': 10284, 'n_tokens_neox': 9579, 'n_words': 2273}