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

Annotation of /FigKernelPackages/ProtSims.pm

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1.11 - (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.10 my $blat_cmd = "blat";
29 :     my $blastall_cmd = "blastall";
30 :     my $formatdb_cmd = "formatdb";
31 :    
32 :     BEGIN {
33 :     eval {
34 :     require FIG_Config;
35 : parrello 1.11 if (defined $FIG_Config::ext_bin) {
36 :     $blat_cmd = "$FIG_Config::ext_bin/blat";
37 :     $blastall_cmd = "$FIG_Config::ext_bin/blastall";
38 :     $formatdb_cmd = "$FIG_Config::ext_bin/formatdb";
39 :     }
40 : olson 1.10 };
41 :     }
42 :    
43 : overbeek 1.1 sub blastP {
44 :     my($q,$db,$min_hits) = @_;
45 :    
46 : overbeek 1.2 my(@q,@db);
47 :    
48 : overbeek 1.1 my $qF;
49 :     if (ref $q)
50 :     {
51 :     $qF = "/tmp/query.$$";
52 :     &gjoseqlib::print_alignment_as_fasta($qF,$q);
53 :     }
54 :     else
55 :     {
56 :     $qF = $q;
57 :     if (-e $qF)
58 :     {
59 : overbeek 1.2 @q = &gjoseqlib::read_fasta($qF);
60 : overbeek 1.1 $q = \@q;
61 :     }
62 :     else
63 :     {
64 :     die "$qF is missing";
65 :     }
66 :     }
67 :    
68 :     my $dbF;
69 :     if (ref $db)
70 :     {
71 :     $dbF = "/tmp/db.$$";
72 :     &gjoseqlib::print_alignment_as_fasta($dbF,$db);
73 : olson 1.10 system "$formatdb_cmd -i $dbF";
74 : overbeek 1.1 }
75 :     else
76 :     {
77 :     $dbF = $db;
78 :     if (-e $dbF)
79 :     {
80 :     if (! ((-e "$dbF.psq") && ((-M "$dbF.psq") < (-M $dbF))))
81 :     {
82 : olson 1.10 system "formatdb_cmd -i $dbF";
83 : overbeek 1.1 }
84 : overbeek 1.2 @db = &gjoseqlib::read_fasta($dbF);
85 : overbeek 1.1 $db = \@db;
86 :     }
87 :     else
88 :     {
89 :     die "$dbF is missing";
90 :     }
91 :     }
92 :    
93 :     my $tmpF = "/tmp/sim.out.$$";
94 : olson 1.10 my $cmd = "$blat_cmd $dbF $qF -prot -out=blast8 $tmpF > /dev/null";
95 : olson 1.4 #print STDERR "$cmd\n";
96 :     my $rc = system $cmd;
97 :     if ($rc != 0)
98 :     {
99 :     die "ProtSims::blastP: blat run failed with rc=$rc: $cmd\n";
100 :     }
101 : overbeek 1.1
102 : overbeek 1.5 my @sims1 = ();
103 :     open(BLAT,"<$tmpF") || die "could not open $tmpF";
104 :     my $sim = <BLAT>;
105 :     while ($sim && ($sim=~ /^(\S+\t\S+)/))
106 :     {
107 :     my $pegs = $1;
108 :     my @pieces = ();
109 :     while ($sim && ($sim=~ /^((\S+\t\S+).*\S)/) && ($2 eq $pegs))
110 :     {
111 :     push(@pieces,[split(/\t/,$1)]);
112 :     $sim = <BLAT>;
113 :     }
114 :     push(@sims1,&condense(\@pieces));
115 :     }
116 :     close(BLAT);
117 : overbeek 1.1 unlink $tmpF;
118 : olson 1.4 #print STDERR &Dumper(sims1 => \@sims1);
119 : overbeek 1.1
120 :     my @rerun = ();
121 :     my @sims = ();
122 :     my $qI = 0;
123 :     my $simI = 0;
124 :     while ($qI < @$q)
125 :     {
126 :     my $peg = $q->[$qI]->[0];
127 : olson 1.4 #print STDERR "processing $peg\n";
128 : overbeek 1.1 my $hits = 0;
129 :     my $simI1 = $simI;
130 :     while (($simI1 < @sims1) && ($sims1[$simI1]->[0] eq $peg))
131 :     {
132 :     if (($simI1 == $#sims1) || ($sims1[$simI1]->[1] ne $sims1[$simI1+1]->[1]))
133 :     {
134 :     $hits++;
135 :     }
136 :     $simI1++;
137 :     }
138 : olson 1.4 #print STDERR "hits=$hits min_hits=$min_hits\n";
139 : overbeek 1.1 if ($hits >= $min_hits)
140 :     {
141 : olson 1.4 push(@sims,@sims1[$simI..$simI1-1]);
142 : overbeek 1.1 }
143 :     else
144 :     {
145 :     push(@rerun,$q->[$qI]);
146 :     }
147 :     $simI = $simI1;
148 :     $qI++;
149 :     }
150 : olson 1.4 #print STDERR &Dumper(rerun => \@rerun);
151 : overbeek 1.1
152 :     if (@rerun > 0)
153 :     {
154 :     my $tmpQ = "/tmp/tmpQ.$$";
155 :     &gjoseqlib::print_alignment_as_fasta($tmpQ,\@rerun);
156 : olson 1.10 my $cmd = "$blastall_cmd -m 8 -i $tmpQ -d $dbF -FF -p blastp -e 1e-5";
157 : olson 1.4 #print STDERR "$cmd\n";
158 :     open(BL, "$cmd|") or die "ProtSims::blastP: pipe to blast failed with $!: $cmd\n";
159 :     my @blastout = map { chomp; [split(/\s+/,$_)] } <BL>;
160 :     close(BL);
161 : overbeek 1.1 push(@sims,@blastout);
162 :     unlink $tmpQ;
163 :     }
164 :    
165 : overbeek 1.3 my %lnQ = map { $_->[0] => length($_->[2]) } @$q;
166 :     my %lnDB = map { $_->[0] => length($_->[2]) } @$db;
167 : overbeek 1.1 @sims = map { push(@$_,$lnQ{$_->[0]},$lnDB{$_->[1]}); bless($_,'Sim') } @sims;
168 :    
169 :     if ($qF eq "/tmp/query.$$") { unlink $qF }
170 : overbeek 1.8 if ($dbF eq "/tmp/db.$$") { unlink($dbF,"$dbF.psq","$dbF.pin","$dbF.phr") }
171 : overbeek 1.9 return sort { ($a->id1 cmp $b->id1) or ($a->psc <=> $b->psc) or ($a->id2 cmp $b->id2) } @sims;
172 : overbeek 1.5 }
173 :    
174 :     sub condense {
175 :     my($pieces) = @_;
176 :    
177 : overbeek 1.6 my $best_sc = $pieces->[0]->[10];
178 : overbeek 1.5 my @sims = sort { ($a->[6] <=> $b->[6]) } @$pieces;
179 : overbeek 1.7 while (@sims > 1)
180 : overbeek 1.5 {
181 : overbeek 1.7 my $gap1 = $sims[1]->[6] - $sims[0]->[7];
182 :     my $gap2 = $sims[1]->[8] - $sims[0]->[9];
183 :     my $diff = abs($gap1 - $gap2);
184 :     if (($gap1 <= 300) && ($gap2 <= 300) && ($diff <= (0.1 * &min($gap1,$gap2))))
185 :     {
186 :     $sims[0]->[7] = $sims[1]->[7];
187 :     $sims[0]->[9] = $sims[1]->[9];
188 :     }
189 : overbeek 1.5 splice(@sims,1,1);
190 :     }
191 : overbeek 1.6 $sims[0]->[10] = $best_sc;
192 : overbeek 1.5 return $sims[0];
193 : overbeek 1.1 }
194 :    
195 : overbeek 1.7 sub min {
196 :     my($x,$y) = @_;
197 :     return ($x <= $y) ? $x : $y;
198 :     }
199 :    
200 : overbeek 1.1 1;

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3