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

Annotation of /FigKernelScripts/trim_sequences.pl

Parent Directory Parent Directory | Revision Log Revision Log


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

1 : overbeek 1.7 ########################################################################
2 : olson 1.6 #
3 :     # Copyright (c) 2003-2006 University of Chicago and Fellowship
4 :     # for Interpretations of Genomes. All Rights Reserved.
5 :     #
6 :     # This file is part of the SEED Toolkit.
7 :     #
8 :     # The SEED Toolkit is free software. You can redistribute
9 :     # it and/or modify it under the terms of the SEED Toolkit
10 :     # Public License.
11 :     #
12 :     # You should have received a copy of the SEED Toolkit Public License
13 :     # along with this program; if not write to the University of Chicago
14 :     # at info@ci.uchicago.edu or the Fellowship for Interpretation of
15 :     # Genomes at veronika@thefig.info or download a copy from
16 :     # http://www.theseed.org/LICENSE.TXT.
17 :     #
18 :    
19 : overbeek 1.7 use strict;
20 : overbeek 1.1
21 : overbeek 1.7 # usage: trim_sequences [-pad_left=nL] [-pad_right=nR] [SimilarityCutoff FracCoverage] < full_seqs > trimmed_seqs
22 :    
23 :     my $sim_cutoff = 1.0e-30;
24 :     my $frac_cov = 0.80;
25 :     my $pad_left = 0;
26 :     my $pad_right = 0;
27 :    
28 :     while ((@ARGV > 0) && ($ARGV[0] =~ /^-/))
29 :     {
30 :     my $arg = shift @ARGV;
31 :     if ($arg =~ /^-pad_left=(\d+)/) { $pad_left = $1 }
32 :     elsif ($arg =~ /^-pad_right=(\d+)/) { $pad_right = $1 }
33 :     else { die "bad argument $arg" }
34 :     }
35 :    
36 :     if (@ARGV == 2)
37 :     {
38 :     ($sim_cutoff,$frac_cov) = @ARGV;
39 :     }
40 : overbeek 1.1
41 :     &make_sure("formatdb","blastall");
42 :    
43 : overbeek 1.7 my $min_sz = 30;
44 :     my $fullF = "/tmp/seqs.$$";
45 :     my $singleF = "/tmp/seq.$$";
46 :    
47 :     my %seq_of;
48 : overbeek 1.9 my %ln;
49 : overbeek 1.3
50 : overbeek 1.1 open(TMP,">$fullF") || die "could not open $fullF";
51 :     $/ = "\n>";
52 :     while (defined($_ = <STDIN>))
53 :     {
54 :     chomp;
55 :     if ($_ =~ /^>?(\S+)[^\n]*\n(.*)/s)
56 :     {
57 : overbeek 1.7 my $id = $1;
58 :     my $seq = $2;
59 : overbeek 1.1 $seq =~ s/\s//gs;
60 :     $seq_of{$id} = $seq;
61 : overbeek 1.9 $ln{$id} = length($seq);
62 : overbeek 1.1 print TMP ">$id\n$seq\n";
63 :     }
64 :     }
65 :     close(TMP);
66 :     $/ = "\n";
67 :    
68 :     &run("formatdb -i $fullF -p T");
69 : overbeek 1.11 my @sims = grep { ($_->[0] ne $_->[1]) && ($_->[6] <= $sim_cutoff)
70 :     # && ((($_->[3] - $_->[2]) / $ln{$_->[0]}) >= $frac_cov)
71 : overbeek 1.7 }
72 :    
73 :     map { $_ =~ /^(\S+)\t(\S+)\t(\S+\t){4}(\S+)\t(\S+)\t(\S+)\t(\S+)\s+(\S+)/;
74 :     [$1,$2,$4,$5,$6,$7,$8]
75 :     }
76 : overbeek 1.11 `blastall -i $fullF -d $fullF -FF -m 8 -p blastp -b 400000`;
77 : overbeek 1.1
78 :     if (@sims == 0)
79 :     {
80 :     print STDERR "No similarities worth processing\n";
81 :     exit(1);
82 :     }
83 :    
84 : overbeek 1.7 my %seen;
85 : overbeek 1.9 my @sims1 = ();
86 : overbeek 1.7 foreach my $sim (@sims)
87 : overbeek 1.1 {
88 :     if (! $seen{"$sim->[0]\t$sim->[1]"}) # pick only the best similarity between two ids
89 :     {
90 : overbeek 1.9 push(@sims1,$sim);
91 : overbeek 1.1 $seen{"$sim->[0]\t$sim->[1]"} = 1;
92 :     }
93 :     }
94 :    
95 : overbeek 1.9 my %counts;
96 :     foreach my $sim (@sims1)
97 :     {
98 :     my($id1,$id2,$b1,$e1,$b2,$e2,$sc) = @$sim;
99 :     if (((($e1 - $b1)+1) / $ln{$id1}) >= $frac_cov) { $counts{$id1}++ }
100 :     }
101 : overbeek 1.1
102 : overbeek 1.9 my @ids = sort { $counts{$b} <=> $counts{$a} } keys(%counts);
103 :     my $i;
104 : overbeek 1.10 for ($i=1; ($i < @ids) && &close_ln($ids[0],$ids[$i],\%ln) &&
105 :     &close_count($ids[0],$ids[1],\%counts); $i++) {}
106 : overbeek 1.1
107 : overbeek 1.9 if ($i < @ids) { $#ids = $i-1 }
108 :     foreach my $id (@ids)
109 : overbeek 1.1 {
110 : overbeek 1.9 # print STDERR "$id\t$counts{$id}\n";
111 : overbeek 1.3 }
112 :    
113 : overbeek 1.9 my %short = map { $_ => 1 } @ids;
114 :     my @short_sims = grep { $short{$_->[1]} } @sims1;
115 : overbeek 1.3
116 : overbeek 1.9 my %starts;
117 :     my %ends;
118 : overbeek 1.3
119 : overbeek 1.9 foreach my $sim (@short_sims)
120 :     {
121 :     my($id1,$id2,$b1,$e1,$b2,$e2,$sc) = @$sim;
122 :     push(@{$starts{$id1}},$b1);
123 :     push(@{$ends{$id1}},$e1);
124 : overbeek 1.1 }
125 :    
126 : overbeek 1.9 foreach my $id (keys(%starts))
127 :     {
128 :     my @x = sort { $a <=> $b } @{$starts{$id}};
129 :     my $b1 = $x[int(@x/2)];
130 :     $b1 = ($b1 > $pad_left) ? $b1 - $pad_left : 1;
131 :    
132 :     @x = sort { $a <=> $b } @{$ends{$id}};
133 :     my $e1 = $x[int(@x/2)];
134 :     $e1 = ($e1 < ($ln{$id} - $pad_right)) ? $e1 + $pad_right : $ln{$id};
135 :     my $s = substr($seq_of{$id},($b1-1),($e1+1-$b1));
136 :     print ">$id $b1-$e1\n$s\n";
137 : overbeek 1.1 }
138 :    
139 : overbeek 1.9 system("rm $fullF\*");
140 : overbeek 1.1
141 :     sub run {
142 :     my($cmd) = @_;
143 :    
144 :     # my @tmp = `date`; chomp @tmp; print STDERR "$tmp[0]: running $cmd\n";
145 :     (system($cmd) == 0) || confess "FAILED: $cmd";
146 :     }
147 :    
148 :     sub make_sure {
149 :     my(@progs) = @_;
150 :    
151 :     my $prog;
152 :     foreach $prog (@progs)
153 :     {
154 :     my @tmp = `which $prog`;
155 :     if ($tmp[0] =~ /^no $prog/)
156 :     {
157 :     print STDERR $tmp[0];
158 :     exit(1);
159 :     }
160 :     }
161 :     }
162 :    
163 : overbeek 1.10 sub close_ln {
164 :     my($id1,$id2,$lns) = @_;
165 :    
166 :     return (abs(($lns->{$id1} - $lns->{$id2})/$lns->{$id1}) < 0.8);
167 :     }
168 :    
169 :     sub close_count {
170 :     my($id1,$id2,$counts) = @_;
171 :    
172 :     return (abs(($counts->{$id1} - $counts->{$id2})/$counts->{$id1}) < 0.8);
173 :     }

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3