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

Annotation of /FigKernelPackages/PartitionSeqs.pm

Parent Directory Parent Directory | Revision Log Revision Log


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

1 : overbeek 1.1 package PartitionSeqs;
2 :    
3 :     use strict;
4 :     use gjoseqlib;
5 :     use FIG;
6 :     use Data::Dumper;
7 :     use FIG_Config;
8 :    
9 :    
10 :     # This module can be used to take a set of similar pegs and partition them,
11 :     # putting those that are quite similar and similar over most of the length
12 :     # of each sequence into the same partition.
13 :     #
14 :     # \@sets = PartitionSeqs:partition( \%options )
15 :     #
16 :     # Options:
17 :     #
18 :     # seqs => [[id1,comment1,seq1],[id2,comment2,seq2]...]
19 :     # pegs => [peg1,peg2,...]
20 :     # use => blast|sims # D = sims
21 :     # expect => maximum e-value for blast
22 :     # coverage => $min_coverage # D = 0.70
23 :     # identity => $min_identity # D = 0.60
24 :     #
25 :     # OUTPUT: ([id1,id2,...],[id3,id4,...]...)
26 :     #
27 :    
28 :     sub partition
29 :     {
30 :     my ($options ) = @_;
31 :     my $min_cover = $options->{ coverage } ||= 0.70; # Minimum fraction of reference covered
32 :     my $max_exp = $options->{ expect } ||= 1.0e-10; # Maximum e-value for blast
33 :     my $min_identity = $options->{ expect } ||= 0.3; # Maximum e-value for blast
34 :     my $use = $options->{ use } ||= 'sims'; # use sims by default
35 :    
36 :     my($pegs,$connections);
37 :     $connections = {};
38 :     my $fig = new FIG;
39 :    
40 :     if (($use eq 'sims') && ($pegs = $options->{pegs}))
41 :     {
42 :     $connections = &connect_by_sims($fig,$pegs,$max_exp,$min_identity,$min_cover);
43 :     }
44 :     else
45 :     {
46 :     my @seqs;
47 :     if ($pegs = $options->{pegs})
48 :     {
49 : overbeek 1.2 @seqs = grep { length($_->[2]) > 10 } map { [$_,'',$fig->get_translation($_)] } @$pegs;
50 : overbeek 1.1 }
51 :     else
52 :     {
53 :     @seqs = @{$options->{seqs}};
54 :     }
55 : overbeek 1.2 (@seqs > 0) || return ();
56 :    
57 : overbeek 1.1 my %length = map { $_->[0] => length($_->[2]) } @seqs;
58 :    
59 :     my $blastF = "$FIG_Config::temp/tmp.$$.fasta";
60 :     &gjoseqlib::print_alignment_as_fasta($blastF,\@seqs);
61 :     system "$FIG_Config::ext_bin/formatdb -p T -i $blastF";
62 :     open(BLAST,"$FIG_Config::ext_bin/blastall -i $blastF -d $blastF -m 8 -e $max_exp -p blastp |")
63 :     || die "blast failed";
64 :     while (defined($_ = <BLAST>))
65 :     {
66 :     my($id1,$id2,$iden,undef,undef,undef,$b1,$e1,$b2,$e2) = split(/\t/,$_);
67 :     if (($id1 ne $id2) && ($min_identity <= $iden))
68 :     {
69 :     my $ln1 = $length{$id1};
70 :     my $ln2 = $length{$id2};
71 :     if (((($e1+1-$b1)/$ln1) >= $min_cover) && ((($e2+1-$b2)/$ln2) >= $min_cover))
72 :     {
73 :     push(@{$connections->{$id1}},$id2);
74 :     push(@{$connections->{$id2}},$id1);
75 :     }
76 :     }
77 :     }
78 :     close(BLAST);
79 :     unlink($blastF);
80 :     }
81 : overbeek 1.3 my @ans = &cluster($connections);
82 :     return wantarray ? @ans : \@ans;
83 : overbeek 1.1 }
84 :    
85 :     sub connect_by_sims {
86 :     my($fig,$pegs,$max_exp,$min_identity,$min_cover) = @_;
87 :    
88 :     my $conn = {};
89 :     my %pegsH = map { $_ => 1 } @$pegs;
90 :     foreach my $peg (@$pegs)
91 :     {
92 :     foreach my $sim ($fig->sims($peg,10000,$max_exp,'fig'))
93 :     {
94 :     if (($pegsH{$sim->id2}) &&
95 :     ($min_identity <= $sim->iden) &&
96 :     ((($sim->e1 + 1 - $sim->b1) / $sim->ln1) >= $min_cover) &&
97 :     ((($sim->e2 + 1 - $sim->b2) / $sim->ln2) >= $min_cover))
98 :     {
99 :     push(@{$conn->{$peg}},$sim->id2);
100 :     push(@{$conn->{$sim->id2}},$peg);
101 :     }
102 :     }
103 :     }
104 :     return $conn;
105 :     }
106 :    
107 :     sub cluster {
108 :     my($conn) = @_;
109 :    
110 :     my %seen;
111 :     my @clusters = ();
112 :     my @pegs = keys(%$conn);
113 :     my($peg);
114 :     while ($peg = shift @pegs)
115 :     {
116 :     if (! $seen{$peg})
117 :     {
118 :     my $cluster = [$peg];
119 :     $seen{$peg} = 1;
120 :     my $i;
121 :     for ($i=0; ($i < @$cluster); $i++)
122 :     {
123 :     my $x = $conn->{$cluster->[$i]};
124 :     foreach my $y (@$x)
125 :     {
126 :     if (! $seen{$y})
127 :     {
128 :     push(@$cluster,$y);
129 :     $seen{$y} = 1;
130 :     }
131 :     }
132 :     }
133 :     push(@clusters,$cluster);
134 :     }
135 :     }
136 :     return @clusters;
137 :     }
138 :    
139 :     1;

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3