[Bio] / FigKernelScripts / build_solid_set.pl Repository:
ViewVC logotype

Annotation of /FigKernelScripts/build_solid_set.pl

Parent Directory Parent Directory | Revision Log Revision Log


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

1 : overbeek 1.1 use strict;
2 :     use Data::Dumper;
3 :     use Carp;
4 :     use SeedUtils;
5 :     use GenerateClusters;
6 : overbeek 1.3 use FIG;
7 :     my $fig = new FIG;
8 : overbeek 1.1
9 :     =head1 svr_build_solid_set -g CloseGenomes < clustersRel
10 :    
11 :     Produces a set of "solid rectangles" separated by "//" lines. Each group of lines
12 :     is the set of rows in one rectangle. Each line in a rectangle begins with a genome ID, and is followed
13 :     by an ordered set of tab-separated Pegs (one
14 :     line per genome. The first line will be an input cluster from the input set
15 :     (from the reference genome)
16 :    
17 :     ------
18 :    
19 :     Example:
20 :    
21 :     svr_build_solid_set -g close.genomes < RefClusters > solid.rectangles
22 :    
23 :     =head2 Command-Line Options
24 :    
25 :     =over 4
26 :    
27 :     =item -g GenomesToTest
28 :    
29 :     =back
30 :    
31 :     =head2 Output Format
32 :    
33 :     The output is a set of lines, one per genome projected to.
34 :     Each line is a set of tab-separated Pegs from a single genome
35 :    
36 :     =cut
37 :    
38 :     use SeedUtils;
39 :     use SAPserver;
40 :     my $sapO = SAPserver->new();
41 :     use Getopt::Long;
42 :     use GenerateClusters;
43 :    
44 :     my $usage = "usage: svr_build_solid_set -g CloseGenomesFile < clusters_from_reference_genome";
45 :    
46 :     my $file;
47 :    
48 :     my $rc = GetOptions('g=s' => \$file);
49 :     if (! $rc) { print STDERR $usage; exit }
50 :    
51 :     if (! $file) { die "you need to give a file of genomes -f argument" }
52 :    
53 :     my @genomes = map { ($_ =~ /^(\S+)/) ? $1 : () } `cat $file`;
54 :     print STDERR scalar @genomes, " genomes to check\n";
55 :    
56 :     my @sets;
57 :     my $last = <STDIN>;
58 :     while ($last && ($last =~ /^(\S+)\t(\S+)/))
59 :     {
60 :     my $curr = $1;
61 :     my $set = [];
62 :     while ($last && ($last =~ /^(\S+)\t(\S+)/) && ($1 eq $curr))
63 :     {
64 :     push(@$set,$2);
65 :     $last = <STDIN>;
66 :     }
67 : overbeek 1.3 if ($set = &ok_set($set,$fig))
68 :     {
69 :     push(@sets,$set);
70 :     }
71 : overbeek 1.1 }
72 :     print STDERR scalar @sets, " sets were read\n";
73 :     my $ref = &SeedUtils::genome_of($sets[0]->[0]);
74 :     print STDERR "reference genome = $ref\n";
75 :     my $corrH = &GenerateClusters::bbhs($ref,\@genomes,$sapO);
76 :     print STDERR "read all correspondence tables\n";
77 :    
78 :     foreach my $set (@sets)
79 :     {
80 :     print join("\t",($ref,@$set)),"\n";
81 :     foreach my $g (@genomes)
82 :     {
83 :     my $instance = &project($set,$ref,$g,$corrH,$sapO);
84 :     if ($instance)
85 :     {
86 :     print join("\t",($g,@$instance)),"\n";
87 :     }
88 :     }
89 :     print "//\n";
90 :     }
91 :    
92 : overbeek 1.3 sub ok_set {
93 :     my($set,$fig) = @_;
94 :    
95 :     my @ok_pegs = grep { my $f = $fig->function_of($_); $f !~ /transport|regulat/ } @$set;
96 :     if (@ok_pegs >= (0.5 * @$set)) { return $set }
97 :     elsif (@ok_pegs > 2) { return \@ok_pegs }
98 :     else { return undef }
99 :     }
100 :    
101 :    
102 : overbeek 1.1 sub project {
103 :     my($set,$ref,$g,$corrH,$sapO) = @_;
104 :    
105 :    
106 :     my @set2 = map { ($corrH->{$g}->{$_}) ? $corrH->{$g}->{$_} : () } @$set;
107 :     if (@set2 == @$set)
108 :     {
109 :     my $sz = &sz_of_cluster($set,$sapO);
110 :     my $sz2 = &sz_of_cluster(\@set2,$sapO);
111 : overbeek 1.2 if ($sz && ((($sz2-$sz) / $sz) < 0.2))
112 : overbeek 1.1 {
113 :     return \@set2;
114 :     }
115 :     }
116 :     return undef;
117 :     }
118 :    
119 :     sub sz_of_cluster {
120 :     my($set,$sapO) = @_;
121 :    
122 :     my $locH = $sapO->fid_locations( -ids => $set, -boundaries => 1);
123 :     my $min;
124 :     my $max;
125 :     foreach my $peg (keys(%$locH))
126 :     {
127 :     my $loc = $locH->{$peg};
128 :     if ($loc && ($loc =~ /^.*:(\S+)_(\d+)[+-](\d+)$/))
129 :     {
130 :     my($contig,$beg,$strand,$len) = ($1,$2,$3,$4);
131 :     if ($strand eq '+')
132 :     {
133 :     if ((! $min) || ($min > $beg)) { $min = $beg }
134 :     if ((! $max) || ($max < ($beg + $len - 1))) { $max = $beg + $len - 1 }
135 :     }
136 :     else
137 :     {
138 :     if ((! $min) || ($min > ($beg - ($len - 1)))) { $min = $beg-($len - 1) }
139 :     if ((! $max) || ($max < $beg)) { $max = $beg }
140 :     }
141 :     }
142 :     }
143 :     return ($max - $min) + 1;
144 :     }
145 :    

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3