[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.5 - (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 : overbeek 1.3 my @all = map { my($id1,$id2,undef,undef,undef,undef,$b1,$e1,$b2,$e2,$psc,undef,$ln1,$ln2) = @$_;
57 :     [$id1,$id2,(($e1+1)-$b1)/$ln1,(($e2+1)-$b2)/$ln2,$psc]
58 :     } &ProtSims::blastP(\@seqs,\@seqs,$min_size);
59 : overbeek 1.1
60 :    
61 : overbeek 1.3 my @sims = grep { my($id1,$id2,$frac1,$frac2,$psc) = @$_;
62 :     ($id1 ne $id2) && ($frac1 >= $frac_cov) && ($frac2 >= $frac_cov) } @all;
63 : overbeek 1.1
64 : overbeek 1.2 my(%seen,%conn);
65 :     foreach my $sim (@sims)
66 : overbeek 1.1 {
67 : overbeek 1.3 my($id1,$id2,$frac1,$frac2,$psc) = @$sim;
68 :     my $key = join("\t",sort ($id1,$id2));
69 :     if (! $seen{$key})
70 :     {
71 :     $seen{$key} = 1;
72 :     push(@{$conn{$id1}},$id2);
73 :     push(@{$conn{$id2}},$id1);
74 : overbeek 1.1 }
75 :     }
76 :     undef %seen;
77 : overbeek 1.3 my @ids = sort { my $x = $conn{$b} ? @{$conn{$b}} : 0;
78 :     my $y = $conn{$a} ? @{$conn{$a}} : 0;
79 :     ($x <=> $y)
80 :     } keys(%seq_of);
81 : overbeek 1.1
82 : overbeek 1.2 my $n = 1;
83 :     my @sets = ();
84 : overbeek 1.1
85 :     open(SZ,">$output/set.sizes") || die "could not open $output/set.sizes";
86 :    
87 : overbeek 1.2 foreach my $central_id (@ids)
88 : overbeek 1.1 {
89 :     if (! $seen{$central_id})
90 :     {
91 :    
92 : overbeek 1.4 my $set = &extract_set($central_id,\%conn,\%seen,\%seq_of);
93 : overbeek 1.1 print SZ join("\t",($n,scalar @$set)),"\n";
94 :    
95 :     push(@sets,[$n,$set]);
96 :    
97 :     open(OUT,">$output/$n") || die "could not open $output/$n";
98 :     $n++;
99 : overbeek 1.2 foreach my $peg (@$set)
100 : overbeek 1.1 {
101 : overbeek 1.2 $seen{$peg} = 1;
102 :     print OUT ">$peg\n$seq_of{$peg}\n";
103 : overbeek 1.1 }
104 :     close(OUT);
105 :     }
106 :     }
107 :     close(SZ);
108 :    
109 : overbeek 1.2 my @distances = &construct_distance_matrix(\@sets,\@all); # entries are [n1,n2,sc]
110 : overbeek 1.1 open(DIST,">$output/distance.matrix\n") || die "could not open $output/distance.matrix";
111 : overbeek 1.2 foreach my $tuple (@distances)
112 : overbeek 1.1 {
113 :     print DIST join("\t",@$tuple),"\n";
114 :     }
115 :     close(DIST);
116 :    
117 :     sub run {
118 :     my($cmd) = @_;
119 :    
120 :     # my @tmp = `date`; chomp @tmp; print STDERR "$tmp[0]: running $cmd\n";
121 :     (system($cmd) == 0) || confess "FAILED: $cmd";
122 :     }
123 :    
124 :     sub extract_set {
125 : overbeek 1.4 my($central_id,$conn,$seen,$seq_of) = @_;
126 : overbeek 1.1
127 : overbeek 1.4 my $center_length = length($seq_of->{$central_id});
128 : overbeek 1.2 my $cluster = [$central_id];
129 :     my %in = ($central_id,1);
130 : overbeek 1.1
131 : overbeek 1.2 my $i;
132 :     while ($i < @$cluster)
133 : overbeek 1.1 {
134 : overbeek 1.2 my $x = $conn->{$cluster->[$i]};
135 :     $i++;
136 : overbeek 1.3 foreach my $peg (@$x)
137 : overbeek 1.1 {
138 : overbeek 1.5 if ((! $in{$peg}) && (! $seen->{$peg}))
139 : overbeek 1.2 {
140 : overbeek 1.4 my $len = length($seq_of->{$peg});
141 :     if ((abs($center_length - $len) / $center_length) <= (1 - $frac_cov))
142 :     {
143 :     push(@$cluster,$peg);
144 :     $in{$peg} = 1;
145 :     }
146 : overbeek 1.2 }
147 : overbeek 1.1 }
148 :     }
149 : overbeek 1.2 return $cluster;
150 : overbeek 1.1 }
151 :    
152 :     sub construct_distance_matrix {
153 :     my($sets,$sims) = @_;
154 : overbeek 1.2 my($tuple,$n,$set,$id1,$id2,$sc,$sim,%best,$key,$n1,$n2);
155 :     my(@distances,$sofar,%in);
156 : overbeek 1.1
157 : overbeek 1.2 foreach my $tuple (@$sets)
158 : overbeek 1.1 {
159 : overbeek 1.2 my($n,$set) = @$tuple;
160 :     foreach my $id (@$set)
161 : overbeek 1.1 {
162 :     $in{$id} = $n;
163 :     }
164 :     }
165 : overbeek 1.2
166 :     my %best;
167 :     foreach my $sim (@$sims)
168 : overbeek 1.1 {
169 : overbeek 1.3 my($id1,$id2,undef,undef,$sc) = @$sim;
170 : overbeek 1.1 if ($in{$id1} && $in{$id2} && ($in{$id1} != $in{$id2}))
171 :     {
172 : overbeek 1.2 my $key = join("\t",($in{$id1},$in{$id2}));
173 : overbeek 1.3 if ((! ($sofar = $best{$key})) || ($sofar < $sc))
174 : overbeek 1.1 {
175 :     $best{$key} = $sc;
176 :     }
177 :    
178 :     $key = join("\t",($in{$id2},$in{$id1}));
179 :     if ((! ($sofar = $best{$key})) || ($sofar > $sc))
180 :     {
181 :     $best{$key} = $sc;
182 :     }
183 :     }
184 :     }
185 :    
186 : overbeek 1.2 my @distances = ();
187 :     foreach my $key (keys(%best))
188 : overbeek 1.1 {
189 : overbeek 1.2 my($n1,$n2) = split(/\t/,$key);
190 : overbeek 1.1 push(@distances,[$n1,$n2,$best{$key}]);
191 :     }
192 :     return sort { ($a->[0] <=> $b->[0]) or ($a->[1] <=> $b->[1]) } @distances;
193 :     }

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3