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

Annotation of /FigKernelPackages/PropagateGBMetadata.pm

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1.2 - (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 :     my $tree = $trees{$ctg};
115 :     $tree->insert($f->{id}, $min, $max);
116 :     }
117 :    
118 :     #
119 :     # Now walk the rna features looking for overlaps with
120 :     for my $f ($gto->features)
121 :     {
122 :     my $fid = $f->{id};
123 :    
124 :     next unless $f->{type} =~ /rna/i;
125 :     next if $seen->{$fid};
126 :     next if exists $f->{genbank_feature};
127 :    
128 :     my($ctg, $min, $max, $dir, $len) = GenomeTypeObject::bounds($f);
129 :     my $tree = $trees{$ctg};
130 :     my $overlap = $tree->fetch($min, $max);
131 :     my @overlap = grep { $_ ne $fid } @$overlap;
132 :     #
133 :     # @overlap contains the genbank feature we are mapping from.
134 :     #
135 :     if (@overlap)
136 :     {
137 :     my @poss;
138 :     for my $fid2 (@overlap)
139 :     {
140 :     my $f2 = $gto->find_feature($fid2);
141 :     my($ctg2, $min2, $max2, $dir2, $len2) = GenomeTypeObject::bounds($f2);
142 :    
143 :     my $olen;
144 :     my $t = "not";
145 :     if ($min > $min2 && $max < $max2)
146 :     {
147 :     # enclosed
148 :     $olen = $len;
149 :     $t = "1 in 2";
150 :     }
151 :     elsif ($min2 > $min && $max2 < $max)
152 :     {
153 :     # enclosed
154 :     $olen = $len2;
155 :     $t = "2 in 1";
156 :     }
157 :     elsif ($min <= $max2 && $min2 <= $max)
158 :     {
159 :     $olen = min($max - $min2, $max2 - $min);
160 :     $t = "overlap";
161 :     }
162 :    
163 :     my $longest = max($max - $min, $max2 - $min2);
164 :     my $opct = (100 * $olen / $longest);
165 :     push(@poss, [$f, $f2, $fid, $fid2, $olen, $longest, $opct, $t, $min2, $max2]);
166 :     }
167 :    
168 :     @poss = sort { $b->[6] <=> $a->[6] } grep { $_->[6] >= $params->{min_rna_pct_coverage} } @poss;
169 :    
170 :     if (@poss)
171 :     {
172 :     my $n = @poss;
173 :    
174 :     for my $ent (@poss)
175 :     {
176 :     my($f, $f2, $fid, $fid2, $olen, $longest, $opct, $t, $min2, $max2) = @$ent;
177 :     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";
178 :     }
179 :    
180 :     propagate_feature_data($gto, $user, $event_id, $poss[0]->[1], $f);
181 :     $seen->{$f->{id}}++;
182 :     }
183 :     else
184 :     {
185 :     print STDERR "No candidates with sufficient overlap for $fid $f->{function}\n";
186 :     }
187 :     }
188 :     else
189 :     {
190 :     print STDERR "No candidates for $fid $f->{function}\n";
191 :     }
192 :     }
193 :    
194 :     }
195 :    
196 :     sub propagate_orfs
197 :     {
198 :     my($gto, $params, $user, $event_id, $orfH, $seen) = @_;
199 :    
200 :     while (my($orf, $features) = each %$orfH)
201 :     {
202 :     next if @$features < 2;
203 :     my @fids = map { $_->{id} } @$features;
204 :     my %types = map { tmap($_->{type}) => 1 } @$features;
205 :     if (keys %types > 1)
206 :     {
207 :     my @types = keys %types;
208 :     warn "Skipping propagation for @fids due to multiple types @types\n";
209 :     next;
210 :     }
211 :     my @gb = grep { exists $_->{genbank_feature} } @$features;
212 :     next if @gb == 0;
213 :     if (@gb > 1)
214 :     {
215 :     warn "Skipping propagation for @fids due to multiple entries with genbank data " . join(" ", map { $_->{id} } @gb) . "\n";
216 :     next;
217 :     }
218 :     my $gb = $gb[0];
219 :    
220 :     my @notgb = grep { $_->{id} ne $gb->{id} } @$features;
221 :    
222 :     for my $f (@notgb)
223 :     {
224 :     $seen->{$f->{id}}++;
225 :     propagate_feature_data($gto, $user, $event_id, $gb, $f);
226 :     }
227 :     }
228 :    
229 :     }
230 :    
231 :     sub propagate_feature_data
232 :     {
233 :     my($gto, $user, $event_id, $gb, $f) = @_;
234 :    
235 :     my $prop;
236 :     my %aliases = map { $_ => 1 } @{$f->{aliases}};
237 :     for my $a (@{$gb->{aliases}})
238 :     {
239 :     if (!$aliases{$a})
240 :     {
241 :     push(@{$f->{aliases}}, $a);
242 :     $prop++;
243 :     # print "Propagated $a to $f->{id}\n";
244 :     }
245 :     }
246 :     my %pairs = map { join($;, @$_) => 1 } @{$f->{alias_pairs}};
247 :     for my $ap (@{$gb->{alias_pairs}})
248 :     {
249 :     if (!$pairs{join($;, @$ap)})
250 :     {
251 :     $prop++;
252 :     push(@{$f->{alias_pairs}}, $ap);
253 :     # print "Propagated @$ap to $f->{id}\n";
254 :     }
255 :     }
256 :     if ($prop)
257 :     {
258 :     $gto->add_annotation($f, "$prop aliases propagated from $gb->{id}",
259 :     $user,
260 :     $event_id);
261 :     }
262 :     }
263 :    
264 :     sub tmap
265 :     {
266 :     my($t) = @_;
267 :     return $t eq 'CDS' ? 'peg' : $t;
268 :     }
269 :    
270 :     1;

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3