[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.3 - (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.2 my $set = &extract_set($central_id,\%conn,\%seen);
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.2 my($central_id,$conn,$seen) = @_;
126 : overbeek 1.1
127 : overbeek 1.2 my $cluster = [$central_id];
128 :     my %in = ($central_id,1);
129 : overbeek 1.1
130 : overbeek 1.2 my $i;
131 :     while ($i < @$cluster)
132 : overbeek 1.1 {
133 : overbeek 1.2 my $x = $conn->{$cluster->[$i]};
134 :     $i++;
135 : overbeek 1.3 foreach my $peg (@$x)
136 : overbeek 1.1 {
137 : overbeek 1.2 if (! $in{$peg})
138 :     {
139 :     push(@$cluster,$peg);
140 :     $in{$peg} = 1;
141 :     }
142 : overbeek 1.1 }
143 :     }
144 : overbeek 1.2 return $cluster;
145 : overbeek 1.1 }
146 :    
147 :     sub construct_distance_matrix {
148 :     my($sets,$sims) = @_;
149 : overbeek 1.2 my($tuple,$n,$set,$id1,$id2,$sc,$sim,%best,$key,$n1,$n2);
150 :     my(@distances,$sofar,%in);
151 : overbeek 1.1
152 : overbeek 1.2 foreach my $tuple (@$sets)
153 : overbeek 1.1 {
154 : overbeek 1.2 my($n,$set) = @$tuple;
155 :     foreach my $id (@$set)
156 : overbeek 1.1 {
157 :     $in{$id} = $n;
158 :     }
159 :     }
160 : overbeek 1.2
161 :     my %best;
162 :     foreach my $sim (@$sims)
163 : overbeek 1.1 {
164 : overbeek 1.3 my($id1,$id2,undef,undef,$sc) = @$sim;
165 : overbeek 1.1 if ($in{$id1} && $in{$id2} && ($in{$id1} != $in{$id2}))
166 :     {
167 : overbeek 1.2 my $key = join("\t",($in{$id1},$in{$id2}));
168 : overbeek 1.3 if ((! ($sofar = $best{$key})) || ($sofar < $sc))
169 : overbeek 1.1 {
170 :     $best{$key} = $sc;
171 :     }
172 :    
173 :     $key = join("\t",($in{$id2},$in{$id1}));
174 :     if ((! ($sofar = $best{$key})) || ($sofar > $sc))
175 :     {
176 :     $best{$key} = $sc;
177 :     }
178 :     }
179 :     }
180 :    
181 : overbeek 1.2 my @distances = ();
182 :     foreach my $key (keys(%best))
183 : overbeek 1.1 {
184 : overbeek 1.2 my($n1,$n2) = split(/\t/,$key);
185 : overbeek 1.1 push(@distances,[$n1,$n2,$best{$key}]);
186 :     }
187 :     return sort { ($a->[0] <=> $b->[0]) or ($a->[1] <=> $b->[1]) } @distances;
188 :     }

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3