<commit_before>=pod

=head1 NAME Bio::EnsEMBL::Analysis::Runnable::TranscriptConsensus

=head1 SYNOPSIS

  my $runnable = Bio::EnsEMBL::Analysis::Runnable::TranscriptConsensus->new
     (
     -query            => $query,
     -analysis        => $analysis,
     -all_genes       => \%biotypes_to_genes , # ref to $hash{biotype_of_gene} = @all_genes_of_this_biotype
     -evidence_sets   => $evidence_sets ,
     -dnadb           => $db ,
     )

=head1 DESCRIPTION

TranscriptConsensus is an extension of TranscriptCoalescer that combines protein coding and est
transcripts to identify the best supported transcript models.
The initial gene sets are clustered, then collapsed into a non-redundant set of exons
and introns which are assigned scores according to the amount of supporting evidence they have.
The similarity genes have UTR added using the est transcripts and the resulting models are assigned
a score by summing the individual exon and intron scores.
The transcripts are sorted by score and the highest scoring models are made into
gene objets and written to the TranscriptCoalescer database.

=head1 CONTACT

Post questions to the EnsEMBL developer list: <<EMAIL>>

=cut


package Bio::EnsEMBL::Analysis::Runnable::TranscriptConsensus;

use strict;
use warnings; 

use Bio::EnsEMBL::Utils::Exception qw(throw warning info);
use Bio::EnsEMBL::Utils::Argument qw( rearrange );
use Bio::EnsEMBL::Analysis::Runnable::TranscriptCoalescer;
use Bio::EnsEMBL::Analysis::Tools::Algorithms::GeneCluster;
use Bio::EnsEMBL::Analysis::Tools::Algorithms::ExonCluster;
use Bio::EnsEMBL::Analysis::Tools::Algorithms::ClusterUtils;
use Bio::EnsEMBL::Analysis::Tools::GeneBuildUtils::TranscriptUtils qw (count_non_canonical_splice_sites) ;
use Bio::EnsEMBL::Analysis::Tools::GeneBuildUtils::EvidenceUtils qw ( clone_Evidence );
use Bio::EnsEMBL::Analysis::Tools::GeneBuildUtils::TranslationUtils qw( clone_Translation );
use Bio::EnsEMBL::Analysis::Config::GeneBuild::TranscriptCoalescer;
use Bio::EnsEMBL::Analysis::Config::GeneBuild::TranscriptConsensus;
use Bio::EnsEMBL::Analysis::Tools::GeneBuildUtils::CollapsedCluster;

use vars qw(@ISA);

@ISA = qw(Bio::EnsEMBL::Analysis::Runnable::TranscriptCoalescer);

=head2 run

   Name       : run
   Args       : none
   Function   : Pre filters genes and clusters them, then collapses each cluster to make non-redundant
              : sets of exons and introns, builds transcripts and adds UTR and assigns them scores,
              : makes genes out of the highest scoring transcripts
   Exceptions : none
   Returnval  : none

=cut

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

  my @allgenes = @{ $self->get_genes_by_evidence_set('simgw') }  ;
  push @allgenes , @{ $self->get_genes_by_evidence_set('est') }  ;

  if ($FILTER_SINGLETONS or   $FILTER_NON_CONSENSUS){
    @allgenes = @{$self->filter_genes(\@allgenes)};
  }
  
  my ($clusters, $non_clusters) = cluster_Genes(\@allgenes, $self->get_all_evidence_sets ) ;

  # create a hash of genes by strand to use when looking for overlapping genes
  # useful to look at ALL genes not just those in the cluster to prevent cluster joining
  foreach my $gene (@allgenes){
    push @{$genes_by_strand->{$gene->strand}}, $gene;
  }

  push @$clusters, @$non_clusters if (scalar(@$non_clusters) > 0 ) ;

  $count = 0;
  foreach my $cluster (@$clusters){
    my $simgw;
    # cluster has to contain at least one similarity gene to be worth continuing with
    next unless $cluster->get_Genes_by_Set('simgw');

    print "\nCluster $count\n" if $VERBOSE;
    print $cluster->start." ".$cluster->end." ".$cluster->strand."\n" if $VERBOSE;
    $count++;
    # collapse the cluster down to make a non redundat set of introns and exons with scores
    my $collapsed_cluster = $self->collapse_cluster($cluster,$genes_by_strand);

    # add UTR to the similarity transcripts and make scores for the transcripts
    my $transcripts = $self->make_transcripts($cluster,$collapsed_cluster);

    # sort the scored transcripts by score and make gene objects out of the best ones
    $self->make_genes($transcripts,$cluster,$collapsed_cluster);
  }
}

