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

Annotation of /FigKernelPackages/CompareMR.pm

Parent Directory Parent Directory | Revision Log Revision Log


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

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3