<commit_before>#
#
# Cared for by EnsEMBL  <<EMAIL>>
#
# Copyright GRL & EBI
#
# You may distribute this module under the same terms as perl itself
#
# POD documentation - main docs before the code

=pod 

=head1 NAME

Bio::EnsEMBL::Pipeline::Runnable::ExonerateMiniEst2Genome

=head1 SYNOPSIS

    my $obj = Bio::EnsEMBL::Pipeline::Runnable::ExonerateMiniEst2Genome->new('-genomic'    => $genseq,
									 '-seqfetcher' => $seqfetcher
									 '-queryseq'   => $queryseq);

    $obj->run

    my @features = $obj->output;


=head1 DESCRIPTION

=head1 CONTACT

Describe contact details here

=head1 APPENDIX

The rest of the documentation details each of the object methods. 
Internal methods are usually preceded with a _

=cut

# Let the code begin...

package Bio::EnsEMBL::Pipeline::Runnable::ExonerateMiniEst2Genome;

use vars qw(@ISA);
use strict;

# Object preamble - inherits from Bio::Root::RootI;
use Bio::EnsEMBL::Pipeline::Runnable::MiniEst2Genome;
use Bio::EnsEMBL::Pipeline::Runnable::Exonerate;
use Bio::EnsEMBL::Pipeline::RunnableI;
#use Bio::EnsEMBL::Analysis::MSPcrunch;
use Bio::PrimarySeqI;
use Bio::SeqIO;
use Bio::DB::RandomAccessI;
use Data::Dumper;

@ISA = qw(Bio::EnsEMBL::Pipeline::RunnableI);

=head2 new

    Title   :   new
    Usage   :   $self->new(-GENOMIC       => $genomicseq,
			   -QUERYSEQ      => $queryseq,
                           -SEQFETCHER    => $sf);
                           
    Function:   creates a 
                Bio::EnsEMBL::Pipeline::Runnable::ExonerateMiniEst2Genome object
    Returns :   A Bio::EnsEMBL::Pipeline::Runnable::ExonerateMiniEst2Genome object
    Args    :   -genomic:    Bio::PrimarySeqI object (genomic sequence)
  -queryseq:   Either path to file containing query seqs or reference to aray of Bio::Seq
                -seqfetcher  Bio::DB::RandomAccessI object
=cut

sub new {
    my ($class,@args) = @_;
    my $self = $class->SUPER::new(@args);

    $self->{'_idlist'} = []; #create key to an array of feature pairs
    
    my( $genomic, $queryseq, $seqfetcher, $id, $length) = $self->_rearrange([qw(GENOMIC
									       QUERYSEQ
									       SEQFETCHER
									       ID_THRESHOLD
									       LENGTH_THRESHOLD)],
									   @args);
	 
    $self->throw("No genomic sequence input")           
      unless defined($genomic);
    $self->throw("[$genomic] is not a Bio::PrimarySeqI") 
      unless $genomic->isa("Bio::PrimarySeqI");
    $self->genomic_sequence($genomic) if defined($genomic);

    $self->throw("No queryseq specified") 
      unless defined($queryseq);
    $self->queryseq($queryseq) if defined($queryseq);

    $self->throw("No seqfetcher provided")           
      unless defined($seqfetcher);
    $self->throw("[$seqfetcher] is not a Bio::DB::RandomAccessI") 
      unless $seqfetcher->isa("Bio::DB::RandomAccessI");
    $self->seqfetcher($seqfetcher) if defined($seqfetcher);

    if(defined $id){
      $self->{'_id_threshold'} = $id;
    }
    else{
      $self->{'_id_threshold'} = 50;
    }
    
    if(defined $length){
      $self->{'_length_threshold'} = $length;
    }
    else{
      $self->{'_length_threshold'} = 50;
    } 

    return $self; 
}

=head2 genomic_sequence

    Title   :   genomic_sequence
    Usage   :   $self->genomic_sequence($seq)
    Function:   Get/set method for genomic sequence
    Returns :   Bio::Seq object
    Args    :   Bio::Seq object

=cut

sub genomic_sequence {
    my( $self, $value ) = @_;    
    if ($value) {
        #need to check if passed sequence is Bio::Seq object
        $value->isa("Bio::PrimarySeqI") || $self->throw("Input isn't a Bio::PrimarySeqI");
        $self->{'_genomic_sequence'} = $value;
    }
    return $self->{'_genomic_sequence'};
}

