Inheritance diagram for Bio::EnsEMBL::IdMapping::ExonScoreBuilder:

Public Member Functions
public	score_exons ()

public	overlap_score ()

public	exonerate_score ()

public	build_overlap_scores ()

public	sort_exons ()

public	compare_exon_containers ()

public	calc_overlap_score ()

public	run_exonerate ()

public	exon_fasta_file ()

public	dump_filtered_exons ()

public	write_filtered_exons ()

public	parse_exonerate_results ()

public	non_mapped_transcript_rescore ()

Public Member Functions inherited from Bio::EnsEMBL::IdMapping::ScoreBuilder
public Bio::EnsEMBL::IdMapping::ScoredMappingMatrix	create_shrinked_matrix ()

public void	internal_id_rescore ()

public void	log_matrix_stats ()

Public Member Functions inherited from Bio::EnsEMBL::IdMapping::BaseObject
public $this	new ()

public Filehandle	get_filehandle ()

public Boolean	file_exists ()

public SQLStatement	fetch_value_from_db ()

public Int	dump_table_to_file ()

public Int	upload_file_into_table ()

public Bio::EnsEMBL::Utils::Logger	logger ()

public Bio::EnsEMBL::Utils::ConfParser	conf ()

public Bio::EnsEMBL::IdMapping::Cache	cache ()

Detailed Description

Synopsis

Description

Combines ExonScoreBuilder, ExonDirectMapper and ExonerateRunner from
Java application.

Definition at line 19 of file ExonScoreBuilder.pm.

Member Function Documentation

◆ build_overlap_scores()

public Bio::EnsEMBL::IdMapping::ExonScoreBuilder::build_overlap_scores ( )

Undocumented method

Code:

click to view

sub build_overlap_scores {
  my ( $self, $matrix ) = @_;
 
  unless ($matrix
      and $matrix->isa('Bio::EnsEMBL::IdMapping::ScoredMappingMatrix') )
  {
    throw('Need a Bio::EnsEMBL::IdMapping::ScoredMappingMatrix.');
  }
 
  # get sorted list of exon containers
  $self->logger->info( "Reading sorted exons from cache...\n",
                       1, 'stamped' );
 
  my @source_exons =
    $self->sort_exons( [
        values %{ $self->cache->get_by_name( 'exons_by_id', 'source' ) }
      ] );
 
  my @target_exons =
    $self->sort_exons( [
        values %{ $self->cache->get_by_name( 'exons_by_id', 'target' ) }
      ] );
 
  $self->logger->info( "Done.\n", 1, 'stamped' );
 
  # get first source and target exon container
  my $source_ec = shift(@source_exons);
  my $target_ec = shift(@target_exons);
 
  my %source_overlap = ();
  my %target_overlap = ();
 
  $self->logger->info( "Scoring...\n", 1, 'stamped' );
 
  while ( defined($source_ec) || defined($target_ec) ) {
    my $add_source = 0;
    my $add_target = 0;
 
    # compare exon containers
    if ( defined($source_ec) && defined($target_ec) ) {
      my $cmp =
        $self->compare_exon_containers( $source_ec, $target_ec );
      if ( $cmp <= 0 ) { $add_source = 1 }
      if ( $cmp >= 0 ) { $add_target = 1 }
    }
    elsif ( defined($source_ec) ) {
      $add_source = 1;
    }
    elsif ( defined($target_ec) ) {
      $add_target = 1;
    }
    else {
      die('The world is a strange place');
    }
 
    if ($add_source) {
      if ( exists( $source_overlap{ $source_ec->[0] } ) ) {
        # remove exon from list of overlapping source exons to score
        # target against
        delete $source_overlap{ $source_ec->[0] };
      }
      else {
        # add exon to list of overlapping source exons to score target
        # against
        $source_overlap{ $source_ec->[0] } = $source_ec->[0];
 
        # score source exon against all target exons in current overlap
        # list
        foreach my $target_exon ( values %target_overlap ) {
          if ( defined( $matrix->get_score(
                                   $source_ec->[0]->id, $target_exon->id
                        ) ) )
          {
            next;
          }
 
          $self->calc_overlap_score( $source_ec->[0], $target_exon,
                                     $matrix );
        }
      }
 
      # get next source exon container
      $source_ec = shift(@source_exons);
    } ## end if ($add_source)
 
    if ($add_target) {
      if ( exists( $target_overlap{ $target_ec->[0] } ) ) {
        # remove exon from list of overlapping target exons to score
        # source against
        delete $target_overlap{ $target_ec->[0] };
      }
      else {
        # add exon to list of overlapping target exons to score source
        # against
        $target_overlap{ $target_ec->[0] } = $target_ec->[0];
 
        # score target exon against all source exons in current overlap
        # list
        foreach my $source_exon ( values %source_overlap ) {
          if ( defined( $matrix->get_score(
                                   $source_exon->id, $target_ec->[0]->id
                        ) ) )
          {
            next;
          }
 
          $self->calc_overlap_score( $source_exon, $target_ec->[0],
                                     $matrix );
        }
      }
 
      # get next target exon container
      $target_ec = shift(@target_exons);
    } ## end if ($add_target)
  } ## end while ( defined($source_ec...))
 
  $self->logger->info( "Done.\n", 1, 'stamped' );
 
  return $matrix;
} ## end sub build_overlap_scores