=head2 collapse_cluster

   Name       : collapse_cluster
   Arg[0]     : Bio::EnsEMBL::Analysis::Tools::Algorithms::GeneCluster
   Arg[1]     : Array ref of Bio::EnsEMBL::Gene objects
   Function   : Cretaes a Bio::EnsEMBL::Analysis::Tools::GeneBuildUtils::CollapsedCluster object and
              : populates it with exons and introns from the GeneCluster
   Exceptions : none
   Returnval  : Bio::EnsEMBL::Analysis::Tools::GeneBuildUtils::CollapsedCluster

=cut

sub collapse_cluster{
  my ($self,$cluster,$genes_by_strand) = @_;
  my @genes = $cluster->get_Genes;
  my @exon_clusters = @{$cluster->get_exon_clustering_from_gene_cluster};
  my $collapsed_cluster = Bio::EnsEMBL::Analysis::Tools::GeneBuildUtils::CollapsedCluster->new();

  foreach my $gene (@genes){
    foreach my $trans (@{$gene->get_all_Transcripts}){
      my @exons = sort {$a->start <=> $b->start } @{$trans->get_all_Exons};
      for (my $i =0 ; $i < scalar(@exons) ; $i++){
	my $exon_cluster;
	my $exon = $exons[$i];
	# figure out which exon cluster the exon belongs to
	# but only once for each non redundant exon
	unless ($collapsed_cluster->contains_exon($exon)){
	EC: foreach my $ec (@exon_clusters){
	    if ($ec->contains_exon($exon)){
	      $exon_cluster = $ec;
	      last EC;
	    }
	  }
	}
	# pass all the exons to the collapsed cluster object
	$collapsed_cluster->add_exon($exon,$exon_cluster);
	# cannot make an intron unless we have 2 exons to look at
	next if $i == 0;
	# make an intron object using exon extended so I can store the evidence sets
	my $intron = $collapsed_cluster->make_intron_from_exons($exons[$i-1],$exons[$i]);
	$intron->ev_set($exon->ev_set);
	$intron->biotype('intron');
	$collapsed_cluster->add_intron($intron);
      }
    }
  }

  # once we have added all of them into a nice non-redundant set
  # we need to assign them a score
  foreach my $exon (@{$collapsed_cluster->get_all_exons}){
    my $score = $self->score($collapsed_cluster,$exon,\@genes,'exon');
    $collapsed_cluster->exon_score($exon,$score);
  }

  foreach my $intron (@{$collapsed_cluster->get_all_introns}){
    # use all the genes on the strand when looking for introns overlapping genes as it helps
    # stop cluster joining
    my $score = $self->score($collapsed_cluster,$intron,$genes_by_strand->{$intron->strand},'intron');
    $collapsed_cluster->intron_score($intron,$score);
  }

  return $collapsed_cluster;
}

=head2 score

   Name       : score
   Arg[0]     : Bio::EnsEMBL::Analysis::Tools::GeneBuildUtils::CollapsedCluster
   Arg[1]     : Bio::EnsEMBL::Analysis:Tools::GeneBuildUtils::ExonExtended
   Arg[2]     : Array ref of Bio::EnsEMBL::Gene objects
   Arg[3]     : String 'exon' or 'intron'
   Function   : Assigns scores to introns and exons essentially the number of identical
              : features divided by the number of overlapping non-identical features.
              : Additional weights are added for end exons and est exons.
   Exceptions : none
   Returnval  : Scalar score

=cut

