Bio/Tools/Promoterwise.pm

   1 # BioPerl module for Bio::Tools::Promoterwise
   2 #
   3 # Please direct questions and support issues to <bioperl-l@bioperl.org>
   4 #
   5 # Cared for by Shawn Hoon <shawnh@fugu-sg.org>
   6 #
   7 # Copyright Shawn Hoon
   8 #
   9 # You may distribute this module under the same terms as perl itself
  10
  11 # POD documentation - main docs before the code
  12
  13 =head1 NAME
  14
  15 Bio::Tools::Promoterwise - parser for Promoterwise tab format output
  16
  17 =head1 SYNOPSIS
  18
  19   use Bio::Tools::Promoterwise;
  20
  21   my $pw = Bio::Tools::Promoterwise->new(-file=>"out",
  22                                          -query1_seq=>$seq1,
  23                                          -query2_seq=>$seq2);
  24   while (my $fp = $pw->next_result){
  25     print "Hit Length: ".$fp->feature1->length."\n";
  26     print "Hit Start: ".$fp->feature1->start."\n";
  27     print "Hit End: ".$fp->feature1->end."\n";
  28     print "Hsps: \n";
  29     my @first_hsp = $fp->feature1->sub_SeqFeature;
  30     my @second_hsp = $fp->feature2->sub_SeqFeature;
  31     foreach my $i (0..$#first_hsp){
  32       print $first_hsp[$i]->start. " ".$first_hsp[$i]->end." ".
  33             $second_hsp[$i]->start. " ".$second_hsp[$i]->end."\n";
  34     }
  35   }
  36
  37 =head1 DESCRIPTION
  38
  39 Promoteriwise is an alignment algorithm that relaxes the constraint
  40 that local alignments have to be co-linear. Otherwise it provides a
  41 similar model to DBA, which is designed for promoter sequence
  42 alignments.  Promoterwise is written by Ewan Birney.  It is part of
  43 the wise2 package available at
  44 L<ftp://ftp.ebi.ac.uk/pub/software/unix/wise2/>
  45
  46 This module is the parser for the Promoterwise output in tab format.
  47
  48 =head1 FEEDBACK
  49
  50 =head2 Mailing Lists
  51
  52 User feedback is an integral part of the evolution of this and other
  53 Bioperl modules. Send your comments and suggestions preferably to
  54 the Bioperl mailing list.  Your participation is much appreciated.
  55
  56   bioperl-l@bioperl.org                  - General discussion
  57   http://bioperl.org/wiki/Mailing_lists  - About the mailing lists
  58
  59 =head2 Support
  60
  61 Please direct usage questions or support issues to the mailing list:
  62
  63 I<bioperl-l@bioperl.org>
  64
  65 rather than to the module maintainer directly. Many experienced and
  66 reponsive experts will be able look at the problem and quickly
  67 address it. Please include a thorough description of the problem
  68 with code and data examples if at all possible.
  69
  70 =head2 Reporting Bugs
  71
  72 Report bugs to the Bioperl bug tracking system to help us keep track
  73 of the bugs and their resolution. Bug reports can be submitted via the
  74 web:
  75
  76   https://github.com/bioperl/bioperl-live/issues
  77
  78 =head1 AUTHOR - Shawn Hoon
  79
  80 Email shawnh@fugu-sg.org
  81
  82 =head1 APPENDIX
  83
  84 The rest of the documentation details each of the object methods.
  85 Internal methods are usually preceded with a _
  86
  87 =cut
  88
  89
  90 # Let the code begin...
  91
  92
  93 package Bio::Tools::Promoterwise;
  94 use strict;
  95
  96 use Bio::SeqFeature::FeaturePair;
  97 use Bio::SeqFeature::Generic;
  98
  99 use base qw(Bio::Root::Root Bio::Root::IO);
 100
 101 =head2 new
 102
 103  Title   : new
 104  Usage   : my $obj = Bio::Tools::Promoterwise->new();
 105  Function: Builds a new Bio::Tools::Promoterwise object
 106  Returns : L<Bio::Tools::Promoterwise>
 107  Args    : -fh/-file => $val, # for initing input, see Bio::Root::IO
 108
 109
 110 =cut
 111
 112 sub new {
 113   my($class,@args) = @_;
 114
 115   my $self = $class->SUPER::new(@args);
 116   $self->_initialize_io(@args);
 117   my ($query1,$query2) = $self->_rearrange([qw(QUERY1_SEQ QUERY2_SEQ)],@args);
 118   $self->query1_seq($query1) if ($query1);
 119   $self->query2_seq($query2) if ($query2);
 120
 121   return $self;
 122 }
 123
 124 =head2 next_result
 125
 126  Title   : next_result
 127  Usage   : my $r = $rpt_masker->next_result
 128  Function: Get the next result set from parser data
 129  Returns : an  L<Bio::SeqFeature::FeaturePair>
 130  Args    : none
 131
 132
 133 =cut
 134
 135 sub next_result {
 136   my ($self) = @_;
 137   $self->_parse unless $self->_parsed;
 138   return $self->_next_result;
 139 }
 140
 141 sub _parse{
 142    my ($self) = @_;
 143    my (%hash,@fp);
 144    while (defined($_ = $self->_readline()) ) {
 145        chomp;
 146        my @array = split;
 147        push @{$hash{$array[-1]}}, \@array;
 148    }
 149    foreach my $key(keys %hash){
 150     my $sf1 = Bio::SeqFeature::Generic->new(-primary=>"conserved_element",
 151                                             -source_tag=>"promoterwise");
 152     $sf1->attach_seq($self->query1_seq) if $self->query1_seq;
 153     my $sf2 = Bio::SeqFeature::Generic->new(-primary=>"conserved_element",
 154                                             -source_tag=>"promoterwise");
 155     $sf2->attach_seq($self->query2_seq) if $self->query2_seq;
 156     foreach my $info(@{$hash{$key}}){
 157
 158       my ($score,$id1,$start_1,$end_1, $strand_1,$s1_len,
 159           $id2,$start_2,$end_2,$strand_2,$s2_len, $group);
 160       if( @{$info} == 12 ) {
 161           ($score,$id1,$start_1,$end_1, $strand_1,$s1_len,
 162            $id2,$start_2,$end_2,$strand_2,$s2_len, $group) = @{$info};
 163       } elsif( @{$info} == 10 ) {
 164           ($score,$id1,$start_1,$end_1, $strand_1,
 165            $id2,$start_2,$end_2,$s2_len, $group) = @{$info};
 166       } else {
 167           $self->throw("unknown promoterwise output, ", scalar @{$info},
 168                        " columns, expected 10 or 12\n");
 169       }
 170       if(!$sf1->strand && !$sf2->strand){
 171         $sf1->strand($strand_1);
 172         $sf2->strand($strand_2);
 173         $sf1->seq_id($id1);
 174         $sf2->seq_id($id2);
 175         $sf1->score($score);
 176         $sf2->score($score);
 177       }
 178
 179       my $sub1 = Bio::SeqFeature::Generic->new(-start=>$start_1,
 180                                               -seq_id=>$id1,
 181                                               -end  =>$end_1,
 182                                               -strand=>$strand_1,
 183                                               -primary=>"conserved_element",
 184                                               -source_tag=>"promoterwise",
 185                                               -score=>$score);
 186       $sub1->attach_seq($self->query1_seq) if $self->query1_seq;
 187
 188       my $sub2 = Bio::SeqFeature::Generic->new(-start=>$start_2,
 189                                               -seq_id=>$id2,
 190                                               -end  =>$end_2,
 191                                               -strand=>$strand_2,
 192                                               -primary=>"conserved_element",
 193                                               -source_tag=>"promoterwise",
 194                                               -score=>$score);
 195       $sub2->attach_seq($self->query2_seq) if $self->query2_seq;
 196       $sf1->add_SeqFeature($sub1,'EXPAND');
 197       $sf2->add_SeqFeature($sub2,'EXPAND');
 198     }
 199
 200     my $fp = Bio::SeqFeature::FeaturePair->new(-feature1=>$sf1,
 201                                                -feature2=>$sf2);
 202     push @fp, $fp;
 203   }
 204     $self->_feature_pairs(\@fp);
 205     $self->_parsed(1);
 206     return;
 207 }
 208
 209 sub _feature_pairs {
 210   my ($self,$fp) = @_;
 211   if($fp){
 212     $self->{'_feature_pairs'} = $fp;
 213   }
 214   return  $self->{'_feature_pairs'};
 215 }
 216
 217 sub _next_result {
 218   my ($self) = @_;
 219   return unless (exists($self->{'_feature_pairs'}) && @{$self->{'_feature_pairs'}});
 220   return shift(@{$self->{'_feature_pairs'}});
 221 }
 222 sub _parsed {
 223   my ($self,$flag) = @_;
 224   if($flag){
 225     $self->{'_flag'} = 1;
 226   }
 227   return $self->{'_flag'};
 228 }
 229
 230 sub query1_seq {
 231   my ($self,$val) = @_;
 232   if($val){
 233     $self->{'query1_seq'} = $val;
 234   }
 235   return $self->{'query1_seq'};
 236 }
 237 sub query2_seq {
 238   my ($self,$val) = @_;
 239   if($val){
 240     $self->{'query2_seq'} = $val;
 241   }
 242   return $self->{'query2_seq'};
 243 }
 244 1;