[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.2 - (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 : overbeek 1.2 use strict;
21 : overbeek 1.1 use FIG;
22 :    
23 : overbeek 1.2 use ProtSims;
24 :     use gjoseqlib;
25 : overbeek 1.1
26 : overbeek 1.2 my $min_sz = 30;
27 :     my $fullF = "/tmp/seqs.$$";
28 :    
29 :     my $usage = "usage: split_sequences_into_sets OutputDirectory [SimilarityCutoff FracCoverage] < full_seqs";
30 : overbeek 1.1
31 : overbeek 1.2 my $output;
32 : overbeek 1.1 (
33 :     ($output = shift @ARGV)
34 :     )
35 :     || die $usage;
36 :    
37 :     (! -e $output) || die "$output already exists -- check it (and delete it, if that is what you really wish)";
38 :     &FIG::verify_dir($output);
39 :    
40 : overbeek 1.2 my @seqs = &gjoseqlib::read_fasta();
41 :     my %seq_of = map { $_->[0] => $_->[2] } @seqs;
42 : overbeek 1.1
43 : overbeek 1.2 my($sim_cutoff,$frac_cov);
44 : overbeek 1.1 if (@ARGV == 2)
45 :     {
46 :     $sim_cutoff = shift @ARGV;
47 :     $frac_cov = shift @ARGV;
48 :     }
49 :     else
50 :     {
51 :     $sim_cutoff = 1.0e-30;
52 :     $frac_cov = 0.70;
53 :     }
54 :    
55 : overbeek 1.2 my $min_size = (@seqs > 20) ? (@seqs / 10) : (@seqs / 2);
56 :     my @all = &ProtSims::blastP(\@seqs,\@seqs,$min_size);
57 :     # &print_sims(\@all);
58 : overbeek 1.1
59 : overbeek 1.2 # $_ = @all; print STDERR "computed $_ similarities\n";
60 : overbeek 1.1
61 : overbeek 1.2 my @sims = grep { ($_->id1 ne $_->id2) &&
62 :     ($_->psc <= $sim_cutoff) &&
63 :     ((($_->e1 - $_->b1) / $_->ln1) >= $frac_cov) &&
64 :     ((($_->e2 - $_->b2) / $_->ln2) >= $frac_cov)
65 :     } @all;
66 : overbeek 1.1
67 : overbeek 1.2 # print STDERR "=======\n";
68 :     # &print_sims(\@sims);
69 : overbeek 1.1
70 : overbeek 1.2 my(%seen,%conn);
71 :     foreach my $sim (@sims)
72 : overbeek 1.1 {
73 : overbeek 1.2 if (! $seen{"$sim->id1\t$sim->id2"}) # pick only the best similarity between two ids
74 : overbeek 1.1 {
75 : overbeek 1.2 push(@{$conn{$sim->id1}},$sim);
76 :     $seen{"$sim->id1\t$sim->id2"} = 1;
77 : overbeek 1.1 }
78 :     }
79 :     undef %seen;
80 : overbeek 1.2 # print STDERR "built connections\n";
81 : overbeek 1.1
82 : overbeek 1.2 my @ids = sort { &cov($conn{$b}) <=> &cov($conn{$a}) }
83 :     keys(%seq_of);
84 : overbeek 1.1 # sort by sum of coverage
85 :    
86 : overbeek 1.2 my $n = 1;
87 :     my @sets = ();
88 : overbeek 1.1
89 :     open(SZ,">$output/set.sizes") || die "could not open $output/set.sizes";
90 :    
91 : overbeek 1.2 foreach my $central_id (@ids)
92 : overbeek 1.1 {
93 :     if (! $seen{$central_id})
94 :     {
95 :    
96 : overbeek 1.2 my $set = &extract_set($central_id,\%conn,\%seen);
97 : overbeek 1.1 print SZ join("\t",($n,scalar @$set)),"\n";
98 :    
99 :     push(@sets,[$n,$set]);
100 :    
101 :     open(OUT,">$output/$n") || die "could not open $output/$n";
102 :     $n++;
103 : overbeek 1.2 foreach my $peg (@$set)
104 : overbeek 1.1 {
105 : overbeek 1.2 $seen{$peg} = 1;
106 :     print OUT ">$peg\n$seq_of{$peg}\n";
107 : overbeek 1.1 }
108 :     close(OUT);
109 :     }
110 :     }
111 :     close(SZ);
112 : overbeek 1.2 # print STDERR "built set sizes\n";
113 : overbeek 1.1
114 : overbeek 1.2 my @distances = &construct_distance_matrix(\@sets,\@all); # entries are [n1,n2,sc]
115 : overbeek 1.1 open(DIST,">$output/distance.matrix\n") || die "could not open $output/distance.matrix";
116 : overbeek 1.2 foreach my $tuple (@distances)
117 : overbeek 1.1 {
118 :     print DIST join("\t",@$tuple),"\n";
119 :     }
120 :     close(DIST);
121 : overbeek 1.2 # print STDERR "built distance matrix\n";
122 : overbeek 1.1
123 : overbeek 1.2 sub print_sims {
124 :     my($sims) = @_;
125 :    
126 :     foreach my $sim (@$sims)
127 : overbeek 1.1 {
128 : overbeek 1.2 print STDERR join("\t",($sim->id1,$sim->id2,$sim->iden,$sim->b1,$sim->e1,$sim->b2,$sim->e2,$sim->psc,$sim->ln1,$sim->ln2)),"\n";
129 : overbeek 1.1 }
130 :     }
131 :    
132 :     sub cov {
133 :     my($sims) = @_;
134 :     my($x,$n);
135 :    
136 :     $n = 0;
137 :     foreach $x (@$sims)
138 :     {
139 : overbeek 1.2 $n += ($x->e1 - $x->b1);
140 : overbeek 1.1 }
141 :     return $n;
142 :     }
143 :    
144 :     sub run {
145 :     my($cmd) = @_;
146 :    
147 :     # my @tmp = `date`; chomp @tmp; print STDERR "$tmp[0]: running $cmd\n";
148 :     (system($cmd) == 0) || confess "FAILED: $cmd";
149 :     }
150 :    
151 :     sub extract_set {
152 : overbeek 1.2 my($central_id,$conn,$seen) = @_;
153 : overbeek 1.1
154 : overbeek 1.2 my $cluster = [$central_id];
155 :     my %in = ($central_id,1);
156 : overbeek 1.1
157 : overbeek 1.2 my $i;
158 :     while ($i < @$cluster)
159 : overbeek 1.1 {
160 : overbeek 1.2 my $x = $conn->{$cluster->[$i]};
161 :     $i++;
162 :     foreach my $peg (map { $_->id2} @$x)
163 : overbeek 1.1 {
164 : overbeek 1.2 if (! $in{$peg})
165 :     {
166 :     push(@$cluster,$peg);
167 :     $in{$peg} = 1;
168 :     }
169 : overbeek 1.1 }
170 :     }
171 : overbeek 1.2 return $cluster;
172 : overbeek 1.1 }
173 :    
174 :     sub construct_distance_matrix {
175 :     my($sets,$sims) = @_;
176 : overbeek 1.2 my($tuple,$n,$set,$id1,$id2,$sc,$sim,%best,$key,$n1,$n2);
177 :     my(@distances,$sofar,%in);
178 : overbeek 1.1
179 : overbeek 1.2 foreach my $tuple (@$sets)
180 : overbeek 1.1 {
181 : overbeek 1.2 my($n,$set) = @$tuple;
182 :     foreach my $id (@$set)
183 : overbeek 1.1 {
184 :     $in{$id} = $n;
185 :     }
186 :     }
187 : overbeek 1.2
188 :     my %best;
189 :     foreach my $sim (@$sims)
190 : overbeek 1.1 {
191 : overbeek 1.2 my $id1 = $sim->id1;
192 :     my $id2 = $sim->id2;
193 :     my $sc = $sim->psc;
194 : overbeek 1.1 if ($in{$id1} && $in{$id2} && ($in{$id1} != $in{$id2}))
195 :     {
196 : overbeek 1.2 my $key = join("\t",($in{$id1},$in{$id2}));
197 : overbeek 1.1 if ((! ($sofar = $best{$key})) || ($sofar > $sc))
198 :     {
199 :     $best{$key} = $sc;
200 :     }
201 :    
202 :     $key = join("\t",($in{$id2},$in{$id1}));
203 :     if ((! ($sofar = $best{$key})) || ($sofar > $sc))
204 :     {
205 :     $best{$key} = $sc;
206 :     }
207 :     }
208 :     }
209 :    
210 : overbeek 1.2 my @distances = ();
211 :     foreach my $key (keys(%best))
212 : overbeek 1.1 {
213 : overbeek 1.2 my($n1,$n2) = split(/\t/,$key);
214 : overbeek 1.1 push(@distances,[$n1,$n2,$best{$key}]);
215 :     }
216 :     return sort { ($a->[0] <=> $b->[0]) or ($a->[1] <=> $b->[1]) } @distances;
217 :     }

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3