# scores are all assigned in the same way
sub score {
  my ($self,$collapsed_cluster,$feature,$genes,$type) = @_;
  my $total_score;
  # scores are calculated as the number of exoct matches / the number of overlaps
  # it is complicated by end exons and evidence sets however.
  $total_score = $collapsed_cluster->get_intron_count($feature) if $type eq 'intron';
  $total_score = $collapsed_cluster->get_exon_count($feature) if $type eq 'exon';
  $collapsed_cluster->exon_info($feature) if $VERBOSE;
  print " " if $VERBOSE;
  my $score = 0 ;
  # this way we penalise anything with less than exons than the minumum number of overlapping exons
  my $overlap = $MIN_CONSENSUS ;
  my $est_trans = $MIN_CONSENSUS;

  # the weighted score is:
  # 50% - proportion of exons in the stack that come from ests / number of overlapping est transcripts
  # 50% - the non est exons over the total number of overlapping transcripts of any type

  # count the number of overlapping exons / introns
  foreach my $gene (@$genes){
    foreach my $tran (@{$gene->get_all_Transcripts}){
      # if its an exon just add up how many genes of each type are in the cluster
      if ($type eq 'exon'){
	$est_trans++ if $tran->ev_set eq 'est';
	$overlap++;
      }
      next if scalar(@{$tran->get_all_Exons}) == 1;
      if ($type eq 'intron'){
	# for introns; count the number of overlapping trans from the entire slice
	# this creates a penatly for introns that span genes that
	# are not in the cluster because they dont share exons
	if ($feature->start < $tran->end() && $feature->end > $tran->start()){
	  $est_trans++ if $tran->ev_set eq 'est';
	  $overlap++;
	}	
      }
    }
  }
  # number of est exons that overlap the feature
  my $est_exons = $collapsed_cluster->ev_count_by_feature($feature,'est');

  # for end exon scores we allow the external boundary to be non identical
  if ($type eq 'exon' && $feature->is_terminal_exon &&  $feature->number_exons > 1){
    print "End exon $est_exons " if $VERBOSE;
    # count the total number of overlapping end exons
    $est_exons = $collapsed_cluster->ev_count_by_terminal_feature($feature,'est');
    $total_score = $collapsed_cluster->count_end_exon($feature);
    $self->throw("No evidence found for exon " . $collapsed_cluster->exon_info($feature) . "\n")
      unless defined($est_exons);
  }
  print "$type " if $VERBOSE;
  print "TOTAL $total_score OVERLAP $overlap " if $VERBOSE;

  if ($est_exons){
    print "WEIGHTED GENES $est_exons EST TRANS $est_trans " if $VERBOSE;
    # number of exons of the weighted type / number of overlapping trans of the weighted type;
    my $weighted_score = 0;
    my $unweighted_score = 0;
    if ($est_trans){
      $weighted_score = $est_exons / $est_trans  ;
      print "weighted score $weighted_score " if $VERBOSE;
    }
    # other componnt of score comes from all the exons
    $unweighted_score =  $total_score / $overlap ;
    print "unweighted score $unweighted_score \t" if $VERBOSE;
    $score = ($weighted_score + $unweighted_score) / 2;
  } else {
    # there are no est exons overlaping this feature score is simple identical / non-identical overlaps
    $score += ($total_score / $overlap) ;
    # if the feature overlaps an est transcript but not est exons, add a penalty
    if ($est_trans &&  $est_trans > 1.5 ){
      $score -= $EST_OVERLAP_PENALTY;
    }
  }
  # favor long end exons
  if ($feature->is_terminal_exon && $type eq 'exon' && $feature->number_exons > 1){
    # longest end exon 
    my $longest_end_exon = $collapsed_cluster->get_end_exon($feature)->length;
    my $multiplyer = $feature->length / $longest_end_exon;
    print "End exon penalty ".$feature->length." / ". $longest_end_exon." -.5\t" if $VERBOSE;
    if ($score >= 0){
      $score *= $multiplyer;
    } else {
      $score /= $multiplyer;
    }
    # end exons are penalised to prevent 'spindly' exons
    $score = $score - $END_EXON_PENALTY;
    print " ENDEXON " if $VERBOSE;
  } else {
    # penalty for  internal short exons / introns
      $score-= 0.5 if $feature->length < $SHORT_INTRON_PENALTY && $type eq 'intron';
      $score-= 0.5 if $feature->length < $SHORT_EXON_PENALTY   && $type eq 'exon';
  }
  print "$score \n" if $VERBOSE;
  return $score;
}

=head2 make_genes

   Name       : make_genes
   Arg[0]     : Array ref of Bio::EnsEMBL::Analysis::Tools::GeneBuildUtils::TranscriptExtended
   Arg[1]     : Bio::EnsEMBL::Analysis::Tools::Algorithms::GeneCluster
   Arg[2]     : Bio::EnsEMBL::Analysis::Tools::GeneBuildUtils::CollapsedCluster
   Function   : Takes an array ref of transcript objects and makes them into gene objects.
              : Deals with translations and orders the genes by score stores a subset
              : of the highest scoring genes on the output array.
   Exceptions : Throws if it cannot work out the translation start or end exons
   Returnval  : none

=cut

