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

Annotation of /FigKernelPackages/PropagateGBMetadata.pm

Parent Directory Parent Directory | Revision Log Revision Log


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

1 : olson 1.1 package PropagateGBMetadata;
2 :    
3 : olson 1.2 #
4 :     # This is a SAS component.
5 :     #
6 :    
7 : olson 1.1 use strict;
8 :     use overlap_resolution;
9 :     use Data::Dumper;
10 :     use Time::HiRes 'gettimeofday';
11 :     use Set::IntervalTree;
12 :     use List::Util qw(min max);
13 :    
14 :     #
15 :     # Given a GTO that has multiple annotations including a GenBank import,
16 :     # propagate the feature metadata from the GB features to the features
17 :     # that have been otherwise called.
18 :     #
19 :    
20 :     sub propagate_gb_metadata
21 :     {
22 :     my($gto, $params, $user) = @_;
23 :    
24 :     $params = {} unless ref($params) eq 'HASH';
25 :     $params->{min_rna_pct_coverage} ||= 90;
26 :    
27 :     $user ||= "PropagateGBMetadata";
28 :    
29 :     my %orf;
30 :    
31 :     my $hostname = `hostname`;
32 :     chomp $hostname;
33 :    
34 :     my $event = { tool_name => "PropagateGBMetadata",
35 :     execution_time => scalar gettimeofday,
36 :     parameters => [],
37 :     hostname => $hostname };
38 :    
39 :     my $event_id = $gto->add_analysis_event($event);
40 :    
41 :     #
42 :     # First propagate to exact matches.
43 :     #
44 :    
45 :     my %seen;
46 :    
47 :     for my $f ($gto->features)
48 :     {
49 :     my($contig, $left, $right, $dir, $size) = overlap_resolution::bounds($f);
50 :    
51 :     my $stop = $dir eq '+' ? $right : $left;
52 :     # print join("\t", $f->{id}, $contig, $stop, $left, $right, $dir), "\n";
53 :     my $orf = join(".", $contig, $stop, $dir, $size);
54 :     push(@{$orf{$orf}}, $f);
55 :     }
56 :    
57 :     propagate_orfs($gto, $params, $user, $event_id, \%orf, \%seen);
58 :    
59 :     #
60 :     # Then propagate fuzzier matches.
61 :     #
62 :    
63 :     %orf = ();
64 :    
65 :     for my $f ($gto->features)
66 :     {
67 :     next if $seen{$f->{id}};
68 :    
69 :     my($contig, $left, $right, $dir, $size) = overlap_resolution::bounds($f);
70 :    
71 :     my $stop = $dir eq '+' ? $right : $left;
72 :     # print join("\t", $f->{id}, $contig, $stop, $left, $right, $dir), "\n";
73 :     my $orf = join(".", $contig, $stop, $dir);
74 :     push(@{$orf{$orf}}, $f);
75 :     }
76 :    
77 :     propagate_orfs($gto, $params, $user, $event_id, \%orf, \%seen);
78 :    
79 :     propagate_trnas($gto, $params, $user, $event_id, \%seen);
80 :    
81 :     return $gto;
82 :     }
83 :    
84 :     #
85 :     # Propagate tRNA metadata based on coverage - we see differing endpoints
86 :     # so the prior orf-based mechanism does not work.
87 :     #
88 :    
89 :     sub propagate_trnas
90 :     {
91 :     my($gto, $params, $user, $event_id, $seen) = @_;
92 :    
93 :     #
94 :     # We set up interval trees to store the GB-annotated
95 :     # features.
96 :     #
97 :    
98 :     my %trees;
99 :    
100 :     my %trees;
101 :     for my $ctg ($gto->contigs())
102 :     {
103 :     my $tree = Set::IntervalTree->new();
104 :     $trees{$ctg->{id}} = $tree;
105 :     }
106 :    
107 :     for my $f ($gto->features)
108 :     {
109 :     my($ctg, $min, $max, $dir, $len) = GenomeTypeObject::bounds($f);
110 :    
111 :     next unless $f->{type} =~ /rna/i;
112 :     next unless exists $f->{genbank_feature};
113 :    
114 : olson 1.3 #
115 :     # From Set::IntervalTree doc:
116 :     # All intervals are half-open, i.e. [1,3), [2,6), etc.
117 :     # Thus we bump the right endpoint.
118 :     #
119 :     $max++;
120 :    
121 : olson 1.1 my $tree = $trees{$ctg};
122 :     $tree->insert($f->{id}, $min, $max);
123 :     }
124 :    
125 :     #
126 :     # Now walk the rna features looking for overlaps with
127 :     for my $f ($gto->features)
128 :     {
129 :     my $fid = $f->{id};
130 :    
131 :     next unless $f->{type} =~ /rna/i;
132 :     next if $seen->{$fid};
133 :     next if exists $f->{genbank_feature};
134 :    
135 :     my($ctg, $min, $max, $dir, $len) = GenomeTypeObject::bounds($f);
136 :     my $tree = $trees{$ctg};
137 : olson 1.3 my $overlap = $tree->fetch($min, $max+1);
138 : olson 1.1 my @overlap = grep { $_ ne $fid } @$overlap;
139 :     #
140 :     # @overlap contains the genbank feature we are mapping from.
141 :     #
142 :     if (@overlap)
143 :     {
144 :     my @poss;
145 :     for my $fid2 (@overlap)
146 :     {
147 :     my $f2 = $gto->find_feature($fid2);
148 :     my($ctg2, $min2, $max2, $dir2, $len2) = GenomeTypeObject::bounds($f2);
149 :    
150 :     my $olen;
151 :     my $t = "not";
152 :     if ($min > $min2 && $max < $max2)
153 :     {
154 :     # enclosed
155 :     $olen = $len;
156 :     $t = "1 in 2";
157 :     }
158 :     elsif ($min2 > $min && $max2 < $max)
159 :     {
160 :     # enclosed
161 :     $olen = $len2;
162 :     $t = "2 in 1";
163 :     }
164 :     elsif ($min <= $max2 && $min2 <= $max)
165 :     {
166 :     $olen = min($max - $min2, $max2 - $min);
167 :     $t = "overlap";
168 :     }
169 :    
170 :     my $longest = max($max - $min, $max2 - $min2);
171 :     my $opct = (100 * $olen / $longest);
172 :     push(@poss, [$f, $f2, $fid, $fid2, $olen, $longest, $opct, $t, $min2, $max2]);
173 :     }
174 :    
175 :     @poss = sort { $b->[6] <=> $a->[6] } grep { $_->[6] >= $params->{min_rna_pct_coverage} } @poss;
176 :    
177 :     if (@poss)
178 :     {
179 :     my $n = @poss;
180 :    
181 :     for my $ent (@poss)
182 :     {
183 :     my($f, $f2, $fid, $fid2, $olen, $longest, $opct, $t, $min2, $max2) = @$ent;
184 :     print STDERR "Overlap: $olen $opct $longest $t\n\t$min\t$max\t$fid\t$f->{function}\n\t$min2\t$max2\t$fid2\t$f2->{function}\n";
185 :     }
186 :    
187 :     propagate_feature_data($gto, $user, $event_id, $poss[0]->[1], $f);
188 :     $seen->{$f->{id}}++;
189 :     }
190 :     else
191 :     {
192 :     print STDERR "No candidates with sufficient overlap for $fid $f->{function}\n";
193 :     }
194 :     }
195 :     else
196 :     {
197 :     print STDERR "No candidates for $fid $f->{function}\n";
198 :     }
199 :     }
200 :    
201 :     }
202 :    
203 :     sub propagate_orfs
204 :     {
205 :     my($gto, $params, $user, $event_id, $orfH, $seen) = @_;
206 :    
207 :     while (my($orf, $features) = each %$orfH)
208 :     {
209 :     next if @$features < 2;
210 :     my @fids = map { $_->{id} } @$features;
211 :     my %types = map { tmap($_->{type}) => 1 } @$features;
212 :     if (keys %types > 1)
213 :     {
214 :     my @types = keys %types;
215 :     warn "Skipping propagation for @fids due to multiple types @types\n";
216 :     next;
217 :     }
218 :     my @gb = grep { exists $_->{genbank_feature} } @$features;
219 :     next if @gb == 0;
220 :     if (@gb > 1)
221 :     {
222 :     warn "Skipping propagation for @fids due to multiple entries with genbank data " . join(" ", map { $_->{id} } @gb) . "\n";
223 :     next;
224 :     }
225 :     my $gb = $gb[0];
226 :    
227 :     my @notgb = grep { $_->{id} ne $gb->{id} } @$features;
228 :    
229 :     for my $f (@notgb)
230 :     {
231 :     $seen->{$f->{id}}++;
232 :     propagate_feature_data($gto, $user, $event_id, $gb, $f);
233 :     }
234 :     }
235 :    
236 :     }
237 :    
238 :     sub propagate_feature_data
239 :     {
240 :     my($gto, $user, $event_id, $gb, $f) = @_;
241 :    
242 :     my $prop;
243 :     my %aliases = map { $_ => 1 } @{$f->{aliases}};
244 :     for my $a (@{$gb->{aliases}})
245 :     {
246 :     if (!$aliases{$a})
247 :     {
248 :     push(@{$f->{aliases}}, $a);
249 :     $prop++;
250 :     # print "Propagated $a to $f->{id}\n";
251 :     }
252 :     }
253 :     my %pairs = map { join($;, @$_) => 1 } @{$f->{alias_pairs}};
254 :     for my $ap (@{$gb->{alias_pairs}})
255 :     {
256 :     if (!$pairs{join($;, @$ap)})
257 :     {
258 :     $prop++;
259 :     push(@{$f->{alias_pairs}}, $ap);
260 :     # print "Propagated @$ap to $f->{id}\n";
261 :     }
262 :     }
263 :     if ($prop)
264 :     {
265 :     $gto->add_annotation($f, "$prop aliases propagated from $gb->{id}",
266 :     $user,
267 :     $event_id);
268 :     }
269 :     }
270 :    
271 :     sub tmap
272 :     {
273 :     my($t) = @_;
274 :     return $t eq 'CDS' ? 'peg' : $t;
275 :     }
276 :    
277 :     1;

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3