<commit_before>package Bio::EnsEMBL::Pipeline::Runnable::ESTDescriminator;

use strict;
use Bio::EnsEMBL::Pipeline::Alignment::EvidenceAlignment;
use Bio::EnsEMBL::Pipeline::GeneDuplication::CodonBasedAlignment;
use Bio::Align::DNAStatistics;
use Bio::AlignIO;
use Bio::LocatableSeq;
use Bio::EnsEMBL::Utils::Exception qw(throw warning info); 
use Bio::EnsEMBL::Utils::Argument qw(rearrange);

sub new {
  my ($class, @args) = @_;

  my $self = bless {},$class;

  my ($genes,
      $est_db,
      $est_coverage_cutoff,
      $seqfetcher,
     ) = rearrange([qw(GENES
                       EST_DB
                       EST_COVERAGE_CUTOFF
		       SEQUENCE_FETCHER
                      )],@args);

  throw("No genes to analyse")
    unless $genes;

  $self->_genes($genes);

  throw("Constructor needs adaptors to the database which " .
	"contains your ests.")
    unless (defined $est_db);

  $self->_est_db($est_db);

  throw("Cant proceed without a sequence fetcher.")
    unless (defined $seqfetcher);  

  $self->_seqfetcher($seqfetcher);

  if ($est_coverage_cutoff) {
    $self->_est_coverage_cutoff($est_coverage_cutoff)
  } else {
    $self->_est_coverage_cutoff(0.8)
  }

  return $self
}

sub run {
  my $self = shift;

  my %hash_of_ests_and_gene_mapping = 
    $self->_cluster_ests_with_genes;

  ### This is were I fell asleep


}

sub print_shared_ESTs {
  my $self = shift;

  my %est_hash;

  foreach my $gene (@{$self->_genes}){
    my $overlapping_ests = $self->_find_overlapping_ests($gene);

    foreach my $est (@$overlapping_ests){
      $est_hash{$est->get_all_Exons->[0]->get_all_supporting_features->[0]->hseqname}++;
    }
  }

  my $shared_ests = 0;

  foreach my $est_name (keys %est_hash){
    if ($est_hash{$est_name} > 1) {
      print STDERR $est_name . "\tmapped to " . $est_hash{$est_name} . " genes.\n";
      $shared_ests++; 
    }
  }

  print STDERR "Had $shared_ests ESTs that mapped to more than one gene.\n";

  return 1
}

