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

Annotation of /FigKernelPackages/CompareMR.pm

Parent Directory Parent Directory | Revision Log Revision Log


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

1 : overbeek 1.1 #
2 :     # Copyright (c) 2003-2006 University of Chicago and Fellowship
3 :     # for Interpretations of Genomes. All Rights Reserved.
4 :     #
5 :     # This file is part of the SEED Toolkit.
6 :     #
7 :     # The SEED Toolkit is free software. You can redistribute
8 :     # it and/or modify it under the terms of the SEED Toolkit
9 :     # Public License.
10 :     #
11 :     # You should have received a copy of the SEED Toolkit Public License
12 :     # along with this program; if not write to the University of Chicago
13 :     # at info@ci.uchicago.edu or the Fellowship for Interpretation of
14 :     # Genomes at veronika@thefig.info or download a copy from
15 :     # http://www.theseed.org/LICENSE.TXT.
16 :     #
17 :    
18 :     package CompareMR;
19 :    
20 :     use Carp;
21 :     use strict;
22 :     use FIG_Config;
23 :     use Tracer;
24 :     use Data::Dumper;
25 :    
26 :     use FIG;
27 :     my $fig = new FIG;
28 :     use FIGV;
29 :    
30 :     sub compare_genomes_MR {
31 : paarmann 1.3 my($arg1,$genome2, $fig_services, $genome1) = @_;
32 : overbeek 1.1
33 : paarmann 1.3 unless ($fig_services) {
34 :     if (-d "$FIG_Config::organisms/$arg1")
35 :     {
36 :     $fig_services = $fig;
37 :     $genome1 = $arg1;
38 :     }
39 :     else
40 : overbeek 1.1 {
41 : paarmann 1.3 if (-d $arg1)
42 : overbeek 1.1 {
43 : paarmann 1.3 my $figV = new FIGV($arg1,$fig);
44 :     $fig_services = $figV;
45 :     if ($arg1 =~ /(\d+\.\d+)$/)
46 :     {
47 :     $genome1 = $1;
48 :     }
49 :     else
50 :     {
51 :     die "invalid first argument: $arg1";
52 :     }
53 : overbeek 1.1 }
54 : paarmann 1.3 else
55 : overbeek 1.1 {
56 : paarmann 1.3 die "Invalid first argument; $arg1";
57 : overbeek 1.1 }
58 :     }
59 :     }
60 :    
61 :     my $data1 = &data_on_subs_roles_pegs($fig_services,$genome1);
62 :     my $data2 = &data_on_subs_roles_pegs($fig_services,$genome2);
63 :    
64 :     my $common = &process_common($fig_services,$genome1,$genome2,$data1,$data2);
65 :     my $in1_not2 = &process_diff($fig_services,$genome1,$genome2,$data1,$data2);
66 :     my $in2_not1 = &process_diff($fig_services,$genome2,$genome1,$data2,$data1);
67 :     return ($common,$in1_not2,$in2_not1);
68 :     }
69 :    
70 :     sub process_common {
71 :     my($fig_services,$genome1,$genome2,$data1,$data2) = @_;
72 :    
73 :     my $common = [];
74 :     my $i1 = 0;
75 :     my $i2 = 0;
76 :     while (($i1 < @$data1) && ($i2 < @$data2))
77 :     {
78 :     if ( ($data1->[$i1]->[0] lt $data2->[$i2]->[0]) ||
79 :     (($data1->[$i1]->[0] eq $data2->[$i2]->[0]) && ($data1->[$i1]->[1] lt $data2->[$i2]->[1])))
80 :     {
81 :     $i1++;
82 :     }
83 :     elsif ( ($data1->[$i1]->[0] gt $data2->[$i2]->[0]) ||
84 :     (($data1->[$i1]->[0] eq $data2->[$i2]->[0]) && ($data1->[$i1]->[1] gt $data2->[$i2]->[1])))
85 :     {
86 :     $i2++;
87 :     }
88 :     else
89 :     {
90 :     my $currS = $data1->[$i1]->[0];
91 :     my $currR = $data1->[$i1]->[1];
92 :     my $pegs1 = [];
93 :     while (($i1 < @$data1) && ($data1->[$i1]->[0] eq $currS) && ($data1->[$i1]->[1] eq $currR))
94 :     {
95 :     push(@$pegs1,$data1->[$i1]->[2]);
96 :     $i1++;
97 :     }
98 :    
99 :     my $pegs2 = [];
100 :     while (($i2 < @$data2) && ($data2->[$i2]->[0] eq $currS) && ($data2->[$i2]->[1] eq $currR))
101 :     {
102 :     push(@$pegs2,$data2->[$i2]->[2]);
103 :     $i2++;
104 :     }
105 :     push(@$common,[$currS,$currR,$pegs1,$pegs2]);
106 :     }
107 :     }
108 :     return $common;
109 :     }
110 :    
111 :     sub process_diff {
112 :     my($fig_services,$genomeA,$genomeB,$dataA,$dataB) = @_;
113 :    
114 :     my $in1_not2 = [];
115 :     my $i1 = 0;
116 :     my $i2 = 0;
117 :     while ($i1 < @$dataA)
118 :     {
119 : overbeek 1.2 if ( ($i2 == @$dataB) || ($dataA->[$i1]->[0] lt $dataB->[$i2]->[0]) ||
120 : overbeek 1.1 (($dataA->[$i1]->[0] eq $dataB->[$i2]->[0]) && ($dataA->[$i1]->[1] lt $dataB->[$i2]->[1])))
121 :     {
122 :     my $currS = $dataA->[$i1]->[0];
123 :     my $currR = $dataA->[$i1]->[1];
124 :     my $pegs1 = [];
125 :     while (($i1 < @$dataA) && ($dataA->[$i1]->[0] eq $currS) && ($dataA->[$i1]->[1] eq $currR))
126 :     {
127 :     push(@$pegs1,$dataA->[$i1]->[2]);
128 :     $i1++;
129 :     }
130 : overbeek 1.5 push(@$in1_not2,[$currS,$currR,$pegs1,[$fig_services->seqs_with_role($currR,undef,$genomeB)]]);
131 : overbeek 1.1 }
132 :     elsif ( ($dataA->[$i1]->[0] gt $dataB->[$i2]->[0]) ||
133 :     (($dataA->[$i1]->[0] eq $dataB->[$i2]->[0]) && ($dataA->[$i1]->[1] gt $dataB->[$i2]->[1])))
134 :     {
135 :     $i2++;
136 :     }
137 :     else
138 :     {
139 : overbeek 1.2 my $currS = $dataA->[$i1]->[0];
140 :     my $currR = $dataA->[$i1]->[1];
141 :     while (($i1 < @$dataA) && ($dataA->[$i1]->[0] eq $currS) && ($dataA->[$i1]->[1] eq $currR))
142 :     {
143 :     $i1++;
144 :     }
145 :    
146 :     while (($i2 < @$dataB) && ($dataB->[$i2]->[0] eq $currS) && ($dataB->[$i2]->[1] eq $currR))
147 :     {
148 :     $i2++;
149 :     }
150 : overbeek 1.1 }
151 :     }
152 :     return $in1_not2;
153 :     }
154 :    
155 :     sub data_on_subs_roles_pegs {
156 :     my($fig_services,$genome) = @_;
157 :    
158 :     my $sub_data = $fig_services->get_genome_subsystem_data($genome);
159 : dsouza 1.4 return [sort { ($a->[0] cmp $b->[0]) or ($a->[1] cmp $b->[1]) or ($a->[3] <=> $b->[3]) } map {$_->[2] =~ /\.(\d+)$/; [@$_, $1]} @$sub_data];
160 :    
161 :     # sort by_fig_id inefficient for metagenomes with large numbers of features
162 :     # return [sort { ($a->[0] cmp $b->[0]) or ($a->[1] cmp $b->[1]) or &FIG::by_fig_id($a->[2],$b->[2]) } @$sub_data];
163 : overbeek 1.1 }
164 :    
165 :     1;
166 :    

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3