[Bio] / FigKernelScripts / get_merged_families_3.pl Repository:
ViewVC logotype

Annotation of /FigKernelScripts/get_merged_families_3.pl

Parent Directory Parent Directory | Revision Log Revision Log


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

1 : olson 1.1
2 :     use strict;
3 :     use Data::Dumper;
4 :     use Getopt::Long::Descriptive;
5 :     use File::Path 'make_path';
6 :     use File::Slurp;
7 :     use Statistics::Descriptive;
8 :     use File::Basename;
9 :     use gjoseqlib;
10 :     use IPC::Run;
11 :     use Proc::ParallelLoop;
12 :     use Graph;
13 :     use DB_File;
14 :     #
15 :     # Use the output of the MCL runs to compute the sets of local families to be merged, and generate
16 :     # the global families.
17 :     #
18 :    
19 :    
20 :     my($opt, $usage) = describe_options("%c %o kmer-dir genus-data-dir merge-dir inflation",
21 :     ["genera-from=s", "Merge genera from the given file"],
22 :     ["help|h", "Show this help message"]);
23 :    
24 :     print($usage->text), exit if $opt->help;
25 :     die($usage->text) if @ARGV != 4;
26 :    
27 :     my $kmer_dir = shift;
28 :     my $genus_data_dir = shift;
29 :     my $merge_dir = shift;
30 :     my $inflation = shift;
31 :    
32 :     my($fn_to_id, $id_to_fn) = read_function_index($kmer_dir);
33 :    
34 :     my @genera;
35 :     if ($opt->genera_from)
36 :     {
37 :     open(F, "<", $opt->genera_from) or die "Cannot open " . $opt->genera_from . ": $!\n";
38 :     while (<F>)
39 :     {
40 :     chomp;
41 :     s/^\s*//;
42 :     s/\s*$//;
43 :     if (! -d "$genus_data_dir/$_")
44 :     {
45 :     die "Genus directory $merge_dir/$_ does not exist\n";
46 :     }
47 :     push(@genera, $_);
48 :     }
49 :     close(F);
50 :     }
51 :     else
52 :     {
53 :     my @fams = <$genus_data_dir/*/families.all>;
54 :     @genera = map { basename(dirname($_)) } @fams;
55 :     }
56 :     my %genera = map { $_ => 1 } @genera;
57 :    
58 :     #
59 :     # Load genus family data.
60 :     #
61 :    
62 :     print STDERR "Read genus families\n";
63 :     my %local_fams;
64 :     for my $genus (@genera)
65 :     {
66 :     open(F, "<", "$genus_data_dir/$genus/families.all") or die "Cannot open $genus_data_dir/$genus/families.all: $!";
67 :     my $last;
68 :     my $lastfun;
69 :     my $pegs = [];
70 :     while (<F>)
71 :     {
72 :     chomp;
73 :     my($fam, $fun, $lfam, $peg, $len) = split(/\t/);
74 :     if ($fam ne $last)
75 :     {
76 :     if ($lastfun)
77 :     {
78 :     # my $fidx = $fn_to_id->{$lastfun};
79 :     # defined($fidx) or die "Cannot map $fun to an index\n";
80 :     $local_fams{$genus, $last} = $pegs;
81 :     }
82 :     $pegs = [];
83 :     $last = $fam;
84 :     $lastfun = $fun;
85 :     }
86 :     push(@$pegs, [$peg, $len, $fun, $fam]);
87 :     }
88 :     if ($lastfun)
89 :     {
90 :     # my $fidx = $fn_to_id->{$lastfun};
91 :     # defined($fidx) or die "Cannot map $fun to an index\n";
92 :     $local_fams{$genus, $last} = $pegs;
93 :     }
94 :     close(F);
95 :     }
96 :    
97 :     print STDERR "Genus families loaded\n";
98 :    
99 :     #
100 :     # Read peg mapping file
101 :     #
102 :    
103 :     my %peginfo;
104 :    
105 :     if (tie %peginfo, 'DB_File', "$merge_dir/peg.map.btree", O_RDONLY, 0, $DB_BTREE)
106 :     {
107 :     print STDERR "Tied to btree\n";
108 :     }
109 :     else
110 :     {
111 :     print STDERR "Reading peg.map\n";
112 :    
113 :     open(PM, "<", "$merge_dir/peg.map") or die "Cannot open $merge_dir/peg.map: $!";
114 :     while (<PM>)
115 :     {
116 :     chomp;
117 :     my @a = split(/\t/);
118 :     $peginfo{$a[0]} = \@a;
119 :     }
120 :     close(PM);
121 :     }
122 :    
123 :     print STDERR "Peg mapping loaded\n";
124 :    
125 :     #
126 :     # For each family in the merge_dir/mcl directory, use the peg map to map to
127 :     # determine the set of local families to merge.
128 :     #
129 :    
130 :     my %merged;
131 :    
132 :     my %merged_in_fam;
133 :    
134 :     my $gfam = 'GF00000000';
135 :     for my $mcl (<$merge_dir/mcl/$inflation/*>)
136 :     {
137 :     my $fam = basename($mcl);
138 :     next unless $fam =~ /^\d+$/;
139 :     open(F, "<", $mcl) or die "Cannot open $mcl: $!";
140 :     print STDERR "$mcl $gfam\n";
141 :    
142 :     #
143 :     # Read MCL families and map to the local family identifiers they join.
144 :     # Construct a graph where vertices are local family ids and the initial paths
145 :     # are the connections formed by MCL.
146 :     # After reading the MCL connections for a function, the connected subgraphs
147 :     # form the merged families. In other words, we're both merging across local
148 :     # families based on MCL, and merging MCL clusters based on the contents of
149 :     # the local families.
150 :     #
151 :    
152 :     my $graph = Graph::Undirected->new();
153 :     while (<F>)
154 :     {
155 :     chomp;
156 :     my @pegs = split(/\t/);
157 :     next unless @pegs > 1;
158 :     my $fam = [];
159 :     my %lf;
160 :     for my $peg (@pegs)
161 :     {
162 :     my $pi = $peginfo{$peg};
163 :     if (!defined($pi))
164 :     {
165 :     warn "No peginfo for $pi\n";
166 :     next;
167 :     }
168 :     my($rep, $genus, $fam, $fidx, $fun);
169 :     if (ref($pi))
170 :     {
171 :     ($rep, $genus, $fam, $fidx, $fun) = @$pi;
172 :     }
173 :     else
174 :     {
175 :     ($rep, $genus, $fam, $fidx, $fun) = split(/\t/, $pi);
176 :     }
177 :    
178 :     my $lfamid = join($;, $genus, $fam);
179 :     $lf{$lfamid} = 1;
180 :     }
181 :     $graph->add_path(keys %lf);
182 :     }
183 :     close(F);
184 :    
185 :     # print STDERR Dumper($graph->connected_components);
186 :     for my $group ($graph->connected_components)
187 :     {
188 :     # print STDERR "merge\t@$group\n";
189 :     my %gmerge;
190 :     my $found;
191 :     for my $lf (@$group)
192 :     {
193 :     my($genus, $kgfam) = split(/$;/, $lf);
194 :     $gmerge{$genus}++;
195 :     }
196 :    
197 :     my $ng = keys %gmerge;
198 :     my $nf = @$group;
199 :    
200 :     for my $lf (@$group)
201 :     {
202 :     my($genus, $kgfam) = split(/$;/, $lf);
203 :     my $pl = delete $local_fams{$lf};
204 :    
205 :     if (ref($pl))
206 :     {
207 :     for my $pi (@$pl)
208 :     {
209 :     my($peg, $len, $fun, $fam) = @$pi;
210 :     print join("\t", $gfam, $nf, $ng, $peg, $len, $fun, $fam, $genus, $kgfam), "\n";
211 :     }
212 :     $merged_in_fam{$lf} = $gfam;
213 :     $found++;
214 :     }
215 :     else
216 :     {
217 :     warn "No local fam for $lf ($merged_in_fam{$lf})\n";
218 :     }
219 :     }
220 :     $gfam++ if $found;
221 :     }
222 :     next;
223 :    
224 :     while (<F>)
225 :     {
226 :     chomp;
227 :     my @pegs = split(/\t/);
228 :     next unless @pegs > 1;
229 :     my %merge_these;
230 :     my %gmerge;
231 :     for my $peg (@pegs)
232 :     {
233 :     my $pi = $peginfo{$peg};
234 :     if (!ref($pi))
235 :     {
236 :     warn "No peginfo for $pi\n";
237 :     next;
238 :     }
239 :     # print $map_fh "$rep\t$genus\t$fam\t$fun_idx\t$fun\t$gname->{$gid}\n";
240 :     my($rep, $genus, $fam, $fidx, $fun) = @$pi;
241 :     next unless $genera{$genus};
242 :     $gmerge{$genus}++;
243 :     $merge_these{$genus, $fam} = 1;
244 :     }
245 :     my $n = keys %merge_these;
246 :     if ($n > 1)
247 :     {
248 :     print STDERR "Merging:\n";
249 :     for my $kk (keys %merge_these)
250 :     {
251 :     my($genus, $fam) = split(/$;/, $kk);
252 :    
253 :     my $lfams = ref($local_fams{$kk}) ? join(" ", map { $_->[0] } @{$local_fams{$kk}}) : 'NA';
254 :    
255 :     print STDERR " $genus\t$fam\t$merged_in_fam{$kk}\t$lfams\n";
256 :     }
257 :     my $found;
258 :     for my $k (keys %merge_these)
259 :     {
260 :     my($genus, $kgfam) = split(/$;/, $k);
261 :     my $pl = delete $local_fams{$k};
262 :     my $ng = keys %gmerge;
263 :     if (ref($pl))
264 :     {
265 :     for my $pi (@$pl)
266 :     {
267 :     my($peg, $len, $fun, $fam) = @$pi;
268 :     print join("\t", $gfam, $n, $ng, $peg, $len, $fun, $fam, $genus, $kgfam), "\n";
269 :     }
270 :     $merged_in_fam{$k} = $gfam;
271 :     $found++;
272 :     }
273 :     else
274 :     {
275 :     warn "No local fam for $k ($merged_in_fam{$k})\n";
276 :     }
277 :     }
278 :     $gfam++ if $found;
279 :     }
280 :     }
281 :     }
282 :    
283 :     for my $lk (sort keys %local_fams)
284 :     {
285 :     my $pl = $local_fams{$lk};
286 :     my($genus, $kgfam) = split(/$;/, $lk);
287 :    
288 :     for my $pi (@$pl)
289 :     {
290 :     my($peg, $len, $fun, $fam) = @$pi;
291 :     print join("\t", $gfam, 1, 1, $peg, $len, $fun, $fam, $genus, $kgfam), "\n";
292 :     }
293 :     $gfam++;
294 :     }
295 :    
296 :     sub read_function_index
297 :     {
298 :     my($dir) = @_;
299 :    
300 :     my $idx = {};
301 :     my $from_id = [];
302 :    
303 :     open(F, "<", "$dir/function.index") or die "Cannot open $dir/function.index:$ !";
304 :     while (<F>)
305 :     {
306 :     chomp;
307 :     my($i, $fun) = split(/\t/);
308 :     $idx->{$fun} = $i;
309 :     $from_id->[$i] = $fun;
310 :     }
311 :     close(F);
312 :    
313 :     return($idx, $from_id);
314 :     }

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3