=head2 seqfetcher

    Title   :   seqfetcher
    Usage   :   $self->seqfetcher($seqfetcher)
    Function:   Get/set method for SeqFetcher
    Returns :   Bio::DB:RandomAccessI
    Args    :   Bio::DB:RandomAccessI

=cut

sub seqfetcher {
    my( $self, $value ) = @_;    
    if ($value) {
      #need to check if passed sequence is Bio::EnsEMBL::Pipeline::SeqFetcherI object
      $self->throw("Input isn't a Bio::DB::RandomAccessI") unless $value->isa("Bio::DB::RandomAccessI");
      $self->{'_seqfetcher'} = $value;
    }
    return $self->{'_seqfetcher'};
}

=head2 queryseq

    Title   :   queryseq
    Usage   :   $self->queryseq($seq)
    Function:   Get/set method for queryseq
    Returns :   
    Args    :   name of a file containing query seq(s), OR reference to an array of Bio::Seq

=cut

sub queryseq {
    my( $self, $queryseq ) = @_;   
    if ($queryseq) { 
      if (ref($queryseq) eq 'ARRAY') {

	my $time = time; chomp($time);
	my $estfile = "/tmp/estfile_.$$.$time.fn";
	$self->queryfilename($estfile);

	foreach my $est(@$queryseq) {
	  $est->isa("Bio::PrimarySeqI") || $self->throw("Input isn't a Bio::PrimarySeqI");
	}

	$self->{'_query_sequences'} = $queryseq;
      }
      else {
	# it's a filename - check the file exists
	$self->throw("[$queryseq] : file does not exist\n") unless -e $queryseq;
	$self->queryfilename($queryseq);
	$self->{'_query_sequences'} = $queryseq;
    }
  }
  
  #NB ref to an array of Bio::Seq
  return $self->{'_query_sequences'};

  }

=head2 queryfilename

    Title   :   queryfilename
    Usage   :   $self->queryfilename($filename)
    Function:   Get/set method for queryfilename
    Returns :   
    Args    :   

=cut

sub queryfilename {
  my ($self, $queryfilename) = @_;
  $self->{'_queryfilename'} = $queryfilename if ($queryfilename);
  return $self->{'_queryfilename'};
}

=head2 get_all_FeatureIds

  Title   : get_all_FeatureIds
  Usage   : my @ids = get_all_FeatureIds
  Function: Returns an array of all distinct feature hids 
  Returns : @string
  Args    : none

=cut

sub get_Ids {
    my ($self) = @_;

    if (!defined($self->{'_idlist'})) {
	$self->{'_idlist'} = [];
    }
    return @{$self->{'_idlist'}};
}


=head2 parse_Header

  Title   : parse_Header
  Usage   : my $newid = $self->parse_Header($id);
  Function: Parses different sequence headers
  Returns : string
  Args    : none

=cut

sub parse_Header {
    my ($self,$id) = @_;

    if (!defined($id)) {
	$self->throw("No id input to parse_Header");
    }

    my $newid = $id;

    if ($id =~ /^(.*)\|(.*)\|(.*)/) {
	$newid = $2;
	$newid =~ s/(.*)\..*/$1/;
	
    } elsif ($id =~ /^..\:(.*)/) {
	$newid = $1;
    }
    $newid =~ s/ //g;
    return $newid;
}


=head2 run

  Title   : run
  Usage   : $self->run()
  Function: Runs exonerate vs input seqs, and runs a MiniEst2Genome runnable for each appropriate set of blast hits
  Returns : none
  Args    : 

=cut

