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   9 use strict;
  2         3  
  2         53  
73              
74 2     2   411 use Text::Wrap;
  2         2049  
  2         88  
75 2     2   284 use Bio::Variation::SeqDiff;
  2         4  
  2         43  
76 2     2   270 use Bio::Variation::DNAMutation;
  2         2  
  2         58  
77 2     2   287 use Bio::Variation::RNAChange;
  2         2  
  2         46  
78 2     2   278 use Bio::Variation::AAChange;
  2         4  
  2         45  
79 2     2   260 use Bio::Variation::Allele;
  2         3  
  2         57  
80              
81              
82 2     2   8 use base qw(Bio::Variation::IO);
  2         4  
  2         5837  
83              
84             sub new {
85 5     5 1 11 my($class, @args) = @_;
86 5         7 my $self = bless {}, $class;
87 5         10 $self->_initialize(@args);
88 5         16 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         40 local $/ = '//';
110 17 50       34 return unless my $entry = $self->_readline;
111            
112 17 100       59 return if $entry =~ /^\s+$/;
113              
114 15 50       46 $entry =~ /\s*ID\s+\S+/ || $self->throw("We do need an ID!");
115              
116 15 50       69 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         55 my $h =Bio::Variation::SeqDiff->new(-id => $id,
120             -offset => $offset,
121             );
122 15 50       20 if ($alphabet) {
123 15 100       28 if ($alphabet eq 'g') {
    50          
124 4         4 $alphabet = 'dna';
125             }
126             elsif ($alphabet eq 'c') {
127 11         11 $alphabet = 'rna';
128             }
129 15         20 $h->alphabet($alphabet);
130             }
131             #
132             # DNA
133             #
134 15         57 my @dna = split ( / DNA;/, $entry );
135 15         12 shift @dna;
136 15         13 my $prevdnaobj;
137 15         19 foreach my $dna (@dna) {
138 17         258 $dna =~ s/Feature[ \t]+//g;
139 17         42 ($dna) = split "RNA; ", $dna;
140             #$self->warn("|$dna|") ;
141             #exit;
142 17         177 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         18 $change =~ s/[ \n]//g;
145 17         33 my ($ori, $mut) = split /[>\|]/, $change;
146 17         17 my ($variation_number, $change_number) = split /\./, $mut_number;
147             #$self->warn("|$mut_number|>|$variation_number|$change_number|");
148 17         9 my $dnamut;
149 17 100 100     38 if ($change_number and $change_number > 1 ) {
150 1         3 my $a3 = Bio::Variation::Allele->new;
151 1 50       4 $a3->seq($mut) if $mut;
152             #$dnamut->add_Allele($a3);
153 1         3 $prevdnaobj->add_Allele($a3);
154             } else {
155 16         16 $upflank =~ s/[ \n]//g;
156 16         24 $dnflank =~ s/[ \n]//g;
157 16         38 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         46 my ($start, $sep, $end) = $location =~ /(-?\d+)(.)?\D?(-?\d+)?/;
164 16 100       31 $end = $start if not defined $end ;
165 16         24 my ($len) = $end - $start +1;
166 16 100 100     34 $len = 0, $start = $end if defined $sep and $sep eq '^';
167 16         12 my $ismut = 0;
168 16 100       31 $ismut = 1 if $change =~ m/>/;
169            
170 16         48 $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         14 $prevdnaobj = $dnamut;
180 16         33 my $a1 = Bio::Variation::Allele->new;
181 16 100       33 $a1->seq($ori) if $ori;
182 16         28 $dnamut->allele_ori($a1);
183 16         29 my $a2 = Bio::Variation::Allele->new;
184 16 100       32 $a2->seq($mut) if $mut;
185 16         30 $dnamut->add_Allele($a2);
186 16 100       21 if ($ismut) {
187 14         21 $dnamut->isMutation(1);
188 14         22 $dnamut->allele_mut($a2);
189             }
190 16 100       29 $dnamut->region($region) if defined $region;
191 16 100       24 $dnamut->region_value($region_value) if defined $region_value;
192 16 100       22 $dnamut->region_dist($region_dist) if defined $region_dist;
193              
194 16         33 $h->add_Variant($dnamut);
195 16         20 $dnamut->SeqDiff($h);
196             }
197             }
198              
199             #
200             # RNA
201             #
202 15         59 my @rna = split ( / RNA;/, $entry );
203 15         15 shift @rna;
204 15         11 my $prevrnaobj;
205 15         21 foreach my $rna (@rna) {
206 17         38 $rna = substr ($rna, 0, index($rna, 'Feature AA'));
207 17         139 $rna =~ s/Feature[ \t]+//g;
208 17         38 ($rna) = split "DNA; ", $rna;
209             #$self->warn("|$rna|") ;
210 17         180 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         75 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         22 $change =~ s/[ \n]//g;
219 17         28 my ($ori, $mut) = split /[>\|]/, $change;
220 17         16 my $rnamut;
221 17         19 my ($variation_number, $change_number) = split /\./, $mut_number;
222 17 100 100     43 if ($change_number and $change_number > 1 ) {
223 1         3 my $a3 = Bio::Variation::Allele->new;
224 1 50       4 $a3->seq($mut) if $mut;
225             #$rnamut->add_Allele($a3);
226 1         3 $prevrnaobj->add_Allele($a3);
227             } else {
228 16         40 my ($start, $sep, $end) = $location =~ /(-?\d+)(.)?\D?(-?\d+)?/;
229 16 100       28 $end = $start if not defined $end ;
230 16         25 my ($len) = $end - $start + 1;
231 16 100 100     36 $len = 0, $start = $end if defined $sep and $sep eq '^';
232 16         10 my $ismut;
233 16 100       28 $ismut = 1 if $change =~ m/>/;
234 16         39 my ($codon_table) = $rna =~ m|.+/codon_table: (\d+)|s;
235 16         32 my ($codon_pos) = $rna =~ m|.+/codon:[^;]+; ([123])|s;
236              
237 16         47 $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         13 $prevrnaobj = $rnamut;
248 16         29 my $a1 = Bio::Variation::Allele->new;
249 16 100       32 $a1->seq($ori) if $ori;
250 16         30 $rnamut->allele_ori($a1);
251 16         29 my $a2 = Bio::Variation::Allele->new;
252 16 100       31 $a2->seq($mut) if $mut;
253 16         29 $rnamut->add_Allele($a2);
254 16 100       22 if ($ismut) {
255 14         19 $rnamut->isMutation(1);
256 14         18 $rnamut->allele_mut($a2);
257             }
258 16 50       39 $rnamut->region($region) if defined $region;
259 16 50       19 $rnamut->region_value($region_value) if defined $region_value;
260 16 100       24 $rnamut->region_dist($region_dist) if defined $region_dist;
261              
262 16 100       33 $rnamut->codon_table($codon_table) if $codon_table;
263 16 100       27 $rnamut->codon_pos($codon_pos) if $codon_pos;
264 16         26 $h->add_Variant($rnamut);
265 16         26 foreach my $mut ($h->each_Variant) {
266 35 100       105 if ($mut->isa('Bio::Variation::DNAMutation') ) {
267 18 100       24 if ($mut->mut_number == $rnamut->mut_number) {
268 16         29 $rnamut->DNAMutation($mut);
269 16         24 $mut->RNAChange($rnamut);
270             }
271             }
272             }
273             }
274             }
275             #
276             # AA
277             #
278 15         80 my @aa = split ( / AA;/, $entry );
279 15         13 shift @aa;
280 15         17 my $prevaaobj;
281 15         13 foreach my $aa (@aa) {
282 13         22 $aa = substr ($aa, 0, index($aa, 'Feature AA'));
283 13         85 $aa =~ s/Feature[ \t]+//g;
284 13         27 ($aa) = split "DNA; ", $aa;
285             #$self->warn("|$aa|") ;
286 13         68 my ($mut_number, $proof, $location, $change) =
287             $aa =~ m|\W+([\d\.]+).+/proof: (\w+).+/location: ([^ \n]+)./change: ([^/;]+)|s;
288 13         25 $change =~ s/[ \n]//g;
289             #my $s = join ("|", $mut_number, $proof, $location, $change);
290             #$self->warn($s);
291             #exit;
292 13         14 $change =~ s/[ \n]//g;
293 13         9 $change =~ s/DNA$//;
294 13         25 my ($ori, $mut) = split /[>\|]/, $change;
295             #print "------$location----$ori-$mut-------------\n";
296 13         14 my ($variation_number, $change_number) = split /\./, $mut_number;
297 13         10 my $aamut;
298 13 100 100     27 if ($change_number and $change_number > 1 ) {
299 1         4 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         35 my ($start, $sep, $end) = $location =~ /(-?\d+)(.)?\D?(-?\d+)?/;
304 12 100       22 $end = $start if not defined $end ;
305 12         19 my ($len) = $end - $start + 1;
306 12 50 66     21 $len = 0, $start = $end if defined $sep and $sep eq '^';
307 12         7 my $ismut;
308 12 100       19 $ismut = 1 if $change =~ m/>/;
309 12         14 my ($region) = $aa =~ m|.+/region: (\w+)|s ;
310 12         30 $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         9 $prevaaobj = $aamut;
318 12         21 my $a1 = Bio::Variation::Allele->new;
319 12 50       20 $a1->seq($ori) if $ori;
320 12         21 $aamut->allele_ori($a1);
321 12         21 my $a2 = Bio::Variation::Allele->new;
322 12 100       25 $a2->seq($mut) if $mut;
323 12         22 $aamut->add_Allele($a2);
324 12 100       17 if ($ismut) {
325 10         15 $aamut->isMutation(1);
326 10         11 $aamut->allele_mut($a2);
327             }
328 12 50       17 $region && $aamut->region($region);
329 12         22 $h->add_Variant($aamut);
330 12         20 foreach my $mut ($h->each_Variant) {
331 41 100       122 if ($mut->isa('Bio::Variation::RNAChange') ) {
332 14 100       21 if ($mut->mut_number == $aamut->mut_number) {
333 12         26 $aamut->RNAChange($mut);
334 12         19 $mut->AAChange($aamut);
335             }
336             }
337             }
338              
339             }
340             }
341 15         60 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 65 my ($self,@h) = @_;
357              
358             #$columns = 75; #default for Text::Wrap
359 16         49 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       31 if( !defined $h[0] ) {
372 0         0 $self->throw("Attempting to write with no information!");
373             }
374              
375 16         16 foreach my $h (@h) {
376            
377 16         16 my @entry =();
378            
379 16         12 my ($text, $tmp, $tmp2, $sep);
380 16         16 my ($count) = 0;
381              
382            
383 16         16 $text = $tag{ID};
384            
385 16         31 $text .= $h->id;
386 16         32 $text .= ":(". $h->offset;
387 16 100       26 $text .= "+1" if $h->sysname =~ /-/;
388 16         25 $text .= ")". $h->sysname;
389 16 100       26 $text .= "; ". $h->trivname if $h->trivname;
390 16         19 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         22 my @allvariants = $h->each_Variant;
397 16         19 my %variants = ();
398 16         24 foreach my $mut ($h->each_Variant) {
399 47         33 push @{$variants{$mut->mut_number} }, $mut;
  47         62  
400             }
401             #my ($variation_number, $change_number) = split /\./, $mut_number;
402 16         36 foreach my $var (sort keys %variants) {
403             #print $var, ": ", join (" ", @{$variants{$var}}), "\n";
404            
405 17         15 foreach my $mut (@{$variants{$var}}) {
  17         23  
406             #
407             # DNA
408             #
409 47 100       214 if ( $mut->isa('Bio::Variation::DNAMutation') ) {
    100          
    50          
410             #collect all non-reference alleles
411 17 50       24 $self->throw("allele_ori needs to be defined in [$mut]")
412             if not $mut->allele_ori;
413 17 100       29 if ($mut->isMutation) {
414 15         15 $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         16 my $count = 0; # two alleles
421 17         17 foreach my $allele (@alleles) {
422 18         14 $count++;
423 18         27 my ($variation_number, $change_number) = split /\./, $mut->mut_number;
424 18 100 100     38 if ($change_number and $change_number != $count){
425 1         3 $mut->mut_number("$change_number.$count");
426             }
427 18         30 $mut->allele_mut($allele);
428             push (@entry,
429 18         44 $tag{FeatureKey}. 'DNA'. "; ". $mut->mut_number
430             );
431             #label
432 18         41 $text=$tag{FeatureQual}. '/label: '. $mut->label;
433 18         19 push (@entry, $text);
434            
435             #proof
436 18 50       34 if ($mut->proof) {
437 18         27 $text = $tag{FeatureQual}. '/proof: '. $mut->proof;
438 18         15 push (@entry, $text) ;
439             }
440             #location
441 18         27 $text = $tag{FeatureQual}. '/location: ';
442             #$mut->id. '; '. $mut->start;
443 18 100       27 if ($mut->length > 1 ) {# if ($mut->end - $mut->start ) {
    100          
444 2         3 my $l = $mut->start + $mut->length -1;
445 2         5 $text .= $mut->start. '..'. $l;
446             }
447             elsif ($mut->length == 0) {
448 2         6 my $tmp_start = $mut->start - 1;
449 2 50       5 $tmp_start-- if $tmp_start == 0;
450 2         5 $text .= $tmp_start. '^'. $mut->end;
451             } else {
452 14         30 $text .= $mut->start;
453             }
454              
455 18 100 66     34 if ($h->alphabet && $h->alphabet eq 'dna') {
456 4         7 $tmp = $mut->start + $h->offset;
457 4 50       8 $tmp-- if $tmp <= 0;
458 4 100       7 $mut->start < 1 && $tmp++;
459             #$text.= ' ('. $h->id. '::'. $tmp;
460 4         6 $tmp2 = $mut->end + $h->offset;
461 4 50       8 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         4 $text .= ')';
473             }
474 18         19 push (@entry, $text);
475             #sequence
476             push (@entry,
477 18         41 $tag{FeatureQual}. '/upflank: '. $mut->upStreamSeq
478             );
479 18         15 $text = '';
480 18 100       29 $text = $mut->allele_ori->seq if $mut->allele_ori->seq;
481 18         20 $text .= $sep;
482 18 100       26 $text .= $mut->allele_mut->seq if $mut->allele_mut->seq;
483             push (@entry,
484             wrap($tag{FeatureQual}. '/change: ', $tag{FeatureWrap},
485 18         53 $text)
486             );
487            
488             push (@entry,
489 18         2204 $tag{FeatureQual}. '/dnflank: '. $mut->dnStreamSeq
490             );
491             #restriction enzyme
492 18 100       43 if ($mut->restriction_changes ne '') {
493 17         39 $text = $mut->restriction_changes;
494 17         56 $text = wrap($tag{FeatureQual}. '/re_site: ', $tag{FeatureWrap}, $text);
495 17         2816 push (@entry,
496             $text
497             );
498             }
499             #region
500 18 100       51 if ($mut->region ) {
501 4         10 $text = $tag{FeatureQual}. '/region: '. $mut->region;
502 4 50 66     5 $text .= ';' if $mut->region_value or $mut->region_dist;
503 4 100       5 $text .= ' '. $mut->region_value if $mut->region_value;
504 4 50       7 if ($mut->region_dist ) {
505 4         3 $tmp = '';
506 4 100       6 $tmp = '+' if $mut->region_dist > 1;
507 4         8 $text .= " (". $tmp. $mut->region_dist. ')';
508             }
509 4         5 push (@entry, $text);
510             }
511             #CpG
512 18 50       34 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       30 $self->throw("allele_ori needs to be defined in [$mut]")
525             if not $mut->allele_ori;
526 17         41 my @alleles = $mut->each_Allele;
527             #push @alleles, $mut->allele_mut if $mut->allele_mut;
528 17 100       33 if ($mut->isMutation) {
529 15         16 $sep = '>';
530             } else {
531 2         2 $sep = '|';
532             }
533              
534 17         15 my $count = 0; # two alleles
535 17         16 foreach my $allele (@alleles) {
536 18         18 $count++;
537 18         31 my ($variation_number, $change_number) = split /\./, $mut->mut_number;
538 18 100 100     40 if ($change_number and $change_number != $count){
539 1         4 $mut->mut_number("$change_number.$count");
540             }
541 18         27 $mut->allele_mut($allele);
542             push (@entry,
543 18         43 $tag{FeatureKey}. 'RNA'. "; ". $mut->mut_number
544             );
545             #label
546 18         45 $text=$tag{FeatureQual}. '/label: '. $mut->label;
547 18         23 push (@entry, $text);
548             #proof
549 18 50       34 if ($mut->proof) {
550 18         41 $text = $tag{FeatureQual}. '/proof: '. $mut->proof;
551 18         19 push (@entry, $text) ;
552             }
553             #location
554 18         25 $text = $tag{FeatureQual}. '/location: ' ;
555 18 100       31 if ($mut->length > 1 ) {
    100          
556 2         5 $text .= $mut->start. '..'. $mut->end;
557 2         4 $tmp2 = $mut->end + $h->offset;
558             }
559             elsif ($mut->length == 0) {
560 2         5 my $tmp_start = $mut->start;
561 2         2 $tmp_start--;
562 2 50       4 $tmp_start-- if $tmp_start == 0;
563 2         5 $text .= $tmp_start. '^'. $mut->end;
564             } else {
565 14         26 $text .= $mut->start;
566             }
567              
568 18 100 66     38 if ($h->alphabet && $h->alphabet eq 'rna') {
569 14         23 $tmp = $mut->start + $h->offset;
570 14 50       23 $tmp-- if $tmp <= 0;
571             #$mut->start < 1 && $tmp++;
572             #$text.= ' ('. $h->id. '::'. $tmp;
573 14         30 $tmp2 = $mut->end + $h->offset;
574             #$mut->end < 1 && $tmp2++;
575 14 100       27 if ( $mut->length > 1 ) {
    100          
576 2         5 $text.= ' ('. $h->id. '::'. $tmp. "..". $tmp2;
577             }
578             elsif ($mut->length == 0) {
579 2         2 $tmp--;
580 2         4 $text .= ' ('. $h->id. '::'. $tmp. '^'. $tmp2;
581             } else {
582 10         16 $text.= ' ('. $h->id. '::'. $tmp;
583             }
584              
585 14         14 $text .= ')';
586             }
587 18         21 push (@entry, $text);
588              
589             #sequence
590             push (@entry,
591 18         35 $tag{FeatureQual}. '/upflank: '. $mut->upStreamSeq
592             );
593 18         16 $text = '';
594 18 100       25 $text = $mut->allele_ori->seq if $mut->allele_ori->seq;
595 18         20 $text .= $sep;
596 18 100       26 $text .= $mut->allele_mut->seq if $mut->allele_mut->seq;
597             push (@entry,
598             wrap($tag{FeatureQual}. '/change: ', $tag{FeatureWrap},
599 18         54 $text)
600             );
601             push (@entry,
602 18         2275 $tag{FeatureQual}. '/dnflank: '. $mut->dnStreamSeq
603             );
604             #restriction
605 18 100       41 if ($mut->restriction_changes ne '') {
606 17         26 $text = $mut->restriction_changes;
607 17         58 $text = wrap($tag{FeatureQual}. '/re_site: ', $tag{FeatureWrap}, $text);
608 17         2643 push (@entry,
609             $text
610             );
611             }
612             #coding
613 18 100       44 if ($mut->region eq 'coding') {
614             #codon table
615 14         19 $text = $tag{FeatureQual}. '/codon_table: ';
616 14         24 $text .= $mut->codon_table;
617 14         16 push (@entry, $text);
618             #codon
619              
620 14         30 $text = $tag{FeatureQual}. '/codon: '. $mut->codon_ori. $sep;
621 14 100       22 if ($mut->DNAMutation->label =~ /.*point/) {
622 9         17 $text .= $mut->codon_mut;
623             }
624             else {
625 5         6 $text .= '-';
626             }
627 14         23 $text .= "; ". $mut->codon_pos;
628 14         19 push (@entry, $text);
629             }
630             #region
631 18 50       27 if ($mut->region ) {
632 18         34 $text = $tag{FeatureQual}. '/region: '. $mut->region;
633 18 100 66     40 $text .= ';' if $mut->region_value or $mut->region_dist;
634 18 50       28 $text .= ' '. $mut->region_value if $mut->region_value;
635 18 100       22 if ($mut->region_dist ) {
636 2         3 $tmp = '';
637 2 100       3 $tmp = '+' if $mut->region_dist > 1;
638 2         3 $text .= " (". $tmp. $mut->region_dist. ')';
639             }
640 18         48 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       21 $self->throw("allele_ori needs to be defined in [$mut]")
650             if not $mut->allele_ori;
651 13 100       30 if ($mut->isMutation) {
652 11         11 $sep = '>';
653             } else {
654 2         3 $sep = '|';
655             }
656 13         20 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         15 foreach my $allele (@alleles) {
660 14         12 $count++;
661 14         23 my ($variation_number, $change_number) = split /\./, $mut->mut_number;
662 14 100 100     36 if ($change_number and $change_number != $count){
663 1         3 $mut->mut_number("$change_number.$count");
664             }
665 14         23 $mut->allele_mut($allele);
666             push (@entry,
667 14         34 $tag{FeatureKey}. 'AA'. "; ". $mut->mut_number
668             );
669             #label
670 14         36 $text=$tag{FeatureQual}. '/label: '. $mut->label;
671 14         16 push (@entry, $text) ;
672             #proof
673 14 50       23 if ($mut->proof) {
674 14         23 $text = $tag{FeatureQual}. '/proof: '. $mut->proof;
675 14         19 push (@entry, $text) ;
676             }
677             #location
678 14         28 $text = $tag{FeatureQual}. '/location: '.
679             #$mut->id. '; '. $mut->start;
680             $mut->start;
681 14 100       27 if ($mut->length > 1 ) {
682 1         3 $tmp = $mut->start + $mut->length -1;
683 1         3 $text .= '..'. $tmp;
684             }
685 14         13 push (@entry, $text);
686             #sequence
687 14         13 $text = '';
688 14 50       19 $text = $mut->allele_ori->seq if $mut->allele_ori->seq;
689 14         17 $text .= $sep;
690 14 100       22 $text .= $mut->allele_mut->seq if $mut->allele_mut->seq;
691             push (@entry,
692             wrap($tag{FeatureQual}. '/change: ', $tag{FeatureWrap},
693 14         42 $text)
694             );
695             #region
696 14 50       1974 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         19 push (@entry,
712             "//"
713             );
714 16         103 my $str = join ("\n", @entry). "\n";
715 16         56 $str =~ s/\t/ /g;
716 16         48 $self->_print($str);
717             }
718 16         48 return 1;
719             }
720              
721             1;