projects/tri/q2.pl

   1 #!/usr/bin/env perl
   2 use strict;
   3 use warnings;
   4 use Bio::AlignIO;
   5 use DBI;
   6
   7 use Data::Dump qw (ddx);
   8 #use IO::Unread;
   9
  10 use Bio::Align::DNAStatistics;
  11 use Bio::EnsEMBL::Registry;
  12
  13 ## Load the registry automatically
  14 my $reg = "Bio::EnsEMBL::Registry";
  15 #my $url = 'mysql://anonymous@ensembldb.ensembl.org';
  16 my $url = 'mysql://root@localhost';
  17 $reg->load_registry_from_url($url);
  18
  19 my $compara_db_adaptor = $reg->get_DBAdaptor('Multi', 'compara');
  20 my $genome_db_adaptor = $compara_db_adaptor->get_GenomeDBAdaptor();
  21 my $dbh = $compara_db_adaptor->dbc->db_handle;
  22
  23 my $sql = qq{
  24 SELECT DISTINCT ss.genome_db_id, gdb.taxon_id, gdb.name, gdb.assembly FROM species_set_tag sst
  25  JOIN species_set ss USING (species_set_id)
  26  JOIN genome_db gdb USING (genome_db_id)
  27  WHERE sst.value = 'mammals';
  28 };
  29
  30 my $sth = $dbh->prepare($sql);
  31 $sth->execute();
  32
  33 my (@IDsGnomeDB, %GnomeDBnfo);
  34 while(my @retdat = $sth->fetchrow()) {
  35         #warn join("\t",@retdat),"\n";
  36         #$GnomeDBnfo{$retdat[0]} = [@retdat[1,2,3]];
  37         my $GnomedbID = shift @retdat;
  38         push @IDsGnomeDB,$GnomedbID;
  39         $GnomeDBnfo{$GnomedbID} = \@retdat;
  40 }
  41
  42 ddx \%GnomeDBnfo;
  43
  44 #my $MammaliaTaxID = 40674; # http://www.ncbi.nlm.nih.gov/Taxonomy/Browser/wwwtax.cgi?mode=Undef&id=40674&lvl=3&keep=1&srchmode=1&unlock
  45 #my $MammaliaGDBs = $genome_db_adaptor->fetch_all_by_ancestral_taxon_id($MammaliaTaxID,1);
  46 #ddx $MammaliaGDBs;
  47
  48 $sql = qq{
  49         SELECT DISTINCT method_link_species_set_id FROM method_link_species_set
  50           JOIN method_link USING (method_link_id)
  51           JOIN species_set ss ON ss.species_set_id=method_link_species_set.species_set_id
  52           JOIN species_set ss0 ON ss.genome_db_id=ss0.genome_db_id
  53           JOIN species_set ss1 ON ss.species_set_id=ss1.species_set_id
  54           JOIN species_set ss2 ON ss1.genome_db_id=ss2.genome_db_id
  55           JOIN species_set_tag sst0 ON sst0.species_set_id=ss0.species_set_id
  56           JOIN species_set_tag sst2 ON sst2.species_set_id=ss2.species_set_id
  57           WHERE sst0.value = 'mammals' AND sst2.value = 'mammals' AND ss1.genome_db_id != ss.genome_db_id AND method_link.type='ENSEMBL_ORTHOLOGUES';
  58 };
  59 $sth = $dbh->prepare($sql);
  60 $sth->execute();
  61
  62 my (@IDsMLSS);
  63 while(my @retdat = $sth->fetchrow()) {
  64         push @IDsMLSS,$retdat[0];
  65 }
  66 warn 'MLSS_ID: ',join(', ',@IDsMLSS),"\n";
  67
  68 ## The BioPerl alignment formatter
  69 my $alignIO = Bio::AlignIO->newFh(-format => "clustalw");
  70
  71 my $HomologyAdaptor = $reg->get_adaptor("Multi", "compara", "Homology");
  72 for my $mlss (@IDsMLSS) {
  73         my $homologies = $HomologyAdaptor->fetch_all_by_MethodLinkSpeciesSet($mlss);
  74         print '-----',$mlss,"-----\n";
  75     ## For each homology
  76     foreach my $this_homology (@{$homologies}) {
  77                 next unless defined $this_homology->dn();
  78           my $description = $this_homology->description;
  79           # next unless ($description =~ /orth/);    # uncomment for orthologs only
  80         my ($a,$b) = @{$this_homology->gene_list};
  81         my $spa = $a->taxon->get_short_name;
  82         my $spb = $b->taxon->get_short_name;
  83         my $labela = $a->stable_id;
  84         $labela .= "(" . $a->display_label . ")" if $a->display_label;
  85         my $labelb = $b->stable_id;
  86         $labelb .= "(" . $b->display_label . ")" if $b->display_label;
  87                  print "$spa,$labela,$spb,$labelb,$description\n";
  88                 #ddx $this_homology;
  89     }
  90     print "\n";
  91 }
  92
  93 __END__