File Coverage

Bio/Variation/IO/flat.pm
Criterion Covered Total %
statement 354 369 95.9
branch 167 214 78.0
condition 34 42 80.9
subroutine 12 12 100.0
pod 3 3 100.0
total 570 640 89.0


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