[Bio] / FigKernelPackages / GenerateClusters.pm Repository:
ViewVC logotype

Annotation of /FigKernelPackages/GenerateClusters.pm

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1.2 - (view) (download) (as text)

1 : overbeek 1.1 # -*- perl -*-
2 :     ########################################################################
3 :     # Copyright (c) 2003-2011 University of Chicago and Fellowship
4 :     # for Interpretations of Genomes. All Rights Reserved.
5 :     #
6 :     # This file is part of the SEED Toolkit.
7 :     #
8 :     # The SEED Toolkit is free software. You can redistribute
9 :     # it and/or modify it under the terms of the SEED Toolkit
10 :     # Public License.
11 :     #
12 :     # You should have received a copy of the SEED Toolkit Public License
13 :     # along with this program; if not write to the University of Chicago
14 :     # at info@ci.uchicago.edu or the Fellowship for Interpretation of
15 :     # Genomes at veronika@thefig.info or download a copy from
16 :     # http://www.theseed.org/LICENSE.TXT.
17 :     ########################################################################
18 :    
19 :     package GenerateClusters;
20 :    
21 :     use strict;
22 :     no warnings 'redefine'; ## prevents spurious warnings due to use recursion
23 :     use SeedUtils;
24 :     use SAPserver;
25 :     use Data::Dumper;
26 :    
27 :     sub generate_clusters {
28 : overbeek 1.2 my($ref,$genomes,$min_iden) = @_;
29 : overbeek 1.1
30 : overbeek 1.2 if (! $min_iden) { $min_iden = 50 }
31 : overbeek 1.1 my $sapO = SAPserver->new();
32 :    
33 :     my %neighH;
34 :     foreach my $g (($ref,@$genomes))
35 :     {
36 :     my $pegH = $sapO->all_features( -ids => [$g], -type => ['peg'] );
37 :     my $pegs = $pegH->{$g};
38 :     my $locH = $sapO->fid_locations( -ids => $pegs, -boundaries => 1);
39 :     $neighH{$g} = &neighbors($locH);
40 :     print STDERR "$g done\n";
41 :     }
42 : overbeek 1.2 my $bbhsH = &bbhs($ref,$genomes,$sapO,$min_iden);
43 : overbeek 1.1 my @connected;
44 :    
45 :     foreach my $peg (sort {&SeedUtils::by_fig_id($a,$b) } keys(%{$neighH{$ref}}))
46 :     {
47 :     my @neigh = keys(%{$neighH{$ref}->{$peg}});
48 :     foreach my $peg2 (sort { &SeedUtils::by_fig_id($a,$b) } @neigh)
49 :     {
50 :     my $preserved = 0;
51 :     my $not_preserved = 0;
52 :    
53 :     foreach my $g (@$genomes)
54 :     {
55 :    
56 :     my $bbh_peg = $bbhsH->{$g}->{$peg};
57 :     my $bbh_neigh = $bbhsH->{$g}->{$peg2};
58 :     if ($bbh_peg && $bbh_neigh)
59 :     {
60 :     if ($neighH{$g}->{$bbh_peg}->{$bbh_neigh})
61 :     {
62 :     $preserved++;
63 :     }
64 :     else
65 :     {
66 :     $not_preserved++;
67 :     }
68 :     }
69 :     }
70 :     push(@connected,[$peg,$peg2,$preserved,$not_preserved]);
71 :     }
72 :     }
73 :     return \@connected;
74 :     }
75 :    
76 :     sub ok_coverage {
77 :     my($b1,$e1,$ln1,$b2,$e2,$ln2) = @_;
78 :    
79 :     my $ln = ($ln1 > $ln2) ? $ln1 : $ln2;
80 :     return ((($ln - $ln1) / $ln) < 0.2) && ((($ln - $ln2) / $ln) < 0.2);
81 :     }
82 :    
83 :     sub bbhs {
84 : overbeek 1.2 my($ref,$genomes,$sapO,$min_iden) = @_;
85 : overbeek 1.1
86 :     my $corrH = {};
87 :     foreach my $g (@$genomes)
88 :     {
89 :     # open(CORR,"<Corr/$g") || die "could not open correspondence for $g";
90 :     # while (defined($_ = <CORR>))
91 :     # {
92 :     # chop;
93 :     # my @x = split(/\t/,$_);
94 :     # if ($x[8] eq "<=>")
95 :     # {
96 :     # $corrH->{$g}->{$x[0]} = $x[1];
97 :     # }
98 :     # }
99 :     # close(CORR);
100 :     print STDERR "bbhs for $ref and $g\n";
101 :     # $corrH->{$g} = $sapO->gene_correspondence_map( -genome1 => $ref, -genome2 => $g );
102 :     my $map = $sapO->gene_correspondence_map( -genome1 => $ref, -genome2 => $g,-fullOutput => 1 );
103 :     foreach my $tuple (@$map)
104 :     {
105 : overbeek 1.2 if (($tuple->[8] eq "<=>") &&
106 :     &ok_coverage($tuple->[11],$tuple->[12],$tuple->[13],$tuple->[14],$tuple->[15],$tuple->[16]) &&
107 :     ($min_iden <= $tuple->[9]))
108 : overbeek 1.1 {
109 :     $corrH->{$g}->{$tuple->[0]} = $tuple->[1];
110 :     }
111 :     }
112 :     }
113 :     return $corrH;
114 :     }
115 :    
116 :     sub neighbors {
117 :     my($locH) = @_;
118 :    
119 :     my @locs;
120 :     foreach my $peg (keys(%$locH))
121 :     {
122 :     my $loc = $locH->{$peg};
123 :     if ($loc && ($loc =~ /^.*:(\S+)_(\d+)[+-](\d+)$/))
124 :     {
125 :     my($contig,$beg,$strand,$len) = ($1,$2,$3,$4);
126 :     my $mid = ($strand eq '+') ? ($beg + ($len/2)) : ($beg - ($len/2));
127 :     push(@locs,[$peg,$contig,$mid]);
128 :     }
129 :     }
130 :     @locs = sort { ($a->[1] cmp $b->[1]) or ($a->[2] <=> $b->[2]) } @locs;
131 :     my $neighH = {};
132 :     my($i,$j);
133 :     for ($i=0; ($i < @locs); $i++)
134 :     {
135 :     $j = ($i >= 5) ? ($i - 5) : 0;
136 :     while ($j < $i)
137 :     {
138 :     if ($locs[$j]->[1] eq $locs[$i]->[1])
139 :     {
140 :     $neighH->{$locs[$i]->[0]}->{$locs[$j]->[0]} = 1
141 :     # push(@$neigh,$locs[$j]->[0]);
142 :     }
143 :     $j++;
144 :     }
145 :     $j = $i+1;
146 :     while (($j < @locs) && ($j <= $i+5))
147 :     {
148 :     if ($locs[$j]->[1] eq $locs[$i]->[1])
149 :     {
150 :     $neighH->{$locs[$i]->[0]}->{$locs[$j]->[0]} = 1
151 :     # push(@$neigh,$locs[$j]->[0]);
152 :     }
153 :     $j++;
154 :     }
155 :     # $neighH->{$locs[$i]->[0]} = $neigh;
156 :     }
157 :     return $neighH;
158 :     }
159 :    
160 :     1;

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3