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

Annotation of /FigKernelPackages/ProtSims.pm

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1.2 - (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 :     use Data::Dumper;
22 :     use Carp;
23 :     use FIG_Config;
24 :     use gjoseqlib;
25 :    
26 :     sub blastP {
27 :     my($q,$db,$min_hits) = @_;
28 :    
29 : overbeek 1.2 my(@q,@db);
30 :    
31 : overbeek 1.1 my $qF;
32 :     if (ref $q)
33 :     {
34 :     $qF = "/tmp/query.$$";
35 :     &gjoseqlib::print_alignment_as_fasta($qF,$q);
36 :     }
37 :     else
38 :     {
39 :     $qF = $q;
40 :     if (-e $qF)
41 :     {
42 :     if (! ((-e "$qF.psq") && ((-M "$qF.psq") < (-M $qF))))
43 :     {
44 :     system "$FIG_Config::ext_bin/formatdb -i $qF";
45 :     }
46 : overbeek 1.2 @q = &gjoseqlib::read_fasta($qF);
47 : overbeek 1.1 $q = \@q;
48 :     }
49 :     else
50 :     {
51 :     die "$qF is missing";
52 :     }
53 :     }
54 :    
55 :     my $dbF;
56 :     if (ref $db)
57 :     {
58 :     $dbF = "/tmp/db.$$";
59 :     &gjoseqlib::print_alignment_as_fasta($dbF,$db);
60 :     system "$FIG_Config::ext_bin/formatdb -i $dbF";
61 :     }
62 :     else
63 :     {
64 :     $dbF = $db;
65 :     if (-e $dbF)
66 :     {
67 :     if (! ((-e "$dbF.psq") && ((-M "$dbF.psq") < (-M $dbF))))
68 :     {
69 :     system "$FIG_Config::ext_bin/formatdb -i $dbF";
70 :     }
71 : overbeek 1.2 @db = &gjoseqlib::read_fasta($dbF);
72 : overbeek 1.1 $db = \@db;
73 :     }
74 :     else
75 :     {
76 :     die "$dbF is missing";
77 :     }
78 :     }
79 :    
80 :     my $tmpF = "/tmp/sim.out.$$";
81 :     system "$FIG_Config::ext_bin/blat $dbF $qF -prot -out=blast8 $tmpF > /dev/null";
82 :    
83 :     my @sims1 = map { chomp; [split(/\s+/,$_)] } `cat $tmpF`;
84 :     unlink $tmpF;
85 :     # print STDERR &Dumper(\@sims1);
86 :    
87 :     my @rerun = ();
88 :     my @sims = ();
89 :     my $qI = 0;
90 :     my $simI = 0;
91 :     while ($qI < @$q)
92 :     {
93 :     my $peg = $q->[$qI]->[0];
94 :     # print STDERR "processing $peg\n";
95 :     my $hits = 0;
96 :     my $simI1 = $simI;
97 :     while (($simI1 < @sims1) && ($sims1[$simI1]->[0] eq $peg))
98 :     {
99 :     if (($simI1 == $#sims1) || ($sims1[$simI1]->[1] ne $sims1[$simI1+1]->[1]))
100 :     {
101 :     $hits++;
102 :     }
103 :     $simI1++;
104 :     }
105 :     # print STDERR "hits=$hits min_hits=$min_hits\n";
106 :     if ($hits >= $min_hits)
107 :     {
108 :     push(@sims,@sims1[$simI..$simI1]);
109 :     }
110 :     else
111 :     {
112 :     push(@rerun,$q->[$qI]);
113 :     }
114 :     $simI = $simI1;
115 :     $qI++;
116 :     }
117 :     # print STDERR &Dumper(\@rerun);
118 :    
119 :     if (@rerun > 0)
120 :     {
121 :     my $tmpQ = "/tmp/tmpQ.$$";
122 :     &gjoseqlib::print_alignment_as_fasta($tmpQ,\@rerun);
123 :     my @blastout = map { chomp; [split(/\s+/,$_)] }
124 :     `$FIG_Config::ext_bin/blastall -m 8 -i $tmpQ -d $dbF -FF -p blastp`;
125 :     push(@sims,@blastout);
126 :     unlink $tmpQ;
127 :     }
128 :    
129 :     my %lnQ = map { $_->[0] => length($_->[2]) } @q;
130 :     my %lnDB = map { $_->[0] => length($_->[2]) } @db;
131 :     @sims = map { push(@$_,$lnQ{$_->[0]},$lnDB{$_->[1]}); bless($_,'Sim') } @sims;
132 :    
133 :     if ($qF eq "/tmp/query.$$") { unlink $qF }
134 :     if ($dbF eq "/tmp/db.$$") { unlink $dbF }
135 :     return @sims;
136 :     }
137 :    
138 :     1;

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3