sub _cluster_ests_with_genes {
  my $self = shift;

  my %ests_clustered_to_genes;

  # Execute the subroutines that will generate the 
  # needed data.

  $self->_calculate_gene_vs_est_distances;
  $self->_calculate_intergene_distances;

  # Derive a list of unique EST ids

  my %all_shared_ests;
  foreach my $gene_id (keys %{$self->_gene_lookup}){
    my $est_ids = $self->_find_shared_ests_by_gene_id;
    foreach my $est_id (@$est_ids){
      $all_shared_ests{$est_id}++;
    }
  }

  # Work through ESTs, grouping each with the most 
  # appropriate gene(s).

  foreach my $est_id (keys %all_shared_ests){

    # Make an array of simple hashes, containing
    # both EST vs gene and inter-gene distances
    # for this particular EST.

    my @est_gene_distances;

    my $genes_that_share_this_est = $self->_find_genes_by_est_id;

      # EST vs gene distances

    foreach my $gene_id (@$genes_that_share_this_est){
      my %pairwise_comparison = 
	('type'     => 'est_vs_gene',
	 'gene_ids'  => [$gene_id],
	 'distance' => $self->_pairwise_distance($gene_id, $est_id));
      push @est_gene_distances, \%pairwise_comparison;
    }

      # Inter-gene distances

    for (my $i = 0; $i < scalar @$genes_that_share_this_est; $i++) {
      for (my $k = $i + 1; $k < scalar @$genes_that_share_this_est; $k++) {
	my %pairwise_comparison = 
	  ('type'     => 'inter_gene',
	   'gene_ids' => [$genes_that_share_this_est->[$i], 
			  $genes_that_share_this_est->[$k]],
	   'distance' => $self->_pairwise_distance($genes_that_share_this_est->[$i],
						   $genes_that_share_this_est->[$k]));
	push @est_gene_distances, \%pairwise_comparison;
      }
    }


    # Sort the array (of hashes) by distance in ascending order.

    @est_gene_distances = 
      sort {$a->{distance} <=> $b->{distance}} @est_gene_distances;

    # Group the EST to specific genes according
    # to following rules:
    # * EST maps to one gene if EST vs gene distance
    #     is smaller than the smallest inter-gene 
    #     distance.
    # * Where genes are more closely related than the
    #     EST vs gene distances, this EST must be mapped
    #     to either both genes or neither gene.
    # * The inherant error attached to the distance
    #     must be considered when saying an EST is more
    #     closely related to one gene rather than
    #     another. ----------------------------------------THIS IS PRESENTLY NOT IMPLEMENTED.


    my $closest_gene_id;
    my @closely_related_genes;
    foreach my $comparison (@est_gene_distances){
      if ($comparison->{type} eq 'est_vs_gene') {
	$closest_gene_id = $comparison->{gene_ids}->[0];
	last
      } elsif ($comparison->{type} eq 'inter_gene') {
	push @closely_related_genes, $comparison;
      }
    }

    throw("Failed to find most closely related gene.")
      unless defined $closest_gene_id;

    if (scalar @closely_related_genes == 0){
      push @{$ests_clustered_to_genes{$est_id}}, $closest_gene_id;
      next
    }

    foreach my $gene_pair (@closely_related_genes){
      if ($gene_pair->{gene_ids}->[0] eq $closest_gene_id) {
	push @{$ests_clustered_to_genes{$est_id}}, $gene_pair->{gene_ids}->[1]
      } elsif ($gene_pair->{gene_ids}->[1] eq $closest_gene_id) {
	push @{$ests_clustered_to_genes{$est_id}}, $gene_pair->{gene_ids}->[0]
      }
    }

    throw("Erroneously didn't map EST to any gene.")
      unless defined $ests_clustered_to_genes{$est_id};
  }

  return \%ests_clustered_to_genes
}

