[Bio] / FigKernelScripts / p3x-update-propagated-family-names.pl Repository:
ViewVC logotype

Annotation of /FigKernelScripts/p3x-update-propagated-family-names.pl

Parent Directory Parent Directory | Revision Log Revision Log


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

1 : olson 1.1 #
2 :     # Read the output logs from p3x-propagate-family-names runs and apply the renames
3 :     # to the input families file.
4 :     #
5 :     # We will read all logs (one for global propagation, one for each genus-specific
6 :     # local propagation) and cache the data, and apply the changes to the input families
7 :     # file(s).
8 :     #
9 :     # Emit a record of the family splits and merges.
10 :     #
11 :    
12 :     use strict;
13 :     use Getopt::Long::Descriptive;
14 :     use Data::Dumper;
15 :     use IO::Handle;
16 :     use File::Basename;
17 :    
18 :     my($opt, $usage) = describe_options("%c %o [log-files] ::: [family-files]",
19 :     ["output-directory|o=s" => "If specified write output files to this directory named using the name of the input family file. Otherwise write to standard output"],
20 :     ["max-id-file=s" => "File containing family max used ID values (required)"],
21 :     ["help|h" => "Show this help message"],
22 :     );
23 :    
24 :     print($usage->text), exit 1 if $opt->help;
25 :     die($usage->text) if @ARGV < 3 || !$opt->max_id_file;
26 :    
27 :     #
28 :     # Find logfiles before the ::: indicator and family files after the ::: indicator.
29 :     #
30 :     my @log_files;
31 :     my @fam_files;
32 :     my $list = \@log_files;
33 :     for my $ent (@ARGV)
34 :     {
35 :     if ($ent eq ':::')
36 :     {
37 :     if ($list == \@fam_files)
38 :     {
39 :     die "Only one instance of ::: allowed in the argument list\n";
40 :     }
41 :     else
42 :     {
43 :     $list = \@fam_files;
44 :     }
45 :     next;
46 :     }
47 :     if (! -f $ent)
48 :     {
49 :     die "File $ent not readable\n";
50 :     }
51 :     push(@$list, $ent);
52 :     }
53 :     @log_files or die "At least one logfile must be provided\n";
54 :     @fam_files or die "At least one family file must be provided\n";
55 :    
56 :     main();
57 :    
58 :     sub main
59 :     {
60 :     my %lf_id_map;
61 :     my %gf_id_map;
62 :     my @splits;
63 :     my @joins;
64 :    
65 :     #
66 :     # Read the max ID file; increment each by one to create the next id to be assigned
67 :     # for our global and local families.
68 :     #
69 :    
70 :     my %next_id;
71 :     my %max_new_id;
72 :     open(MID, "<", $opt->max_id_file) or die "Max id file " . $opt->max_id_file . " cannot be opened: $!";
73 :     while (<MID>)
74 :     {
75 :     chomp;
76 :     my($fam, $val) = split(/\t/);
77 :     $next_id{$fam} = 1 + $val;
78 :     }
79 :     close(MID);
80 :    
81 :     for my $lf (@log_files)
82 :     {
83 :     read_log_file($lf, \%lf_id_map, \%gf_id_map, \@splits, \@joins, \%next_id, \%max_new_id);
84 :     }
85 :    
86 :     # print Dumper(\@splits, \@joins);
87 :     # print Dumper(\%lf_id_map, \%gf_id_map);
88 :    
89 :     my %fam_fun;
90 :     for my $ff (@fam_files)
91 :     {
92 :     process_fam_file($ff, $opt->output_directory, \%lf_id_map, \%gf_id_map, \%fam_fun, \%next_id);
93 :     }
94 :     }
95 :    
96 :     sub process_fam_file
97 :     {
98 :     my($file, $out_dir, $lf_id_map, $gf_id_map, $fam_fun, $next_id) = @_;
99 :    
100 :     open(F, "<", $file) or die "Cannot open $file: $!";
101 :     my $out_fh = \*STDOUT;
102 :     if ($out_dir)
103 :     {
104 :     my $out_file = $out_dir . "/" . basename($file);
105 :     $out_fh = IO::Handle->new();
106 :     open($out_fh, ">", $out_file) or die "Cannot write $out_file: $!";
107 :     }
108 :    
109 :     while (<F>)
110 :     {
111 :     chomp;
112 :     my($gf, $n, $ng, $peg, $len, $fun, $fam, $genus) = split(/\t/);
113 :     my $lf = "$genus.$fam";
114 :    
115 :     my $ngf = $gf_id_map->{$gf};
116 :     my $nlf = $lf_id_map->{$lf};
117 :     $ngf or die "Cannot map $gf\n";
118 :    
119 :     #
120 :     # If there is no local family mapping, see if this is a new genus.
121 :     # A new genus will not have an entry in the $next_id hash that maintains the
122 :     # next id to be assigned for a local family. In this case we just carry across
123 :     # the newly created ID from the incoming family file.
124 :     #
125 :    
126 :     if (!$nlf)
127 :     {
128 :     if (exists($next_id->{$genus}))
129 :     {
130 :     die "Cannot map $lf\n";
131 :     }
132 :     else
133 :     {
134 :     $nlf = $lf;
135 :     }
136 :     }
137 :     $fam_fun->{$ngf} = $fam_fun->{$nlf} = $fun;
138 :     my($ngenus, $nfam) = $nlf =~ /^(.*?)\.(\d+)$/;
139 :     print $out_fh join("\t", $ngf, $n, $ng, $peg, $len, $fun, $nfam, $ngenus, $nfam), "\n";
140 :     }
141 :    
142 :     if ($out_dir)
143 :     {
144 :     close($out_fh);
145 :     }
146 :     }
147 :    
148 :     sub read_log_file
149 :     {
150 :     my($file, $lf_id_map, $gf_id_map, $splits, $joins, $next_id, $max_new_id) = @_;
151 :    
152 :     open(F, "<", $file) or die "Cannot open $file: $!";
153 :    
154 :     my $find_unmapped;
155 :     while (<F>)
156 :     {
157 :     chomp;
158 :     if (/^Unmapped\s+new/)
159 :     {
160 :     $find_unmapped++;
161 :     last;
162 :     }
163 :     elsif (/^(.*?)\s+NOW\s+(.*?)(\s+weight=.*)?$/)
164 :     {
165 :     my($new, $old) = ($1, $2);
166 :     my $map;
167 :     my $genus;
168 :     if ($new =~ /^(.*?)\.(\d+)/)
169 :     {
170 :     $genus = $1;
171 :     $map = $lf_id_map;
172 :     }
173 :     else
174 :     {
175 :     $genus = "GLOBAL";
176 :     $map = $gf_id_map;
177 :     }
178 :     # print "'$_' '$new' '$old' '$genus'\n";
179 :     my $old_upd = $old;
180 :     if ($old =~ /NEW_(\d+)/)
181 :     {
182 :     my $val = $next_id->{$genus} + $1;
183 :     $max_new_id->{$genus} = $1 if $1 > $max_new_id->{$genus};
184 :     if ($genus eq 'GLOBAL')
185 :     {
186 :     $old_upd = sprintf("GF%08d", $val);
187 :     }
188 :     else
189 :     {
190 :     $old_upd = "$genus.$val";
191 :     }
192 :     }
193 :    
194 :     defined($map->{$new}) && die "map entry for $new already exists with $old";
195 :     # print "MAP '$new' '$old_upd'\t'$_'\n";
196 :     $map->{$new} = $old_upd;
197 :     }
198 :     elsif (/^SPLIT\s+O\s+(.*?)\s+=>\s+N\s+(.*)$/)
199 :     {
200 :     my @newids;
201 :     my $oid = $1;
202 :     my $newids = $2;
203 :     if ($oid =~ /^(.*?)\.(\d+)$/)
204 :     {
205 :     my $g = $1;
206 :     @newids = $newids =~ /($g\.\d+)/mg;
207 :     }
208 :     else
209 :     {
210 :     @newids = split(/\s+/, $2);
211 :     }
212 :     push(@$splits, [$oid, [@newids]]);
213 :     }
214 :     elsif (/^JOIN\s+(.*?)\s+=>\s+(.*)$/)
215 :     {
216 :     my @oids = split(/\s+/, $1);
217 :     my $nid = $2;
218 :     push(@$joins, [[@oids], $nid]);
219 :     }
220 :     else
221 :     {
222 :     warn "Unmatched: $_\n";
223 :     }
224 :     }
225 :    
226 :     if ($find_unmapped)
227 :     {
228 :     my %next_new;
229 :     for my $genus (keys %$max_new_id)
230 :     {
231 :     $next_new{$genus} = 1 + $max_new_id->{$genus};
232 :     }
233 :     while (<F>)
234 :     {
235 :     if (/^\t([^\t]+)\t/)
236 :     {
237 :     my $new = $1;
238 :     #
239 :     # This is a new id that was not mapped to an old ID; we need to assign a new
240 :     # family id.
241 :     #
242 :     my $map;
243 :     my $old;
244 :     my $genus;
245 :     if ($new =~ /^(.*)\.(\d+)$/)
246 :     {
247 :     $genus = $1;
248 :     my $id = $next_new{$genus}++;
249 :     $map = $lf_id_map;
250 :     $old = "$genus.$id";
251 :     }
252 :     else
253 :     {
254 :     $genus = "GLOBAL";
255 :     my $id = $next_new{$genus}++;
256 :     $map = $gf_id_map;
257 :     $old = sprintf("GF%08d", $id);
258 :     }
259 :     defined($map->{$new}) && die "map entry for $new already exists with $old";
260 :     $map->{$new} = $old;
261 :     }
262 :     }
263 :     }
264 :     }
265 :    
266 :    
267 :    
268 :     sub is_localfam_id
269 :     {
270 :     my($id) = @_;
271 :     return $id =~ /^\S+\.\d+$/;
272 :     }

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3