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

Annotation of /FigKernelScripts/package_sims.pl

Parent Directory Parent Directory | Revision Log Revision Log


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

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

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3