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