◆ calc_overlap_score()

public Bio::EnsEMBL::IdMapping::ExonScoreBuilder::calc_overlap_score ( )

Undocumented method

Code:

click to view

◆ compare_exon_containers()

public Bio::EnsEMBL::IdMapping::ExonScoreBuilder::compare_exon_containers ( )

Undocumented method

Code:

click to view

◆ dump_filtered_exons()

public Bio::EnsEMBL::IdMapping::ExonScoreBuilder::dump_filtered_exons ( )

Undocumented method

Code:

click to view

◆ exon_fasta_file()

public Bio::EnsEMBL::IdMapping::ExonScoreBuilder::exon_fasta_file ( )

Undocumented method

Code:

click to view

◆ exonerate_score()

public Bio::EnsEMBL::IdMapping::ExonScoreBuilder::exonerate_score ( )

Undocumented method

Code:

click to view

◆ non_mapped_transcript_rescore()

public Bio::EnsEMBL::IdMapping::ExonScoreBuilder::non_mapped_transcript_rescore ( )

Undocumented method

Code:

click to view

◆ overlap_score()

public Bio::EnsEMBL::IdMapping::ExonScoreBuilder::overlap_score ( )

Undocumented method

Code:

click to view

◆ parse_exonerate_results()

public Bio::EnsEMBL::IdMapping::ExonScoreBuilder::parse_exonerate_results ( )

Undocumented method

Code:

click to view

sub parse_exonerate_results {
  my ( $self, $exonerate_matrix ) = @_;
 
  unless ( defined($exonerate_matrix)
           &&
           $exonerate_matrix->isa(
                         'Bio::EnsEMBL::IdMapping::ScoredMappingMatrix')
    )
  {
    throw('You must provide a ScoredMappingMatrix.');
  }
 
  $self->logger->info( "Parsing exonerate results...\n", 0, 'stamped' );
 
  # loop over all result files
  my $dump_path = $self->cache->dump_path;
  my $num_files = 0;
  my $num_lines = 0;
 
  opendir( my $dumpdir, $dump_path ) or
    $self->logger->error("Can't open $dump_path for reading: $!");
 
  my $penalised = 0;
  my $killed    = 0;
 
  while ( defined( my $file = readdir($dumpdir) ) ) {
    unless ( $file =~ /exonerate_map\.\d+/ ) { next }
 
    $num_files++;
 
    open( my $fh, '<', "$dump_path/$file" ); 
 
    my $threshold = $self->conf->param('exonerate_threshold') || 0.5;
 
    while (<$fh>) {
      $num_lines++;
      chomp;
 
  # line format:
  # myinfo: source_id target_id match_length source_length target_length
      my ( undef, $source_id, $target_id, $match_length, $source_length,
           $target_length )
        = split;
 
      my $score = 0;
 
      if ( $source_length == 0 or $target_length == 0 ) {
        $self->logger->warning(
               "Alignment length is 0 for $source_id or $target_id.\n");
      }
      else {
        $score = 2*$match_length/( $source_length + $target_length );
 
      }
 
      if ( $score > $threshold ) {
        my $source_sr =
          $self->cache()
          ->get_by_key( 'exons_by_id', 'source', $source_id )
          ->seq_region_name();
        my $target_sr =
          $self->cache()
          ->get_by_key( 'exons_by_id', 'target', $target_id )
          ->seq_region_name();
 
        if ( $source_sr ne $target_sr ) {
          # PENALTY: The target and source are not on the same
          # seq_region.
          $score *= 0.70;
 
          # With a penalty of 0.7, exonerate scores need to be above
          # approximately 0.714 to make the 0.5 threshold.
 
          ++$penalised;
        }
 
        if ( $score > $threshold ) {
          $exonerate_matrix->add_score( $source_id, $target_id,
                                        $score );
        }
        else {
          ++$killed;
        }
      } ## end if ( $score > $threshold)
 
    } ## end while (<$fh>)
 
    close($fh);
  } ## end while ( defined( my $file...))
 
  closedir($dumpdir);
 
  $self->logger->info(
        "Done parsing $num_lines lines from $num_files result files.\n",
        0, 'stamped' );
  $self->logger->info( "Penalised $penalised exon alignments " .
                         "for not being on the same seq_region " .
                         "($killed killed).\n",
                       0,
                       'stamped' );
 
  return $exonerate_matrix;
} ## end sub parse_exonerate_results

