#!/usr/bin/env perl ############################################################ # # $Id: get-ensembl-genome-only-objects.pl,v 1.15 2011/02/17 05:07:46 rsat Exp $ # # Time-stamp: <2003-07-04 12:48:55 jvanheld> # ############################################################ #use strict; use DBI(); ## use Devel::Peek; ## This requires a specific compilation of Perl BEGIN { if ($0 =~ /([^(\/)]+)$/) { push (@INC, "$`lib/"); } } require "RSA.lib"; require RSAT::util; ## parsing libraries require "RSA.lib"; push @INC, $ENV{RSAT}."/perl-scripts/parsers/" if ($ENV{RSAT}); require "lib/load_classes.pl"; require "lib/parsing_util.pl"; ## EnsEMBL libraries use Bio::EnsEMBL::DBSQL::DBAdaptor; use Bio::EnsEMBL::DBSQL::SliceAdaptor; ## TO DO ## Introns and exons as transcript attributes ## Get CDS and check start/stop codons ## Test with other genomes (Anopheles) ## Check the start and stop codons of Anopheles. Half of them are false. ## Add cross-references to the RSAT objects ($rsat_gene, $rsat_transcript, $rsat_cds, ...) ################################################################ #### Class for EMBL feature package EMBL::Feature; { @ISA = qw ( classes::DatabaseObject ); ### class attributes $_count = 0; $_prefix = "ft_"; @_objects = (); %_name_index = (); %_id_index = (); %_attribute_count = (); %_attribute_cardinality = (id=>"SCALAR", names=>"ARRAY", organism=>"SCALAR", type=>"SCALAR", description=>"SCALAR", position=>"SCALAR", contig=>"SCALAR", strand=>"SCALAR", start_pos=>"SCALAR", end_pos=>"SCALAR", xrefs=>"EXPANDED", gene=>"SCALAR" ); } ################################################################ #### Class for transcript package EMBL::Transcript; { @ISA = qw ( classes::DatabaseObject ); ### class attributes $_count = 0; $_prefix = "trans_"; @_objects = (); %_name_index = (); %_id_index = (); %_attribute_count = (); %_attribute_cardinality = (id=>"SCALAR", names=>"ARRAY", organism=>"SCALAR", type=>"SCALAR", description=>"SCALAR", position=>"SCALAR", contig=>"SCALAR", strand=>"SCALAR", start_pos=>"SCALAR", end_pos=>"SCALAR", xrefs=>"EXPANDED", gene=>"SCALAR" ); } ################################################################ #### Class for mRNA package EMBL::mRNA; { @ISA = qw ( classes::DatabaseObject ); ### class attributes $_count = 0; $_prefix = "mRNA_"; @_objects = (); %_name_index = (); %_id_index = (); %_attribute_count = (); %_attribute_cardinality = (id=>"SCALAR", names=>"ARRAY", organism=>"SCALAR", type=>"SCALAR", description=>"SCALAR", position=>"SCALAR", contig=>"SCALAR", strand=>"SCALAR", start_pos=>"SCALAR", end_pos=>"SCALAR", xrefs=>"EXPANDED", gene=>"SCALAR" ); } ################################################################ #### Class for scRNA package EMBL::scRNA; { @ISA = qw ( classes::DatabaseObject ); ### class attributes $_count = 0; $_prefix = "scRNA_"; @_objects = (); %_name_index = (); %_id_index = (); %_attribute_count = (); %_attribute_cardinality = (id=>"SCALAR", names=>"ARRAY", organism=>"SCALAR", type=>"SCALAR", description=>"SCALAR", position=>"SCALAR", contig=>"SCALAR", strand=>"SCALAR", start_pos=>"SCALAR", end_pos=>"SCALAR", xrefs=>"EXPANDED", gene=>"SCALAR" ); } ################################################################ #### Class for tRNA package EMBL::tRNA; { @ISA = qw ( classes::DatabaseObject ); ### class attributes $_count = 0; $_prefix = "tRNA_"; @_objects = (); %_name_index = (); %_id_index = (); %_attribute_count = (); %_attribute_cardinality = (id=>"SCALAR", names=>"ARRAY", organism=>"SCALAR", type=>"SCALAR", description=>"SCALAR", position=>"SCALAR", contig=>"SCALAR", strand=>"SCALAR", start_pos=>"SCALAR", end_pos=>"SCALAR", xrefs=>"EXPANDED", gene=>"SCALAR" ); } ################################################################ #### Class for rRNA package EMBL::rRNA; { @ISA = qw ( classes::DatabaseObject ); ### class attributes $_count = 0; $_prefix = "rRNA_"; @_objects = (); %_name_index = (); %_id_index = (); %_attribute_count = (); %_attribute_cardinality = (id=>"SCALAR", names=>"ARRAY", organism=>"SCALAR", type=>"SCALAR", description=>"SCALAR", position=>"SCALAR", contig=>"SCALAR", strand=>"SCALAR", start_pos=>"SCALAR", end_pos=>"SCALAR", xrefs=>"EXPANDED", gene=>"SCALAR" ); } ################################################################ #### Class for miscellaneous RNA package EMBL::misc_RNA; { @ISA = qw ( classes::DatabaseObject ); ### class attributes $_count = 0; $_prefix = "misc_RNA_"; @_objects = (); %_name_index = (); %_id_index = (); %_attribute_count = (); %_attribute_cardinality = (id=>"SCALAR", names=>"ARRAY", organism=>"SCALAR", type=>"SCALAR", description=>"SCALAR", position=>"SCALAR", contig=>"SCALAR", strand=>"SCALAR", start_pos=>"SCALAR", end_pos=>"SCALAR", xrefs=>"EXPANDED", gene=>"SCALAR" ); } ################################################################ #### Class for repeat regions package EMBL::repeat_region; { @ISA = qw ( classes::DatabaseObject ); ### class attributes $_count = 0; $_prefix = "rep_"; @_objects = (); %_name_index = (); %_id_index = (); %_attribute_count = (); %_attribute_cardinality = (id=>"SCALAR", names=>"ARRAY", organism=>"SCALAR", type=>"SCALAR", description=>"SCALAR", position=>"SCALAR", contig=>"SCALAR", strand=>"SCALAR", start_pos=>"SCALAR", end_pos=>"SCALAR", xrefs=>"EXPANDED" ); } ################################################################ #### Class for CDS package EMBL::CDS; { @ISA = qw ( classes::DatabaseObject ); ### class attributes $_count = 0; $_prefix = "cds_"; @_objects = (); %_name_index = (); %_id_index = (); %_attribute_count = (); %_attribute_cardinality = (id=>"SCALAR", names=>"ARRAY", organism=>"SCALAR", type=>"SCALAR", description=>"SCALAR", position=>"SCALAR", contig=>"SCALAR", strand=>"SCALAR", start_pos=>"SCALAR", end_pos=>"SCALAR", xrefs=>"EXPANDED", transcript=>"SCALAR" ); } ################################################################ #### Class for sources package EMBL::Source; { @ISA = qw ( classes::DatabaseObject ); ### class attributes $_count = 0; $_prefix = "src_"; @_objects = (); %_name_index = (); %_id_index = (); %_attribute_count = (); %_attribute_cardinality = (id=>"SCALAR", names=>"ARRAY", organism=>"SCALAR", type=>"SCALAR", description=>"SCALAR", position=>"SCALAR", contig=>"SCALAR", strand=>"SCALAR", start_pos=>"SCALAR", end_pos=>"SCALAR", xrefs=>"EXPANDED" ); } ################################################################ #### Class for miscellaneous features package EMBL::misc_feature; { @ISA = qw ( classes::DatabaseObject ); ### class attributes $_count = 0; $_prefix = "misc_feature_"; @_objects = (); %_name_index = (); %_id_index = (); %_attribute_count = (); %_attribute_cardinality = (id=>"SCALAR", names=>"ARRAY", organism=>"SCALAR", type=>"SCALAR", description=>"SCALAR", position=>"SCALAR", contig=>"SCALAR", strand=>"SCALAR", start_pos=>"SCALAR", end_pos=>"SCALAR", xrefs=>"EXPANDED" ); } ################################################################ #### Class for Gene package EMBL::Gene; { @ISA = qw ( classes::DatabaseObject ); ### class attributes $_count = 0; $_prefix = "gn_"; @_objects = (); %_name_index = (); %_id_index = (); %_attribute_count = (); %_attribute_cardinality = (id=>"SCALAR", names=>"ARRAY", ); } ################################################################ #### Class for EMBL organism package EMBL::Organism; { @ISA = qw ( classes::DatabaseObject ); ### class attributes $_count = 0; $_prefix = "org_"; @_objects = (); %_name_index = (); %_id_index = (); %_attribute_count = (); %_attribute_cardinality = (id=>"SCALAR", names=>"ARRAY", taxonomy=>"SCALAR", ); } ################################################################ #### Class for EMBL contig package EMBL::Contig; { @ISA = qw ( classes::DatabaseObject ); ### class attributes $_count = 0; $_prefix = "ctg_"; @_objects = (); %_name_index = (); %_id_index = (); %_attribute_count = (); %_attribute_cardinality = (id=>"SCALAR", accession=>"SCALAR", version=>"SCALAR", type=>"SCALAR", length=>"SCALAR", description=>"SCALAR", ); } ################################################################ #### main package package main; { ################################################################ #### initialise parameters local $start_time = &RSAT::util::StartScript(); local $slice_type = "chromosome"; local @chromnames = (); local $no_seq = 0; local $no_masked = 0; local $test = 0; local $test_number=20; %dir = (); ### Working directory $dir{main} = `pwd`; #### remember working directory chomp($dir{main}); local %outfile = (); $outfile{stats}="get-ensembl-genome_stats.txt"; $outfile{log}="get-ensembl-genome_log.txt"; $outfile{err}="get-ensembl-genome_err.txt"; $outfile{contigs} = "contigs.txt"; local $verbose = 0; ## Connection to the EnsEMBL MYSQL database $ensembl_host = 'ensembldb.ensembl.org'; $ensembl_user = "anonymous"; $dbname = ''; $org = ''; #### Options for the exported SQL database $host= "localhost"; $schema="ensembl"; $user="ensembl"; $password="ensembl"; ################################################################ ## Class factories #### Organism factory my $organisms = classes::ClassFactory->new_class(object_type=>"EMBL::Organism",prefix=>"org_"); $organisms->set_out_fields(qw (id taxonomy name names )); #### Contig factory my $contigs = classes::ClassFactory->new_class(object_type=>"EMBL::Contig",prefix=>"ctg_"); $contigs->set_out_fields(qw (id accession version form type file length description )); ## Gene factory my $genes = classes::ClassFactory->new_class(object_type=>"EMBL::Gene",prefix=>"gene_"); $genes->set_out_fields(qw(id type name contig start_pos end_pos strand description names db_xref introns exons )); ## Transcript factory my $transcripts = classes::ClassFactory->new_class(object_type=>"EMBL::Transcript",prefix=>"trans_"); $transcripts->set_out_fields(qw( id type name contig start_pos end_pos strand description names db_xref introns exons gene )); ## mRNA factory my $mRNAs = classes::ClassFactory->new_class(object_type=>"EMBL::mRNA",prefix=>"mRNA_"); $mRNAs->set_out_fields(qw( id type name contig start_pos end_pos strand description names db_xref introns exons gene )); ## scRNA factory my $scRNAs = classes::ClassFactory->new_class(object_type=>"EMBL::scRNA",prefix=>"scRNA_"); $scRNAs->set_out_fields(qw( id type name contig start_pos end_pos strand description names db_xref introns exons gene )); ## tRNA factory my $tRNAs = classes::ClassFactory->new_class(object_type=>"EMBL::tRNA",prefix=>"tRNA_"); $tRNAs->set_out_fields(qw( id type name contig start_pos end_pos strand description names db_xref introns exons gene )); ## rRNA factory my $rRNAs = classes::ClassFactory->new_class(object_type=>"EMBL::rRNA",prefix=>"rRNA_"); $rRNAs->set_out_fields(qw( id type name contig start_pos end_pos strand description names db_xref introns exons gene )); ## Repeated regions factory my $repeat_regions = classes::ClassFactory->new_class(object_type=>"EMBL::repeat_region",prefix=>"rep_"); $repeat_regions->set_out_fields(qw( id type name contig start_pos end_pos strand )); ## Miscellaneous RNA factory my $misc_RNAs = classes::ClassFactory->new_class(object_type=>"EMBL::misc_RNA",prefix=>"misc_RNA_"); $misc_RNAs->set_out_fields(qw( id type name contig start_pos end_pos strand description names db_xref introns exons gene )); ## CDS factory my $CDSs = classes::ClassFactory->new_class(object_type=>"EMBL::CDS",prefix=>"cds_"); $CDSs->set_out_fields(qw( id type name contig start_pos end_pos strand description names db_xref introns exons transcript )); ## Miscellaneous features factory my $misc_features = classes::ClassFactory->new_class(object_type=>"EMBL::misc_feature",prefix=>"misc_feature_"); $misc_features->set_out_fields(qw( id type name contig start_pos end_pos strand description names db_xref introns exons )); ## Sources factory my $sources = classes::ClassFactory->new_class(object_type=>"EMBL::Source",prefix=>"src_"); $sources->set_out_fields(qw( id type name contig start_pos end_pos strand description names db_xref introns exons )); #### classes @classes = qw (EMBL::Organism EMBL::Contig EMBL::Gene EMBL::mRNA EMBL::scRNA EMBL::tRNA EMBL::rRNA EMBL::misc_RNA EMBL::repeat_region EMBL::CDS EMBL::Feature EMBL::misc_feature EMBL::Source); ################################################################ ## Read arguments &ReadArguments(); ################################################################ ## Check argument values # Either an organism or a database name must be provided unless ($org || $dbname) { die "; You must provide either an organism name (-org) or a database name (-dbname)\n"; } # Database name has priority over organism if both provided if ($org && $dbname) { $org = ''; } ### verbose ### if ($verbose >= 1) { print "; get-ensembl-genome "; &PrintArguments(); } ################################################################ ## If option -org is used, connect to ensembldb to get list of ## databases and pick the latest one corresponding to chosen organism if ($org) { &RSAT::message::TimeWarn (join("\t", "Connecting EnsEMBL to get the dbname for organism ", $org, "host=".$ensembl_host, "user=".$ensembl_user )) if ($main::verbose >= 1); my $dbh = DBI->connect("DBI:mysql:host=$ensembl_host", "$ensembl_user", "", {'RaiseError' => 1}); my $sth = $dbh->prepare("SHOW DATABASES"); $sth->execute(); while (my $ref = $sth->fetchrow_hashref()) { if ($ref->{'Database'} =~ /($org)_core_\d+/) { $dbname = $ref->{'Database'}; } } &RSAT::message::Info (join("\t", "dbname = ", $dbname)) if ($main::verbose >= 1); $sth->finish(); $dbh->disconnect(); } ################################################################ ### open output streams unless ($dir{output}) { $dir{output} = join("_", "EnsEMBL", $dbname); if (scalar(@chromnames) > 0) { $dir{output} .= "_chrom_"; $dir{output} .= join("_", @chromnames); } if ($test) { $dir{output} .= "_test".$test_number; } } &RSAT::util::CheckOutDir($dir{output}); chdir($dir{output}); ## log file open $log, ">".$outfile{log} || die "cannot open error log file".$outfile{log}."\n"; ## error file open ERR, ">".$outfile{err} || die "cannot open error log file".$outfile{err}."\n"; ## List of contig sequence files open CTG, ">", $outfile{contigs} || die "cannot open contig file".$outfile{contigs}."\n"; # file with contig IDs ## cds sequences in fasta format unless ($seqonly) { @dbsplit = split /_core_/, $dbname; $org = $dbsplit[0]; $outfile{pp} = $org."_aa.fasta"; open PP, ">", $outfile{pp} || die "cannot open sequence file".$outfile{pp}."\n"; } #### print verbose &Verbose() if ($verbose); ################################################################ ## Open a new connection to EnsEMBL database, but this time we specify the DB name &RSAT::message::TimeWarn("Connecting EnsEMBL to retrieve the organism", "host=".$ensembl_host, "user=".$ensembl_user, "dbname=".$dbname, ) if ($main::verbose >= 1); my $db = new Bio::EnsEMBL::DBSQL::DBAdaptor(-host => $ensembl_host, -user => $ensembl_user, -dbname => $dbname); ## Get Species object my $meta_container = $db->get_MetaContainer(); my $tax_id = $meta_container->get_taxonomy_id(); my $species = $meta_container->get_Species(); my $organism_name = $species->binomial(); my $common_name = $species->common_name(); my @classification = $species->classification(); &RSAT::message::Info (join (":", "; Classification ", reverse(@classification))) if ($main::verbose >= 1); &RSAT::message::Info(join ("\t", "; Organism = ", $organism_name, "common name = ", $common_name, "NCBI taxid = ", $tax_id)) if ($main::verbose >= 1); ## Instantiate an object for the organism my $rsat_organism = $organisms->new_object(name=>$organism_name); $rsat_organism->push_attribute("names", $organism_name); $rsat_organism->set_attribute("taxonomy", join (";", reverse(@classification))); &RSAT::message::Debug("Db", $db) if ($main::verbose >= 3); my $slice_adaptor = $db->get_SliceAdaptor(); &RSAT::message::Debug("Adaptor", $slice_adaptor) if ($main::verbose >= 3); ## Collect the list of chromosomes my @slices; if (scalar(@chromnames) > 0) { ## Get selected chromosomes foreach $chromname (@chromnames) { &RSAT::message::Info(join("\t", "Getting slice", $slice_type, $chromname)) if ($main::verbose >= 1); push @slices, $slice_adaptor->fetch_by_region($slice_type, $chromname); } } else { ## get all the chromosomes @slices = @{$slice_adaptor->fetch_all($slice_type)}; } my $slices_number = scalar(@slices); my $s=0; while (my $slice = shift(@slices)) { # &RSAT::message::Debug("MEMORY USAGE", &Devel::Peek::mstat()) if ($main::verbose >= 0); $s++; my $slice_id = $slice->id(); my $slice_name = $slice->name(); ################################################################ ## TEMPORARY : a tricky fix for human genome, which contains 109 chromosome slices ! ## We discard slices which do not correspond to full chromosomes. ## These slices have a different name (they contain _NT_) if ($slice_name =~ /_NT_/) { &RSAT::message::Warning(join "\t", "Skipping slice (fragment of chromosome)", $s."/".$slices_number, $slice_id, $slice_name) if ($main::verbose >= 1); next; } &RSAT::message::psWarn("BEFORE GETTING SLICE: ".$s."/".$slices_number) if ($main::verbose >= 0); &RSAT::message::TimeWarn("Collecting data for slice", $slice_type, $s."/".$slices_number, "name:".$slice->name(), "seq_region_name:".$slice->seq_region_name(), ) if ($main::verbose >= 1); my $rsat_contig = $contigs->new_object(); $rsat_contig->force_attribute("id", $slice_id); $rsat_contig->set_attribute("name", $slice_name); $rsat_contig->set_attribute("type", $slice_type); $rsat_contig->set_attribute("accession", $slice->accession_number()); $rsat_contig->set_attribute("length", $slice->length()); $rsat_contig->set_attribute("description", $slice->desc()); $rsat_contig->set_attribute("chromosome", $slice->seq_region_name()); # $rsat_contig->set_attribute("is_circular", $slice->is_circular()); ## Get all repeat regions (features) unless (($no_rep) || ($seq_only)) { &RSAT::message::psWarn("Before collecting repeats for slice", $slice_name) if ($main::verbose >= 0); &RSAT::message::TimeWarn("\tGetting repeats for slice", $slice_name) if ($main::verbose >= 1); my $rep = 0; my @ensembl_repeats = @{$slice->get_all_RepeatFeatures()}; my $repeats_number = scalar(@ensembl_repeats); &RSAT::message::psWarn("After collecting repeats for slice", $slice_name) if ($main::verbose >= 0); while (my $ensembl_repeat = shift(@ensembl_repeats)){ $rep++; if (($test) && ($rep > $test_number)) { &RSAT::message::Info(join ("\t","TEST", $test_number, "skipping next repeats for contig", $rsat_contig->get_attribute("name"))) if ($main::verbose >= 1); last; } my $repeat_name = $ensembl_repeat->display_id(); $repeat_id = $repeat_name.".".$rep; &RSAT::message::Info(join("\t", ";", $slice_type, $slice->seq_region_name(), $s."/".$slices_number, "repeat", $rep."/".$repeats_number, $repeat_name)) if ($main::verbose >= 3); my $rsat_repeat = $repeat_regions->new_object(source=>"ensembl", name=>$repeat_name); $rsat_repeat->force_attribute("id", $repeat_id); $rsat_repeat->force_attribute("type", "repeat_region"); $rsat_repeat->set_attribute("contig", $slice_id); $rsat_repeat->set_attribute("start_pos", $ensembl_repeat->hstart); $rsat_repeat->set_attribute("end_pos", $ensembl_repeat->hend); my $ensembl_repeat_strand = $ensembl_repeat->hstrand; if ($ensembl_repeat_strand == 1){ $rsat_repeat->set_attribute("strand", "D"); } else { $rsat_repeat->set_attribute("strand", "R"); } } @ensembl_repeats = undef; &RSAT::message::psWarn("After undefining repeats for slice", $slice_name) if ($main::verbose >= 0); } ## Get all Gene objects unless ($seq_only) { &RSAT::message::psWarn("Before collecting genes for slice", $slice_name) if ($main::verbose >= 0); &RSAT::message::TimeWarn("\tGetting genes for slice", $slice_name) if ($main::verbose >= 1); my $g = 0; my @ensembl_genes = @{$slice->get_all_Genes()}; my $genes_number = scalar(@ensembl_genes); while (my $ensembl_gene = shift(@ensembl_genes)) { $g++; if (($test) && ($g > $test_number)) { &RSAT::message::Info(join ("\t","TEST", $test_number, "skipping next genes for contig", $rsat_contig->get_attribute("name"))) if ($main::verbose >= 1); last; } ## Create a new gene object my $gene_name = $ensembl_gene->external_name() || $ensembl_gene->stable_id(); warn join("\t", ";", $slice_type, $slice->seq_region_name(), $s."/".$slices_number, "gene", $g."/".$genes_number, $gene_name), "\n" if ($main::verbose >= 3); my $rsat_gene = $genes->new_object(source=>"ensembl", name=>$gene_name); $rsat_gene->force_attribute("id", $ensembl_gene->stable_id()); $rsat_gene->force_attribute("name", $gene_name); $rsat_gene->push_attribute("names", $gene_name); if ($ensembl_gene->stable_id() ne $gene_name) { $rsat_gene->push_attribute("names", $ensembl_gene->stable_id()); } $rsat_gene->set_attribute("contig", $ensembl_gene->slice->id()); $rsat_gene->set_attribute("start_pos", $ensembl_gene->start()); $rsat_gene->set_attribute("end_pos", $ensembl_gene->end()); if ($ensembl_gene->strand() == 1) { $rsat_gene->set_attribute("strand", "D"); } else { $rsat_gene->set_attribute("strand", "R"); } if ($ensembl_gene->description()) { $rsat_gene->set_attribute("description", $ensembl_gene->description()); } else { $rsat_gene->set_attribute("description",""); print ERR join ("\t", "No description for gene", $ensembl_gene->stable_id()), "\n"; } $rsat_gene->set_attribute("type", "gene"); $rsat_gene->force_attribute("organism", $organism_name); # print_DBEntries($ensembl_gene->get_all_DBLinks()); ## Get all Transcript objects for the current gene my $tr = 0; my @ensembl_transcript = @{$ensembl_gene->get_all_Transcripts()}; while (my $trans = shift(@ensembl_transcript)) { $tr++; warn join("\t", "transcript", $trans), "\n" if ($main::verbose >= 5); my $rsat_transcript = $transcripts->new_object(source=>"ensembl"); $rsat_transcript->force_attribute("id", $trans->stable_id()); $rsat_transcript->set_attribute("type", $trans->biotype()); # $rsat_transcript->set_attribute("contig", $trans->slice->id()); ## =CIRCULAR REFERENCE? $rsat_transcript->set_attribute("contig", $slice_id); ## NO CIRCULAR REFERENCE? $rsat_transcript->set_attribute("start_pos", $trans->start()); $rsat_transcript->set_attribute("end_pos", $trans->end()); if ($trans->strand() == 1) { $rsat_transcript->set_attribute("strand", "D"); } else { $rsat_transcript->set_attribute("strand", "R"); } if ($trans->description()) { $rsat_transcript->set_attribute("description", $trans->description()); } else { $rsat_transcript->set_attribute("description",""); print ERR join ("\t", "No description for gene", $trans->stable_id()), "\n"; } $rsat_transcript->force_attribute("organism", $organism_name); my $transcript_name = $trans->external_name() || $trans->stable_id(); # $transcript_name .= ".".$tr; $rsat_transcript->force_attribute("name", $transcript_name); $rsat_transcript->push_attribute("names", $transcript_name); if ($trans->stable_id() ne $transcript_name) { $rsat_transcript->push_attribute("names", $trans->stable_id()); } unless (($rsat_transcript->get_attribute("id") eq $rsat_gene->get_attribute("id")) || ($rsat_transcript->get_attribute("name") eq $rsat_gene->get_attribute("id"))){ $rsat_transcript->push_attribute("names", $rsat_gene->get_attribute("id")); } unless (($rsat_transcript->get_attribute("id") eq $rsat_gene->get_attribute("name")) || ($rsat_transcript->get_attribute("name") eq $rsat_gene->get_attribute("name"))){ $rsat_transcript->push_attribute("names", $rsat_gene->get_attribute("name")); } $rsat_transcript->set_attribute("gene", $rsat_gene->get_attribute("id")); if ($rsat_transcript->get_attribute("description") eq "") { $rsat_transcript->force_attribute("description", $rsat_gene->get_attribute("description")); } if ($rsat_transcript->get_attribute("type") eq "protein_coding"){ my $rsat_mrna = $mRNAs->new_object(source=>"ensembl"); $rsat_transcript->force_attribute("type", "mRNA"); $rsat_mrna->force_attribute("id", $rsat_transcript->get_attribute("id")); $rsat_mrna->set_attribute("type", $rsat_transcript->get_attribute("type")); $rsat_mrna->set_attribute("organism", $organism_name); $rsat_mrna->set_attribute("name", $rsat_transcript->get_attribute("name")); $rsat_mrna->push_attribute("names", $rsat_transcript->get_attribute("names")); $rsat_mrna->set_attribute("gene", $rsat_transcript->get_attribute("gene")); $rsat_mrna->set_attribute("description", $rsat_transcript->get_attribute("description")); $rsat_mrna->set_attribute("start_pos", $rsat_transcript->get_attribute("start_pos")); $rsat_mrna->set_attribute("end_pos", $rsat_transcript->get_attribute("end_pos")); $rsat_mrna->set_attribute("strand", $rsat_transcript->get_attribute("strand")); $rsat_mrna->set_attribute("contig", $rsat_transcript->get_attribute("contig")); } if ($rsat_transcript->get_attribute("type") eq "tRNA"){ my $rsat_trna = $tRNAs->new_object(source=>"ensembl"); $rsat_trna->force_attribute("id", $rsat_transcript->get_attribute("id")); $rsat_trna->set_attribute("type", $rsat_transcript->get_attribute("type")); $rsat_trna->set_attribute("organism", $organism_name); $rsat_trna->set_attribute("name", $rsat_transcript->get_attribute("name")); $rsat_trna->push_attribute("names", $rsat_transcript->get_attribute("names")); $rsat_trna->set_attribute("gene", $rsat_transcript->get_attribute("gene")); $rsat_trna->set_attribute("description", $rsat_transcript->get_attribute("description")); $rsat_trna->set_attribute("start_pos", $rsat_transcript->get_attribute("start_pos")); $rsat_trna->set_attribute("end_pos", $rsat_transcript->get_attribute("end_pos")); $rsat_trna->set_attribute("strand", $rsat_transcript->get_attribute("strand")); $rsat_trna->set_attribute("contig", $rsat_transcript->get_attribute("contig")); } if ($rsat_transcript->get_attribute("type") eq "rRNA"){ my $rsat_rrna = $rRNAs->new_object(source=>"ensembl"); $rsat_rrna->force_attribute("id", $rsat_transcript->get_attribute("id")); $rsat_rrna->set_attribute("type", $rsat_transcript->get_attribute("type")); $rsat_rrna->set_attribute("organism", $organism_name); $rsat_rrna->set_attribute("name", $rsat_transcript->get_attribute("name")); $rsat_rrna->push_attribute("names", $rsat_transcript->get_attribute("names")); $rsat_rrna->set_attribute("gene", $rsat_transcript->get_attribute("gene")); $rsat_rrna->set_attribute("description", $rsat_transcript->get_attribute("description")); $rsat_rrna->set_attribute("start_pos", $rsat_transcript->get_attribute("start_pos")); $rsat_rrna->set_attribute("end_pos", $rsat_transcript->get_attribute("end_pos")); $rsat_rrna->set_attribute("strand", $rsat_transcript->get_attribute("strand")); $rsat_rrna->set_attribute("contig", $rsat_transcript->get_attribute("contig")); } if ($rsat_transcript->get_attribute("type") eq "scRNA"){ my $rsat_scrna = $scRNAs->new_object(source=>"ensembl"); $rsat_scrna->force_attribute("id", $rsat_transcript->get_attribute("id")); $rsat_scrna->set_attribute("type", $rsat_transcript->get_attribute("type")); $rsat_scrna->set_attribute("organism", $organism_name); $rsat_scrna->set_attribute("name", $rsat_transcript->get_attribute("name")); $rsat_scrna->push_attribute("names", $rsat_transcript->get_attribute("names")); $rsat_scrna->set_attribute("gene", $rsat_transcript->get_attribute("gene")); $rsat_scrna->set_attribute("description", $rsat_transcript->get_attribute("description")); $rsat_scrna->set_attribute("start_pos", $rsat_transcript->get_attribute("start_pos")); $rsat_scrna->set_attribute("end_pos", $rsat_transcript->get_attribute("end_pos")); $rsat_scrna->set_attribute("strand", $rsat_transcript->get_attribute("strand")); $rsat_scrna->set_attribute("contig", $rsat_transcript->get_attribute("contig")); } if ($rsat_transcript->get_attribute("type") eq "misc_RNA"){ my $rsat_miscrna = $misc_RNAs->new_object(source=>"ensembl"); $rsat_miscrna->force_attribute("id", $rsat_transcript->get_attribute("id")); $rsat_miscrna->set_attribute("type", $rsat_transcript->get_attribute("type")); $rsat_miscrna->set_attribute("organism", $organism_name); $rsat_miscrna->set_attribute("name", $rsat_transcript->get_attribute("name")); $rsat_miscrna->push_attribute("names", $rsat_transcript->get_attribute("names")); $rsat_miscrna->set_attribute("gene", $rsat_transcript->get_attribute("gene")); $rsat_miscrna->set_attribute("description", $rsat_transcript->get_attribute("description")); $rsat_miscrna->set_attribute("start_pos", $rsat_transcript->get_attribute("start_pos")); $rsat_miscrna->set_attribute("end_pos", $rsat_transcript->get_attribute("end_pos")); $rsat_miscrna->set_attribute("strand", $rsat_transcript->get_attribute("strand")); $rsat_miscrna->set_attribute("contig", $rsat_transcript->get_attribute("contig")); } ## Get CDS ID and coordinates (relative to chromosome) - there is a strand trick (see API doc) ## Problem: coordinates are strange if($trans->translation()) { my $rsat_cds = $CDSs->new_object(source=>"ensembl"); $rsat_cds->force_attribute("id", $trans->translation()->stable_id()); $rsat_cds->set_attribute("type", "CDS"); $rsat_cds->set_attribute("name", $rsat_transcript->get_attribute("name")); $rsat_cds->set_attribute("contig", $rsat_transcript->get_attribute("contig")); $rsat_cds->set_attribute("start_pos", $trans->coding_region_start()); $rsat_cds->set_attribute("end_pos", $trans->coding_region_end()); $rsat_cds->set_attribute("strand", $rsat_transcript->get_attribute("strand")); $rsat_cds->set_attribute("description", $rsat_transcript->get_attribute("description")); $rsat_cds->set_attribute("organism", $rsat_transcript->get_attribute("organism")); $rsat_cds->set_attribute("transcript", $rsat_transcript->get_attribute("id")); $rsat_cds->push_attribute("names", $rsat_gene->get_attribute("name")); unless ($rsat_gene->get_attribute("name") eq $rsat_gene->get_attribute("id")){ $rsat_cds->push_attribute("names", $rsat_gene->get_attribute("id")); } unless (($rsat_gene->get_attribute("id") eq $rsat_transcript->get_attribute("name")) || ($rsat_gene->get_attribute("name") eq $rsat_transcript->get_attribute("name"))){ $rsat_cds->push_attribute("names", $rsat_transcript->get_attribute("name")); } ## Print the translated sequence for the current CDS &PrintNextSequence(PP,"fasta",60,$trans->translation()->seq(),$trans->translation()->stable_id(),$rsat_gene->get_attribute("description")); } # unless ($no_ons) { # ## Get all Exon objects # foreach my $exon (@{$trans->get_all_Exons()}) { # my @exonfeature = &get_exonfeature($exon); # } # ## Get all Intron objects # foreach my $intron (@{$trans->get_all_Introns()}) { # my @intronfeature = &get_intronfeature($intron); # $intronfeature[0] = "Trnscrpt - ".$transcriptID; # } # } } @ensembl_transcript = undef; } &RSAT::message::psWarn("After collecting genes for slice", $s, $slice_name) if ($main::verbose >= 0); @ensembl_genes = undef; &RSAT::message::psWarn("After undefining genes for slice", $s, $slice_name) if ($main::verbose >= 0); } ################################################################ ## Export sequence unless otherwise specified my $seq_file = $rsat_contig->get_attribute("id"); $seq_file =~ s/\:/_/g; $seq_file .= ".raw"; my $masked_seq_file = $seq_file; $masked_seq_file =~ s/\.raw$/_repeat_masked.raw/; # if ($rsat_contig->get_attribute("is_circular")) { # $contig_shape = 'circular'; # } else { # $contig_shape = 'linear'; # } print CTG join ("\t", $seq_file, $rsat_contig->get_attribute("id")), "\n"; unless ($no_seq) { ## Export slice sequence (unmasked) &RSAT::message::TimeWarn("Saving sequence for slice", $s."/".$slices_number, $slice_type, $slice->seq_region_name(), $slice_name) if ($main::verbose >= 1); open SEQ, ">".$seq_file || die "cannot open error log file".$seq_file."\n"; print SEQ $slice->seq(); close SEQ; ## not sure this is useful, but to try improving garbage collection &RSAT::message::psWarn("After saving sequence for slice", $s, $slice_name) if ($main::verbose >= 0); ## Export slice sequence (hard masked) unless ($no_masked) { &RSAT::message::TimeWarn("Getting masked sequence for slice", $s."/".$slices_number, $slice_type, $slice->seq_region_name(), $slice_name) if ($main::verbose >= 1); &RSAT::message::psWarn("Before collecting repeatmasked sequence for slice", $s, $slice_name) if ($main::verbose >= 0); my $masked_sequence_slice = $slice->get_repeatmasked_seq(); open MASKED_SEQ, ">".$masked_seq_file || die "cannot open error log file".$masked_seq_file."\n"; print MASKED_SEQ $masked_sequence_slice->seq(); close MASKED_SEQ; ## not sure this is useful, but to try improving garbage collection &RSAT::message::psWarn("After collecting repeatmasked sequence for slice", $s, $slice_name) if ($main::verbose >= 0); $masked_sequence_slice = undef; &RSAT::message::psWarn("After destroying repeatmasked sequence for slice", $s, $slice_name) if ($main::verbose >= 0); } } &RSAT::message::psWarn("AFTER COLLECTING INFO FOR SLICE: ".$s."/".$slices_number) if ($main::verbose >= 0); $slice = undef; &RSAT::message::psWarn("AFTER UNDEFINING SLICE: ".$s."/".$slices_number) if ($main::verbose >= 0); } ################################################################ ### Save result in tab files unless ($seq_only) { &RSAT::message::TimeWarn("Exporting collected objects to directory",$dir{output}) if ($main::verbose >= 1); chdir $dir{main}; warn "; Main directory\t", $dir{main}, "\n" if ($main::verbose >= 3); foreach $class_factory ($organisms, $contigs, $genes, $repeat_regions, $CDSs, $mRNAs, $tRNAs, $rRNAs, $scRNAs, $misc_RNAs) { $class_factory->dump_tables(); $class_factory->generate_sql(dir=>"$dir{output}/sql_scripts", schema=>$schema, host=>$host, user=>$user, password=>$password ); } &ExportMakefile(@classes); chdir($dir{output}); &PrintStats($outfile{stats}, @classes); } ################################################################ ## Report execution time and close output stream my $exec_time = &RSAT::util::ReportExecutionTime($start_time); ## This has to be exectuted by all scripts print $log $exec_time if ($main::verbose >= 1); ## only report exec time if verbosity is specified close $log if ($outfile{log}); close ERR if ($outfile{err}); close CTG if ($outfile{contigs}); close PP if ($outfile{pp}); ################################################################ ## Report the output directory warn "; Results stored in directory\t", $dir{output}, "\n" if ($main::verbose >= 1); exit(0); } ################################################################ ################### SUBROUTINE DEFINITION ###################### ################################################################ ################################################################ #### display full help message sub PrintHelp { open HELP, "| more"; print HELP <