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

Annotation of /FigKernelPackages/PartitionSeqs.pm

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1.1 - (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 :     @seqs = map { [$_,'',$fig->get_translation($_)] } @$pegs;
50 :     }
51 :     else
52 :     {
53 :     @seqs = @{$options->{seqs}};
54 :     }
55 :     my %length = map { $_->[0] => length($_->[2]) } @seqs;
56 :    
57 :     my $blastF = "$FIG_Config::temp/tmp.$$.fasta";
58 :     &gjoseqlib::print_alignment_as_fasta($blastF,\@seqs);
59 :     system "$FIG_Config::ext_bin/formatdb -p T -i $blastF";
60 :     open(BLAST,"$FIG_Config::ext_bin/blastall -i $blastF -d $blastF -m 8 -e $max_exp -p blastp |")
61 :     || die "blast failed";
62 :     while (defined($_ = <BLAST>))
63 :     {
64 :     my($id1,$id2,$iden,undef,undef,undef,$b1,$e1,$b2,$e2) = split(/\t/,$_);
65 :     if (($id1 ne $id2) && ($min_identity <= $iden))
66 :     {
67 :     my $ln1 = $length{$id1};
68 :     my $ln2 = $length{$id2};
69 :     if (((($e1+1-$b1)/$ln1) >= $min_cover) && ((($e2+1-$b2)/$ln2) >= $min_cover))
70 :     {
71 :     push(@{$connections->{$id1}},$id2);
72 :     push(@{$connections->{$id2}},$id1);
73 :     }
74 :     }
75 :     }
76 :     close(BLAST);
77 :     unlink($blastF);
78 :     }
79 :     return &cluster($connections);
80 :     }
81 :    
82 :     sub connect_by_sims {
83 :     my($fig,$pegs,$max_exp,$min_identity,$min_cover) = @_;
84 :    
85 :     my $conn = {};
86 :     my %pegsH = map { $_ => 1 } @$pegs;
87 :     foreach my $peg (@$pegs)
88 :     {
89 :     foreach my $sim ($fig->sims($peg,10000,$max_exp,'fig'))
90 :     {
91 :     if (($pegsH{$sim->id2}) &&
92 :     ($min_identity <= $sim->iden) &&
93 :     ((($sim->e1 + 1 - $sim->b1) / $sim->ln1) >= $min_cover) &&
94 :     ((($sim->e2 + 1 - $sim->b2) / $sim->ln2) >= $min_cover))
95 :     {
96 :     push(@{$conn->{$peg}},$sim->id2);
97 :     push(@{$conn->{$sim->id2}},$peg);
98 :     }
99 :     }
100 :     }
101 :     return $conn;
102 :     }
103 :    
104 :     sub cluster {
105 :     my($conn) = @_;
106 :    
107 :     my %seen;
108 :     my @clusters = ();
109 :     my @pegs = keys(%$conn);
110 :     my($peg);
111 :     while ($peg = shift @pegs)
112 :     {
113 :     if (! $seen{$peg})
114 :     {
115 :     my $cluster = [$peg];
116 :     $seen{$peg} = 1;
117 :     my $i;
118 :     for ($i=0; ($i < @$cluster); $i++)
119 :     {
120 :     my $x = $conn->{$cluster->[$i]};
121 :     foreach my $y (@$x)
122 :     {
123 :     if (! $seen{$y})
124 :     {
125 :     push(@$cluster,$y);
126 :     $seen{$y} = 1;
127 :     }
128 :     }
129 :     }
130 :     push(@clusters,$cluster);
131 :     }
132 :     }
133 :     return @clusters;
134 :     }
135 :    
136 :     1;

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3