<commit_before>#!/usr/local/ensembl/bin/perl -w

=head1 Synopsis

load_mitochondria.pl 

=head1 Description

Parses mitochondrial genes out of genbank format NC_ file  and writes them to the database specified.
Can also load the sequence of the chromosome, the corresponding assembly entries and the appropriate attribute types if needed.

=head1 Config

All configuration is done through MitConf.pm

=cut

use strict;
use Carp;

use Bio::EnsEMBL::DBSQL::DBAdaptor;
use Bio::EnsEMBL::Exon;
use Bio::EnsEMBL::Transcript;
use Bio::EnsEMBL::Translation;
use Bio::EnsEMBL::Gene;
use Bio::EnsEMBL::DBEntry;
use Bio::EnsEMBL::Analysis;
use Bio::SeqIO;
use MitConf;
use Getopt::Long;
use Bio::EnsEMBL::Utils::Exception qw(stack_trace);
my $help;
my @genes;;
my $translation; # JUST INCASE
my $count;

GetOptions(
           '-h|help'    => \$help,
           ) or &usage();

unless ($MIT_DBHOST && $MIT_DBUSER && $MIT_DBNAME && $MIT_GENBANK_FILE && !$help){
  warn("Can't run without MitConf.pm values:
MIT_DBHOST $MIT_DBHOST 
MIT_DBUSER $MIT_DBUSER 
MIT_DBNAME $MIT_DBNAME
MIT_DBPASS $MIT_DBPASS
MIT_GENBANK_FILE $MIT_GENBANK_FILE
MIT_LOGIC_NAME $MIT_LOGIC_NAME
MIT_NAME $MIT_NAME
MIT_TOPLEVEL $MIT_TOPLEVEL
MIT_CODON_TABLE  $MIT_CODON_TABLE
MIT_GENE_TYPE $MIT_GENE_TYPE
MIT_TRNA_TYPE $MIT_TRNA_TYPE
MIT_RRNA_TYPE $MIT_RRNA_TYPE
");
  $help = 1;
}

if ($help) {
    exec('perldoc', $0);
}
############################
#PARSE GENBANK FILE FIRST

my ($genbank_ref,$assembly_ref)  = &_parse_coordinates;
my @genbank = @{$genbank_ref};
my %assembly = %{$assembly_ref};

##########################
# Get chromosome sequence 

my $seq_file = Bio::SeqIO->new(
			       -file   => $MIT_GENBANK_FILE,
			       -format => 'genbank',
			      );
my $chromosome_seq = $seq_file->next_seq;

####################
#Open db for writing

my $output_db = new Bio::EnsEMBL::DBSQL::DBAdaptor(
						   '-host'   => $MIT_DBHOST,
						   '-user'   => $MIT_DBUSER,
						   '-pass'   => $MIT_DBPASS,
						   '-dbname' => $MIT_DBNAME,
						   '-port'   => $MIT_DBPORT,	
						  );

########################
# Check taxon is correct

my $meta_container = $output_db->get_MetaContainer();
my $db_taxon_id = $meta_container->get_taxonomy_id();
my $gb_taxon_id = $genbank[1]{'db_xref'};
$gb_taxon_id =~ s/^taxon://;
if ($db_taxon_id  == $gb_taxon_id){
print "Genbank taxon id $gb_taxon_id matches the taxon id in the meta table\n";
}else {
  croak "taxonomy id in the genbank file $gb_taxon_id is not equal to the id found in the meta table $db_taxon_id \n";
}

###########################################
# write chromsome if it is not already there

my $dbe_adaptor = $output_db->get_DBEntryAdaptor;
my $slice_adaptor = $output_db->get_SliceAdaptor;
my $slice = $slice_adaptor->fetch_by_region('toplevel',$MIT_NAME);
my $slices_ref =  &get_chromosomes($assembly_ref,$genbank_ref, $output_db);
my %slices = %{$slices_ref};

if ($slice){
  print "Found chromosome ".$slice->name."\n" ;
}  else {
  print "There is no chromosome in $MIT_DBNAME called $MIT_NAME\nHave found the following slices to load:\n";
  foreach my $cs (keys %slices){
    print "$cs    \t".$slices{$cs}->name."\n";
  }
  print "Do you want to load the chromosome and the associated assembly entries
from the genbank file?(Y/N) "; 
  my $answer = <>;
  chomp $answer;
  if ($answer eq "y" or $answer eq "Y"){
    &load_chromosomes(\%slices, $output_db, $chromosome_seq );
    $slice = $slice_adaptor->fetch_by_region('toplevel',$MIT_NAME);
  }
  else {
    print "Ok gene load aborted.\n";
    exit 0;
  }
}

#########################################
#Fetch analysis object

my $logic_name = 'ncRNA';
my $ncRNA_analysis = $output_db->get_AnalysisAdaptor->fetch_by_logic_name($logic_name);
if(!defined $ncRNA_analysis){
#croak "analysis $logic_name not found\n";
  $ncRNA_analysis = Bio::EnsEMBL::Analysis->new(
                                                -logic_name => 'ncRNA',
                                               );
  print "You have no ".$logic_name." defined creating new object\n";
}

if ($MIT_LOGIC_NAME){
  $logic_name = $MIT_LOGIC_NAME;
} else {
  $logic_name = 'ensembl';
  print "Cannot find MIT_LOGIC_NAME - using standard logic name ensembl\n" if $MIT_DEBUG;
}
my $ensembl_analysis = $output_db->get_AnalysisAdaptor->fetch_by_logic_name($logic_name);
if(!defined $ensembl_analysis){
#croak "analysis $logic_name not found\n";
  $ensembl_analysis = Bio::EnsEMBL::Analysis->new(
                                                  -logic_name => 'ensembl',
                                                 );
  print "You have no ".$logic_name." defined creating new object\n";
}

 #########################################################
 # Create gene objects from array of hashes
 # contining parsed coords
 # index 0 contains the accession 
 # index 1 contains source data (taxon and sequence file)
 # subsequent indecies contain all the other annotations

for (my $i=2; $i <= $#genbank; $i++){
  my $desc = $genbank[$i]{'product'};
  my $type = $genbank[$i]{'type'};

  my $acc;
  if ($genbank[$i]{'protein_id'}){
    $acc =  $genbank[$i]{'protein_id'};
  }
  else{
    $acc = $genbank[0]{'accession'}."  ".$genbank[$i]{'start'}."  ".$genbank[$i]{'end'};
  }

  my $strand =  $genbank[$i]{'strand'};
  if ($type eq 'tRNA'){
    $type = $MIT_TRNA_TYPE;
  }
  if ($type eq 'rRNA'){
    $type = $MIT_RRNA_TYPE;
  }
  foreach my $key (keys %{$genbank[$i]}){
    print "$key ".$genbank[$i]{$key}."\n" if $MIT_DEBUG;
  }
  print "\n***********************************************************\n" if $MIT_DEBUG;
  # By default analysis is set to ncRNA
  my $analysis = $ncRNA_analysis;
  my $transcript = new Bio::EnsEMBL::Transcript;
  my $start_exon;
  my $end_exon;
  my $total_exons = scalar(@{$genbank[$i]{'start'}}-1);
  # exons are stored in an array within the hash incase of multi exon genes

  for (my $exon_number=0; $exon_number <= $total_exons ; $exon_number ++){
    my $start = $genbank[$i]{'start'}[$exon_number];
    my $end   = $genbank[$i]{ 'end' }[$exon_number];

    ###########################################################
    # Dodgy hack to solve problem of frameshifts in pseudogenes

    if ($exon_number > 0 && $start <= $genbank[$i]{ 'end' }[$exon_number-1]){
      $start = $genbank[$i]{ 'end' }[$exon_number-1]+1;
    }

    #######
    # EXONS

    my $exon = new Bio::EnsEMBL::Exon;
    $exon->start($start);
    $exon->end($end);
    $exon->strand($strand);
    $exon->slice($slice);
    $exon->phase(0);
    $exon->end_phase(($exon->end - $exon->start + 1)%3);

    #############
    # TRANSCRIPTS

    $transcript->add_Exon($exon);
    $start_exon = $exon if($exon_number == 0);
    $end_exon = $exon if($exon_number == $total_exons);
  }

  $transcript->start_Exon($start_exon);
  $transcript->end_Exon  ( $end_exon );

  #############
  # TRANSLATION

  if ($type eq 'CDS'){
    $translation = new  Bio::EnsEMBL::Translation(
						 -START_EXON => $start_exon,
						 -END_EXON   => $end_exon,
						 -SEQ_START  => 1,
						 -SEQ_END    => $transcript->length,
						 );
    $transcript->translation($translation);
    $analysis = $ensembl_analysis;
    $type = $MIT_GENE_TYPE;
  }

  #############
  # Make  genes

  my $gene = new Bio::EnsEMBL::Gene;
  eval {
  $gene->type($type);
  $gene->analysis($analysis);
  $gene->add_Transcript($transcript);
  $count++;
};
  if ($@){<commit_msg>Gene status KNOWN added when loading gene into database<commit_after>#!/usr/local/ensembl/bin/perl -w

=head1 Synopsis

load_mitochondria.pl 

=head1 Description

Parses mitochondrial genes out of genbank format NC_ file  and writes them to the database specified.
Can also load the sequence of the chromosome, the corresponding assembly entries and the appropriate attribute types if needed.

=head1 Config

All configuration is done through MitConf.pm

=cut

use strict;
use Carp;

use Bio::EnsEMBL::DBSQL::DBAdaptor;
use Bio::EnsEMBL::Exon;
use Bio::EnsEMBL::Transcript;
use Bio::EnsEMBL::Translation;
use Bio::EnsEMBL::Gene;
use Bio::EnsEMBL::DBEntry;
use Bio::EnsEMBL::Analysis;
use Bio::SeqIO;
use MitConf;
use Getopt::Long;
use Bio::EnsEMBL::Utils::Exception qw(stack_trace);
my $help;
my @genes;;
my $translation; # JUST INCASE
my $count;

GetOptions(
           '-h|help'    => \$help,
           ) or &usage();

unless ($MIT_DBHOST && $MIT_DBUSER && $MIT_DBNAME && $MIT_GENBANK_FILE && !$help){
  warn("Can't run without MitConf.pm values:
MIT_DBHOST $MIT_DBHOST 
MIT_DBUSER $MIT_DBUSER 
MIT_DBNAME $MIT_DBNAME
MIT_DBPASS $MIT_DBPASS
MIT_GENBANK_FILE $MIT_GENBANK_FILE
MIT_LOGIC_NAME $MIT_LOGIC_NAME
MIT_NAME $MIT_NAME
MIT_TOPLEVEL $MIT_TOPLEVEL
MIT_CODON_TABLE  $MIT_CODON_TABLE
MIT_GENE_TYPE $MIT_GENE_TYPE
MIT_TRNA_TYPE $MIT_TRNA_TYPE
MIT_RRNA_TYPE $MIT_RRNA_TYPE
");
  $help = 1;
}

if ($help) {
    exec('perldoc', $0);
}
############################
#PARSE GENBANK FILE FIRST

my ($genbank_ref,$assembly_ref)  = &_parse_coordinates;
my @genbank = @{$genbank_ref};
my %assembly = %{$assembly_ref};

##########################
# Get chromosome sequence 

my $seq_file = Bio::SeqIO->new(
			       -file   => $MIT_GENBANK_FILE,
			       -format => 'genbank',
			      );
my $chromosome_seq = $seq_file->next_seq;

####################
#Open db for writing

my $output_db = new Bio::EnsEMBL::DBSQL::DBAdaptor(
						   '-host'   => $MIT_DBHOST,
						   '-user'   => $MIT_DBUSER,
						   '-pass'   => $<PASSWORD>,
						   '-dbname' => $MIT_DBNAME,
						   '-port'   => $MIT_DBPORT,	
						  );

########################
# Check taxon is correct

my $meta_container = $output_db->get_MetaContainer();
my $db_taxon_id = $meta_container->get_taxonomy_id();
my $gb_taxon_id = $genbank[1]{'db_xref'};
$gb_taxon_id =~ s/^taxon://;
if ($db_taxon_id  == $gb_taxon_id){
print "Genbank taxon id $gb_taxon_id matches the taxon id in the meta table\n";
}else {
  croak "taxonomy id in the genbank file $gb_taxon_id is not equal to the id found in the meta table $db_taxon_id \n";
}

###########################################
# write chromsome if it is not already there

my $dbe_adaptor = $output_db->get_DBEntryAdaptor;
my $slice_adaptor = $output_db->get_SliceAdaptor;
my $slice = $slice_adaptor->fetch_by_region('toplevel',$MIT_NAME);
my $slices_ref =  &get_chromosomes($assembly_ref,$genbank_ref, $output_db);
my %slices = %{$slices_ref};

if ($slice){
  print "Found chromosome ".$slice->name."\n" ;
}  else {
  print "There is no chromosome in $MIT_DBNAME called $MIT_NAME\nHave found the following slices to load:\n";
  foreach my $cs (keys %slices){
    print "$cs    \t".$slices{$cs}->name."\n";
  }
  print "Do you want to load the chromosome and the associated assembly entries
from the genbank file?(Y/N) "; 
  my $answer = <>;
  chomp $answer;
  if ($answer eq "y" or $answer eq "Y"){
    &load_chromosomes(\%slices, $output_db, $chromosome_seq );
    $slice = $slice_adaptor->fetch_by_region('toplevel',$MIT_NAME);
  }
  else {
    print "Ok gene load aborted.\n";
    exit 0;
  }
}

#########################################
#Fetch analysis object

my $logic_name = 'ncRNA';
my $ncRNA_analysis = $output_db->get_AnalysisAdaptor->fetch_by_logic_name($logic_name);
if(!defined $ncRNA_analysis){
#croak "analysis $logic_name not found\n";
  $ncRNA_analysis = Bio::EnsEMBL::Analysis->new(
                                                -logic_name => 'ncRNA',
                                               );
  print "You have no ".$logic_name." defined creating new object\n";
}

if ($MIT_LOGIC_NAME){
  $logic_name = $MIT_LOGIC_NAME;
} else {
  $logic_name = 'ensembl';
  print "Cannot find MIT_LOGIC_NAME - using standard logic name ensembl\n" if $MIT_DEBUG;
}
my $ensembl_analysis = $output_db->get_AnalysisAdaptor->fetch_by_logic_name($logic_name);
if(!defined $ensembl_analysis){
#croak "analysis $logic_name not found\n";
  $ensembl_analysis = Bio::EnsEMBL::Analysis->new(
                                                  -logic_name => 'ensembl',
                                                 );
  print "You have no ".$logic_name." defined creating new object\n";
}

 #########################################################
 # Create gene objects from array of hashes
 # contining parsed coords
 # index 0 contains the accession 
 # index 1 contains source data (taxon and sequence file)
 # subsequent indecies contain all the other annotations

for (my $i=2; $i <= $#genbank; $i++){
  my $desc = $genbank[$i]{'product'};
  my $type = $genbank[$i]{'type'};
  my $status = "KNOWN";

  my $acc;
  if ($genbank[$i]{'protein_id'}){
    $acc =  $genbank[$i]{'protein_id'};
  }
  else{
    $acc = $genbank[0]{'accession'}."  ".$genbank[$i]{'start'}."  ".$genbank[$i]{'end'};
  }

  my $strand =  $genbank[$i]{'strand'};
  if ($type eq 'tRNA'){
    $type = $MIT_TRNA_TYPE;
  }
  if ($type eq 'rRNA'){
    $type = $MIT_RRNA_TYPE;
  }
  foreach my $key (keys %{$genbank[$i]}){
    print "$key ".$genbank[$i]{$key}."\n" if $MIT_DEBUG;
  }
  print "\n***********************************************************\n" if $MIT_DEBUG;
  # By default analysis is set to ncRNA
  my $analysis = $ncRNA_analysis;
  my $transcript = new Bio::EnsEMBL::Transcript;
  my $start_exon;
  my $end_exon;
  my $total_exons = scalar(@{$genbank[$i]{'start'}}-1);
  # exons are stored in an array within the hash incase of multi exon genes

  for (my $exon_number=0; $exon_number <= $total_exons ; $exon_number ++){
    my $start = $genbank[$i]{'start'}[$exon_number];
    my $end   = $genbank[$i]{ 'end' }[$exon_number];

    ###########################################################
    # Dodgy hack to solve problem of frameshifts in pseudogenes

    if ($exon_number > 0 && $start <= $genbank[$i]{ 'end' }[$exon_number-1]){
      $start = $genbank[$i]{ 'end' }[$exon_number-1]+1;
    }

    #######
    # EXONS

    my $exon = new Bio::EnsEMBL::Exon;
    $exon->start($start);
    $exon->end($end);
    $exon->strand($strand);
    $exon->slice($slice);
    $exon->phase(0);
    $exon->end_phase(($exon->end - $exon->start + 1)%3);

    #############
    # TRANSCRIPTS

    $transcript->add_Exon($exon);
    $start_exon = $exon if($exon_number == 0);
    $end_exon = $exon if($exon_number == $total_exons);
  }

  $transcript->start_Exon($start_exon);
  $transcript->end_Exon  ( $end_exon );

  #############
  # TRANSLATION

  if ($type eq 'CDS'){
    $translation = new  Bio::EnsEMBL::Translation(
						 -START_EXON => $start_exon,
						 -END_EXON   => $end_exon,
						 -SEQ_START  => 1,
						 -SEQ_END    => $transcript->length,
						 );
    $transcript->translation($translation);
    $analysis = $ensembl_analysis;
    $type = $MIT_GENE_TYPE;
  }

  #############
  # Make  genes

  my $gene = new Bio::EnsEMBL::Gene;
  eval {
  $gene->type($type);
  $gene->analysis($analysis);
  $gene->status($status);
  $gene->description($desc);
  $gene->add_Transcript($transcript);
  $count++;
};
  if ($@){
{'fraction_non_alphanumeric': 0.1979872525997987, 'fraction_numerical': 0.002952029520295203, 'mean_word_length': 3.800644122383253, 'pattern_counts': {'":': 0, '<': 12, '<?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': 16, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '2557086', 'n_tokens_mistral': 5696, 'n_tokens_neox': 5276, 'n_words': 1231}