<commit_before>

=pod

=head1 NAME
  
Bio::EnsEMBL::Pipeline::Runnable::GeneWise - aligns protein hits to genomic sequence

=head1 SYNOPSIS

    # Initialise the genewise module  
    my $genewise = new  Bio::EnsEMBL::Pipeline::Runnable::Genewise  ('genomic'  => $genomic,
								       'protein'  => $protein,
								       'memory'   => $memory);


   $genewise->run;
    
   my @genes = $blastwise->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

package Bio::EnsEMBL::Pipeline::Runnable::Genewise;

use strict;  
use vars   qw(@ISA);

# Object preamble - inherits from Bio::Root::RootI

use Bio::Root::RootI;
use Bio::EnsEMBL::Analysis::Programs qw(genewise); 
use Bio::EnsEMBL::SeqFeature;
use Bio::EnsEMBL::FeaturePair;
use Bio::EnsEMBL::Pipeline::RunnableI;

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

sub new {
  my($class,@args) = @_;
  my $self = $class->SUPER::new(@args);
  
  my ($genomic, $protein, $memory,$reverse,$endbias,$genewise) = 
      $self->_rearrange([qw(GENOMIC PROTEIN MEMORY REVERSE ENDBIAS GENEWISE)], @args);
  
#  print $genomic . "\n";

  $genewise = 'genewise' unless defined($genewise);

  $self->genewise($self->find_executable($genewise));

  $self->genomic($genomic) || $self->throw("No genomic sequence entered for blastwise");
  $self->protein($protein) || $self->throw("No protein sequence entered for blastwise");

  $self->reverse($reverse)   if (defined($reverse));             
  $self->endbias($endbias)   if (defined($endbias));             
  $self->memory ($memory)    if (defined($memory));

  return $self;
}

sub genewise {
  my ($self,$arg) = @_;

  if (defined($arg)) {
    $self->{_genewise} = $arg;
  }
  return $self->{_genewise};
}

# RunnableI methods

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

    $self->_wise_init;
    $self->_align_protein;
}

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

#    return $self->eachGene;
    return $self->eachExon;
}

# Now the blastwise stuff

sub set_environment {
  my ($self) = @_;
  
  if (! -d $ENV{WISECONFIGDIR}) {
    $self->throw("No WISECONFIGDIR ["  . $ENV{WISECONFIGDIR} . "]");
  }
}

sub _wise_init {
  my($self,$genefile) = @_;

#  print("Initializing genewise\n");

  $self->set_environment;

}


=pod

=head2 _align_protein

  Title   : _align_protein
  Usage   : $bw->align_protein($pro,$start,$end)
  Function: aligns a protein to a genomic sequence between the genomic coords $start,$end
  Returns : nothing
  Args    : Bio:Seq,int,int

=cut


