Archive Ensembl HomeArchive Ensembl Home
conservation_score.pm
Go to the documentation of this file.
00001 #
00002 # Bio::Das::ProServer::SourceAdaptor::conservation_score
00003 #
00004 # Copyright EnsEMBL Team
00005 #
00006 # You may distribute this module under the same terms as perl itself
00007 #
00008 # pod documentation - main docs before the code
00009 
00010 =head1 NAME
00011 
00012 Bio::Das::ProServer::SourceAdaptor::conservation_score - Extension of the ProServer for e! conservation scores
00013 
00014 =head1 INHERITANCE
00015 
00016 This module inherits attributes and methods from Bio::Das::ProServer::SourceAdaptor
00017 
00018 =head1 DAS CONFIGURATION FILE
00019 
00020 There are some specific parameters for this module you can use in the DAS server configuration file
00021 
00022 =head2 registry
00023 
00024 Your registry configuration file to connect to the compara database
00025 
00026 =head2 database
00027 
00028 The species name in your Registry configuration file.
00029 
00030 =head2 this_species
00031 
00032 The main species. Features will be shown for this species.
00033 
00034 =head2 other_species
00035 
00036 The other species. This DAS track will show alignments between this_species and other_species.
00037 You can add more than one other species separated by commas.
00038 
00039 =head2 analysis
00040 
00041 The method_link_type. This defines the type of score. E.g. GERP_CONSERVATION_SCORE
00042 See perldoc Bio::EnsEMBL::Compara::MethodLinkSpeciesSet for more details about the
00043 method_link_type
00044 
00045 =head2 Example
00046 
00047 =head3 registry configuration file
00048 
00049   use strict;
00050   use Bio::EnsEMBL::Utils::ConfigRegistry;
00051   use Bio::EnsEMBL::Compara::DBSQL::DBAdaptor;
00052 
00053   new Bio::EnsEMBL::Compara::DBSQL::DBAdaptor(
00054       -host => 'ensembldb.ensembl.org',
00055       -user => 'anonymous',
00056       -port => 3306,
00057       -species => 'ensembl-compara-41',
00058       -dbname => 'ensembl_compara_41');
00059 
00060 
00061 =head3 DAS server configuration file
00062 
00063   [general]
00064   hostname    = ecs4b.internal.sanger.ac.uk
00065   prefork     = 6
00066   maxclients  = 100
00067   port        = 9013
00068 
00069   [conservation_score]
00070   registry        = /home/foo/ProServer/eg/reg.pl
00071   state           = on
00072   adaptor         = compara
00073   database        = ensembl-compara-41
00074   this_species    = Homo sapiens
00075   other_species   = Mus musculus, Rattus norvegicus, Canis familiaris, Gallus gallus, Bos taurus, Monodelphis domestica
00076   analysis        = GERP_CONSTRAINED_ELEMENT
00077   description     = 7 way mlagan alignment
00078   group_type      = default
00079 
00080 =cut
00081 
00082 package Bio::Das::ProServer::SourceAdaptor::conservation_score;
00083 
00084 use strict;
00085 use Bio::EnsEMBL::Registry;
00086 use Bio::EnsEMBL::Utils::Exception;
00087 
00088 use base qw( Bio::Das::ProServer::SourceAdaptor );
00089 
00090 sub init
00091 {
00092     my ($self) = @_;
00093 
00094     $self->{'capabilities'} = { 'features' => '1.0',
00095                                 'stylesheet' => '1.0' };
00096 
00097     my $registry = $self->config()->{'registry'};
00098     unless (defined $registry) {
00099      throw("registry not defined\n");
00100     }
00101 
00102     if (not $Bio::EnsEMBL::Registry::registry_register->{'seen'}) {
00103         Bio::EnsEMBL::Registry->load_all($registry);
00104     }
00105 
00106 }
00107 
00108 sub build_features
00109 {
00110     my ($self, $opts) = @_;
00111 
00112     my $db      = "Bio::EnsEMBL::Registry";
00113     $db->no_version_check(1);
00114     my $dbname  = $self->config()->{'database'};
00115 
00116     #need to put adaptors here and not in init 
00117     my $meta_con =
00118         $db->get_adaptor($dbname, 'compara', 'MetaContainer') or
00119             die "no metadbadaptor:$dbname, 'compara','MetaContainer' \n";
00120 
00121     my $mlss_adaptor =
00122         $db->get_adaptor($dbname, 'compara', 'MethodLinkSpeciesSet') or
00123             die "can't get $dbname, 'compara', 'MethodLinkSpeciesSet'\n";
00124 
00125     my $cs_adaptor =
00126         $db->get_adaptor($dbname, 'compara', 'ConservationScore') or
00127             die "can't get $dbname, 'compara', 'ConservationScore'\n";
00128 
00129     my $species  = $self->config()->{'this_species'};
00130     my $slice_adaptor =
00131         $db->get_adaptor($species, 'core', 'Slice') or
00132             die "can't get $species, 'core', 'Slice'\n";  
00133 
00134     my $genome_db_adaptor =
00135         $db->get_adaptor($dbname, 'compara', 'GenomeDB') or
00136             die "can't get $dbname, 'compara', 'GenomeDB'\n";
00137 
00138     my $genomedbs = $genome_db_adaptor->fetch_all();
00139 
00140     my $daschr      = $opts->{'segment'} || return ( );
00141     my $dasstart    = $opts->{'start'} || return ( );
00142     my $dasend      = $opts->{'end'} || return ( );
00143 
00144     my $species1    = $self->config()->{'this_species'};
00145     my @other_species = split(/\s*\,\s*/, $self->config()->{'other_species'});
00146     my $chr1        = $daschr;
00147     my $start1      = $dasstart;
00148     my $end1        = $dasend;
00149 
00150     my $method_link = $self->config()->{'analysis'};
00151 
00152     my $stored_max_alignment_length;
00153 
00154     my $values = $meta_con->list_value_by_key("max_alignment_length");
00155 
00156     if(@$values) {
00157         $stored_max_alignment_length = $values->[0];
00158     }
00159 
00160     my $species1_genome_db;
00161     my @other_species_genome_dbs;
00162 
00163     ## Get the Bio::EnsEMBL::Compara::GenomeDB object for the primary species
00164     foreach my $this_genome_db (@$genomedbs){
00165       if ($this_genome_db->name eq $species1) {
00166         $species1_genome_db = $this_genome_db;
00167       }
00168     }
00169     if (!defined($species1_genome_db)) {
00170       die "No species called $species1 in the database -- check spelling\n";
00171     }
00172 
00173     ## Get the Bio::EnsEMBL::Compara::GenomeDB objects for the remaining species
00174     foreach my $this_other_species (@other_species) {
00175       my $this_other_genome_db;
00176       foreach my $this_genome_db (@$genomedbs){
00177         if ($this_genome_db->name eq $this_other_species) {
00178           $this_other_genome_db = $this_genome_db;
00179           last;
00180         }
00181       }
00182       if (!defined($this_other_genome_db)) {
00183         die "No species called $this_other_species in the database -- check spelling\n";
00184       }
00185       push(@other_species_genome_dbs, $this_other_genome_db);
00186     }
00187 
00188     ## Fetch the Bio::EnsEMBL::Compara::MethodLinkSpeciesSet object
00189     my $method_link_species_set;
00190     $method_link_species_set =
00191         $mlss_adaptor->fetch_by_method_link_type_GenomeDBs(
00192             $method_link, [$species1_genome_db, @other_species_genome_dbs]);
00193 
00194     ##Fetch the Bio::EnsEMBL::Slice object
00195     my $slice = $slice_adaptor->fetch_by_region(undef, $chr1, $start1, $end1);
00196 
00197     #Fetch conservation scores
00198     my $conservation_scores = $cs_adaptor->fetch_all_by_MethodLinkSpeciesSet_Slice($method_link_species_set, $slice);
00199 
00200     ## Build the results array
00201     my @results = ();
00202 
00203     foreach my $score (@$conservation_scores) {
00204 
00205     unless (defined $score->diff_score) {
00206         next;
00207     }
00208 
00209     #my $id = $score->genomic_align_block_id;
00210 
00211     my $id = $score->genomic_align_block_id ."_" . $score->position;
00212     my $label = "Conservation scores";
00213 
00214     # note will contain expected and observed scores and window size 
00215     my $note = sprintf("Expected %.3f Observed %.3f Diff %.3f Window Size %d Max", $score->expected_score, $score->observed_score, $score->diff_score, $score->window_size); 
00216 
00217     #my $start_pos = $start1 + $score->position;
00218     #my $end_pos = $start1 + $score->position + $score->window_size;
00219     my $start_pos = $start1 + $score->position - 1;
00220     my $end_pos = $start_pos + $score->window_size - 1;
00221 
00222     my $new_score = $score->diff_score()*-1;
00223     push @results, {
00224         'id'    => $id,
00225         'label' => $label,
00226         'method'=> $method_link,
00227         'start' => $start_pos,
00228         'end'   => $end_pos,
00229         'ori'   => '+',
00230         'score' => $new_score,
00231         'note'  => $note,
00232         'typecategory' => 'Conservation scores',
00233         'type'  => 'histogram'
00234         };
00235     }
00236 
00237     return @results;
00238 }
00239 
00240 sub das_stylesheet
00241 {
00242     my $self = shift;
00243 
00244     return <<EOT;
00245 <!DOCTYPE DASSTYLE SYSTEM "http://www.biodas.org/dtd/dasstyle.dtd">
00246 <DASSTYLE>
00247     <STYLESHEET version="1.0">
00248         <CATEGORY id="Conservation scores">
00249            <TYPE id="histogram"><GLYPH><HISTOGRAM>
00250               <MIN>-3</MIN>
00251               <MAX>3</MAX>
00252               <HEIGHT>100</HEIGHT>
00253               <STEPS>50</STEPS>
00254               <COLOR1>red</COLOR1>
00255               <COLOR2>yellow</COLOR2>
00256               <COLOR3>blue</COLOR3>
00257             </HISTOGRAM></GLYPH></TYPE>
00258         </CATEGORY>
00259     </STYLESHEET>
00260 </DASSTYLE>
00261 EOT
00262 }
00263 
00264 1;