<commit_before>=head1 LICENSE

Copyright [1999-2015] Wellcome Trust Sanger Institute and the EMBL-European Bioinformatics Institute
Copyright [2016-2017] EMBL-European Bioinformatics Institute

Licensed under the Apache License, Version 2.0 (the "License");
you may not use this file except in compliance with the License.
You may obtain a copy of the License at

     http://www.apache.org/licenses/LICENSE-2.0

Unless required by applicable law or agreed to in writing, software
distributed under the License is distributed on an "AS IS" BASIS,
WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
See the License for the specific language governing permissions and
limitations under the License.

=head1 CONTACT

Please email comments or questions to the public Ensembl
developers list at <http://lists.ensembl.org/mailman/listinfo/dev>.

Questions may also be sent to the Ensembl help desk at
<http://www.ensembl.org/Help/Contact>.

=head1 NAME

Bio::EnsEMBL::Analysis::Hive::RunnableDB::HiveCesar2

=cut

=head1 DESCRIPTION

HiveCesar2 fetches the gene corresponding to gene_id in source_db, projects its exons based on the given Compara lastz alignment and the CESAR2.0 aligner and builds single-transcript genes from these projections to be written to target_db while filtering out the specified transcripts by applying the filter in TRANSCRIPT_FILTER.
seqEdits and selenocysteine attributes are inserted in order to deal with alignment gaps and seleno-like TGA stops which will be converted into NNN triplets to make the best possible aligment.

=head1 OPTIONS

-iid                  gene_id from the source_db corresponding to the gene to be projected from the source_dna_db to the target_dna_db.
-output_path          Path where the output files will be stored.
-source_dna_db        Ensembl database containing the DNA sequences that correspond to the input gene_id from the source_db.
-target_dna_db        Ensembl database containing the DNA sequences corresponding to the target_db species where the gene_id will be projected to.
-source_db            Ensembl database containing the genes whose transcripts will be projected to the target species db target_db.
-target_db            Ensembl database containing the DNA sequences corresponding to the target species the input gene_id will be projected to.
-compara_db           Compara database containing the lastz alignments of the source and target species.
-method_link_type     Default to 'LASTZ_NET' so it works with the Compara lastz alignments.
-exon_region_padding  Default to 50 to add 50 bases at each side of the exons to be projected.
-cesar_path           Path to the directory containing the CESAR2.0 binary to be run (excluding the binary filename).
-TRANSCRIPT_FILTER    Hash containing the parameters required to apply to the projected transcript to exclude some of them. Default to ExonerateTranscriptFilter pid,cov 50,50 although note that the actual implementation of this filter allows pid,cov below 50,50 in some cases.

=cut

package Bio::EnsEMBL::Analysis::Hive::RunnableDB::HiveCesar2;

use warnings;
use strict;
use feature 'say';
use Data::Dumper;
#use File::chdir;

use Bio::EnsEMBL::Analysis::Tools::GeneBuildUtils::GeneUtils qw(empty_Gene);

use Bio::SeqIO;
use Bio::EnsEMBL::DBSQL::DBAdaptor;
use Bio::EnsEMBL::Utils::Exception qw(throw warning);
use Bio::EnsEMBL::Compara::DBSQL::DBAdaptor;
use Bio::EnsEMBL::Analysis::Tools::WGA2Genes::GeneScaffold;
use Bio::EnsEMBL::Analysis::Tools::ClusterFilter;
use Bio::EnsEMBL::Analysis::Tools::GeneBuildUtils::TranscriptUtils qw(replace_stops_with_introns
                                                                      calculate_exon_phases
                                                                      set_alignment_supporting_features
                                                                      features_overlap);
use Bio::EnsEMBL::Analysis::Tools::Utilities qw(align_proteins);
use parent ('Bio::EnsEMBL::Analysis::Hive::RunnableDB::HiveBaseRunnableDB');

