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

Annotation of /FigKernelScripts/split_sequences_into_sets.pl

Parent Directory Parent Directory | Revision Log Revision Log


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

1 : overbeek 1.1 ########################################################################
2 :     # -*- perl -*-
3 :     #
4 :     # Copyright (c) 2003-2006 University of Chicago and Fellowship
5 :     # for Interpretations of Genomes. All Rights Reserved.
6 :     #
7 :     # This file is part of the SEED Toolkit.
8 :     #
9 :     # The SEED Toolkit is free software. You can redistribute
10 :     # it and/or modify it under the terms of the SEED Toolkit
11 :     # Public License.
12 :     #
13 :     # You should have received a copy of the SEED Toolkit Public License
14 :     # along with this program; if not write to the University of Chicago
15 :     # at info@ci.uchicago.edu or the Fellowship for Interpretation of
16 :     # Genomes at veronika@thefig.info or download a copy from
17 :     # http://www.theseed.org/LICENSE.TXT.
18 :     #
19 :    
20 :     use FIG;
21 :     $min_sz = 30;
22 :     $fullF = "/tmp/seqs.$$";
23 :    
24 :     $usage = "usage: split_sequences_into_sets OutputDirectory [SimilarityCutoff FracCoverage] [blastout=BlastFile] < full_seqs";
25 :    
26 :     for ($i= $#ARGV; ($i >= 0); $i--)
27 :     {
28 :     if ($ARGV[$i] =~ /blastout=(\S+)$/)
29 :     {
30 :     $blastout = $1;
31 :     splice(@ARGV,$i,1);
32 :     }
33 :     }
34 :    
35 :     (
36 :     ($output = shift @ARGV)
37 :     )
38 :     || die $usage;
39 :    
40 :     (! -e $output) || die "$output already exists -- check it (and delete it, if that is what you really wish)";
41 :     &FIG::verify_dir($output);
42 :    
43 :     open(TMP,">$fullF") || die "could not open $fullF";
44 :     $/ = "\n>";
45 :     while (defined($_ = <STDIN>))
46 :     {
47 :     chomp;
48 :     if ($_ =~ /^>?(\S+)[^\n]*\n(.*)/s)
49 :     {
50 :     $id = $1;
51 :     $seq = $2;
52 :     $seq =~ s/\s//gs;
53 :     $seq =~ s/[uU]/x/g;
54 :     $seq_of{$id} = $seq;
55 :     print TMP ">$id\n$seq\n";
56 :     }
57 :     }
58 :     close(TMP);
59 :     $/ = "\n";
60 :    
61 :     if (@ARGV == 2)
62 :     {
63 :     $sim_cutoff = shift @ARGV;
64 :     $frac_cov = shift @ARGV;
65 :     }
66 :     else
67 :     {
68 :     $sim_cutoff = 1.0e-30;
69 :     $frac_cov = 0.70;
70 :     }
71 :    
72 :     &run("$FIG_Config::ext_bin/formatdb -i $fullF -p T");
73 :    
74 :     open(BLASTOUT,"$FIG_Config::ext_bin/blastall -i $fullF -d $fullF -FF -m 8 -p blastp -b 20000 -v 20000 |")
75 :     || die "could not run blastall -i $fullF -d $fullF -FF -m 8 -p blastp -b 20000 -v 20000";
76 :     my @all = ();
77 :     while (defined($_ = <BLASTOUT>))
78 :     {
79 :     if (($_ =~ /^(\S+)\t(\S+)\t(\S+\t){4}(\S+)\t(\S+)\t(\S+)\t(\S+)\s+(\S+)/) &&
80 :     ($1 ne $2) &&
81 :     ($8 <= $sim_cutoff) &&
82 :     ((@all == 0) || (($1 ne $all[$#all]->[0]) || ($2 ne $all[$#all]->[1]))))
83 :     {
84 :     push(@all,[$1,$2,$4,$5,$6,$7,$8]);
85 :     }
86 :     }
87 :     close(BLASTOUT);
88 :    
89 :     if (defined($blastout))
90 :     {
91 :     open(BLASTOUT,">$blastout") || die "could not open $blastout";
92 :     foreach $_ (@all)
93 :     {
94 :     print BLASTOUT join("\t",(@$_,length($seq_of{$_->[0]}),length($seq_of{$_->[1]}))),"\n";
95 :     }
96 :     close(BLASTOUT);
97 :     }
98 :    
99 :     @sims = grep { ($_->[6] <= $sim_cutoff) }
100 :     grep { ($_->[0] ne $_->[1]) &&
101 :     ((($_->[3] - $_->[2]) / length($seq_of{$_->[0]})) >= $frac_cov) &&
102 :     ((($_->[5] - $_->[4]) / length($seq_of{$_->[1]})) >= $frac_cov)
103 :     }
104 :     @all;
105 :    
106 :     if (@sims == 0)
107 :     {
108 :     print STDERR "No similarities worth processing:\n"
109 :     , join("\n", map { join(", ", @ { $_ }) } @all), "\n"
110 :     , "--- deleting $output\n";
111 :     system("rm -fR $output") && warn "Could not delete directory $output";
112 :     exit(0);
113 :     }
114 :    
115 :     foreach $sim (@sims)
116 :     {
117 :     if (! $seen{"$sim->[0]\t$sim->[1]"}) # pick only the best similarity between two ids
118 :     {
119 :     push(@{$conn{$sim->[0]}},$sim);
120 :     $seen{"$sim->[0]\t$sim->[1]"} = 1;
121 :     }
122 :     }
123 :     undef %seen;
124 :    
125 :     if (0) # This is just for debugging -- display the similarities
126 :     {
127 :     foreach $id (sort keys(%seq_of))
128 :     {
129 :     &print_sims($id,\@sims);
130 :     }
131 :     }
132 :    
133 :     # @ids = sort { (@{$conn{$b}} <=> @{$conn{$a}}) or (&cov($conn{$b}) <=> &cov($conn{$a})) }
134 :     # sort by number of connections; within that on sum of coverage
135 :    
136 :     @ids = sort { &cov($conn{$b}) <=> &cov($conn{$a}) }
137 :     keys(%seq_of);
138 :     # sort by sum of coverage
139 :    
140 :     @tmp = map { $x = $_; [$x,scalar @{$conn{$x}},&cov($conn{$x})] } @ids;
141 :     $n = 1;
142 :     @sets = ();
143 :    
144 :     open(SZ,">$output/set.sizes") || die "could not open $output/set.sizes";
145 :    
146 :     foreach $central_id (@ids)
147 :     {
148 :     if (! $seen{$central_id})
149 :     {
150 :     # print STDERR "========\n";
151 :     # print STDERR "picking central_id = $central_id\n";
152 :    
153 :     $set = &extract_set($central_id,$conn{$central_id},\%seq_of,\%seen);
154 :     # print STDERR "========\n";
155 :     # print STDERR "calculated set:", &Dumper($et);
156 :    
157 :     print SZ join("\t",($n,scalar @$set)),"\n";
158 :    
159 :     push(@sets,[$n,$set]);
160 :    
161 :     open(OUT,">$output/$n") || die "could not open $output/$n";
162 :     $n++;
163 :     foreach $tuple (@$set)
164 :     {
165 :     ($id,$b,$e,$seq) = @$tuple;
166 :     $seen{$id} = 1;
167 :     print OUT ">$id\n$seq\n";
168 :     }
169 :     close(OUT);
170 :     }
171 :     }
172 :     close(SZ);
173 :    
174 :     @distances = &construct_distance_matrix(\@sets,\@all); # entries are [n1,n2,sc]
175 :     open(DIST,">$output/distance.matrix\n") || die "could not open $output/distance.matrix";
176 :     foreach $tuple (@distances)
177 :     {
178 :     print DIST join("\t",@$tuple),"\n";
179 :     }
180 :     close(DIST);
181 :    
182 :    
183 :     system("rm $fullF\*");
184 :    
185 :     sub possible_region {
186 :     my($r1,$r2) = @_;
187 :    
188 :     $r2 = $r2 ? $r2 : [1,100000];
189 :    
190 :     my($b1,$e1) = @$r1;
191 :     my($b2,$e2) = @$r2;
192 :    
193 :     my $b = ($b1 > $b2) ? $b1 : $b2;
194 :     my $e = ($e1 < $e2) ? $e1 : $e2;
195 :    
196 :     if ((($e-$b)+1) >= $min_sz)
197 :     {
198 :     return [$b,$e];
199 :     }
200 :     return undef;
201 :     }
202 :    
203 :     sub cov {
204 :     my($sims) = @_;
205 :     my($x,$n);
206 :    
207 :     $n = 0;
208 :     foreach $x (@$sims)
209 :     {
210 :     $n += ($x->[3] - $x->[2]);
211 :     }
212 :     return $n;
213 :     }
214 :    
215 :     sub run {
216 :     my($cmd) = @_;
217 :    
218 :     # my @tmp = `date`; chomp @tmp; print STDERR "$tmp[0]: running $cmd\n";
219 :     (system($cmd) == 0) || confess "FAILED: $cmd";
220 :     }
221 :    
222 :     sub print_sims {
223 :     my($id,$sims) = @_;
224 :     my($sim,@tmp);
225 :    
226 :     print STDERR "$id\n";
227 :     foreach $sim (@$sims)
228 :     {
229 :     if ($sim->[0] eq $id)
230 :     {
231 :     @tmp = @$sim;
232 :     $tmp[0] = "";
233 :     print STDERR join("\t",@tmp),"\n";
234 :     }
235 :     }
236 :     }
237 :    
238 :     sub extract_set {
239 :     my($central_id,$sims,$seq_of,$seen) = @_;
240 :     my($sim,$r1,$id1,$id2,$b1,$e1,$b2,$e2);
241 :    
242 :     my(@tuples) = ();
243 :    
244 :     my @sims = sort { (($b->[3] - $b->[2]) <=> ($a->[3] - $a->[2])) } @$sims;
245 :     my $region = [1,length($seq_of->{$central_id})];
246 :    
247 :     foreach $sim (@sims)
248 :     {
249 :     ($id1,$id2,$b1,$e1,$b2,$e2) = @$sim;
250 :     if ((! $seen->{$id2}) && ($r1 = &possible_region([$b1,$e1],$region)))
251 :     {
252 :     $region = $r1;
253 :     }
254 :     }
255 :    
256 :     my($b1F,$e1F) = @$region;
257 :     my $seq = substr($seq_of->{$central_id},$b1F-1,($e1F-$b1F)+1);
258 :     my $core_ln = length($seq);
259 :    
260 :     my $singleF = "/tmp/seq.$$";
261 :    
262 :     open(TMP,">$singleF") || die "could not open $singleF";
263 :     print TMP ">$central_id\n$seq\n";
264 :     close(TMP);
265 :    
266 :     @sims = grep { ($_->[6] <= 1.0e-5) && (($_->[3] - $_->[2]) > (0.8 * $core_ln)) }
267 :     map { $_ =~ /^(\S+)\t(\S+)\t(\S+\t){4}(\S+)\t(\S+)\t(\S+)\t(\S+)\s+(\S+)/;
268 :     [$1,$2,$4,$5,$6,$7,$8] }
269 :     `$FIG_Config::ext_bin/blastall -i $singleF -d $fullF -FF -m 8 -p blastp -b 20000 -v 20000`;
270 :     system "rm $singleF\*";
271 :     undef %seen1;
272 :    
273 :     foreach $sim (@sims)
274 :     {
275 :     ($id1,$id2,$b1,$e1,$b2,$e2) = @$sim;
276 :     if ((! $seen->{$id2}) && (! $seen1{$id2}) && ((($e2-$b2)+1) >= $min_sz))
277 :     {
278 :     $seq = $seq_of{$id2};
279 :     ## $seq = substr($seq_of{$id2},$b2-1,($e2-$b2)+1); <<<
280 :     push(@tuples,[$id2,$b2,$e2,$seq]);
281 :     $seen1{$id2} = 1;
282 :     }
283 :     }
284 :     return [@tuples];
285 :     }
286 :    
287 :     sub construct_distance_matrix {
288 :     my($sets,$sims) = @_;
289 :     my($tuple1,$tuple2,$n,$set,$id1,$id2,$sc,$sim,%best,$key,$n1,$n2);
290 :     my(@distances,$sofar);
291 :    
292 :     foreach $tuple1 (@$sets)
293 :     {
294 :     ($n,$set) = @$tuple1;
295 :     foreach $tuple2 (@$set)
296 :     {
297 :     ($id,$b,$e,$seq) = @$tuple2;
298 :     $in{$id} = $n;
299 :     }
300 :     }
301 :     foreach $sim (@$sims)
302 :     {
303 :     ($id1,$id2,undef,undef,undef,undef,$sc) = @$sim;
304 :     if ($in{$id1} && $in{$id2} && ($in{$id1} != $in{$id2}))
305 :     {
306 :     $key = join("\t",($in{$id1},$in{$id2}));
307 :     if ((! ($sofar = $best{$key})) || ($sofar > $sc))
308 :     {
309 :     $best{$key} = $sc;
310 :     }
311 :    
312 :     $key = join("\t",($in{$id2},$in{$id1}));
313 :     if ((! ($sofar = $best{$key})) || ($sofar > $sc))
314 :     {
315 :     $best{$key} = $sc;
316 :     }
317 :     }
318 :     }
319 :    
320 :     @distances = ();
321 :     foreach $key (keys(%best))
322 :     {
323 :     ($n1,$n2) = split(/\t/,$key);
324 :     push(@distances,[$n1,$n2,$best{$key}]);
325 :     }
326 :     return sort { ($a->[0] <=> $b->[0]) or ($a->[1] <=> $b->[1]) } @distances;
327 :     }

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3