◆ run_exonerate()

public Bio::EnsEMBL::IdMapping::ExonScoreBuilder::run_exonerate ( )

Undocumented method

Code:

click to view

sub run_exonerate {
  my $self = shift;
 
  my $source_file = $self->exon_fasta_file('source');
  my $target_file = $self->exon_fasta_file('target');
  my $source_size = -s $source_file;
  my $target_size = -s $target_file;
 
  # check if fasta files exist and are not empty
  unless ($source_size and $target_size) {
    throw("Can't find exon fasta files.");
  }
 
  # create an empty lsf log directory
  my $logpath = path_append($self->logger->logpath, 'exonerate');
  system("rm -rf $logpath") == 0 or
    $self->logger->error("Unable to delete lsf log dir $logpath: $!\n");
  system("mkdir -p $logpath") == 0 or
    $self->logger->error("Can't create lsf log dir $logpath: $!\n");
 
  # delete exonerate output from previous runs
  my $dump_path = $self->cache->dump_path;
 
  opendir(my $dumpdir, $dump_path) or
    $self->logger->error("Can't open $dump_path for reading: $!");
 
  while (defined(my $file = readdir($dumpdir))) {
    next unless /exonerate_map\.\d+/;
 
    unlink("$dump_path/$file") or
      $self->logger->error("Can't delete $dump_path/$file: $!");
  }
  
  closedir($dumpdir);
 
  # determine number of jobs to split task into
  my $bytes_per_job = $self->conf->param('exonerate_bytes_per_job')
    || 250000;
  my $num_jobs = $self->conf->param('exonerate_jobs');
  $num_jobs ||= int( $source_size/$bytes_per_job + 1 );
 
  my $percent =
    int( ( $self->conf->param('exonerate_threshold') || 0.5 )*100 );
  my $lsf_name       = 'idmapping_exonerate_' . time;
  my $exonerate_path = $self->conf->param('exonerate_path');
  my $exonerate_extra_params =
    $self->conf->param('exonerate_extra_params');
 
  #
  # run exonerate jobs using lsf
  #
  my $exonerate_job =
    qq{$exonerate_path } .
    qq{--query $source_file --target $target_file } .
    q{--querychunkid $LSB_JOBINDEX } .
    qq{--querychunktotal $num_jobs } .
    q{--model ungapped -M 1000 -D 100 } .
    q{--showalignment FALSE --subopt no } . qq{--percent $percent } .
    $self->conf->param('exonerate_extra_params') . " " .
    q{--ryo 'myinfo: %qi %ti %et %ql %tl\n' } .
    qq{| grep '^myinfo:' > $dump_path/exonerate_map.\$LSB_JOBINDEX} .
    "\n";
 
  $self->logger->info("Submitting $num_jobs exonerate jobs to lsf.\n");
  $self->logger->debug("$exonerate_job\n\n");
 
  my $bsub_cmd = sprintf(
               "|bsub -J '%s[1-%d]%%%d' -o %s/exonerate.%%I.out %s",
               $lsf_name,
               $num_jobs,
               $self->conf()->param('exonerate_concurrent_jobs') || 200,
               $logpath,
               $self->conf()->param('lsf_opt_exonerate') );
 
  local *BSUB;
  open( BSUB, $bsub_cmd ) ## no critic
    or $self->logger->error("Could not open open pipe to bsub: $!\n");
 
  print BSUB $exonerate_job;
  $self->logger->error("Error submitting exonerate jobs: $!\n")
    unless ($? == 0); 
  close BSUB;
 
  # submit dependent job to monitor finishing of exonerate jobs
  $self->logger->info("Waiting for exonerate jobs to finish...\n", 0, 'stamped');
 
  my $dependent_job =
    qq{bsub -K -w "ended($lsf_name)" -q production } .
    qq{-M 1000 -R 'select[mem>1000]' -R 'rusage[mem=1000]' } .
    qq{-o $logpath/exonerate_depend.out /bin/true};
 
  system($dependent_job) == 0 or
    $self->logger->error("Error submitting dependent job: $!\n");
 
  $self->logger->info("All exonerate jobs finished.\n", 0, 'stamped');
 
  #
  # check results
  #
  my @missing;
  my @error;
  
  for (my $i = 1; $i <= $num_jobs; $i++) {
  
    # check that output file exists
    my $outfile = "$dump_path/exonerate_map.$i";
    push @missing, $outfile unless (-f "$outfile");
 
    # check no errors occurred
    my $errfile = "$logpath/exonerate.$i.err";
    push @error, $errfile if (-s "$errfile");
  }
 
  if (@missing) {
    $self->logger->info("Couldn't find all exonerate output files. These are missing:\n");
    foreach (@missing) {
      $self->logger->info("$_\n", 1);
    }
 
    exit(1);
  }
 
  if (@error) {
    $self->logger->info("One or more exonerate jobs failed. Check these error files:\n");
    foreach (@error) {
      $self->logger->info("$_\n", 1);
    }
 
    exit(1);
  }
 
}

