modified: makesrc3.mk
[GalaxyCodeBases.git] / projects / tri / q2.pl
blob1232c86e5b02dc80d4728e4bb6055cbe210fb515
1 #!/usr/bin/env perl
2 use strict;
3 use warnings;
4 use Bio::AlignIO;
5 use DBI;
7 use Data::Dump qw (ddx);
8 #use IO::Unread;
10 use Bio::Align::DNAStatistics;
11 use Bio::EnsEMBL::Registry;
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);
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;
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';
30 my $sth = $dbh->prepare($sql);
31 $sth->execute();
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;
42 ddx \%GnomeDBnfo;
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;
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';
59 $sth = $dbh->prepare($sql);
60 $sth->execute();
62 my (@IDsMLSS);
63 while(my @retdat = $sth->fetchrow()) {
64 push @IDsMLSS,$retdat[0];
66 warn 'MLSS_ID: ',join(', ',@IDsMLSS),"\n";
68 ## The BioPerl alignment formatter
69 my $alignIO = Bio::AlignIO->newFh(-format => "clustalw");
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;
90 print "\n";
93 __END__