sub run {
    my ($self) = @_;

    # filter ESTs using exonerate
    my @exonerate_res = $self->run_exonerate();

    print STDERR "**got " . scalar(@exonerate_res) . " hits\n";
    
    my %exonerate_ests;
    my @feat;
    foreach my $res(@exonerate_res ) {
      my $seqname = $res->hseqname;       #gb|AA429061.1|AA429061
      $seqname =~ s/\S+\|(\S+)\|\S+/$1/;
      $res->hseqname($seqname);

      # may move this out of here.

      # score cutoff 500 for exonerate ... take all features for a sequence as long as one of them gets over this threshold
      if($res->score > 500 || defined $exonerate_ests{$seqname}) {
	push(@{$exonerate_ests{$seqname}}, $res);
	push (@feat, $res);
      }
    }

    # filter features
    my %filtered_ests;
    my $filter = new Bio::EnsEMBL::Pipeline::Runnable::FeatureFilter( '-coverage' => 10,
								      '-minscore' => 500);
    my @filteredfeats = $filter->run(@feat);

    foreach my $f(@filteredfeats){
      push(@{$filtered_ests{$f->hseqname}}, $f);
    }

  ID:    
    foreach my $id(keys %filtered_ests) {
      my @features = @{$filtered_ests{$id}};
 
      # only use ESTs that have >1 exonerate hit to cut down on how many e2gs we run.
#      next ID unless scalar(@features) > 1; # ??? too strict?

      # make MiniEst2Genome runnables
      
      my $e2g = new Bio::EnsEMBL::Pipeline::Runnable::MiniEst2Genome('-genomic'  => $self->genomic_sequence,
								     '-features' => \@features,
								     '-seqfetcher' => $self->seqfetcher);

      # run runnable
      $e2g->run;
      
      # sort out output
      my @f = $e2g->output;
      
      push(@{$self->{'_output'}},@f);
    }
}

