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

Annotation of /FigKernelPackages/CompareMR.pm

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1.6 - (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 :     use FIGV;
28 :    
29 :     sub compare_genomes_MR {
30 : paarmann 1.3 my($arg1,$genome2, $fig_services, $genome1) = @_;
31 : overbeek 1.1
32 : paarmann 1.3 unless ($fig_services) {
33 :     if (-d "$FIG_Config::organisms/$arg1")
34 :     {
35 : parrello 1.6 $fig_services = FIG->new();
36 : paarmann 1.3 $genome1 = $arg1;
37 :     }
38 :     else
39 : overbeek 1.1 {
40 : paarmann 1.3 if (-d $arg1)
41 : overbeek 1.1 {
42 : parrello 1.6 my $figV = new FIGV($arg1,FIG->new());
43 : paarmann 1.3 $fig_services = $figV;
44 :     if ($arg1 =~ /(\d+\.\d+)$/)
45 :     {
46 :     $genome1 = $1;
47 :     }
48 :     else
49 :     {
50 :     die "invalid first argument: $arg1";
51 :     }
52 : overbeek 1.1 }
53 : paarmann 1.3 else
54 : overbeek 1.1 {
55 : paarmann 1.3 die "Invalid first argument; $arg1";
56 : overbeek 1.1 }
57 :     }
58 :     }
59 :    
60 :     my $data1 = &data_on_subs_roles_pegs($fig_services,$genome1);
61 :     my $data2 = &data_on_subs_roles_pegs($fig_services,$genome2);
62 :    
63 :     my $common = &process_common($fig_services,$genome1,$genome2,$data1,$data2);
64 :     my $in1_not2 = &process_diff($fig_services,$genome1,$genome2,$data1,$data2);
65 :     my $in2_not1 = &process_diff($fig_services,$genome2,$genome1,$data2,$data1);
66 :     return ($common,$in1_not2,$in2_not1);
67 :     }
68 :    
69 :     sub process_common {
70 :     my($fig_services,$genome1,$genome2,$data1,$data2) = @_;
71 :    
72 :     my $common = [];
73 :     my $i1 = 0;
74 :     my $i2 = 0;
75 :     while (($i1 < @$data1) && ($i2 < @$data2))
76 :     {
77 :     if ( ($data1->[$i1]->[0] lt $data2->[$i2]->[0]) ||
78 :     (($data1->[$i1]->[0] eq $data2->[$i2]->[0]) && ($data1->[$i1]->[1] lt $data2->[$i2]->[1])))
79 :     {
80 :     $i1++;
81 :     }
82 :     elsif ( ($data1->[$i1]->[0] gt $data2->[$i2]->[0]) ||
83 :     (($data1->[$i1]->[0] eq $data2->[$i2]->[0]) && ($data1->[$i1]->[1] gt $data2->[$i2]->[1])))
84 :     {
85 :     $i2++;
86 :     }
87 :     else
88 :     {
89 :     my $currS = $data1->[$i1]->[0];
90 :     my $currR = $data1->[$i1]->[1];
91 :     my $pegs1 = [];
92 :     while (($i1 < @$data1) && ($data1->[$i1]->[0] eq $currS) && ($data1->[$i1]->[1] eq $currR))
93 :     {
94 :     push(@$pegs1,$data1->[$i1]->[2]);
95 :     $i1++;
96 :     }
97 :    
98 :     my $pegs2 = [];
99 :     while (($i2 < @$data2) && ($data2->[$i2]->[0] eq $currS) && ($data2->[$i2]->[1] eq $currR))
100 :     {
101 :     push(@$pegs2,$data2->[$i2]->[2]);
102 :     $i2++;
103 :     }
104 :     push(@$common,[$currS,$currR,$pegs1,$pegs2]);
105 :     }
106 :     }
107 :     return $common;
108 :     }
109 :    
110 :     sub process_diff {
111 :     my($fig_services,$genomeA,$genomeB,$dataA,$dataB) = @_;
112 :    
113 :     my $in1_not2 = [];
114 :     my $i1 = 0;
115 :     my $i2 = 0;
116 :     while ($i1 < @$dataA)
117 :     {
118 : overbeek 1.2 if ( ($i2 == @$dataB) || ($dataA->[$i1]->[0] lt $dataB->[$i2]->[0]) ||
119 : overbeek 1.1 (($dataA->[$i1]->[0] eq $dataB->[$i2]->[0]) && ($dataA->[$i1]->[1] lt $dataB->[$i2]->[1])))
120 :     {
121 :     my $currS = $dataA->[$i1]->[0];
122 :     my $currR = $dataA->[$i1]->[1];
123 :     my $pegs1 = [];
124 :     while (($i1 < @$dataA) && ($dataA->[$i1]->[0] eq $currS) && ($dataA->[$i1]->[1] eq $currR))
125 :     {
126 :     push(@$pegs1,$dataA->[$i1]->[2]);
127 :     $i1++;
128 :     }
129 : overbeek 1.5 push(@$in1_not2,[$currS,$currR,$pegs1,[$fig_services->seqs_with_role($currR,undef,$genomeB)]]);
130 : overbeek 1.1 }
131 :     elsif ( ($dataA->[$i1]->[0] gt $dataB->[$i2]->[0]) ||
132 :     (($dataA->[$i1]->[0] eq $dataB->[$i2]->[0]) && ($dataA->[$i1]->[1] gt $dataB->[$i2]->[1])))
133 :     {
134 :     $i2++;
135 :     }
136 :     else
137 :     {
138 : overbeek 1.2 my $currS = $dataA->[$i1]->[0];
139 :     my $currR = $dataA->[$i1]->[1];
140 :     while (($i1 < @$dataA) && ($dataA->[$i1]->[0] eq $currS) && ($dataA->[$i1]->[1] eq $currR))
141 :     {
142 :     $i1++;
143 :     }
144 :    
145 :     while (($i2 < @$dataB) && ($dataB->[$i2]->[0] eq $currS) && ($dataB->[$i2]->[1] eq $currR))
146 :     {
147 :     $i2++;
148 :     }
149 : overbeek 1.1 }
150 :     }
151 :     return $in1_not2;
152 :     }
153 :    
154 :     sub data_on_subs_roles_pegs {
155 :     my($fig_services,$genome) = @_;
156 :    
157 :     my $sub_data = $fig_services->get_genome_subsystem_data($genome);
158 : 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];
159 :    
160 :     # sort by_fig_id inefficient for metagenomes with large numbers of features
161 :     # return [sort { ($a->[0] cmp $b->[0]) or ($a->[1] cmp $b->[1]) or &FIG::by_fig_id($a->[2],$b->[2]) } @$sub_data];
162 : overbeek 1.1 }
163 :    
164 :     1;
165 :    

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3