◆ score_exons()

public Bio::EnsEMBL::IdMapping::ExonScoreBuilder::score_exons ( )

Undocumented method

Code:

click to view

◆ sort_exons()

public Bio::EnsEMBL::IdMapping::ExonScoreBuilder::sort_exons ( )

Undocumented method

Code:

click to view

◆ write_filtered_exons()

public Bio::EnsEMBL::IdMapping::ExonScoreBuilder::write_filtered_exons ( )

Undocumented method

Code:

click to view

sub write_filtered_exons {
  my $self = shift;
  my $type = shift;
  my $matrix = shift;
 
  throw("You must provide a type.") unless $type;
  unless ($matrix and
          $matrix->isa('Bio::EnsEMBL::IdMapping::ScoredMappingMatrix')) {
    throw('You must provide a ScoredMappingMatrix.');
  }
 
  $self->logger->info("\nDumping $type exons to fasta file...\n", 0, 'stamped');
 
  # don't dump exons shorter than this
  my $min_exon_length = $self->conf->param('min_exon_length') || 15;
 
  # counters
  my $total_exons = 0;
  my $dumped_exons = 0;
 
  # filehandle for fasta files
  my $fh;
  my $file = $self->exon_fasta_file($type);
  open($fh, '>', $file) or throw("Unable to open $file for writing: $!");
 
  # loop over exons, dump sequence to fasta file if longer than threshold and
  # score < 1
  EXON:
  foreach my $eid (sort { $b <=> $a }
                   keys %{ $self->cache->get_by_name('exons_by_id', $type) }) {
 
    my $exon = $self->cache->get_by_key('exons_by_id', $type, $eid);
 
    $total_exons++;
 
    # skip if exon shorter than threshold
    next EXON if ($exon->length < $min_exon_length);
 
    # skip if overlap score with any other exon is 1
    if ( $type eq 'source' ) {
      foreach my $target ( @{ $matrix->get_targets_for_source($eid) } )
      {
        if ( $matrix->get_score( $eid, $target ) > 0.9999 ) {
          next EXON;
        }
      }
    } else {
      foreach my $source ( @{ $matrix->get_sources_for_target($eid) } )
      {
        if ( $matrix->get_score( $source, $eid ) > 0.9999 ) {
          next EXON;
        }
      }
    }
 
    # write exon to fasta file
    print $fh '>', $eid, "\n", $exon->seq, "\n";
 
    $dumped_exons++;
 
  }
 
  close($fh);
 
  # log
  my $fmt = "%-30s%10s\n";
  my $size = -s $file;
  $self->logger->info(sprintf($fmt, 'Total exons:', $total_exons), 1);
  $self->logger->info(sprintf($fmt, 'Dumped exons:', $dumped_exons), 1);
  $self->logger->info(sprintf($fmt, 'Dump file size:', parse_bytes($size)), 1);
  $self->logger->info("Done.\n\n", 0, 'stamped');
 
  return $dumped_exons;
}

The documentation for this class was generated from the following file:

ensembl/modules/Bio/EnsEMBL/IdMapping/ExonScoreBuilder.pm

Public Member Functions

Detailed Description

Synopsis

Description

Member Function Documentation

◆ build_overlap_scores()

◆ calc_overlap_score()

◆ compare_exon_containers()

◆ dump_filtered_exons()

◆ exon_fasta_file()

◆ exonerate_score()

◆ non_mapped_transcript_rescore()

◆ overlap_score()

◆ parse_exonerate_results()

◆ run_exonerate()

◆ score_exons()

◆ sort_exons()

◆ write_filtered_exons()