[Bio] / FigKernelPackages / ProtSims.pm Repository:
ViewVC logotype

Annotation of /FigKernelPackages/ProtSims.pm

Parent Directory Parent Directory | Revision Log Revision Log


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

1 : overbeek 1.1
2 : olson 1.10 # This is a SAS component
3 :    
4 : overbeek 1.1 #
5 :     # Copyright (c) 2003-2006 University of Chicago and Fellowship
6 :     # for Interpretations of Genomes. All Rights Reserved.
7 :     #
8 :     # This file is part of the SEED Toolkit.
9 :     #
10 :     # The SEED Toolkit is free software. You can redistribute
11 :     # it and/or modify it under the terms of the SEED Toolkit
12 :     # Public License.
13 :     #
14 :     # You should have received a copy of the SEED Toolkit Public License
15 :     # along with this program; if not write to the University of Chicago
16 :     # at info@ci.uchicago.edu or the Fellowship for Interpretation of
17 :     # Genomes at veronika@thefig.info or download a copy from
18 :     # http://www.theseed.org/LICENSE.TXT.
19 :     #
20 :     package ProtSims;
21 :    
22 :     use strict;
23 : olson 1.4 use Sim;
24 : overbeek 1.1 use Data::Dumper;
25 :     use Carp;
26 :     use gjoseqlib;
27 :    
28 : olson 1.13 use SeedAware;
29 :    
30 :     my $blat_cmd = SeedAware::executable_for("blat");
31 :     my $blastall_cmd = SeedAware::executable_for("blastall");
32 :     my $formatdb_cmd = SeedAware::executable_for("formatdb");
33 : olson 1.10
34 : overbeek 1.1 sub blastP {
35 : overbeek 1.12 my($q,$db,$min_hits,$use_blast) = @_;
36 : overbeek 1.1
37 : overbeek 1.2 my(@q,@db);
38 :    
39 : olson 1.13 my $tmp_dir = SeedAware::location_of_tmp();
40 :    
41 : overbeek 1.1 my $qF;
42 :     if (ref $q)
43 :     {
44 : olson 1.13 $qF = "$tmp_dir/query.$$";
45 : overbeek 1.1 &gjoseqlib::print_alignment_as_fasta($qF,$q);
46 :     }
47 :     else
48 :     {
49 :     $qF = $q;
50 :     if (-e $qF)
51 :     {
52 : overbeek 1.2 @q = &gjoseqlib::read_fasta($qF);
53 : overbeek 1.1 $q = \@q;
54 :     }
55 :     else
56 :     {
57 :     die "$qF is missing";
58 :     }
59 :     }
60 :    
61 :     my $dbF;
62 :     if (ref $db)
63 :     {
64 : olson 1.13 $dbF = "$tmp_dir/db.$$";
65 : overbeek 1.1 &gjoseqlib::print_alignment_as_fasta($dbF,$db);
66 : olson 1.13 system($formatdb_cmd, '-i', $dbF);
67 : overbeek 1.1 }
68 :     else
69 :     {
70 :     $dbF = $db;
71 :     if (-e $dbF)
72 :     {
73 :     if (! ((-e "$dbF.psq") && ((-M "$dbF.psq") < (-M $dbF))))
74 :     {
75 : olson 1.13 system($formatdb_cmd, '-i', $dbF);;
76 : overbeek 1.1 }
77 : overbeek 1.2 @db = &gjoseqlib::read_fasta($dbF);
78 : overbeek 1.1 $db = \@db;
79 :     }
80 :     else
81 :     {
82 :     die "$dbF is missing";
83 :     }
84 :     }
85 :    
86 : olson 1.13 my $tmpF = "$tmp_dir/sim.out.$$";
87 :    
88 :     if ($use_blast)
89 :     {
90 :     open(my $fh, ">", $tmpF);
91 :     close($fh);
92 :     }
93 :     else
94 : olson 1.4 {
95 : olson 1.13 my @cmd = ($blat_cmd, $dbF, $qF, "-prot", "-out=blast8", $tmpF);
96 :    
97 :     my $rc = system_with_redirect(\@cmd, { stdout => '/dev/null' } );
98 :     print STDERR "Blat returns $rc: @cmd\n";
99 :     if ($rc != 0)
100 :     {
101 :     warn "Blat run failed with rc=$rc\n";
102 :     open(my $fh, ">", $tmpF);
103 :     close($fh);
104 :     }
105 : olson 1.4 }
106 : overbeek 1.1
107 : olson 1.13 # my $cmd = $use_blast ? "touch $tmpF" : "$blat_cmd $dbF $qF -prot -out=blast8 $tmpF > /dev/null";
108 :     # #print STDERR "$cmd\n";
109 :     # my $rc = system $cmd;
110 :     # if ($rc != 0)
111 :     # {
112 :     # die "ProtSims::blastP: blat run failed with rc=$rc: $cmd\n";
113 :     # }
114 :    
115 : overbeek 1.5 my @sims1 = ();
116 :     open(BLAT,"<$tmpF") || die "could not open $tmpF";
117 :     my $sim = <BLAT>;
118 :     while ($sim && ($sim=~ /^(\S+\t\S+)/))
119 :     {
120 :     my $pegs = $1;
121 :     my @pieces = ();
122 :     while ($sim && ($sim=~ /^((\S+\t\S+).*\S)/) && ($2 eq $pegs))
123 :     {
124 :     push(@pieces,[split(/\t/,$1)]);
125 :     $sim = <BLAT>;
126 :     }
127 :     push(@sims1,&condense(\@pieces));
128 :     }
129 :     close(BLAT);
130 : overbeek 1.1 unlink $tmpF;
131 : olson 1.4 #print STDERR &Dumper(sims1 => \@sims1);
132 : overbeek 1.1
133 :     my @rerun = ();
134 :     my @sims = ();
135 :     my $qI = 0;
136 :     my $simI = 0;
137 :     while ($qI < @$q)
138 :     {
139 :     my $peg = $q->[$qI]->[0];
140 : olson 1.4 #print STDERR "processing $peg\n";
141 : overbeek 1.1 my $hits = 0;
142 :     my $simI1 = $simI;
143 :     while (($simI1 < @sims1) && ($sims1[$simI1]->[0] eq $peg))
144 :     {
145 :     if (($simI1 == $#sims1) || ($sims1[$simI1]->[1] ne $sims1[$simI1+1]->[1]))
146 :     {
147 :     $hits++;
148 :     }
149 :     $simI1++;
150 :     }
151 : olson 1.4 #print STDERR "hits=$hits min_hits=$min_hits\n";
152 : overbeek 1.1 if ($hits >= $min_hits)
153 :     {
154 : olson 1.4 push(@sims,@sims1[$simI..$simI1-1]);
155 : overbeek 1.1 }
156 :     else
157 :     {
158 :     push(@rerun,$q->[$qI]);
159 :     }
160 :     $simI = $simI1;
161 :     $qI++;
162 :     }
163 : olson 1.4 #print STDERR &Dumper(rerun => \@rerun);
164 : overbeek 1.1
165 :     if (@rerun > 0)
166 :     {
167 : olson 1.13 my $tmpQ = "$tmp_dir/tmpQ.$$";
168 : overbeek 1.1 &gjoseqlib::print_alignment_as_fasta($tmpQ,\@rerun);
169 : olson 1.13 # my $cmd = "$blastall_cmd -m 8 -i $tmpQ -d $dbF -FF -p blastp -e 1e-5";
170 :     my @cmd = ($blastall_cmd, '-m', 8, '-i', $tmpQ, '-d', $dbF, '-FF', '-p', 'blastp', '-e', '1e-5');
171 : olson 1.4 #print STDERR "$cmd\n";
172 : olson 1.13 #open(BL, "$cmd|") or die "ProtSims::blastP: pipe to blast failed with $!: $cmd\n";
173 :    
174 : olson 1.14 #
175 :     # It'd be nice to do this but windows doesn't support it.
176 :     #
177 :     #open(BL, "-|", @cmd) or die "ProtSims::blastP: pipe to blast failed with $!: @cmd\n";
178 :     my $out_tmp = "$tmp_dir/blast_out.$$";
179 :     push(@cmd, "-o", $out_tmp);
180 :     print STDERR "@cmd\n";
181 :     my $rc = system(@cmd);
182 :     if ($rc != 0)
183 :     {
184 :     warn "Error rc=$rc running blast: @cmd\n";
185 :     }
186 :     else
187 :     {
188 :     if (open(BL, "<", $out_tmp))
189 :     {
190 :     while (<BL>)
191 :     {
192 :     chomp;
193 :     push @sims, [ split(/\s+/, $_) ];
194 :     }
195 :     #my @blastout = map { chomp; [split(/\s+/,$_)] } <BL>;
196 :     # push(@sims,@blastout);
197 :     close(BL);
198 :     }
199 :     }
200 :     unlink $out_tmp;
201 : overbeek 1.1 unlink $tmpQ;
202 :     }
203 :    
204 : overbeek 1.3 my %lnQ = map { $_->[0] => length($_->[2]) } @$q;
205 :     my %lnDB = map { $_->[0] => length($_->[2]) } @$db;
206 : overbeek 1.1 @sims = map { push(@$_,$lnQ{$_->[0]},$lnDB{$_->[1]}); bless($_,'Sim') } @sims;
207 :    
208 : olson 1.13 if ($qF eq "$tmp_dir/query.$$") { unlink $qF }
209 :     if ($dbF eq "$tmp_dir/db.$$") { unlink($dbF,"$dbF.psq","$dbF.pin","$dbF.phr") }
210 : overbeek 1.9 return sort { ($a->id1 cmp $b->id1) or ($a->psc <=> $b->psc) or ($a->id2 cmp $b->id2) } @sims;
211 : overbeek 1.5 }
212 :    
213 :     sub condense {
214 :     my($pieces) = @_;
215 :    
216 : overbeek 1.6 my $best_sc = $pieces->[0]->[10];
217 : overbeek 1.5 my @sims = sort { ($a->[6] <=> $b->[6]) } @$pieces;
218 : overbeek 1.7 while (@sims > 1)
219 : overbeek 1.5 {
220 : overbeek 1.7 my $gap1 = $sims[1]->[6] - $sims[0]->[7];
221 :     my $gap2 = $sims[1]->[8] - $sims[0]->[9];
222 :     my $diff = abs($gap1 - $gap2);
223 :     if (($gap1 <= 300) && ($gap2 <= 300) && ($diff <= (0.1 * &min($gap1,$gap2))))
224 :     {
225 :     $sims[0]->[7] = $sims[1]->[7];
226 :     $sims[0]->[9] = $sims[1]->[9];
227 :     }
228 : overbeek 1.5 splice(@sims,1,1);
229 :     }
230 : overbeek 1.6 $sims[0]->[10] = $best_sc;
231 : overbeek 1.5 return $sims[0];
232 : overbeek 1.1 }
233 :    
234 : overbeek 1.7 sub min {
235 :     my($x,$y) = @_;
236 :     return ($x <= $y) ? $x : $y;
237 :     }
238 :    
239 : overbeek 1.1 1;

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3