sub _calculate_gene_vs_est_distances {
  my $self = shift;

  my $genes_and_shared_ests = $self->_identify_shared_ESTs;
  my %gene_vs_est_distance_matrix;

  foreach my $gene_id (keys %$genes_and_shared_ests) {
    # THERE IS A SMARTER WAY TO DO THIS-------------------
    my $transcript = $self->_gene_lookup->{$gene_id}->get_all_Transcripts->[0];

    # Gather all EST features
    my @est_supporting_features;
    foreach my $est_id (keys %{$genes_and_shared_ests->{$gene_id}}){
      push @est_supporting_features, 
           @{$genes_and_shared_ests->{$gene_id}->{$est_id}};    
    }

    # Obtain an alignment of our transcript with the shared ESTs
    my $evidence_align = 
      Bio::EnsEMBL::Pipeline::Alignment::EvidenceAlignment->new(
        -transcript          => $transcript,
        -dbadaptor           => $self->_gene_lookup->{$gene_id}->adaptor->db,
        -seqfetcher          => $self->_seqfetcher,
        -padding             => 15,
        -supporting_features => \@est_supporting_features 
        );

    my $ens_align = $evidence_align->retrieve_alignment( 
                   -type            => 'nucleotide',
                   -remove_introns  => 0,
                   -merge_sequences => 1);

    # A horrible and hopefully temporary hack to 
    # ensure all sequences are of the same length.
    # This is due to a failing in Bio::EnsEMBL::Pipeline::Alignmenti::EvidenceAlignment.

    my @ens_align_seqs = @{$ens_align->fetch_AlignmentSeqs};
    my $longest_seq = 0;
    foreach my $ens_align_seq (@ens_align_seqs) {
      my $length = length($ens_align_seq->seq);
      if ($length > $longest_seq){
        $longest_seq = length($ens_align_seq->seq)
      } 
    }  
    for (my $i = 0; $i < scalar @ens_align_seqs; $i++) {
      my $length = length($ens_align_seqs[$i]->seq);
      if ($length < $longest_seq) {
        my $shortfall = $longest_seq - $length;
        my $tack_on = '-' x $shortfall;
        my $augmented_seq = $ens_align_seqs[$i]->seq . $tack_on;
        $ens_align_seqs[$i]->seq($augmented_seq); 
      }  
    } 

    # Convert this alignment into a bioperl AlignI object.
    my $bp_align = Bio::SimpleAlign->new();

    my @sequence_order;
#    foreach my $align_seq (@{$ens_align->fetch_AlignmentSeqs}){
    foreach my $align_seq (@ens_align_seqs){

      # Skip the genomic sequence in our alignment - use the exon
      # sequence for calculating distance
      next if $align_seq->name eq 'genomic_sequence';

      my $seq = Bio::Seq->new(-display_id => $align_seq->name,
	  		      -seq        => $align_seq->seq);
      my $locatable_seq = Bio::LocatableSeq->new(
                            -id       => $align_seq->name,
                            -seq      => $align_seq->seq,
                            -start    => 1,
                            -end      => length($align_seq->seq),
                            -alphabet => 'dna');

      $bp_align->add_seq($locatable_seq);

      push @sequence_order, $align_seq->name;
    }

    # Use bioperl object to calculate our pairwise sequence distances.
    my $dna_stat = Bio::Align::DNAStatistics->new;

    my $distances = $dna_stat->distance(-align  => $bp_align,
					-method => 'Kimura');

    # Store result of first column of result matrix (ie. gene vs ests)
    for (my $row = 1; $row < scalar @$distances; $row++) {
      next if $sequence_order[$row-1] eq 'exon_sequence';
      next unless defined $distances->[$row];

      $self->_pairwise_distance($gene_id, 
				$sequence_order[$row-1], 
				$distances->[$row]->[1])
    }
  }
<commit_msg>Fixed small error where appropriate argument was not passed to subroutine.<commit_after>package Bio::EnsEMBL::Pipeline::Runnable::ESTDescriminator;

use strict;
use Bio::EnsEMBL::Pipeline::Alignment::EvidenceAlignment;
use Bio::EnsEMBL::Pipeline::GeneDuplication::CodonBasedAlignment;
use Bio::Align::DNAStatistics;
use Bio::AlignIO;
use Bio::LocatableSeq;
use Bio::EnsEMBL::Utils::Exception qw(throw warning info); 
use Bio::EnsEMBL::Utils::Argument qw(rearrange);

sub new {
  my ($class, @args) = @_;

  my $self = bless {},$class;

  my ($genes,
      $est_db,
      $est_coverage_cutoff,
      $seqfetcher,
     ) = rearrange([qw(GENES
                       EST_DB
                       EST_COVERAGE_CUTOFF
		       SEQUENCE_FETCHER
                      )],@args);

  throw("No genes to analyse")
    unless $genes;

  $self->_genes($genes);

  throw("Constructor needs adaptors to the database which " .
	"contains your ests.")
    unless (defined $est_db);

  $self->_est_db($est_db);

  throw("Cant proceed without a sequence fetcher.")
    unless (defined $seqfetcher);  

  $self->_seqfetcher($seqfetcher);

  if ($est_coverage_cutoff) {
    $self->_est_coverage_cutoff($est_coverage_cutoff)
  } else {
    $self->_est_coverage_cutoff(0.8)
  }

  return $self
}

sub run {
  my $self = shift;

  my %hash_of_ests_and_gene_mapping = 
    %{$self->_cluster_ests_with_genes};

  ### This is were I fell asleep


}

sub print_shared_ESTs {
  my $self = shift;

  my %est_hash;

  foreach my $gene (@{$self->_genes}){
    my $overlapping_ests = $self->_find_overlapping_ests($gene);

    foreach my $est (@$overlapping_ests){
      $est_hash{$est->get_all_Exons->[0]->get_all_supporting_features->[0]->hseqname}++;
    }
  }

  my $shared_ests = 0;

  foreach my $est_name (keys %est_hash){
    if ($est_hash{$est_name} > 1) {
      print STDERR $est_name . "\tmapped to " . $est_hash{$est_name} . " genes.\n";
      $shared_ests++; 
    }
  }

  print STDERR "Had $shared_ests ESTs that mapped to more than one gene.\n";

  return 1
}

sub _cluster_ests_with_genes {
  my $self = shift;

  my %ests_clustered_to_genes;

  # Execute the subroutines that will generate the 
  # needed data.

  $self->_calculate_gene_vs_est_distances;
  $self->_calculate_intergene_distances;

  # Derive a list of unique EST ids

  my %all_shared_ests;
  foreach my $gene_id (keys %{$self->_gene_lookup}){
    my $est_ids = $self->_find_shared_ests_by_gene_id($gene_id);
    foreach my $est_id (@$est_ids){
      $all_shared_ests{$est_id}++;
    }
  }

  # Work through ESTs, grouping each with the most 
  # appropriate gene(s).

  foreach my $est_id (keys %all_shared_ests){

    # Make an array of simple hashes, containing
    # both EST vs gene and inter-gene distances
    # for this particular EST.

    my @est_gene_distances;

    my $genes_that_share_this_est = $self->_find_genes_by_est_id($est_id);

      # EST vs gene distances

    foreach my $gene_id (@$genes_that_share_this_est){
      my %pairwise_comparison = 
	('type'     => 'est_vs_gene',
	 'gene_ids'  => [$gene_id],
	 'distance' => $self->_pairwise_distance($gene_id, $est_id));
      push @est_gene_distances, \%pairwise_comparison;
    }

      # Inter-gene distances

    for (my $i = 0; $i < scalar @$genes_that_share_this_est; $i++) {
      for (my $k = $i + 1; $k < scalar @$genes_that_share_this_est; $k++) {
	my %pairwise_comparison = 
	  ('type'     => 'inter_gene',
	   'gene_ids' => [$genes_that_share_this_est->[$i], 
			  $genes_that_share_this_est->[$k]],
	   'distance' => $self->_pairwise_distance($genes_that_share_this_est->[$i],
						   $genes_that_share_this_est->[$k]));
	push @est_gene_distances, \%pairwise_comparison;
      }
    }


    # Sort the array (of hashes) by distance in ascending order.

    @est_gene_distances = 
      sort {$a->{distance} <=> $b->{distance}} @est_gene_distances;

    # Group the EST to specific genes according
    # to following rules:
    # * EST maps to one gene if EST vs gene distance
    #     is smaller than the smallest inter-gene 
    #     distance.
    # * Where genes are more closely related than the
    #     EST vs gene distances, this EST must be mapped
    #     to either both genes or neither gene.
    # * The inherant error attached to the distance
    #     must be considered when saying an EST is more
    #     closely related to one gene rather than
    #     another. ----------------------------------------THIS IS PRESENTLY NOT IMPLEMENTED.


    my $closest_gene_id;
    my @closely_related_genes;
    foreach my $comparison (@est_gene_distances){
      if ($comparison->{type} eq 'est_vs_gene') {
	$closest_gene_id = $comparison->{gene_ids}->[0];
	last
      } elsif ($comparison->{type} eq 'inter_gene') {
	push @closely_related_genes, $comparison;
      }
    }

    throw("Failed to find most closely related gene.")
      unless defined $closest_gene_id;

    if (scalar @closely_related_genes == 0){
      push @{$ests_clustered_to_genes{$est_id}}, $closest_gene_id;
      next
    }

    foreach my $gene_pair (@closely_related_genes){
      if ($gene_pair->{gene_ids}->[0] eq $closest_gene_id) {
	push @{$ests_clustered_to_genes{$est_id}}, $gene_pair->{gene_ids}->[1]
      } elsif ($gene_pair->{gene_ids}->[1] eq $closest_gene_id) {
	push @{$ests_clustered_to_genes{$est_id}}, $gene_pair->{gene_ids}->[0]
      }
    }

    throw("Erroneously didn't map EST to any gene.")
      unless defined $ests_clustered_to_genes{$est_id};
  }

  return \%ests_clustered_to_genes
}

sub _calculate_gene_vs_est_distances {
  my $self = shift;

  my $genes_and_shared_ests = $self->_identify_shared_ESTs;
  my %gene_vs_est_distance_matrix;

  foreach my $gene_id (keys %$genes_and_shared_ests) {
    # THERE IS A SMARTER WAY TO DO THIS-------------------
    my $transcript = $self->_gene_lookup->{$gene_id}->get_all_Transcripts->[0];

    # Gather all EST features
    my @est_supporting_features;
    foreach my $est_id (keys %{$genes_and_shared_ests->{$gene_id}}){
      push @est_supporting_features, 
           @{$genes_and_shared_ests->{$gene_id}->{$est_id}};    
    }

    # Obtain an alignment of our transcript with the shared ESTs
    my $evidence_align = 
      Bio::EnsEMBL::Pipeline::Alignment::EvidenceAlignment->new(
        -transcript          => $transcript,
        -dbadaptor           => $self->_gene_lookup->{$gene_id}->adaptor->db,
        -seqfetcher          => $self->_seqfetcher,
        -padding             => 15,
        -supporting_features => \@est_supporting_features 
        );

    my $ens_align = $evidence_align->retrieve_alignment( 
                   -type            => 'nucleotide',
                   -remove_introns  => 0,
                   -merge_sequences => 1);

    # A horrible and hopefully temporary hack to 
    # ensure all sequences are of the same length.
    # This is due to a failing in Bio::EnsEMBL::Pipeline::Alignmenti::EvidenceAlignment.

    my @ens_align_seqs = @{$ens_align->fetch_AlignmentSeqs};
    my $longest_seq = 0;
    foreach my $ens_align_seq (@ens_align_seqs) {
      my $length = length($ens_align_seq->seq);
      if ($length > $longest_seq){
        $longest_seq = length($ens_align_seq->seq)
      } 
    }  
    for (my $i = 0; $i < scalar @ens_align_seqs; $i++) {
      my $length = length($ens_align_seqs[$i]->seq);
      if ($length < $longest_seq) {
        my $shortfall = $longest_seq - $length;
        my $tack_on = '-' x $shortfall;
        my $augmented_seq = $ens_align_seqs[$i]->seq . $tack_on;
        $ens_align_seqs[$i]->seq($augmented_seq); 
      }  
    } 

    # Convert this alignment into a bioperl AlignI object.
    my $bp_align = Bio::SimpleAlign->new();

    my @sequence_order;
#    foreach my $align_seq (@{$ens_align->fetch_AlignmentSeqs}){
    foreach my $align_seq (@ens_align_seqs){

      # Skip the genomic sequence in our alignment - use the exon
      # sequence for calculating distance
      next if $align_seq->name eq 'genomic_sequence';

      my $locatable_seq = Bio::LocatableSeq->new(
                            -id       => $align_seq->name,
                            -seq      => $align_seq->seq,
                            -start    => 1,
                            -end      => length($align_seq->seq),
                            -alphabet => 'dna');

      $bp_align->add_seq($locatable_seq);

      push @sequence_order, $align_seq->name;
    }

    # Use bioperl object to calculate our pairwise sequence distances.
    my $dna_stat = Bio::Align::DNAStatistics->new;

    my $distances = $dna_stat->distance(-align  => $bp_align,
					-method => 'Kimura');

    # Store result of first column of result matrix (ie. gene vs ests)
print STDERR "Pairwise distances between gene sequence [$gene_id] and ESTs : \n";
    for (my $row = 1; $row < scalar @$distances; $row++) {
      next if $sequence_order[$row-1] eq 'exon_sequence';
      next unless defined $distances->[$row];
print STDERR "Pairwise distance : [$row] : " . $distances->[$row]->[1] . "\n";
      $self->_pairwise_distance($gene_id, 
				$sequence_order[$row-1], 
				$distances->[$row]->[1])
    }
  }

{'fraction_non_alphanumeric': 0.13484815140845072, 'fraction_numerical': 0.003025968309859155, 'mean_word_length': 2.8437301755127935, 'pattern_counts': {'":': 0, '<': 15, '<?xml version=': 0, '>': 206, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '6034138', 'n_tokens_mistral': 6582, 'n_tokens_neox': 6013, 'n_words': 1408}