sub _align_protein {
  my ($self) = @_;
  my $memory   = $self->memory;
  
  my $gwfile   = "/tmp/gw." . $$ . ".out";
  my $genfile  = "/tmp/gw." . $$ . "." . $self->protein->id . ".gen.fa";
  my $protfile = "/tmp/gw." . $$ . "." . $self->protein->id . ".pro.fa";
  
  my $genio  = new Bio::SeqIO(-file   => ">$genfile",
			      '-format' => 'fasta');
  
  $genio->write_seq($self->genomic);
  $genio = undef;
  
  my $proio  = new Bio::SeqIO(-file   => ">$protfile",
			      '-format' => 'fasta');
  
  $proio->write_seq($self->protein);
  $proio = undef;

  my $genewise = $self->genewise;
#  my $command = "/nfs/acari/birney/wise2-bin-snapshot/genewise $protfile $genfile -genesf -kbyte $memory -ext 2 -gap 12 -subs 0.0000001 -quiet";
  my $command = "$genewise $protfile $genfile -genesf -kbyte $memory -ext 2 -gap 12 -subs 0.0000001 -quiet";
    
  if ($self->endbias == 1) {
    $command .= " -init endbias -splice flat ";
  }
  
  if ($self->reverse == 1) {
    $command .= " -trev ";
  }
  
#  print STDERR "Command is $command\n";
  
  open(GW, "$command |") or $self->throw("error piping to genewise: $!\n");

  # for genesf parsing
  my @genesf_exons;
  my $curr_exon;
  my $curr_gene;

  # making assumption of only 1 gene prediction ... this will change once we start using hmms
 GENESF: while (<GW>) {
    chomp;
    my @f = split;

    next unless (defined $f[0] && ($f[0] eq 'Gene' || $f[0] eq 'Exon' || $f[0] eq 'Supporting'));
    
    if($f[0] eq 'Gene'){
      # flag a frameshift - ultimately we will do something clever here but for now ...
      if (/^(Gene\s+\d+)$/){
	if(!defined($curr_gene)){
	  $curr_gene = $1;
	}
	elsif($1 ne $curr_gene){
	  $self->warn("frameshift!\n");
	}
      }
    }

    elsif($f[0] eq 'Exon'){
      #      make a new "exon"
      $curr_exon = new Bio::EnsEMBL::SeqFeature;
      $curr_exon->seqname  ($self->genomic->id);
      $curr_exon->id        ($self->protein->id);
      $curr_exon->source_tag('genewise');
      $curr_exon->primary_tag('similarity');<commit_msg>removed the sourece tag and primary tag from exon creation as now depreciated<commit_after>

=pod

=head1 NAME
  
Bio::EnsEMBL::Pipeline::Runnable::GeneWise - aligns protein hits to genomic sequence

=head1 SYNOPSIS

    # Initialise the genewise module  
    my $genewise = new  Bio::EnsEMBL::Pipeline::Runnable::Genewise  ('genomic'  => $genomic,
								       'protein'  => $protein,
								       'memory'   => $memory);


   $genewise->run;
    
   my @genes = $blastwise->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

package Bio::EnsEMBL::Pipeline::Runnable::Genewise;

use strict;  
use vars   qw(@ISA);

# Object preamble - inherits from Bio::Root::RootI

use Bio::Root::RootI;
use Bio::EnsEMBL::Analysis::Programs qw(genewise); 
use Bio::EnsEMBL::SeqFeature;
use Bio::EnsEMBL::FeaturePair;
use Bio::EnsEMBL::Pipeline::RunnableI;

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

sub new {
  my($class,@args) = @_;
  my $self = $class->SUPER::new(@args);
  
  my ($genomic, $protein, $memory,$reverse,$endbias,$genewise) = 
      $self->_rearrange([qw(GENOMIC PROTEIN MEMORY REVERSE ENDBIAS GENEWISE)], @args);
  
#  print $genomic . "\n";

  $genewise = 'genewise' unless defined($genewise);

  $self->genewise($self->find_executable($genewise));

  $self->genomic($genomic) || $self->throw("No genomic sequence entered for blastwise");
  $self->protein($protein) || $self->throw("No protein sequence entered for blastwise");

  $self->reverse($reverse)   if (defined($reverse));             
  $self->endbias($endbias)   if (defined($endbias));             
  $self->memory ($memory)    if (defined($memory));

  return $self;
}

sub genewise {
  my ($self,$arg) = @_;

  if (defined($arg)) {
    $self->{_genewise} = $arg;
  }
  return $self->{_genewise};
}

# RunnableI methods

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

    $self->_wise_init;
    $self->_align_protein;
}

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

#    return $self->eachGene;
    return $self->eachExon;
}

# Now the blastwise stuff

sub set_environment {
  my ($self) = @_;
  
  if (! -d $ENV{WISECONFIGDIR}) {
    $self->throw("No WISECONFIGDIR ["  . $ENV{WISECONFIGDIR} . "]");
  }
}

sub _wise_init {
  my($self,$genefile) = @_;

#  print("Initializing genewise\n");

  $self->set_environment;

}


=pod

=head2 _align_protein

  Title   : _align_protein
  Usage   : $bw->align_protein($pro,$start,$end)
  Function: aligns a protein to a genomic sequence between the genomic coords $start,$end
  Returns : nothing
  Args    : Bio:Seq,int,int

=cut


sub _align_protein {
  my ($self) = @_;
  my $memory   = $self->memory;
  
  my $gwfile   = "/tmp/gw." . $$ . ".out";
  my $genfile  = "/tmp/gw." . $$ . "." . $self->protein->id . ".gen.fa";
  my $protfile = "/tmp/gw." . $$ . "." . $self->protein->id . ".pro.fa";
  
  my $genio  = new Bio::SeqIO(-file   => ">$genfile",
			      '-format' => 'fasta');
  
  $genio->write_seq($self->genomic);
  $genio = undef;
  
  my $proio  = new Bio::SeqIO(-file   => ">$protfile",
			      '-format' => 'fasta');
  
  $proio->write_seq($self->protein);
  $proio = undef;

  my $genewise = $self->genewise;
#  my $command = "/nfs/acari/birney/wise2-bin-snapshot/genewise $protfile $genfile -genesf -kbyte $memory -ext 2 -gap 12 -subs 0.0000001 -quiet";
  my $command = "$genewise $protfile $genfile -genesf -kbyte $memory -ext 2 -gap 12 -subs 0.0000001 -quiet";
    
  if ($self->endbias == 1) {
    $command .= " -init endbias -splice flat ";
  }
  
  if ($self->reverse == 1) {
    $command .= " -trev ";
  }
  
#  print STDERR "Command is $command\n";
  
  open(GW, "$command |") or $self->throw("error piping to genewise: $!\n");

  # for genesf parsing
  my @genesf_exons;
  my $curr_exon;
  my $curr_gene;

  # making assumption of only 1 gene prediction ... this will change once we start using hmms
 GENESF: while (<GW>) {
    chomp;
    my @f = split;

    next unless (defined $f[0] && ($f[0] eq 'Gene' || $f[0] eq 'Exon' || $f[0] eq 'Supporting'));
    
    if($f[0] eq 'Gene'){
      # flag a frameshift - ultimately we will do something clever here but for now ...
      if (/^(Gene\s+\d+)$/){
	if(!defined($curr_gene)){
	  $curr_gene = $1;
	}
	elsif($1 ne $curr_gene){
	  $self->warn("frameshift!\n");
	}
      }
    }

    elsif($f[0] eq 'Exon'){
      #      make a new "exon"
      $curr_exon = new Bio::EnsEMBL::SeqFeature;
      $curr_exon->seqname  ($self->genomic->id);
      $curr_exon->id        ($self->protein->id);
      
{'fraction_non_alphanumeric': 0.1763180639585134, 'fraction_numerical': 0.00864304235090752, 'mean_word_length': 3.1436884512085945, 'pattern_counts': {'":': 0, '<': 5, '<?xml version=': 0, '>': 109, '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': '5824578', 'n_tokens_mistral': 3709, 'n_tokens_neox': 3431, 'n_words': 793}