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

Annotation of /FigKernelPackages/ProtSims.pm

Parent Directory Parent Directory | Revision Log Revision Log


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

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

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3