sub run_exonerate {
  my ($self) = @_;
  my @res;

  my $estseq  = $self->queryseq;
  my $estfile = $self->queryfilename;

  # do we need to write out the est sequences?
  if(ref($estseq) eq 'ARRAY'){
    eval{
      if (-e $estfile) { $self->throw("alreayd using $estfile\n"); }
      my $estOutput = Bio::SeqIO->new(-file => ">$estfile" , '-format' => 'Fasta')
	or $self->throw("Can't create new Bio::SeqIO from $estfile '$' : $!");
      
      foreach my $eseq(@$estseq) {
	$estOutput->write_seq($eseq);<commit_msg>getting rid of an undef warning.<commit_after>#
#
# Cared for by EnsEMBL  <<EMAIL>>
#
# Copyright GRL & EBI
#
# You may distribute this module under the same terms as perl itself
#
# POD documentation - main docs before the code

=pod 

=head1 NAME

Bio::EnsEMBL::Pipeline::Runnable::ExonerateMiniEst2Genome

=head1 SYNOPSIS

    my $obj = Bio::EnsEMBL::Pipeline::Runnable::ExonerateMiniEst2Genome->new('-genomic'    => $genseq,
									 '-seqfetcher' => $seqfetcher
									 '-queryseq'   => $queryseq);

    $obj->run

    my @features = $obj->output;


=head1 DESCRIPTION

=head1 CONTACT

Describe contact details here

=head1 APPENDIX

The rest of the documentation details each of the object methods. 
Internal methods are usually preceded with a _

=cut

# Let the code begin...

package Bio::EnsEMBL::Pipeline::Runnable::ExonerateMiniEst2Genome;

use vars qw(@ISA);
use strict;

# Object preamble - inherits from Bio::Root::RootI;
use Bio::EnsEMBL::Pipeline::Runnable::MiniEst2Genome;
use Bio::EnsEMBL::Pipeline::Runnable::Exonerate;
use Bio::EnsEMBL::Pipeline::RunnableI;
#use Bio::EnsEMBL::Analysis::MSPcrunch;
use Bio::PrimarySeqI;
use Bio::SeqIO;
use Bio::DB::RandomAccessI;
use Data::Dumper;

@ISA = qw(Bio::EnsEMBL::Pipeline::RunnableI);

=head2 new

    Title   :   new
    Usage   :   $self->new(-GENOMIC       => $genomicseq,
			   -QUERYSEQ      => $queryseq,
                           -SEQFETCHER    => $sf);
                           
    Function:   creates a 
                Bio::EnsEMBL::Pipeline::Runnable::ExonerateMiniEst2Genome object
    Returns :   A Bio::EnsEMBL::Pipeline::Runnable::ExonerateMiniEst2Genome object
    Args    :   -genomic:    Bio::PrimarySeqI object (genomic sequence)
  -queryseq:   Either path to file containing query seqs or reference to aray of Bio::Seq
                -seqfetcher  Bio::DB::RandomAccessI object
=cut

sub new {
    my ($class,@args) = @_;
    my $self = $class->SUPER::new(@args);

    $self->{'_idlist'} = []; #create key to an array of feature pairs
    
    my( $genomic, $queryseq, $seqfetcher, $id, $length) = $self->_rearrange([qw(GENOMIC
									       QUERYSEQ
									       SEQFETCHER
									       ID_THRESHOLD
									       LENGTH_THRESHOLD)],
									   @args);
	 
    $self->throw("No genomic sequence input")           
      unless defined($genomic);
    $self->throw("[$genomic] is not a Bio::PrimarySeqI") 
      unless $genomic->isa("Bio::PrimarySeqI");
    $self->genomic_sequence($genomic) if defined($genomic);

    $self->throw("No queryseq specified") 
      unless defined($queryseq);
    $self->queryseq($queryseq) if defined($queryseq);

    $self->throw("No seqfetcher provided")           
      unless defined($seqfetcher);
    $self->throw("[$seqfetcher] is not a Bio::DB::RandomAccessI") 
      unless $seqfetcher->isa("Bio::DB::RandomAccessI");
    $self->seqfetcher($seqfetcher) if defined($seqfetcher);

    if(defined $id){
      $self->{'_id_threshold'} = $id;
    }
    else{
      $self->{'_id_threshold'} = 50;
    }
    
    if(defined $length){
      $self->{'_length_threshold'} = $length;
    }
    else{
      $self->{'_length_threshold'} = 50;
    } 

    return $self; 
}

=head2 genomic_sequence

    Title   :   genomic_sequence
    Usage   :   $self->genomic_sequence($seq)
    Function:   Get/set method for genomic sequence
    Returns :   Bio::Seq object
    Args    :   Bio::Seq object

=cut

sub genomic_sequence {
    my( $self, $value ) = @_;    
    if ($value) {
        #need to check if passed sequence is Bio::Seq object
        $value->isa("Bio::PrimarySeqI") || $self->throw("Input isn't a Bio::PrimarySeqI");
        $self->{'_genomic_sequence'} = $value;
    }
    return $self->{'_genomic_sequence'};
}

=head2 seqfetcher

    Title   :   seqfetcher
    Usage   :   $self->seqfetcher($seqfetcher)
    Function:   Get/set method for SeqFetcher
    Returns :   Bio::DB:RandomAccessI
    Args    :   Bio::DB:RandomAccessI

=cut

sub seqfetcher {
    my( $self, $value ) = @_;    
    if ($value) {
      #need to check if passed sequence is Bio::EnsEMBL::Pipeline::SeqFetcherI object
      $self->throw("Input isn't a Bio::DB::RandomAccessI") unless $value->isa("Bio::DB::RandomAccessI");
      $self->{'_seqfetcher'} = $value;
    }
    return $self->{'_seqfetcher'};
}

=head2 queryseq

    Title   :   queryseq
    Usage   :   $self->queryseq($seq)
    Function:   Get/set method for queryseq
    Returns :   
    Args    :   name of a file containing query seq(s), OR reference to an array of Bio::Seq

=cut

sub queryseq {
    my( $self, $queryseq ) = @_;   
    if ($queryseq) { 
      if (ref($queryseq) eq 'ARRAY') {

	my $time = time; chomp($time);
	my $estfile = "/tmp/estfile_.$$.$time.fn";
	$self->queryfilename($estfile);

	foreach my $est(@$queryseq) {
	  $est->isa("Bio::PrimarySeqI") || $self->throw("Input isn't a Bio::PrimarySeqI");
	}

	$self->{'_query_sequences'} = $queryseq;
      }
      else {
	# it's a filename - check the file exists
	$self->throw("[$queryseq] : file does not exist\n") unless -e $queryseq;
	$self->queryfilename($queryseq);
	$self->{'_query_sequences'} = $queryseq;
    }
  }
  
  #NB ref to an array of Bio::Seq
  return $self->{'_query_sequences'};

  }

=head2 queryfilename

    Title   :   queryfilename
    Usage   :   $self->queryfilename($filename)
    Function:   Get/set method for queryfilename
    Returns :   
    Args    :   

=cut

sub queryfilename {
  my ($self, $queryfilename) = @_;
  $self->{'_queryfilename'} = $queryfilename if ($queryfilename);
  return $self->{'_queryfilename'};
}

=head2 get_all_FeatureIds

  Title   : get_all_FeatureIds
  Usage   : my @ids = get_all_FeatureIds
  Function: Returns an array of all distinct feature hids 
  Returns : @string
  Args    : none

=cut

sub get_Ids {
    my ($self) = @_;

    if (!defined($self->{'_idlist'})) {
	$self->{'_idlist'} = [];
    }
    return @{$self->{'_idlist'}};
}


=head2 parse_Header

  Title   : parse_Header
  Usage   : my $newid = $self->parse_Header($id);
  Function: Parses different sequence headers
  Returns : string
  Args    : none

=cut

sub parse_Header {
    my ($self,$id) = @_;

    if (!defined($id)) {
	$self->throw("No id input to parse_Header");
    }

    my $newid = $id;

    if ($id =~ /^(.*)\|(.*)\|(.*)/) {
	$newid = $2;
	$newid =~ s/(.*)\..*/$1/;
	
    } elsif ($id =~ /^..\:(.*)/) {
	$newid = $1;
    }
    $newid =~ s/ //g;
    return $newid;
}


=head2 run

  Title   : run
  Usage   : $self->run()
  Function: Runs exonerate vs input seqs, and runs a MiniEst2Genome runnable for each appropriate set of blast hits
  Returns : none
  Args    : 

=cut

sub run {
    my ($self) = @_;

    # filter ESTs using exonerate
    my @exonerate_res = $self->run_exonerate();

    print STDERR "**got " . scalar(@exonerate_res) . " hits\n";
    
    my %exonerate_ests;
    my @feat;
    foreach my $res(@exonerate_res ) {
      my $seqname = $res->hseqname;       #gb|AA429061.1|AA429061
      $seqname =~ s/\S+\|(\S+)\|\S+/$1/;
      $res->hseqname($seqname);

      # may move this out of here.

      # score cutoff 500 for exonerate ... take all features for a sequence as long as one of them gets over this threshold
      if($res->score > 500 || defined $exonerate_ests{$seqname}) {
	push(@{$exonerate_ests{$seqname}}, $res);
	push (@feat, $res);
      }
    }

    # filter features
    my %filtered_ests;
    my $filter = new Bio::EnsEMBL::Pipeline::Runnable::FeatureFilter( '-coverage' => 10,
								      '-minscore' => 500);
    my @filteredfeats = $filter->run(@feat);

    foreach my $f(@filteredfeats){
      push(@{$filtered_ests{$f->hseqname}}, $f);
    }

  ID:    
    foreach my $id(keys %filtered_ests) {
      my @features = @{$filtered_ests{$id}};
 
      # only use ESTs that have >1 exonerate hit to cut down on how many e2gs we run.
#      next ID unless scalar(@features) > 1; # ??? too strict?

      # make MiniEst2Genome runnables
      
      my $e2g = new Bio::EnsEMBL::Pipeline::Runnable::MiniEst2Genome('-genomic'  => $self->genomic_sequence,
								     '-features' => \@features,
								     '-seqfetcher' => $self->seqfetcher);

      # run runnable
      $e2g->run;
      
      # sort out output
      my @f = $e2g->output;
      
      push(@{$self->{'_output'}},@f);
    }
}

sub run_exonerate {
  my ($self) = @_;

  my $estseq  = $self->queryseq;
  my $estfile = $self->queryfilename;

  # do we need to write out the est sequences?
  if(ref($estseq) eq 'ARRAY'){
    eval{
      if (-e $estfile) { $self->throw("alreayd using $estfile\n"); }
      my $estOutput = Bio::SeqIO->new(-file => ">$estfile" , '-format' => 'Fasta')
	or $self->throw("Can't create new Bio::SeqIO from $estfile '$' : $!");
      
      foreach my $eseq(@$estseq) {
	$estOutput->write_seq($eseq);
{'fraction_non_alphanumeric': 0.1544876886417792, 'fraction_numerical': 0.006808124361738341, 'mean_word_length': 2.8261341436943783, 'pattern_counts': {'":': 0, '<': 7, '<?xml version=': 0, '>': 173, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '2275858', 'n_tokens_mistral': 6598, 'n_tokens_neox': 6047, 'n_words': 1454}