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

Annotation of /FigKernelScripts/package_sims.pl

Parent Directory Parent Directory | Revision Log Revision Log


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

1 : olson 1.1 # -*- perl -*-
2 : olson 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 :    
19 : olson 1.1 use Carp 'cluck';
20 :     use Data::Dumper;
21 :     use P2Pupdate;
22 :     use strict;
23 :     use Sim;
24 :    
25 :     use FIG;
26 :     my $fig = new FIG;
27 :     my $dbh = $fig->db_handle();
28 :    
29 :     my %file_map;
30 :    
31 :     # usage: package_sims SimsPackage genome
32 :    
33 :     sub usage
34 :     {
35 :     die "Usage: $0 SimsPackage genome\n";
36 :     }
37 :    
38 :     my $package = shift @ARGV;
39 :     my $genome = shift @ARGV;
40 :    
41 :     $genome =~ /^\d+\.\d+$/ or die "Invalid genome $genome\n";
42 :    
43 :     open(PKG, ">$package") or die "Could not open package file $package for writing: $!\n";
44 :    
45 :     my $warning_count = 0;
46 :     my $max_warnings = 20;
47 :    
48 :     #
49 :     # Gather pegs.
50 :     #
51 :    
52 :     my @pegs = sort { &FIG::by_fig_id($a, $b) } $fig->pegs_of($genome);
53 :    
54 :     my %pegs;
55 :    
56 :     map { $pegs{$_}++ } @pegs;
57 :    
58 :     my %pegs_left = %pegs;
59 :    
60 :     my $dbi = $dbh->{_dbh};
61 :    
62 :     my $sim_seek_stmt = $dbi->prepare(qq(SELECT filen, seek, len
63 :     FROM sim_seeks
64 :     where id = ?));
65 :    
66 :     for my $peg (@pegs)
67 :     {
68 :     #
69 :     # Determine if this id is the representative ID in the peg.synonyms.
70 :     #
71 :    
72 :     my @mapped = $fig->mapped_prot_ids($peg);
73 :    
74 :     my $rep = $mapped[0]->[0];
75 :    
76 :     if (not defined($rep))
77 :     {
78 :     print "Couldn't map $peg\n";
79 :     next;
80 :     }
81 :    
82 :     #
83 :     # Retrieve sims based on $rep
84 :     #
85 :    
86 :     $sim_seek_stmt->execute($rep);
87 :     my $sim_entries = $sim_seek_stmt->fetchall_arrayref();
88 :    
89 :     my $delta = 0;
90 :     my @entry;
91 :     if ($rep ne $peg)
92 :     {
93 :     #
94 :     # Determine the offsets we need to use in adjusting sim output.
95 :     #
96 :    
97 :    
98 :     @entry = grep { $_->[0] eq $peg } @mapped;
99 :     if (not (@entry == 1 and defined($entry[0]->[1])))
100 :     {
101 :     die "bad entry from mapped list\n";
102 :     }
103 :    
104 :     $delta = $mapped[0]->[1] - $entry[0]->[1];
105 :     }
106 :    
107 :     #
108 :     # Read the sims, and map names/lengths if necessary.
109 :     #
110 :    
111 :     if (not $sim_entries or @$sim_entries == 0)
112 :     {
113 :     print STDERR "No sims found for peg=$peg rep=$rep ", Dumper(\@mapped, \@entry);
114 :     next;
115 :     }
116 :    
117 :     for my $ent (@$sim_entries)
118 :     {
119 :     my($filen, $seek, $len) = @$ent;
120 :    
121 :     my $file = map_file($filen);
122 :     my $fh = $fig->openF($file);
123 :     $fh or die "Cannot opensims file for $file $filen\n";
124 :    
125 :     my $block = &FIG::read_block($fh, $seek, $len);
126 :    
127 :     #
128 :     # Walk thru the sims in the block.
129 :     #
130 :     # If this peg isn't the representative synonym, adjust the id1
131 :     # name to our peg and shift the match position in id1. (since
132 :     # the tails are the same, but the positions are reported with
133 :     # respect to the head).
134 :     #
135 :    
136 :     for my $line (@$block)
137 :     {
138 :     my $fixed = ($line =~ s/\t\t/\t/g);
139 :    
140 :     my $sim = Sim->new_from_line($line);
141 :    
142 :     if (!$sim->validate)
143 :     {
144 :     $warning_count++;
145 :    
146 :     cluck "Invalid similarity for $sim, skipping" if $warning_count < $max_warnings;
147 :     next;
148 :     }
149 :    
150 :     if ($fixed)
151 :     {
152 :     $warning_count++;
153 :     cluck "Fixed double-tab in similarity $sim" if $warning_count < $max_warnings;
154 :     }
155 :    
156 :     if ($rep ne $peg)
157 :     {
158 :     $sim->[0] = $peg;
159 :     $sim->[6] -= $delta;
160 :     $sim->[7] -= $delta;
161 :     }
162 :    
163 :     print PKG $sim->as_line;
164 :     }
165 :     }
166 :     }
167 :    
168 :     print "Package complete with $warning_count warnings\n";
169 :     close(PKG);
170 :    
171 :     sub map_file
172 :     {
173 :     my($filen) = @_;
174 :    
175 :     my $mapped = $file_map{$filen};
176 :     if (not defined($mapped))
177 :     {
178 :     $mapped = $fig->N2file($filen);
179 :     $file_map{$filen} = $mapped;
180 :     }
181 :     return $mapped;
182 :     }

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3