<commit_before>=head1 NAME

ExonPair

=head1 SYNOPSIS

my $exon_pair = Bio::EnsEMBL::Pipeline::GeneComparison::ExonPair->new();

my $score = $exon_pair->blast_Exons( $exon1, $exon2 );

$gene_pair->compare_isoforms();

=head1 DESCRIPTION

Class to operate on potential orthologous exons.

=head1 CONTACT

<EMAIL>

=cut

# Let the code begin ...

package Bio::EnsEMBL::Pipeline::GeneComparison::ExonPair;

use vars qw(@ISA);
use strict;

use Bio::EnsEMBL::Root;
use Bio::EnsEMBL::Pipeline::Runnable::Blast;
use Bio::SeqIO;

@ISA = qw(Bio::EnsEMBL::Root);

=head1 METHODS

=cut

#########################################################################

sub new {
    my ($class,@args) = @_;
    my $self = $class->SUPER::new(@args);
  
    my( $gap_penalty ) = $self->_rearrange([qw(
					       GAP_PENALTY
					       )], 
					   @args);
 
    if ( $gap_penalty ){
	$self->gap_penalty( $gap_penalty );
    }
    else{
	$self->gap_penalty( -100 );
    }
    
    return $self;
}

############################################################


sub blast_Exons{
  my ($self,$exon1, $exon2) =@_;
  
  my $id1;
  if ( $exon1->dbID ){
    $id1 = $exon1->stable_id || $exon2->dbID;
  }
  else{
    $id1 = $exon1;
  }
  
  my $id2;
  if ( $exon2->dbID ){
    $id2 = $exon2->stable_id || $exon2->dbID;
  }
  else{
    $id2 = $exon2;
  }
  
  my $seq1    = $exon1->seq;
  my $length1 = $seq1->length;
  unless ( $seq1->display_id ){
    $seq1->display_id($id1);
  }

  my $seq2    = $exon2->seq;
  my $length2 = $seq2->length;
  unless ( $seq2->display_id ){
    $seq2->display_id($id2);
  }
  
  my $min_length = $length1;
  if ( $length2 < $length1 ){
    $min_length = $length2;
  }
  my $word = 5;
  if ( 3*$word >= $min_length ){
    $word = int($min_length/3) - 1;
  }
  if ( $word < 1 ){
    return 0;
  }

  #print STDERR "word: $word\n";
  #print STDERR "comparing $id1: ".$seq1->seq."\n";
  #print STDERR "      and $id2: ".$seq2->seq."\n";


  ############################################################
  # create database
  my $file = 'seq_'.$$.'.fa';
  my $database = "/tmp/".$file;
  open( DB_SEQ,">$database") || die("Could not open $database $!");
  
  my $seqout = Bio::SeqIO->new('-format' => 'Fasta',
			       '-fh'     => \*DB_SEQ);
  
  $seqout->write_seq($seq2);
  close( DB_SEQ );
  
  system("pressdb $database > /dev/null 2>&1");
  
  ############################################################
  # Ian's parameters:
  #my $options = "W=5 M=1 N=-1 Q=3 R=3";
  
  my $options = "W=5";

  # tblastx options:
  #my $options = 'altscore="* any na" altscore="any * na" S2=13';
  #$options .= " V=200 B=200 ";
  #$options .= " -nogap ";
  #$options .= " W=$word ";
  #print STDERR "options: $options\n";
  #my $options = 'V=200 B=200 altscore="* any na" altscore="any * na" W=4 E=0.01 E2=0.01 -nogap';
  #my $options = 'V=200 B=200 W=9 E=0.01 E2=0.01';
  my $blast =  
    Bio::EnsEMBL::Pipeline::Runnable::Blast->new ('-query'          => $seq1,
						  #'-program'        => 'wutblastx',
						  '-program'        => 'wublastn',
						  '-database'       => $database,
						  '-threshold_type' => "PVALUE",
						  #'-threshold'      => 1e-10,
						  '-options'        => $options,
						 );
  
  
  $blast->add_regex($file,'(\S+)');
  $blast->run();
  
  unlink( $database );
  
  my @featurepairs = $blast->output();
  
  if ( @featurepairs ){
      #my @pos_strand = grep { $_->strand == 1} @featurepairs;  
      #my @neg_strand = grep { $_->strand == -1} @featurepairs;  
      #foreach my $fp (sort{ $a->hstart <=> $b->hstart} @pos_strand) {
      #  print $fp->gffstring . "\n";<commit_msg>changed to blastn to compare Exons<commit_after>=head1 NAME

ExonPair

=head1 SYNOPSIS

my $exon_pair = Bio::EnsEMBL::Pipeline::GeneComparison::ExonPair->new();

my $score = $exon_pair->blast_Exons( $exon1, $exon2 );

$gene_pair->compare_isoforms();

=head1 DESCRIPTION

Class to operate on potential orthologous exons.

=head1 CONTACT

<EMAIL>

=cut

# Let the code begin ...

package Bio::EnsEMBL::Pipeline::GeneComparison::ExonPair;

use vars qw(@ISA);
use strict;

use Bio::EnsEMBL::Root;
use Bio::EnsEMBL::Pipeline::Runnable::Blast;
use Bio::SeqIO;

@ISA = qw(Bio::EnsEMBL::Root);

=head1 METHODS

=cut

#########################################################################

sub new {
    my ($class,@args) = @_;
    my $self = $class->SUPER::new(@args);
  
    my( $gap_penalty ) = $self->_rearrange([qw(
					       GAP_PENALTY
					       )], 
					   @args);
 
    if ( $gap_penalty ){
	$self->gap_penalty( $gap_penalty );
    }
    else{
	$self->gap_penalty( -100 );
    }
    
    return $self;
}

############################################################


sub blast_Exons{
  my ($self,$exon1, $exon2) =@_;
  
  my $id1;
  if ( $exon1->dbID ){
    $id1 = $exon1->stable_id || $exon2->dbID;
  }
  else{
    $id1 = $exon1;
  }
  
  my $id2;
  if ( $exon2->dbID ){
    $id2 = $exon2->stable_id || $exon2->dbID;
  }
  else{
    $id2 = $exon2;
  }
  
  my $seq1    = $exon1->seq;
  my $length1 = $seq1->length;
  unless ( $seq1->display_id ){
    $seq1->display_id($id1);
  }

  my $seq2    = $exon2->seq;
  my $length2 = $seq2->length;
  unless ( $seq2->display_id ){
    $seq2->display_id($id2);
  }
  
  my $min_length = $length1;
  if ( $length2 < $length1 ){
    $min_length = $length2;
  }
  my $word = 5;
  #if ( $word*3 >= $min_length ){
  #  $word = int($min_length/3) - 1;
  #}
  if ( $word >= $min_length ){
      $word = $min_length - 1;
  }
  if ( $word < 1 ){
    return 0;
  }

  #print STDERR "word: $word\n";
  #print STDERR "comparing $id1: ".$seq1->seq."\n";
  #print STDERR "      and $id2: ".$seq2->seq."\n";


  ############################################################
  # create database
  my $file = 'seq_'.$$.'.fa';
  my $database = "/tmp/".$file;
  open( DB_SEQ,">$database") || die("Could not open $database $!");
  
  my $seqout = Bio::SeqIO->new('-format' => 'Fasta',
			       '-fh'     => \*DB_SEQ);
  
  $seqout->write_seq($seq2);
  close( DB_SEQ );
  
  system("pressdb $database > /dev/null 2>&1");
  
  ############################################################
  # Ian's parameters:
  my $options = "W=$word M=1 N=-1 Q=3 R=3"; 
  #my $options = "W=5";
  
  # tblastx options:
  #my $options = 'altscore="* any na" altscore="any * na" S2=12';
  #$options .= " V=200 B=200 ";
  #$options .= " -nogap ";
  #$options .= " W=$word ";
  
  #print STDERR "options: $options\n";
  
 #my $options = 'V=200 B=200 altscore="* any na" altscore="any * na" W=4 E=0.01 E2=0.01 -nogap';
  #my $options = 'V=200 B=200 W=9 E=0.01 E2=0.01';
  my $blast =  
    Bio::EnsEMBL::Pipeline::Runnable::Blast->new ('-query'          => $seq1,
						  #'-program'        => 'wutblastx',
						  '-program'        => 'wublastn',
						  '-database'       => $database,
						  '-threshold_type' => "PVALUE",
						  #'-threshold'      => 1e-10,
						  '-options'        => $options,
						 );
  
  
  $blast->add_regex($file,'(\S+)');
  $blast->run();
  
  unlink( $database );
  
  my @featurepairs = $blast->output();
  
  if ( @featurepairs ){
      #my @pos_strand = grep { $_->strand == 1} @featurepairs;  
      #my @neg_strand = grep { $_->strand == -1} @featurepairs;  
      #foreach my $fp (sort{ $a->hstart <=> $b->hstart} @pos_strand) {
      #  print $fp->gffstring . "\n";
{'fraction_non_alphanumeric': 0.2588347659980898, 'fraction_numerical': 0.02892618365397735, 'mean_word_length': 2.878306878306878, 'pattern_counts': {'":': 0, '<': 11, '<?xml version=': 0, '>': 102, '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': '976495', 'n_tokens_mistral': 3198, 'n_tokens_neox': 2885, 'n_words': 564}