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

Diff of /FigKernelPackages/ProtSims.pm

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 1.2, Tue Aug 18 15:16:21 2009 UTC revision 1.15, Tue Sep 14 19:09:17 2010 UTC
# Line 1  Line 1 
1    
2    # This is a SAS component
3    
4  #  #
5  # Copyright (c) 2003-2006 University of Chicago and Fellowship  # Copyright (c) 2003-2006 University of Chicago and Fellowship
6  # for Interpretations of Genomes. All Rights Reserved.  # for Interpretations of Genomes. All Rights Reserved.
# Line 18  Line 20 
20  package ProtSims;  package ProtSims;
21    
22  use strict;  use strict;
23    use Sim;
24  use Data::Dumper;  use Data::Dumper;
25  use Carp;  use Carp;
 use FIG_Config;  
26  use gjoseqlib;  use gjoseqlib;
27    use DB_File;
28    
29    use SeedAware;
30    
31    my $blat_cmd = SeedAware::executable_for("blat");
32    my $blastall_cmd = SeedAware::executable_for("blastall");
33    my $formatdb_cmd = SeedAware::executable_for("formatdb");
34    
35  sub blastP {  sub blastP {
36      my($q,$db,$min_hits) = @_;      my($q,$db,$min_hits,$use_blast) = @_;
37    
38      my(@q,@db);      my(@q,@db);
39    
40        my $tmp_dir = SeedAware::location_of_tmp();
41    
42      my $qF;      my $qF;
43      if (ref $q)      if (ref $q)
44      {      {
45          $qF = "/tmp/query.$$";          $qF = "$tmp_dir/query.$$";
46          &gjoseqlib::print_alignment_as_fasta($qF,$q);          &gjoseqlib::print_alignment_as_fasta($qF,$q);
47      }      }
48      else      else
# Line 39  Line 50 
50          $qF = $q;          $qF = $q;
51          if (-e $qF)          if (-e $qF)
52          {          {
             if (! ((-e "$qF.psq") && ((-M "$qF.psq") < (-M $qF))))  
             {  
                 system "$FIG_Config::ext_bin/formatdb -i $qF";  
             }  
53              @q = &gjoseqlib::read_fasta($qF);              @q = &gjoseqlib::read_fasta($qF);
54              $q = \@q;              $q = \@q;
55          }          }
# Line 52  Line 59 
59          }          }
60      }      }
61    
62        my $db_lengths = {};
63    
64      my $dbF;      my $dbF;
65      if (ref $db)      if (ref $db)
66      {      {
67          $dbF = "/tmp/db.$$";          $dbF = "$tmp_dir/db.$$";
68          &gjoseqlib::print_alignment_as_fasta($dbF,$db);          &gjoseqlib::print_alignment_as_fasta($dbF,$db);
69          system "$FIG_Config::ext_bin/formatdb -i $dbF";          system($formatdb_cmd, '-i', $dbF);
70            $db_lengths->{$_->[0]} = length($_->[2]) for @$db;
71      }      }
72      else      else
73      {      {
# Line 66  Line 76 
76          {          {
77              if (! ((-e "$dbF.psq") && ((-M "$dbF.psq") < (-M $dbF))))              if (! ((-e "$dbF.psq") && ((-M "$dbF.psq") < (-M $dbF))))
78              {              {
79                  system "$FIG_Config::ext_bin/formatdb -i $dbF";                  system($formatdb_cmd, '-i', $dbF);;
80                }
81    
82                my $db_tie;
83                if (-f "$dbF.lengths")
84                {
85                    $db_tie = tie %$db_lengths, 'DB_File', "$dbF.lengths", O_RDONLY, 0, $DB_BTREE;
86                    $db_tie or warn "Cannot tie $dbF.lengths: $!";
87                }
88                if (!defined($db_tie))
89                {
90                    if (open(my $dbfh, "<", $dbF))
91                    {
92                        while (my($id, $def, $seq) = gjoseqlib::read_next_fasta_seq($dbfh))
93                        {
94                            $db_lengths->{$id} = length($seq);
95                        }
96                        close($dbfh);
97                    }
98                    else
99                    {
100                        warn "Cannot open $dbF: $!";
101                    }
102              }              }
             @db = &gjoseqlib::read_fasta($dbF);  
             $db = \@db;  
103          }          }
104          else          else
105          {          {
# Line 77  Line 107 
107          }          }
108      }      }
109    
110      my $tmpF = "/tmp/sim.out.$$";      my $tmpF = "$tmp_dir/sim.out.$$";
111      system "$FIG_Config::ext_bin/blat $dbF $qF -prot -out=blast8 $tmpF > /dev/null";  
112        if ($use_blast)
113        {
114            open(my $fh, ">", $tmpF);
115            close($fh);
116        }
117        else
118        {
119            my @cmd = ($blat_cmd, $dbF, $qF, "-prot", "-out=blast8", $tmpF);
120    
121            my $rc = system_with_redirect(\@cmd, { stdout => '/dev/null' } );
122            print STDERR "Blat returns $rc: @cmd\n";
123            if ($rc != 0)
124            {
125                warn "Blat run failed with rc=$rc\n";
126                open(my $fh, ">", $tmpF);
127                close($fh);
128            }
129        }
130    
131    #     my $cmd = $use_blast ? "touch $tmpF" : "$blat_cmd $dbF $qF -prot -out=blast8 $tmpF > /dev/null";
132    #     #print STDERR "$cmd\n";
133    #     my $rc = system $cmd;
134    #     if ($rc != 0)
135    #     {
136    #       die "ProtSims::blastP: blat run failed with rc=$rc: $cmd\n";
137    #     }
138    
139      my @sims1 = map { chomp; [split(/\s+/,$_)] } `cat $tmpF`;      my @sims1 = ();
140        open(BLAT,"<$tmpF") || die "could not open $tmpF";
141        my $sim = <BLAT>;
142        while ($sim && ($sim=~ /^(\S+\t\S+)/))
143        {
144            my $pegs = $1;
145            my @pieces = ();
146            while ($sim && ($sim=~ /^((\S+\t\S+).*\S)/) && ($2 eq $pegs))
147            {
148                push(@pieces,[split(/\t/,$1)]);
149                $sim = <BLAT>;
150            }
151            push(@sims1,&condense(\@pieces));
152        }
153        close(BLAT);
154      unlink $tmpF;      unlink $tmpF;
155  #    print STDERR &Dumper(\@sims1);      #print STDERR &Dumper(sims1 => \@sims1);
156    
157      my @rerun = ();      my @rerun = ();
158      my @sims  = ();      my @sims  = ();
# Line 105  Line 175 
175  #       print STDERR "hits=$hits min_hits=$min_hits\n";  #       print STDERR "hits=$hits min_hits=$min_hits\n";
176          if ($hits >= $min_hits)          if ($hits >= $min_hits)
177          {          {
178              push(@sims,@sims1[$simI..$simI1]);              push(@sims,@sims1[$simI..$simI1-1]);
179          }          }
180          else          else
181          {          {
# Line 114  Line 184 
184          $simI = $simI1;          $simI = $simI1;
185          $qI++;          $qI++;
186      }      }
187  #    print STDERR &Dumper(\@rerun);      #print STDERR &Dumper(rerun => \@rerun);
188    
189      if (@rerun > 0)      if (@rerun > 0)
190      {      {
191          my $tmpQ = "/tmp/tmpQ.$$";          my $tmpQ = "$tmp_dir/tmpQ.$$";
192          &gjoseqlib::print_alignment_as_fasta($tmpQ,\@rerun);          &gjoseqlib::print_alignment_as_fasta($tmpQ,\@rerun);
193          my @blastout = map { chomp; [split(/\s+/,$_)] }          # my $cmd = "$blastall_cmd -m 8 -i $tmpQ -d $dbF -FF -p blastp -e 1e-5";
194                         `$FIG_Config::ext_bin/blastall -m 8 -i $tmpQ -d $dbF -FF -p blastp`;          my @cmd = ($blastall_cmd,  '-m', 8, '-i', $tmpQ, '-d', $dbF, '-FF', '-p', 'blastp', '-e', '1e-5');
195          push(@sims,@blastout);          #print STDERR "$cmd\n";
196            #open(BL, "$cmd|") or die "ProtSims::blastP: pipe to blast failed with $!: $cmd\n";
197    
198            #
199            # It'd be nice to do this but windows doesn't support it.
200            #
201            #open(BL, "-|", @cmd) or die "ProtSims::blastP: pipe to blast failed with $!: @cmd\n";
202            my $out_tmp = "$tmp_dir/blast_out.$$";
203            push(@cmd, "-o", $out_tmp);
204            print STDERR "@cmd\n";
205            my $rc = system(@cmd);
206            if ($rc != 0)
207            {
208                warn "Error rc=$rc running blast: @cmd\n";
209            }
210            else
211            {
212                if (open(BL, "<", $out_tmp))
213                {
214                    while (<BL>)
215                    {
216                        chomp;
217                        push @sims, [ split(/\s+/, $_) ];
218                    }
219                    #my @blastout = map { chomp; [split(/\s+/,$_)] } <BL>;
220                    # push(@sims,@blastout);
221                    close(BL);
222                }
223            }
224            unlink $out_tmp;
225          unlink $tmpQ;          unlink $tmpQ;
226      }      }
227    
228      my %lnQ   = map { $_->[0] => length($_->[2]) } @q;      my %lnQ   = map { $_->[0] => length($_->[2]) } @$q;
229      my %lnDB  = map { $_->[0] => length($_->[2]) } @db;  
230      @sims = map { push(@$_,$lnQ{$_->[0]},$lnDB{$_->[1]}); bless($_,'Sim') } @sims;      @sims = map { push(@$_, $lnQ{$_->[0]}, $db_lengths->{$_->[1]}); bless($_,'Sim') } @sims;
231    
232      if ($qF      eq "/tmp/query.$$")   { unlink $qF  }      if ($qF      eq "$tmp_dir/query.$$")   { unlink $qF  }
233      if ($dbF     eq "/tmp/db.$$")      { unlink $dbF }      if ($dbF     eq "$tmp_dir/db.$$")      { unlink($dbF,"$dbF.psq","$dbF.pin","$dbF.phr") }
234      return @sims;      return sort { ($a->id1 cmp $b->id1) or ($a->psc <=> $b->psc) or ($a->id2 cmp $b->id2) } @sims;
235    }
236    
237    sub condense {
238        my($pieces) = @_;
239    
240        my $best_sc = $pieces->[0]->[10];
241        my @sims = sort { ($a->[6] <=> $b->[6]) } @$pieces;
242        while (@sims > 1)
243        {
244            my $gap1 = $sims[1]->[6] - $sims[0]->[7];
245            my $gap2 = $sims[1]->[8] - $sims[0]->[9];
246            my $diff = abs($gap1 - $gap2);
247            if (($gap1 <= 300) && ($gap2 <= 300) && ($diff <= (0.1 * &min($gap1,$gap2))))
248            {
249                $sims[0]->[7] = $sims[1]->[7];
250                $sims[0]->[9] = $sims[1]->[9];
251            }
252            splice(@sims,1,1);
253        }
254        $sims[0]->[10] = $best_sc;
255        return $sims[0];
256    }
257    
258    sub min {
259        my($x,$y) = @_;
260        return ($x <= $y) ? $x : $y;
261  }  }
262    
263  1;  1;

Legend:
Removed from v.1.2  
changed lines
  Added in v.1.15

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3