sub make_genes{
  my ($self,$transcripts,$cluster,$collapsed_cluster)= @_;
  my @genes;
  my @final_genes;
  $self->throw("No Transcripts to use\n") unless defined($transcripts);
 SCORE:  foreach my $t (sort {$b->score <=> $a->score } @$transcripts){
    my $transcript = $self->clone_transcript_extended($t);
    my $cds_start_exon;
    my $cds_end_exon;
    my $similarity_tran = $transcript->similarity;
    # each smilarity gene only allowed once, just the highest scoring one

    my $translation = Bio::EnsEMBL::Translation->new();
    my $strand= 1;
    my $slice =  $similarity_tran->slice;
    my $weighted_score;
    # dont compute the translation just set it
    # need the similarity gene that it was built from
    # so I can set the cds start and end on the transcript
    # if the start is in the same exon in the new transcript
    # just transfer it, if it is over the end of the new exon
    # trim back so that it is in the exon if you see what I mean
    # print stuff
    print "Transcript ".$transcript->start.":".$transcript->end.":".$transcript->strand." ".$transcript->biotype." " if $VERBOSE;
    print scalar(@{$transcript->get_all_Exons})." Exons\n" if $VERBOSE;
    my $last_exon;
    foreach my $exon (@{$transcript->get_all_Exons}){
      if ($last_exon){
	my $intron = $collapsed_cluster->make_intron_from_exons($last_exon,$exon);
      }
      $last_exon = $exon;
    }
    if ($transcript->strand == -1){
      my $invertedslice = $similarity_tran->slice->invert;
      $strand = -1;
      $similarity_tran = $similarity_tran->transfer($invertedslice);
      $transcript = $transcript->transfer($invertedslice);
    }
    my $cds_start = $similarity_tran->coding_region_start;
    my $cds_end = $similarity_tran->coding_region_end;
    my @exon_array = @{$transcript->get_all_Exons};
    my @similarity_exons = sort { $a->start <=> $b->start } @{$similarity_tran->get_all_Exons};


    my $number = scalar(@exon_array);
    # translation start and end is upside down when on the opposite strand
    for ( my $i = 0 ; $i < scalar(@exon_array) ; $i ++ ){
      # check for overlaps 
      my $exon = $exon_array[$i];	  
      # set phases to non coding by default
      $exon->phase( -1 );
      $exon->end_phase( -1 );
      foreach my $se (@similarity_exons){
	# set coding exons to have same phase as similarity exons
	if ($exon->start == $se->start){
	  $exon->phase($se->phase);
	}
	if ( $exon->end == $se->end ){
	  $exon->end_phase($se->end_phase);
	}
      }
      if ($exon->start <= $cds_start && $exon->end >= $cds_start){
	# translation start is in this exon
	$translation->start_Exon($exon);
	$translation->start($cds_start - $exon->start+1);
	$exon->phase( -1 ) 	if $exon->start < $cds_start;
      }
      if ($exon->start <= $cds_end && $exon->end >= $cds_end){
	# translation stop is in this exon
	$translation->end_Exon($exon);
	$translation->end($cds_end - $exon->start+1);
	$exon->end_phase( -1 ) 	if $exon->end > $cds_end;
      }
      # find the utr exon that shares the internal boundary with the first simgw exon
      if ($exon->end == $similarity_exons[0]->end){
	$cds_start_exon = $exon;
      }
      # find the utr exon that shares the internal boundary with the last simgw exon
      if ($exon->start == $similarity_exons[-1]->start){
	$cds_end_exon = $exon;
      }
    }
    unless ( $translation->start_Exon && $translation->start_Exon eq $cds_start_exon ) {
      # the start is not inside an exon
      throw("Cannot find the coding sequence start exon\n") unless $cds_start_exon;
      # trim back to the exon that matches the internal boundary
      # stay in frame though
      $translation->start_Exon($cds_start_exon);
      my $trim = $cds_start_exon->start - $cds_start;
      if ( $trim%3 == 0 ){
	$translation->start(1);
      } else {	
	$translation->start(1+(3-$trim%3));
      }
      $cds_start_exon->phase( -1 ) 	if $cds_start_exon->start < $cds_start;
    }
    unless ($translation->end_Exon && $translation->end_Exon eq $cds_end_exon ){
      # the start is not inside an exon
      throw("Cannot find the coding sequence end exon\n") unless $cds_end_exon;
      # trim back to the exon that matches the internal boundary
      # stay in frame 
      $translation->end_Exon($cds_end_exon);
      my $trim = $cds_end - $cds_end_exon->end;
      if ( $trim%3 == 0 ){
	$translation->end($cds_end_exon->length);
      } else {
	$translation->end($cds_end_exon->length-(3-$trim%3));
      }
      $cds_end_exon->end_phase( -1 ) 	if $cds_end_exon->end > $cds_end;
    }
    $transcript->translation($translation);
    if ($strand == -1){
      # flip it back onto the - strand
      my $slice = $transcript->slice->invert;
      $transcript = $transcript->transfer($slice);
    }	

    # penalise transcripts where a utr exon overlaps another exon that was previously coding
    $weighted_score = $self->compare_translation($transcript,$collapsed_cluster);
    print $transcript->score."\t" if $VERBOSE;
    print "\npeptide weighted score $weighted_score\n\n" if $VERBOSE;
    $transcript->score( $weighted_score );
    $transcript->analysis($similarity_tran->analysis);
    $transcript->slice($similarity_tran->slice);
    my $gene = Bio::EnsEMBL::Gene->new;<commit_msg>Fix translations when not adding UTR<commit_after>=pod

=head1 NAME Bio::EnsEMBL::Analysis::Runnable::TranscriptConsensus

=head1 SYNOPSIS

  my $runnable = Bio::EnsEMBL::Analysis::Runnable::TranscriptConsensus->new
     (
     -query            => $query,
     -analysis        => $analysis,
     -all_genes       => \%biotypes_to_genes , # ref to $hash{biotype_of_gene} = @all_genes_of_this_biotype
     -evidence_sets   => $evidence_sets ,
     -dnadb           => $db ,
     )

=head1 DESCRIPTION

TranscriptConsensus is an extension of TranscriptCoalescer that combines protein coding and est
transcripts to identify the best supported transcript models.
The initial gene sets are clustered, then collapsed into a non-redundant set of exons
and introns which are assigned scores according to the amount of supporting evidence they have.
The similarity genes have UTR added using the est transcripts and the resulting models are assigned
a score by summing the individual exon and intron scores.
The transcripts are sorted by score and the highest scoring models are made into
gene objets and written to the TranscriptCoalescer database.

=head1 CONTACT

Post questions to the EnsEMBL developer list: <<EMAIL>>

=cut


package Bio::EnsEMBL::Analysis::Runnable::TranscriptConsensus;

use strict;
use warnings; 

use Bio::EnsEMBL::Utils::Exception qw(throw warning info);
use Bio::EnsEMBL::Utils::Argument qw( rearrange );
use Bio::EnsEMBL::Analysis::Runnable::TranscriptCoalescer;
use Bio::EnsEMBL::Analysis::Tools::Algorithms::GeneCluster;
use Bio::EnsEMBL::Analysis::Tools::Algorithms::ExonCluster;
use Bio::EnsEMBL::Analysis::Tools::Algorithms::ClusterUtils;
use Bio::EnsEMBL::Analysis::Tools::GeneBuildUtils::TranscriptUtils qw (count_non_canonical_splice_sites) ;
use Bio::EnsEMBL::Analysis::Tools::GeneBuildUtils::EvidenceUtils qw ( clone_Evidence );
use Bio::EnsEMBL::Analysis::Tools::GeneBuildUtils::TranslationUtils qw( clone_Translation );
use Bio::EnsEMBL::Analysis::Config::GeneBuild::TranscriptCoalescer;
use Bio::EnsEMBL::Analysis::Config::GeneBuild::TranscriptConsensus;
use Bio::EnsEMBL::Analysis::Tools::GeneBuildUtils::CollapsedCluster;

use vars qw(@ISA);

@ISA = qw(Bio::EnsEMBL::Analysis::Runnable::TranscriptCoalescer);

=head2 run

   Name       : run
   Args       : none
   Function   : Pre filters genes and clusters them, then collapses each cluster to make non-redundant
              : sets of exons and introns, builds transcripts and adds UTR and assigns them scores,
              : makes genes out of the highest scoring transcripts
   Exceptions : none
   Returnval  : none

=cut

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

  my @allgenes = @{ $self->get_genes_by_evidence_set('simgw') }  ;
  push @allgenes , @{ $self->get_genes_by_evidence_set('est') }  ;

  if ($FILTER_SINGLETONS or   $FILTER_NON_CONSENSUS){
    @allgenes = @{$self->filter_genes(\@allgenes)};
  }
  
  my ($clusters, $non_clusters) = cluster_Genes(\@allgenes, $self->get_all_evidence_sets ) ;

  # create a hash of genes by strand to use when looking for overlapping genes
  # useful to look at ALL genes not just those in the cluster to prevent cluster joining
  foreach my $gene (@allgenes){
    push @{$genes_by_strand->{$gene->strand}}, $gene;
  }

  push @$clusters, @$non_clusters if (scalar(@$non_clusters) > 0 ) ;

  $count = 0;
  foreach my $cluster (@$clusters){
    my $simgw;
    # cluster has to contain at least one similarity gene to be worth continuing with
    next unless $cluster->get_Genes_by_Set('simgw');

    print "\nCluster $count\n" if $VERBOSE;
    print $cluster->start." ".$cluster->end." ".$cluster->strand."\n" if $VERBOSE;
    $count++;
    # collapse the cluster down to make a non redundat set of introns and exons with scores
    my $collapsed_cluster = $self->collapse_cluster($cluster,$genes_by_strand);

    # add UTR to the similarity transcripts and make scores for the transcripts
    my $transcripts = $self->make_transcripts($cluster,$collapsed_cluster);

    # sort the scored transcripts by score and make gene objects out of the best ones
    $self->make_genes($transcripts,$cluster,$collapsed_cluster);
  }
}

=head2 collapse_cluster

   Name       : collapse_cluster
   Arg[0]     : Bio::EnsEMBL::Analysis::Tools::Algorithms::GeneCluster
   Arg[1]     : Array ref of Bio::EnsEMBL::Gene objects
   Function   : Cretaes a Bio::EnsEMBL::Analysis::Tools::GeneBuildUtils::CollapsedCluster object and
              : populates it with exons and introns from the GeneCluster
   Exceptions : none
   Returnval  : Bio::EnsEMBL::Analysis::Tools::GeneBuildUtils::CollapsedCluster

=cut

sub collapse_cluster{
  my ($self,$cluster,$genes_by_strand) = @_;
  my @genes = $cluster->get_Genes;
  my @exon_clusters = @{$cluster->get_exon_clustering_from_gene_cluster};
  my $collapsed_cluster = Bio::EnsEMBL::Analysis::Tools::GeneBuildUtils::CollapsedCluster->new();

  foreach my $gene (@genes){
    foreach my $trans (@{$gene->get_all_Transcripts}){
      my @exons = sort {$a->start <=> $b->start } @{$trans->get_all_Exons};
      for (my $i =0 ; $i < scalar(@exons) ; $i++){
	my $exon_cluster;
	my $exon = $exons[$i];
	# figure out which exon cluster the exon belongs to
	# but only once for each non redundant exon
	unless ($collapsed_cluster->contains_exon($exon)){
	EC: foreach my $ec (@exon_clusters){
	    if ($ec->contains_exon($exon)){
	      $exon_cluster = $ec;
	      last EC;
	    }
	  }
	}
	# pass all the exons to the collapsed cluster object
	$collapsed_cluster->add_exon($exon,$exon_cluster);
	# cannot make an intron unless we have 2 exons to look at
	next if $i == 0;
	# make an intron object using exon extended so I can store the evidence sets
	my $intron = $collapsed_cluster->make_intron_from_exons($exons[$i-1],$exons[$i]);
	$intron->ev_set($exon->ev_set);
	$intron->biotype('intron');
	$collapsed_cluster->add_intron($intron);
      }
    }
  }

  # once we have added all of them into a nice non-redundant set
  # we need to assign them a score
  foreach my $exon (@{$collapsed_cluster->get_all_exons}){
    my $score = $self->score($collapsed_cluster,$exon,\@genes,'exon');
    $collapsed_cluster->exon_score($exon,$score);
  }

  foreach my $intron (@{$collapsed_cluster->get_all_introns}){
    # use all the genes on the strand when looking for introns overlapping genes as it helps
    # stop cluster joining
    my $score = $self->score($collapsed_cluster,$intron,$genes_by_strand->{$intron->strand},'intron');
    $collapsed_cluster->intron_score($intron,$score);
  }

  return $collapsed_cluster;
}

=head2 score

   Name       : score
   Arg[0]     : Bio::EnsEMBL::Analysis::Tools::GeneBuildUtils::CollapsedCluster
   Arg[1]     : Bio::EnsEMBL::Analysis:Tools::GeneBuildUtils::ExonExtended
   Arg[2]     : Array ref of Bio::EnsEMBL::Gene objects
   Arg[3]     : String 'exon' or 'intron'
   Function   : Assigns scores to introns and exons essentially the number of identical
              : features divided by the number of overlapping non-identical features.
              : Additional weights are added for end exons and est exons.
   Exceptions : none
   Returnval  : Scalar score

=cut

# scores are all assigned in the same way
sub score {
  my ($self,$collapsed_cluster,$feature,$genes,$type) = @_;
  my $total_score;
  # scores are calculated as the number of exoct matches / the number of overlaps
  # it is complicated by end exons and evidence sets however.
  $total_score = $collapsed_cluster->get_intron_count($feature) if $type eq 'intron';
  $total_score = $collapsed_cluster->get_exon_count($feature) if $type eq 'exon';
  $collapsed_cluster->exon_info($feature) if $VERBOSE;
  print " " if $VERBOSE;
  my $score = 0 ;
  # this way we penalise anything with less than exons than the minumum number of overlapping exons
  my $overlap = $MIN_CONSENSUS ;
  my $est_trans = $MIN_CONSENSUS;

  # the weighted score is:
  # 50% - proportion of exons in the stack that come from ests / number of overlapping est transcripts
  # 50% - the non est exons over the total number of overlapping transcripts of any type

  # count the number of overlapping exons / introns
  foreach my $gene (@$genes){
    foreach my $tran (@{$gene->get_all_Transcripts}){
      # if its an exon just add up how many genes of each type are in the cluster
      if ($type eq 'exon'){
	$est_trans++ if $tran->ev_set eq 'est';
	$overlap++;
      }
      next if scalar(@{$tran->get_all_Exons}) == 1;
      if ($type eq 'intron'){
	# for introns; count the number of overlapping trans from the entire slice
	# this creates a penatly for introns that span genes that
	# are not in the cluster because they dont share exons
	if ($feature->start < $tran->end() && $feature->end > $tran->start()){
	  $est_trans++ if $tran->ev_set eq 'est';
	  $overlap++;
	}	
      }
    }
  }
  # number of est exons that overlap the feature
  my $est_exons = $collapsed_cluster->ev_count_by_feature($feature,'est');

  # for end exon scores we allow the external boundary to be non identical
  if ($type eq 'exon' && $feature->is_terminal_exon &&  $feature->number_exons > 1){
    print "End exon $est_exons " if $VERBOSE;
    # count the total number of overlapping end exons
    $est_exons = $collapsed_cluster->ev_count_by_terminal_feature($feature,'est');
    $total_score = $collapsed_cluster->count_end_exon($feature);
    $self->throw("No evidence found for exon " . $collapsed_cluster->exon_info($feature) . "\n")
      unless defined($est_exons);
  }
  print "$type " if $VERBOSE;
  print "TOTAL $total_score OVERLAP $overlap " if $VERBOSE;

  if ($est_exons){
    print "WEIGHTED GENES $est_exons EST TRANS $est_trans " if $VERBOSE;
    # number of exons of the weighted type / number of overlapping trans of the weighted type;
    my $weighted_score = 0;
    my $unweighted_score = 0;
    if ($est_trans){
      $weighted_score = $est_exons / $est_trans  ;
      print "weighted score $weighted_score " if $VERBOSE;
    }
    # other componnt of score comes from all the exons
    $unweighted_score =  $total_score / $overlap ;
    print "unweighted score $unweighted_score \t" if $VERBOSE;
    $score = ($weighted_score + $unweighted_score) / 2;
  } else {
    # there are no est exons overlaping this feature score is simple identical / non-identical overlaps
    $score += ($total_score / $overlap) ;
    # if the feature overlaps an est transcript but not est exons, add a penalty
    if ($est_trans &&  $est_trans > 1.5 ){
      $score -= $EST_OVERLAP_PENALTY;
    }
  }
  # favor long end exons
  if ($feature->is_terminal_exon && $type eq 'exon' && $feature->number_exons > 1){
    # longest end exon 
    my $longest_end_exon = $collapsed_cluster->get_end_exon($feature)->length;
    my $multiplyer = $feature->length / $longest_end_exon;
    print "End exon penalty ".$feature->length." / ". $longest_end_exon." -.5\t" if $VERBOSE;
    if ($score >= 0){
      $score *= $multiplyer;
    } else {
      $score /= $multiplyer;
    }
    # end exons are penalised to prevent 'spindly' exons
    $score = $score - $END_EXON_PENALTY;
    print " ENDEXON " if $VERBOSE;
  } else {
    # penalty for  internal short exons / introns
      $score-= 0.5 if $feature->length < $SHORT_INTRON_PENALTY && $type eq 'intron';
      $score-= 0.5 if $feature->length < $SHORT_EXON_PENALTY   && $type eq 'exon';
  }
  print "$score \n" if $VERBOSE;
  return $score;
}

=head2 make_genes

   Name       : make_genes
   Arg[0]     : Array ref of Bio::EnsEMBL::Analysis::Tools::GeneBuildUtils::TranscriptExtended
   Arg[1]     : Bio::EnsEMBL::Analysis::Tools::Algorithms::GeneCluster
   Arg[2]     : Bio::EnsEMBL::Analysis::Tools::GeneBuildUtils::CollapsedCluster
   Function   : Takes an array ref of transcript objects and makes them into gene objects.
              : Deals with translations and orders the genes by score stores a subset
              : of the highest scoring genes on the output array.
   Exceptions : Throws if it cannot work out the translation start or end exons
   Returnval  : none

=cut

sub make_genes{
  my ($self,$transcripts,$cluster,$collapsed_cluster)= @_;
  my @genes;
  my @final_genes;
  $self->throw("No Transcripts to use\n") unless defined($transcripts);
 SCORE:  foreach my $t (sort {$b->score <=> $a->score } @$transcripts){
    my $transcript = $self->clone_transcript_extended($t);
    my $cds_start_exon;
    my $cds_end_exon;
    my $similarity_tran = $transcript->similarity;
    # each smilarity gene only allowed once, just the highest scoring one

    my $translation = Bio::EnsEMBL::Translation->new();
    my $strand= 1;
    my $slice =  $similarity_tran->slice;
    my $weighted_score;
    if ($ADD_UTR){
      # dont compute the translation just set it
      # need the similarity gene that it was built from
      # so I can set the cds start and end on the transcript
      # if the start is in the same exon in the new transcript
      # just transfer it, if it is over the end of the new exon
      # trim back so that it is in the exon if you see what I mean
      # print stuff
      print "Transcript ".$transcript->start.":".$transcript->end.":".$transcript->strand." ".$transcript->biotype." " if $VERBOSE;
      print scalar(@{$transcript->get_all_Exons})." Exons\n" if $VERBOSE;
      my $last_exon;
      foreach my $exon (@{$transcript->get_all_Exons}){
	if ($last_exon){
	  my $intron = $collapsed_cluster->make_intron_from_exons($last_exon,$exon);
	}
	$last_exon = $exon;
      }
      if ($transcript->strand == -1){
	my $invertedslice = $similarity_tran->slice->invert;
	$strand = -1;
	$similarity_tran = $similarity_tran->transfer($invertedslice);
	$transcript = $transcript->transfer($invertedslice);
      }
      my $cds_start = $similarity_tran->coding_region_start;
      my $cds_end = $similarity_tran->coding_region_end;
      my @exon_array = @{$transcript->get_all_Exons};
      my @similarity_exons = sort { $a->start <=> $b->start } @{$similarity_tran->get_all_Exons};
      
      
      my $number = scalar(@exon_array);
      # translation start and end is upside down when on the opposite strand
      for ( my $i = 0 ; $i < scalar(@exon_array) ; $i ++ ){
	# check for overlaps 
	my $exon = $exon_array[$i];	  
	# set phases to non coding by default
	$exon->phase( -1 );
	$exon->end_phase( -1 );
	foreach my $se (@similarity_exons){
	  # set coding exons to have same phase as similarity exons
	  if ($exon->start == $se->start){
	    $exon->phase($se->phase);
	}
	  if ( $exon->end == $se->end ){
	    $exon->end_phase($se->end_phase);
	  }
	}
	if ($exon->start <= $cds_start && $exon->end >= $cds_start){
	  # translation start is in this exon
	  $translation->start_Exon($exon);
	  $translation->start($cds_start - $exon->start+1);
	  $exon->phase( -1 ) 	if $exon->start < $cds_start;
	}
	if ($exon->start <= $cds_end && $exon->end >= $cds_end){
	  # translation stop is in this exon
	  $translation->end_Exon($exon);
	  $translation->end($cds_end - $exon->start+1);
	  $exon->end_phase( -1 ) 	if $exon->end > $cds_end;
	}
	# find the utr exon that shares the internal boundary with the first simgw exon
	if ($exon->end == $similarity_exons[0]->end){
	  $cds_start_exon = $exon;
	}
	# find the utr exon that shares the internal boundary with the last simgw exon
	if ($exon->start == $similarity_exons[-1]->start){
	  $cds_end_exon = $exon;
	}
      }
      unless ( $translation->start_Exon && $translation->start_Exon eq $cds_start_exon ) {
	# the start is not inside an exon
	throw("Cannot find the coding sequence start exon\n") unless $cds_start_exon;
	# trim back to the exon that matches the internal boundary
	# stay in frame though
	$translation->start_Exon($cds_start_exon);
	my $trim = $cds_start_exon->start - $cds_start;
	if ( $trim%3 == 0 ){
	  $translation->start(1);
	} else {	
	  $translation->start(1+(3-$trim%3));
	}
	$cds_start_exon->phase( -1 ) 	if $cds_start_exon->start < $cds_start;
      }
      unless ($translation->end_Exon && $translation->end_Exon eq $cds_end_exon ){
	# the start is not inside an exon
	throw("Cannot find the coding sequence end exon\n") unless $cds_end_exon;
	# trim back to the exon that matches the internal boundary
	# stay in frame 
      $translation->end_Exon($cds_end_exon);
	my $trim = $cds_end - $cds_end_exon->end;
	if ( $trim%3 == 0 ){
	  $translation->end($cds_end_exon->length);
	} else {
	  $translation->end($cds_end_exon->length-(3-$trim%3));
	}
	$cds_end_exon->end_phase( -1 ) 	if $cds_end_exon->end > $cds_end;
      }
      $transcript->translation($translation);
      if ($strand == -1){
	# flip it back onto the - strand
	my $slice = $transcript->slice->invert;
	$transcript = $transcript->transfer($slice);
      }	
    } else {
      # dont need to mess with translations if I am not adding UTR...
      my $newtranslation = clone_Translation($similarity_tran,$transcript);
      $transcript->translation($newtranslation);
    }
	
    # penalise transcripts where a utr exon overlaps another exon that was previously coding
    $weighted_score = $self->compare_translation($transcript,$collapsed_cluster);
    print $transcript->score."\t" if $VERBOSE;
    print "\npeptide weighted score $weighted_score\n\n" if $VERBOSE;
    $transcript->score( $weighted_score );
    $transcript->analysis($similarity_tran->analysis);
    $transcript->slice($similarity_tran->slice);
    my $gene = Bio::EnsEMBL::Gene->new;
{'fraction_non_alphanumeric': 0.12303883467196361, 'fraction_numerical': 0.003915133719089156, 'mean_word_length': 3.793431764868221, 'pattern_counts': {'":': 4, '<': 31, '<?xml version=': 0, '>': 336, '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': '7658337', 'n_tokens_mistral': 12011, 'n_tokens_neox': 11195, 'n_words': 3287}