sub param_defaults {
    return {
      iid => '',
      output_path => '',
      source_dna_db => '',
      target_dna_db => '',
      source_db => '',
      target_db => '',
      compara_db => '',
      method_link_type => 'LASTZ_NET',
      exon_region_padding => 50,
      cesar_path => '',
      TRANSCRIPT_FILTER => {
                             OBJECT     => 'Bio::EnsEMBL::Analysis::Tools::ExonerateTranscriptFilter',
                             PARAMETERS => {
                               -coverage => 50,
                               -percent_id => 50,
                             },
                           }
   }
}

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

  unless(-e $self->param('output_path')) {
    system("mkdir -p ".$self->param('output_path'));
  }

  # if CESAR2.0 fails this will get -1 value so the module does not write anything and exits.
  #$self->param('himem_required') = 0;

  my $input_id = $self->param('iid');

  # Define the dna dbs
  my $source_dna_dba = $self->hrdb_get_dba($self->param('source_dna_db'));
  my $target_dna_dba = $self->hrdb_get_dba($self->param('target_dna_db'));
  $self->hrdb_set_con($source_dna_dba,'source_dna_db');
  $self->hrdb_set_con($target_dna_dba,'target_dna_db');

  # Define the source transcript and target transcript dbs
  my $source_transcript_dba = $self->hrdb_get_dba($self->param('source_db'));
  my $target_transcript_dba = $self->hrdb_get_dba($self->param('target_db'));
  $self->hrdb_set_con($source_transcript_dba,'source_transcript_db');
  $self->hrdb_set_con($target_transcript_dba,'target_transcript_db');

  # Define the compara db
  my $compara_dba = $self->hrdb_get_dba($self->param('compara_db'),undef,'Compara');
  $self->hrdb_set_con($compara_dba,'compara_db');

  # Get the genome db adpator
  my $genome_dba = $compara_dba->get_GenomeDBAdaptor;

  # Retrieve the production names for the query and target species
  my $source_species = $source_transcript_dba->get_MetaContainerAdaptor->get_production_name();
  my $target_species = $target_transcript_dba->get_MetaContainerAdaptor->get_production_name();

  my $source_genome_db = $genome_dba->fetch_by_core_DBAdaptor($source_transcript_dba);
  my $target_genome_db = $genome_dba->fetch_by_core_DBAdaptor($target_transcript_dba);

  ########
  # check that the default assembly for the query and target agrees
  # with that for the method_link_species_set GenomeDBs
  ########

  my $source_assembly = $source_genome_db->assembly;
  my $target_assembly = $target_genome_db->assembly;

  my ($source_assembly_version, $target_assembly_version);
  eval {
    $source_assembly_version = $source_transcript_dba->get_CoordSystemAdaptor->fetch_by_name('toplevel',$source_genome_db->assembly);
    $target_assembly_version = $target_transcript_dba->get_CoordSystemAdaptor->fetch_by_name('toplevel',$target_genome_db->assembly);
  };
  if($@) {
    $self->throw("Had trouble fetching coord systems for ".$source_genome_db->assembly . " and " .$target_genome_db->assembly . " from core dbs:\n".$@);
  }

  my $gene = $source_transcript_dba->get_GeneAdaptor->fetch_by_dbID($input_id);
  $self->parent_gene($gene);

  my @unique_translateable_exons = $self->get_unique_translateable_exons($gene);
  my $exon_align_slices;

  #########
  # get the compara data: MethodLinkSpeciesSet, reference DnaFrag,
  # and all GenomicAlignBlocks
  #########
  my $mlss = $compara_dba->get_MethodLinkSpeciesSetAdaptor->fetch_by_method_link_type_GenomeDBs($self->param('method_link_type'),
                                                                                                [$source_genome_db,
                                                                                                 $target_genome_db]);
  unless($mlss) {
    $self->throw("No MethodLinkSpeciesSet for :\n" .$self->param('method_link_type') . "\n" .$source_species . "\n" .$target_species);
  }

  my $dna_fragment = $compara_dba->get_DnaFragAdaptor->fetch_by_GenomeDB_and_name($source_genome_db,$gene->slice->seq_region_name);
  my $genomic_align_block_adaptor = $compara_dba->get_GenomicAlignBlockAdaptor;
  my $exon_region_padding = $self->param('exon_region_padding');
  foreach my $exon (@unique_translateable_exons) {
    my $exon_padded_start = $exon->start - $exon_region_padding;
    if($exon_padded_start < 0) {
      $exon_padded_start = 0;
    }
<commit_msg>Removed unused dna_fragment variable.<commit_after>=head1 LICENSE

Copyright [1999-2015] Wellcome Trust Sanger Institute and the EMBL-European Bioinformatics Institute
Copyright [2016-2017] EMBL-European Bioinformatics Institute

Licensed under the Apache License, Version 2.0 (the "License");
you may not use this file except in compliance with the License.
You may obtain a copy of the License at

     http://www.apache.org/licenses/LICENSE-2.0

Unless required by applicable law or agreed to in writing, software
distributed under the License is distributed on an "AS IS" BASIS,
WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
See the License for the specific language governing permissions and
limitations under the License.

=head1 CONTACT

Please email comments or questions to the public Ensembl
developers list at <http://lists.ensembl.org/mailman/listinfo/dev>.

Questions may also be sent to the Ensembl help desk at
<http://www.ensembl.org/Help/Contact>.

=head1 NAME

Bio::EnsEMBL::Analysis::Hive::RunnableDB::HiveCesar2

=cut

=head1 DESCRIPTION

HiveCesar2 fetches the gene corresponding to gene_id in source_db, projects its exons based on the given Compara lastz alignment and the CESAR2.0 aligner and builds single-transcript genes from these projections to be written to target_db while filtering out the specified transcripts by applying the filter in TRANSCRIPT_FILTER.
seqEdits and selenocysteine attributes are inserted in order to deal with alignment gaps and seleno-like TGA stops which will be converted into NNN triplets to make the best possible aligment.

=head1 OPTIONS

-iid                  gene_id from the source_db corresponding to the gene to be projected from the source_dna_db to the target_dna_db.
-output_path          Path where the output files will be stored.
-source_dna_db        Ensembl database containing the DNA sequences that correspond to the input gene_id from the source_db.
-target_dna_db        Ensembl database containing the DNA sequences corresponding to the target_db species where the gene_id will be projected to.
-source_db            Ensembl database containing the genes whose transcripts will be projected to the target species db target_db.
-target_db            Ensembl database containing the DNA sequences corresponding to the target species the input gene_id will be projected to.
-compara_db           Compara database containing the lastz alignments of the source and target species.
-method_link_type     Default to 'LASTZ_NET' so it works with the Compara lastz alignments.
-exon_region_padding  Default to 50 to add 50 bases at each side of the exons to be projected.
-cesar_path           Path to the directory containing the CESAR2.0 binary to be run (excluding the binary filename).
-TRANSCRIPT_FILTER    Hash containing the parameters required to apply to the projected transcript to exclude some of them. Default to ExonerateTranscriptFilter pid,cov 50,50 although note that the actual implementation of this filter allows pid,cov below 50,50 in some cases.

=cut

package Bio::EnsEMBL::Analysis::Hive::RunnableDB::HiveCesar2;

use warnings;
use strict;
use feature 'say';
use Data::Dumper;
#use File::chdir;

use Bio::EnsEMBL::Analysis::Tools::GeneBuildUtils::GeneUtils qw(empty_Gene);

use Bio::SeqIO;
use Bio::EnsEMBL::DBSQL::DBAdaptor;
use Bio::EnsEMBL::Utils::Exception qw(throw warning);
use Bio::EnsEMBL::Compara::DBSQL::DBAdaptor;
use Bio::EnsEMBL::Analysis::Tools::WGA2Genes::GeneScaffold;
use Bio::EnsEMBL::Analysis::Tools::ClusterFilter;
use Bio::EnsEMBL::Analysis::Tools::GeneBuildUtils::TranscriptUtils qw(replace_stops_with_introns
                                                                      calculate_exon_phases
                                                                      set_alignment_supporting_features
                                                                      features_overlap);
use Bio::EnsEMBL::Analysis::Tools::Utilities qw(align_proteins);
use parent ('Bio::EnsEMBL::Analysis::Hive::RunnableDB::HiveBaseRunnableDB');

sub param_defaults {
    return {
      iid => '',
      output_path => '',
      source_dna_db => '',
      target_dna_db => '',
      source_db => '',
      target_db => '',
      compara_db => '',
      method_link_type => 'LASTZ_NET',
      exon_region_padding => 50,
      cesar_path => '',
      TRANSCRIPT_FILTER => {
                             OBJECT     => 'Bio::EnsEMBL::Analysis::Tools::ExonerateTranscriptFilter',
                             PARAMETERS => {
                               -coverage => 50,
                               -percent_id => 50,
                             },
                           }
   }
}

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

  unless(-e $self->param('output_path')) {
    system("mkdir -p ".$self->param('output_path'));
  }

  # if CESAR2.0 fails this will get -1 value so the module does not write anything and exits.
  #$self->param('himem_required') = 0;

  my $input_id = $self->param('iid');

  # Define the dna dbs
  my $source_dna_dba = $self->hrdb_get_dba($self->param('source_dna_db'));
  my $target_dna_dba = $self->hrdb_get_dba($self->param('target_dna_db'));
  $self->hrdb_set_con($source_dna_dba,'source_dna_db');
  $self->hrdb_set_con($target_dna_dba,'target_dna_db');

  # Define the source transcript and target transcript dbs
  my $source_transcript_dba = $self->hrdb_get_dba($self->param('source_db'));
  my $target_transcript_dba = $self->hrdb_get_dba($self->param('target_db'));
  $self->hrdb_set_con($source_transcript_dba,'source_transcript_db');
  $self->hrdb_set_con($target_transcript_dba,'target_transcript_db');

  # Define the compara db
  my $compara_dba = $self->hrdb_get_dba($self->param('compara_db'),undef,'Compara');
  $self->hrdb_set_con($compara_dba,'compara_db');

  # Get the genome db adpator
  my $genome_dba = $compara_dba->get_GenomeDBAdaptor;

  # Retrieve the production names for the query and target species
  my $source_species = $source_transcript_dba->get_MetaContainerAdaptor->get_production_name();
  my $target_species = $target_transcript_dba->get_MetaContainerAdaptor->get_production_name();

  my $source_genome_db = $genome_dba->fetch_by_core_DBAdaptor($source_transcript_dba);
  my $target_genome_db = $genome_dba->fetch_by_core_DBAdaptor($target_transcript_dba);

  ########
  # check that the default assembly for the query and target agrees
  # with that for the method_link_species_set GenomeDBs
  ########

  my $source_assembly = $source_genome_db->assembly;
  my $target_assembly = $target_genome_db->assembly;

  my ($source_assembly_version, $target_assembly_version);
  eval {
    $source_assembly_version = $source_transcript_dba->get_CoordSystemAdaptor->fetch_by_name('toplevel',$source_genome_db->assembly);
    $target_assembly_version = $target_transcript_dba->get_CoordSystemAdaptor->fetch_by_name('toplevel',$target_genome_db->assembly);
  };
  if($@) {
    $self->throw("Had trouble fetching coord systems for ".$source_genome_db->assembly . " and " .$target_genome_db->assembly . " from core dbs:\n".$@);
  }

  my $gene = $source_transcript_dba->get_GeneAdaptor->fetch_by_dbID($input_id);
  $self->parent_gene($gene);

  my @unique_translateable_exons = $self->get_unique_translateable_exons($gene);
  my $exon_align_slices;

  #########
  # get the compara data: MethodLinkSpeciesSet, reference DnaFrag,
  # and all GenomicAlignBlocks
  #########
  my $mlss = $compara_dba->get_MethodLinkSpeciesSetAdaptor->fetch_by_method_link_type_GenomeDBs($self->param('method_link_type'),
                                                                                                [$source_genome_db,
                                                                                                 $target_genome_db]);
  unless($mlss) {
    $self->throw("No MethodLinkSpeciesSet for :\n" .$self->param('method_link_type') . "\n" .$source_species . "\n" .$target_species);
  }

  my $genomic_align_block_adaptor = $compara_dba->get_GenomicAlignBlockAdaptor;
  my $exon_region_padding = $self->param('exon_region_padding');
  foreach my $exon (@unique_translateable_exons) {
    my $exon_padded_start = $exon->start - $exon_region_padding;
    if($exon_padded_start < 0) {
      $exon_padded_start = 0;
    }

{'fraction_non_alphanumeric': 0.09709434413487983, 'fraction_numerical': 0.006815735979911515, 'mean_word_length': 3.599120153973055, 'pattern_counts': {'":': 0, '<': 9, '<?xml version=': 0, '>': 139, 'https://': 0, 'lorem ipsum': 0, 'www.': 4, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 34, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '5214740', 'n_tokens_mistral': 5266, 'n_tokens_neox': 4947, 'n_words': 1326}