|
0
|
1 # $Id: flat.pm,v 1.12 2002/10/22 07:38:50 lapp Exp $
|
|
|
2 # BioPerl module for Bio::Variation::IO::flat
|
|
|
3 #
|
|
|
4 # Cared for by Heikki Lehvaslaiho <Heikki@ebi.ac.uk>
|
|
|
5 #
|
|
|
6 # Copyright Heikki Lehvaslaiho
|
|
|
7 #
|
|
|
8 # You may distribute this module under the same terms as perl itself
|
|
|
9 #
|
|
|
10
|
|
|
11 # POD documentation - main docs before the code
|
|
|
12
|
|
|
13 =head1 NAME
|
|
|
14
|
|
|
15 Bio::Variation::IO::flat - flat file sequence variation input/output stream
|
|
|
16
|
|
|
17 =head1 SYNOPSIS
|
|
|
18
|
|
|
19 Do not use this module directly. Use it via the Bio::Variation::IO class.
|
|
|
20
|
|
|
21 =head1 DESCRIPTION
|
|
|
22
|
|
|
23 This object can transform Bio::Variation::SeqDiff objects to and from
|
|
|
24 flat file databases. The format used is EMBL like extension of what is
|
|
|
25 used by the "EBI Mutation Checker" at
|
|
|
26 http://www.ebi.ac.uk/cgi-bin/mutations/check.cgi and will eventually
|
|
|
27 replace it.
|
|
|
28
|
|
|
29 More information of the attributes and values use can be found at
|
|
|
30 http://www.ebi.ac.uk/mutations/recommendations/.
|
|
|
31
|
|
|
32 =head1 FEEDBACK
|
|
|
33
|
|
|
34 =head2 Mailing Lists
|
|
|
35
|
|
|
36 User feedback is an integral part of the evolution of this and other
|
|
|
37 Bioperl modules. Send your comments and suggestions preferably to the
|
|
|
38 Bioperl mailing lists Your participation is much appreciated.
|
|
|
39
|
|
|
40 bioperl-l@bioperl.org - General discussion
|
|
|
41 http://bio.perl.org/MailList.html - About the mailing lists
|
|
|
42
|
|
|
43 =head2 Reporting Bugs
|
|
|
44
|
|
|
45 report bugs to the Bioperl bug tracking system to help us keep track
|
|
|
46 the bugs and their resolution. Bug reports can be submitted via
|
|
|
47 email or the web:
|
|
|
48
|
|
|
49 bioperl-bugs@bio.perl.org
|
|
|
50 http://bugzilla.bioperl.org/
|
|
|
51
|
|
|
52 =head1 AUTHOR - Heikki Lehvaslaiho
|
|
|
53
|
|
|
54 Email: heikki@ebi.ac.uk
|
|
|
55 Address:
|
|
|
56
|
|
|
57 EMBL Outstation, European Bioinformatics Institute
|
|
|
58 Wellcome Trust Genome Campus, Hinxton
|
|
|
59 Cambs. CB10 1SD, United Kingdom
|
|
|
60
|
|
|
61
|
|
|
62
|
|
|
63 =head1 APPENDIX
|
|
|
64
|
|
|
65 The rest of the documentation details each of the object
|
|
|
66 methods. Internal methods are usually preceded with a _
|
|
|
67
|
|
|
68 =cut
|
|
|
69
|
|
|
70 # Let the code begin...
|
|
|
71
|
|
|
72 package Bio::Variation::IO::flat;
|
|
|
73 my $VERSION=1.0;
|
|
|
74 use vars qw(@ISA);
|
|
|
75 use strict;
|
|
|
76
|
|
|
77 use Text::Wrap;
|
|
|
78 use Bio::Variation::IO;
|
|
|
79 use Bio::Variation::SeqDiff;
|
|
|
80 use Bio::Variation::DNAMutation;
|
|
|
81 use Bio::Variation::RNAChange;
|
|
|
82 use Bio::Variation::AAChange;
|
|
|
83 use Bio::Variation::Allele;
|
|
|
84
|
|
|
85
|
|
|
86 @ISA = qw(Bio::Variation::IO);
|
|
|
87
|
|
|
88 sub new {
|
|
|
89 my($class, @args) = @_;
|
|
|
90 my $self = bless {}, $class;
|
|
|
91 $self->_initialize(@args);
|
|
|
92 return $self;
|
|
|
93 }
|
|
|
94
|
|
|
95 sub _initialize {
|
|
|
96 my($self,@args) = @_;
|
|
|
97 return unless $self->SUPER::_initialize(@args);
|
|
|
98 }
|
|
|
99
|
|
|
100 =head2 next
|
|
|
101
|
|
|
102
|
|
|
103 Title : next
|
|
|
104 Usage : $haplo = $stream->next()
|
|
|
105 Function: returns the next seqDiff in the stream
|
|
|
106 Returns : Bio::Variation::SeqDiff object
|
|
|
107 Args : NONE
|
|
|
108
|
|
|
109 =cut
|
|
|
110
|
|
|
111 sub next {
|
|
|
112 my( $self ) = @_;
|
|
|
113 local $/ = '//';
|
|
|
114 return unless my $entry = $self->_readline;
|
|
|
115
|
|
|
116 return if $entry =~ /^\s+$/;
|
|
|
117
|
|
|
118 $entry =~ /\s*ID\s+\S+/ || $self->throw("We do need an ID!");
|
|
|
119
|
|
|
120 my ($id, $offset, $alphabet) = $entry =~ /\s*ID +([^:]+)..(\d+)[^\)]*.\[?([cg])?/
|
|
|
121 or $self->throw("Can't parse ID line");
|
|
|
122 # $self->throw("$1|$2|$3");
|
|
|
123 my $h =Bio::Variation::SeqDiff->new(-id => $id,
|
|
|
124 -offset => $offset,
|
|
|
125 );
|
|
|
126 if ($alphabet) {
|
|
|
127 if ($alphabet eq 'g') {
|
|
|
128 $alphabet = 'dna';
|
|
|
129 }
|
|
|
130 elsif ($alphabet eq 'c') {
|
|
|
131 $alphabet = 'rna';
|
|
|
132 }
|
|
|
133 $h->alphabet($alphabet);
|
|
|
134 }
|
|
|
135 #
|
|
|
136 # DNA
|
|
|
137 #
|
|
|
138 my @dna = split ( / DNA;/, $entry );
|
|
|
139 shift @dna;
|
|
|
140 my $prevdnaobj;
|
|
|
141 foreach my $dna (@dna) {
|
|
|
142 $dna =~ s/Feature[ \t]+//g;
|
|
|
143 ($dna) = split "RNA; ", $dna;
|
|
|
144 #$self->warn("|$dna|") ;
|
|
|
145 #exit;
|
|
|
146 my ($mut_number, $proof, $location, $upflank, $change, $dnflank) =
|
|
|
147 $dna =~ m|\W+([\d\.]+).+/proof: (\w+).+/location: ([^ \n]+).+/upflank: ([ \n\w]+).+/change: ([^ /]+).+/dnflank: ([ \n\w]+)|s;
|
|
|
148 $change =~ s/[ \n]//g;
|
|
|
149 my ($ori, $mut) = split /[>\|]/, $change;
|
|
|
150 my ($variation_number, $change_number) = split /\./, $mut_number;
|
|
|
151 #$self->warn("|$mut_number|>|$variation_number|$change_number|");
|
|
|
152 my $dnamut;
|
|
|
153 if ($change_number and $change_number > 1 ) {
|
|
|
154 my $a3 = Bio::Variation::Allele->new;
|
|
|
155 $a3->seq($mut) if $mut;
|
|
|
156 #$dnamut->add_Allele($a3);
|
|
|
157 $prevdnaobj->add_Allele($a3);
|
|
|
158 } else {
|
|
|
159 $upflank =~ s/[ \n]//g;
|
|
|
160 $dnflank =~ s/[ \n]//g;
|
|
|
161 my ($region, $junk, $region_value, $junk2, $region_dist) =
|
|
|
162 $dna =~ m|.+/region: ([\w\']+)(; )?(\w+)?( ?\(\+?)?(-?\d+)?|s;
|
|
|
163 #my $s = join ("|", $mut_number, $proof, $location, $upflank,
|
|
|
164 # $change, $dnflank, $region, $region_value, $region_dist, $1,$2,$3,$4,$5);
|
|
|
165 #$self->warn($s);
|
|
|
166 #exit;
|
|
|
167 my ($start, $sep, $end) = $location =~ /(-?\d+)(.)?\D?(-?\d+)?/;
|
|
|
168 $end = $start if not $end ;
|
|
|
169 my ($len) = $end - $start +1;
|
|
|
170 $len = 0, $start = $end if defined $sep and $sep eq '^';
|
|
|
171 my $ismut = 0;
|
|
|
172 $ismut = 1 if $change =~ m/>/;
|
|
|
173
|
|
|
174 $dnamut = Bio::Variation::DNAMutation->new
|
|
|
175 ('-start' => $start,
|
|
|
176 '-end' => $end,
|
|
|
177 '-length' => $len,
|
|
|
178 '-upStreamSeq' => $upflank,
|
|
|
179 '-dnStreamSeq' => $dnflank,
|
|
|
180 '-proof' => $proof,
|
|
|
181 '-mut_number' => $mut_number
|
|
|
182 );
|
|
|
183 $prevdnaobj = $dnamut;
|
|
|
184 my $a1 = Bio::Variation::Allele->new;
|
|
|
185 $a1->seq($ori) if $ori;
|
|
|
186 $dnamut->allele_ori($a1);
|
|
|
187 my $a2 = Bio::Variation::Allele->new;
|
|
|
188 $a2->seq($mut) if $mut;
|
|
|
189 $dnamut->add_Allele($a2);
|
|
|
190 if ($ismut) {
|
|
|
191 $dnamut->isMutation(1);
|
|
|
192 $dnamut->allele_mut($a2);
|
|
|
193 }
|
|
|
194 $dnamut->region($region) if defined $region;
|
|
|
195 $dnamut->region_value($region_value) if defined $region_value;
|
|
|
196 $dnamut->region_dist($region_dist) if defined $region_dist;
|
|
|
197
|
|
|
198 $h->add_Variant($dnamut);
|
|
|
199 $dnamut->SeqDiff($h);
|
|
|
200 }
|
|
|
201 }
|
|
|
202
|
|
|
203 #
|
|
|
204 # RNA
|
|
|
205 #
|
|
|
206 my @rna = split ( / RNA;/, $entry );
|
|
|
207 shift @rna;
|
|
|
208 my $prevrnaobj;
|
|
|
209 foreach my $rna (@rna) {
|
|
|
210 $rna = substr ($rna, 0, index($rna, 'Feature AA'));
|
|
|
211 $rna =~ s/Feature[ \t]+//g;
|
|
|
212 ($rna) = split "DNA; ", $rna;
|
|
|
213 #$self->warn("|$rna|") ;
|
|
|
214 my ($mut_number, $proof, $location, $upflank, $change, $dnflank) =
|
|
|
215 $rna =~ m|\W+([\d\.]+).+/proof: (\w+).+/location: ([^ \n]+).+/upflank: (\w+).+/change: ([^/]+).+/dnflank: (\w+)|s ;#'
|
|
|
216 my ($region, $junk, $region_value, $junk2, $region_dist) =
|
|
|
217 $rna =~ m|.+/region: ([\w\']+)(; )?(\w+)?( ?\(\+?)?(-?\d+)?|s;
|
|
|
218 #my $s = join ("|", $mut_number, $proof, $location, $upflank,
|
|
|
219 # $change, $dnflank, $region, $region_value, $region_dist, $1,$2,$3,$4,$5);
|
|
|
220 #$self->warn($s);
|
|
|
221 #exit;
|
|
|
222 $change =~ s/[ \n]//g;
|
|
|
223 my ($ori, $mut) = split /[>\|]/, $change;
|
|
|
224 my $rnamut;
|
|
|
225 my ($variation_number, $change_number) = split /\./, $mut_number;
|
|
|
226 if ($change_number and $change_number > 1 ) {
|
|
|
227 my $a3 = Bio::Variation::Allele->new;
|
|
|
228 $a3->seq($mut) if $mut;
|
|
|
229 #$rnamut->add_Allele($a3);
|
|
|
230 $prevrnaobj->add_Allele($a3);
|
|
|
231 } else {
|
|
|
232 my ($start, $sep, $end) = $location =~ /(-?\d+)(.)?\D?(-?\d+)?/;
|
|
|
233 $end = $start if not $end ;
|
|
|
234 my ($len) = $end - $start + 1;
|
|
|
235 $len = 0, $start = $end if defined $sep and $sep eq '^';
|
|
|
236 my $ismut;
|
|
|
237 $ismut = 1 if $change =~ m/>/;
|
|
|
238 my ($codon_table) = $rna =~ m|.+/codon_table: (\d+)|s;
|
|
|
239 my ($codon_pos) = $rna =~ m|.+/codon:[^;]+; ([123])|s;
|
|
|
240
|
|
|
241 $rnamut = Bio::Variation::RNAChange->new
|
|
|
242 ('-start' => $start,
|
|
|
243 '-end' => $end,
|
|
|
244 '-length' => $len,
|
|
|
245 '-upStreamSeq' => $upflank,
|
|
|
246 '-dnStreamSeq' => $dnflank,
|
|
|
247 '-proof' => $proof,
|
|
|
248 '-mut_number' => $mut_number
|
|
|
249
|
|
|
250 );
|
|
|
251 $prevrnaobj = $rnamut;
|
|
|
252 my $a1 = Bio::Variation::Allele->new;
|
|
|
253 $a1->seq($ori) if $ori;
|
|
|
254 $rnamut->allele_ori($a1);
|
|
|
255 my $a2 = Bio::Variation::Allele->new;
|
|
|
256 $a2->seq($mut) if $mut;
|
|
|
257 $rnamut->add_Allele($a2);
|
|
|
258 if ($ismut) {
|
|
|
259 $rnamut->isMutation(1);
|
|
|
260 $rnamut->allele_mut($a2);
|
|
|
261 }
|
|
|
262 $rnamut->region($region) if defined $region;
|
|
|
263 $rnamut->region_value($region_value) if defined $region_value;
|
|
|
264 $rnamut->region_dist($region_dist) if defined $region_dist;
|
|
|
265
|
|
|
266 $rnamut->codon_table($codon_table) if $codon_table;
|
|
|
267 $rnamut->codon_pos($codon_pos) if $codon_pos;
|
|
|
268 $h->add_Variant($rnamut);
|
|
|
269 foreach my $mut ($h->each_Variant) {
|
|
|
270 if ($mut->isa('Bio::Variation::DNAMutation') ) {
|
|
|
271 if ($mut->mut_number == $rnamut->mut_number) {
|
|
|
272 $rnamut->DNAMutation($mut);
|
|
|
273 $mut->RNAChange($rnamut);
|
|
|
274 }
|
|
|
275 }
|
|
|
276 }
|
|
|
277 }
|
|
|
278 }
|
|
|
279 #
|
|
|
280 # AA
|
|
|
281 #
|
|
|
282 my @aa = split ( / AA;/, $entry );
|
|
|
283 shift @aa;
|
|
|
284 my $prevaaobj;
|
|
|
285 foreach my $aa (@aa) {
|
|
|
286 $aa = substr ($aa, 0, index($aa, 'Feature AA'));
|
|
|
287 $aa =~ s/Feature[ \t]+//g;
|
|
|
288 ($aa) = split "DNA; ", $aa;
|
|
|
289 #$self->warn("|$aa|") ;
|
|
|
290 my ($mut_number, $proof, $location, $change) =
|
|
|
291 $aa =~ m|\W+([\d\.]+).+/proof: (\w+).+/location: ([^ \n]+)./change: ([^/;]+)|s;
|
|
|
292 $change =~ s/[ \n]//g;
|
|
|
293 #my $s = join ("|", $mut_number, $proof, $location, $change);
|
|
|
294 #$self->warn($s);
|
|
|
295 #exit;
|
|
|
296 $change =~ s/[ \n]//g;
|
|
|
297 $change =~ s/DNA$//;
|
|
|
298 my ($ori, $mut) = split /[>\|]/, $change;
|
|
|
299 #print "------$location----$ori-$mut-------------\n";
|
|
|
300 my ($variation_number, $change_number) = split /\./, $mut_number;
|
|
|
301 my $aamut;
|
|
|
302 if ($change_number and $change_number > 1 ) {
|
|
|
303 my $a3 = Bio::Variation::Allele->new;
|
|
|
304 $a3->seq($mut) if $mut;
|
|
|
305 $prevaaobj->add_Allele($a3);
|
|
|
306 } else {
|
|
|
307 my ($start, $sep, $end) = $location =~ /(-?\d+)(.)?\D?(-?\d+)?/;
|
|
|
308 $end = $start if not $end ;
|
|
|
309 my ($len) = $end - $start + 1;
|
|
|
310 $len = 0, $start = $end if defined $sep and $sep eq '^';
|
|
|
311 my $ismut;
|
|
|
312 $ismut = 1 if $change =~ m/>/;
|
|
|
313 my ($region) = $aa =~ m|.+/region: (\w+)|s ;
|
|
|
314 $aamut = Bio::Variation::AAChange->new
|
|
|
315 ('-start' => $start,
|
|
|
316 '-end' => $end,
|
|
|
317 '-length' => $len,
|
|
|
318 '-proof' => $proof,
|
|
|
319 '-mut_number' => $mut_number
|
|
|
320 );
|
|
|
321 $prevaaobj = $aamut;
|
|
|
322 my $a1 = Bio::Variation::Allele->new;
|
|
|
323 $a1->seq($ori) if $ori;
|
|
|
324 $aamut->allele_ori($a1);
|
|
|
325 my $a2 = Bio::Variation::Allele->new;
|
|
|
326 $a2->seq($mut) if $mut;
|
|
|
327 $aamut->add_Allele($a2);
|
|
|
328 if ($ismut) {
|
|
|
329 $aamut->isMutation(1);
|
|
|
330 $aamut->allele_mut($a2);
|
|
|
331 }
|
|
|
332 $region && $aamut->region($region);
|
|
|
333 $h->add_Variant($aamut);
|
|
|
334 foreach my $mut ($h->each_Variant) {
|
|
|
335 if ($mut->isa('Bio::Variation::RNAChange') ) {
|
|
|
336 if ($mut->mut_number == $aamut->mut_number) {
|
|
|
337 $aamut->RNAChange($mut);
|
|
|
338 $mut->AAChange($aamut);
|
|
|
339 }
|
|
|
340 }
|
|
|
341 }
|
|
|
342
|
|
|
343 }
|
|
|
344 }
|
|
|
345 return $h;
|
|
|
346 }
|
|
|
347
|
|
|
348 =head2 write
|
|
|
349
|
|
|
350 Title : write
|
|
|
351 Usage : $stream->write(@seqDiffs)
|
|
|
352 Function: writes the $seqDiff object into the stream
|
|
|
353 Returns : 1 for success and 0 for error
|
|
|
354 Args : Bio::Variation::SeqDiff object
|
|
|
355
|
|
|
356
|
|
|
357 =cut
|
|
|
358
|
|
|
359 sub write {
|
|
|
360 my ($self,@h) = @_;
|
|
|
361
|
|
|
362 #$columns = 75; #default for Text::Wrap
|
|
|
363 my %tag =
|
|
|
364 (
|
|
|
365 'ID' => 'ID ',
|
|
|
366 'Description' => 'Description ',
|
|
|
367 'FeatureKey' => 'Feature ',
|
|
|
368 'FeatureQual' => "Feature ",
|
|
|
369 'FeatureWrap' => "Feature ",
|
|
|
370 'ErrorComment' => 'Comment '
|
|
|
371 #'Comment' => 'Comment -!-',
|
|
|
372 #'CommentLine' => 'Comment ',
|
|
|
373 );
|
|
|
374
|
|
|
375 if( !defined $h[0] ) {
|
|
|
376 $self->throw("Attempting to write with no information!");
|
|
|
377 }
|
|
|
378
|
|
|
379 foreach my $h (@h) {
|
|
|
380
|
|
|
381 my @entry =();
|
|
|
382
|
|
|
383 my ($text, $tmp, $tmp2, $sep);
|
|
|
384 my ($count) = 0;
|
|
|
385
|
|
|
386
|
|
|
387 $text = $tag{ID};
|
|
|
388
|
|
|
389 $text .= $h->id;
|
|
|
390 $text .= ":(". $h->offset;
|
|
|
391 $text .= "+1" if $h->sysname =~ /-/;
|
|
|
392 $text .= ")". $h->sysname;
|
|
|
393 $text .= "; ". $h->trivname if $h->trivname;
|
|
|
394 push (@entry, $text);
|
|
|
395
|
|
|
396 #Variants need to be ordered accoding to mutation_number attribute
|
|
|
397 #put them into a hash of arrays holding the Variant objects
|
|
|
398 #This is necessary for cases like several distict mutations present
|
|
|
399 # in the same sequence.
|
|
|
400 my @allvariants = $h->each_Variant;
|
|
|
401 my %variants = ();
|
|
|
402 foreach my $mut ($h->each_Variant) {
|
|
|
403 push @{$variants{$mut->mut_number} }, $mut;
|
|
|
404 }
|
|
|
405 #my ($variation_number, $change_number) = split /\./, $mut_number;
|
|
|
406 foreach my $var (sort keys %variants) {
|
|
|
407 #print $var, ": ", join (" ", @{$variants{$var}}), "\n";
|
|
|
408
|
|
|
409 foreach my $mut (@{$variants{$var}}) {
|
|
|
410 #
|
|
|
411 # DNA
|
|
|
412 #
|
|
|
413 if ( $mut->isa('Bio::Variation::DNAMutation') ) {
|
|
|
414 #collect all non-reference alleles
|
|
|
415 $self->throw("allele_ori needs to be defined in [$mut]")
|
|
|
416 if not $mut->allele_ori;
|
|
|
417 if ($mut->isMutation) {
|
|
|
418 $sep = '>';
|
|
|
419 } else {
|
|
|
420 $sep = '|';
|
|
|
421 }
|
|
|
422 my @alleles = $mut->each_Allele;
|
|
|
423 #push @alleles, $mut->allele_mut if $mut->allele_mut;
|
|
|
424 my $count = 0; # two alleles
|
|
|
425 foreach my $allele (@alleles) {
|
|
|
426 $count++;
|
|
|
427 my ($variation_number, $change_number) = split /\./, $mut->mut_number;
|
|
|
428 if ($change_number and $change_number != $count){
|
|
|
429 $mut->mut_number("$change_number.$count");
|
|
|
430 }
|
|
|
431 $mut->allele_mut($allele);
|
|
|
432 push (@entry,
|
|
|
433 $tag{FeatureKey}. 'DNA'. "; ". $mut->mut_number
|
|
|
434 );
|
|
|
435 #label
|
|
|
436 $text=$tag{FeatureQual}. '/label: '. $mut->label;
|
|
|
437 push (@entry, $text);
|
|
|
438
|
|
|
439 #proof
|
|
|
440 if ($mut->proof) {
|
|
|
441 $text = $tag{FeatureQual}. '/proof: '. $mut->proof;
|
|
|
442 push (@entry, $text) ;
|
|
|
443 }
|
|
|
444 #location
|
|
|
445 $text = $tag{FeatureQual}. '/location: ';
|
|
|
446 #$mut->id. '; '. $mut->start;
|
|
|
447 if ($mut->length > 1 ) {# if ($mut->end - $mut->start ) {
|
|
|
448 my $l = $mut->start + $mut->length -1;
|
|
|
449 $text .= $mut->start. '..'. $l;
|
|
|
450 }
|
|
|
451 elsif ($mut->length == 0) {
|
|
|
452 my $tmp_start = $mut->start - 1;
|
|
|
453 $tmp_start-- if $tmp_start == 0;
|
|
|
454 $text .= $tmp_start. '^'. $mut->end;
|
|
|
455 } else {
|
|
|
456 $text .= $mut->start;
|
|
|
457 }
|
|
|
458
|
|
|
459 if ($h->alphabet && $h->alphabet eq 'dna') {
|
|
|
460 $tmp = $mut->start + $h->offset;
|
|
|
461 $tmp-- if $tmp <= 0;
|
|
|
462 $mut->start < 1 && $tmp++;
|
|
|
463 #$text.= ' ('. $h->id. '::'. $tmp;
|
|
|
464 $tmp2 = $mut->end + $h->offset;
|
|
|
465 if ( $mut->length > 1 ) {
|
|
|
466 $mut->end < 1 && $tmp2++;
|
|
|
467 $text.= ' ('. $h->id. '::'. $tmp. "..". $tmp2;
|
|
|
468 }
|
|
|
469 elsif ($mut->length == 0) {
|
|
|
470 $tmp--;
|
|
|
471 $tmp-- if $tmp == 0;
|
|
|
472 $text .= ' ('. $h->id. '::'. $tmp. '^'. $tmp2;
|
|
|
473 } else {
|
|
|
474 $text.= ' ('. $h->id. '::'. $tmp;
|
|
|
475 }
|
|
|
476 $text .= ')';
|
|
|
477 }
|
|
|
478 push (@entry, $text);
|
|
|
479 #sequence
|
|
|
480 push (@entry,
|
|
|
481 $tag{FeatureQual}. '/upflank: '. $mut->upStreamSeq
|
|
|
482 );
|
|
|
483 $text = '';
|
|
|
484 $text = $mut->allele_ori->seq if $mut->allele_ori->seq;
|
|
|
485 $text .= $sep;
|
|
|
486 $text .= $mut->allele_mut->seq if $mut->allele_mut->seq;
|
|
|
487 push (@entry,
|
|
|
488 wrap($tag{FeatureQual}. '/change: ', $tag{FeatureWrap},
|
|
|
489 $text)
|
|
|
490 );
|
|
|
491
|
|
|
492 push (@entry,
|
|
|
493 $tag{FeatureQual}. '/dnflank: '. $mut->dnStreamSeq
|
|
|
494 );
|
|
|
495 #restriction enzyme
|
|
|
496 if ($mut->restriction_changes ne '') {
|
|
|
497 $text = $mut->restriction_changes;
|
|
|
498 $text = wrap($tag{FeatureQual}. '/re_site: ', $tag{FeatureWrap}, $text);
|
|
|
499 push (@entry,
|
|
|
500 $text
|
|
|
501 );
|
|
|
502 }
|
|
|
503 #region
|
|
|
504 if ($mut->region ) {
|
|
|
505 $text = $tag{FeatureQual}. '/region: '. $mut->region;
|
|
|
506 $text .= ';' if $mut->region_value or $mut->region_dist;
|
|
|
507 $text .= ' '. $mut->region_value if $mut->region_value;
|
|
|
508 if ($mut->region_dist ) {
|
|
|
509 $tmp = '';
|
|
|
510 $tmp = '+' if $mut->region_dist > 1;
|
|
|
511 $text .= " (". $tmp. $mut->region_dist. ')';
|
|
|
512 }
|
|
|
513 push (@entry, $text);
|
|
|
514 }
|
|
|
515 #CpG
|
|
|
516 if ($mut->CpG) {
|
|
|
517 push (@entry,
|
|
|
518 $tag{FeatureQual}. "/CpG"
|
|
|
519 );
|
|
|
520 }
|
|
|
521 }
|
|
|
522 }
|
|
|
523 #
|
|
|
524 # RNA
|
|
|
525 #
|
|
|
526 elsif ($mut->isa('Bio::Variation::RNAChange') ) {
|
|
|
527 #collect all non-reference alleles
|
|
|
528 $self->throw("allele_ori needs to be defined in [$mut]")
|
|
|
529 if not $mut->allele_ori;
|
|
|
530 my @alleles = $mut->each_Allele;
|
|
|
531 #push @alleles, $mut->allele_mut if $mut->allele_mut;
|
|
|
532 if ($mut->isMutation) {
|
|
|
533 $sep = '>';
|
|
|
534 } else {
|
|
|
535 $sep = '|';
|
|
|
536 }
|
|
|
537
|
|
|
538 my $count = 0; # two alleles
|
|
|
539 foreach my $allele (@alleles) {
|
|
|
540 $count++;
|
|
|
541 my ($variation_number, $change_number) = split /\./, $mut->mut_number;
|
|
|
542 if ($change_number and $change_number != $count){
|
|
|
543 $mut->mut_number("$change_number.$count");
|
|
|
544 }
|
|
|
545 $mut->allele_mut($allele);
|
|
|
546 push (@entry,
|
|
|
547 $tag{FeatureKey}. 'RNA'. "; ". $mut->mut_number
|
|
|
548 );
|
|
|
549 #label
|
|
|
550 $text=$tag{FeatureQual}. '/label: '. $mut->label;
|
|
|
551 push (@entry, $text);
|
|
|
552 #proof
|
|
|
553 if ($mut->proof) {
|
|
|
554 $text = $tag{FeatureQual}. '/proof: '. $mut->proof;
|
|
|
555 push (@entry, $text) ;
|
|
|
556 }
|
|
|
557 #location
|
|
|
558 $text = $tag{FeatureQual}. '/location: ' ;
|
|
|
559 if ($mut->length > 1 ) {
|
|
|
560 $text .= $mut->start. '..'. $mut->end;
|
|
|
561 $tmp2 = $mut->end + $h->offset;
|
|
|
562 }
|
|
|
563 elsif ($mut->length == 0) {
|
|
|
564 my $tmp_start = $mut->start;
|
|
|
565 $tmp_start--;
|
|
|
566 $tmp_start-- if $tmp_start == 0;
|
|
|
567 $text .= $tmp_start. '^'. $mut->end;
|
|
|
568 } else {
|
|
|
569 $text .= $mut->start;
|
|
|
570 }
|
|
|
571
|
|
|
572 if ($h->alphabet && $h->alphabet eq 'rna') {
|
|
|
573 $tmp = $mut->start + $h->offset;
|
|
|
574 $tmp-- if $tmp <= 0;
|
|
|
575 #$mut->start < 1 && $tmp++;
|
|
|
576 #$text.= ' ('. $h->id. '::'. $tmp;
|
|
|
577 $tmp2 = $mut->end + $h->offset;
|
|
|
578 #$mut->end < 1 && $tmp2++;
|
|
|
579 if ( $mut->length > 1 ) {
|
|
|
580 $text.= ' ('. $h->id. '::'. $tmp. "..". $tmp2;
|
|
|
581 }
|
|
|
582 elsif ($mut->length == 0) {
|
|
|
583 $tmp--;
|
|
|
584 $text .= ' ('. $h->id. '::'. $tmp. '^'. $tmp2;
|
|
|
585 } else {
|
|
|
586 $text.= ' ('. $h->id. '::'. $tmp;
|
|
|
587 }
|
|
|
588
|
|
|
589 $text .= ')';
|
|
|
590 }
|
|
|
591 push (@entry, $text);
|
|
|
592
|
|
|
593 #sequence
|
|
|
594 push (@entry,
|
|
|
595 $tag{FeatureQual}. '/upflank: '. $mut->upStreamSeq
|
|
|
596 );
|
|
|
597 $text = '';
|
|
|
598 $text = $mut->allele_ori->seq if $mut->allele_ori->seq;
|
|
|
599 $text .= $sep;
|
|
|
600 $text .= $mut->allele_mut->seq if $mut->allele_mut->seq;
|
|
|
601 push (@entry,
|
|
|
602 wrap($tag{FeatureQual}. '/change: ', $tag{FeatureWrap},
|
|
|
603 $text)
|
|
|
604 );
|
|
|
605 push (@entry,
|
|
|
606 $tag{FeatureQual}. '/dnflank: '. $mut->dnStreamSeq
|
|
|
607 );
|
|
|
608 #restriction
|
|
|
609 if ($mut->restriction_changes ne '') {
|
|
|
610 $text = $mut->restriction_changes;
|
|
|
611 $text = wrap($tag{FeatureQual}. '/re_site: ', $tag{FeatureWrap}, $text);
|
|
|
612 push (@entry,
|
|
|
613 $text
|
|
|
614 );
|
|
|
615 }
|
|
|
616 #coding
|
|
|
617 if ($mut->region eq 'coding') {
|
|
|
618 #codon table
|
|
|
619 $text = $tag{FeatureQual}. '/codon_table: ';
|
|
|
620 $text .= $mut->codon_table;
|
|
|
621 push (@entry, $text);
|
|
|
622 #codon
|
|
|
623
|
|
|
624 $text = $tag{FeatureQual}. '/codon: '. $mut->codon_ori. $sep;
|
|
|
625 if ($mut->DNAMutation->label =~ /.*point/) {
|
|
|
626 $text .= $mut->codon_mut;
|
|
|
627 }
|
|
|
628 else {
|
|
|
629 $text .= '-';
|
|
|
630 }
|
|
|
631 $text .= "; ". $mut->codon_pos;
|
|
|
632 push (@entry, $text);
|
|
|
633 }
|
|
|
634 #region
|
|
|
635 if ($mut->region ) {
|
|
|
636 $text = $tag{FeatureQual}. '/region: '. $mut->region;
|
|
|
637 $text .= ';' if $mut->region_value or $mut->region_dist;
|
|
|
638 $text .= ' '. $mut->region_value if $mut->region_value;
|
|
|
639 if ($mut->region_dist ) {
|
|
|
640 $tmp = '';
|
|
|
641 $tmp = '+' if $mut->region_dist > 1;
|
|
|
642 $text .= " (". $tmp. $mut->region_dist. ')';
|
|
|
643 }
|
|
|
644 push (@entry, $text);
|
|
|
645 }
|
|
|
646 }
|
|
|
647 }
|
|
|
648 #
|
|
|
649 # AA
|
|
|
650 #
|
|
|
651 elsif ($mut->isa('Bio::Variation::AAChange')) {
|
|
|
652 #collect all non-reference alleles
|
|
|
653 $self->throw("allele_ori needs to be defined in [$mut]")
|
|
|
654 if not $mut->allele_ori;
|
|
|
655 if ($mut->isMutation) {
|
|
|
656 $sep = '>';
|
|
|
657 } else {
|
|
|
658 $sep = '|';
|
|
|
659 }
|
|
|
660 my @alleles = $mut->each_Allele;
|
|
|
661 #push @alleles, $mut->allele_mut if $mut->allele_mut;
|
|
|
662 my $count = 0; # two alleles
|
|
|
663 foreach my $allele (@alleles) {
|
|
|
664 $count++;
|
|
|
665 my ($variation_number, $change_number) = split /\./, $mut->mut_number;
|
|
|
666 if ($change_number and $change_number != $count){
|
|
|
667 $mut->mut_number("$change_number.$count");
|
|
|
668 }
|
|
|
669 $mut->allele_mut($allele);
|
|
|
670 push (@entry,
|
|
|
671 $tag{FeatureKey}. 'AA'. "; ". $mut->mut_number
|
|
|
672 );
|
|
|
673 #label
|
|
|
674 $text=$tag{FeatureQual}. '/label: '. $mut->label;
|
|
|
675 push (@entry, $text) ;
|
|
|
676 #proof
|
|
|
677 if ($mut->proof) {
|
|
|
678 $text = $tag{FeatureQual}. '/proof: '. $mut->proof;
|
|
|
679 push (@entry, $text) ;
|
|
|
680 }
|
|
|
681 #location
|
|
|
682 $text = $tag{FeatureQual}. '/location: '.
|
|
|
683 #$mut->id. '; '. $mut->start;
|
|
|
684 $mut->start;
|
|
|
685 if ($mut->length > 1 ) {
|
|
|
686 $tmp = $mut->start + $mut->length -1;
|
|
|
687 $text .= '..'. $tmp;
|
|
|
688 }
|
|
|
689 push (@entry, $text);
|
|
|
690 #sequence
|
|
|
691 $text = '';
|
|
|
692 $text = $mut->allele_ori->seq if $mut->allele_ori->seq;
|
|
|
693 $text .= $sep;
|
|
|
694 $text .= $mut->allele_mut->seq if $mut->allele_mut->seq;
|
|
|
695 push (@entry,
|
|
|
696 wrap($tag{FeatureQual}. '/change: ', $tag{FeatureWrap},
|
|
|
697 $text)
|
|
|
698 );
|
|
|
699 #region
|
|
|
700 if ($mut->region ) {
|
|
|
701 $text = $tag{FeatureQual}. '/region: '. $mut->region;
|
|
|
702 $text .= ';' if $mut->region_value or $mut->region_dist;
|
|
|
703 $text .= ' '. $mut->region_value if $mut->region_value;
|
|
|
704 if ($mut->region_dist ) {
|
|
|
705 $tmp = '';
|
|
|
706 $tmp = '+' if $mut->region_dist > 1;
|
|
|
707 $text .= " (". $tmp. $mut->region_dist. ')';
|
|
|
708 }
|
|
|
709 push (@entry, $text);
|
|
|
710 }
|
|
|
711 }
|
|
|
712 }
|
|
|
713 }
|
|
|
714 }
|
|
|
715 push (@entry,
|
|
|
716 "//"
|
|
|
717 );
|
|
|
718 my $str = join ("\n", @entry). "\n";
|
|
|
719 $str =~ s/\t/ /g;
|
|
|
720 $self->_print($str);
|
|
|
721 }
|
|
|
722 return 1;
|
|
|
723 }
|
|
|
724
|
|
|
725 1;
|
