lib/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
  95 use strict;
  96
  97 use Bio::SeqFeature::FeaturePair;
  98 use Bio::SeqFeature::Generic;
  99
 100 use base qw(Bio::Root::Root Bio::Root::IO);
 101
 102 =head2 new
 103
 104  Title   : new
 105  Usage   : my $obj = Bio::Tools::Promoterwise->new();
 106  Function: Builds a new Bio::Tools::Promoterwise object
 107  Returns : L<Bio::Tools::Promoterwise>
 108  Args    : -fh/-file => $val, # for initing input, see Bio::Root::IO
 109
 110
 111 =cut
 112
 113 sub new {
 114   my($class,@args) = @_;
 115
 116   my $self = $class->SUPER::new(@args);
 117   $self->_initialize_io(@args);
 118   my ($query1,$query2) = $self->_rearrange([qw(QUERY1_SEQ QUERY2_SEQ)],@args);
 119   $self->query1_seq($query1) if ($query1);
 120   $self->query2_seq($query2) if ($query2);
 121
 122   return $self;
 123 }
 124
 125 =head2 next_result
 126
 127  Title   : next_result
 128  Usage   : my $r = $rpt_masker->next_result
 129  Function: Get the next result set from parser data
 130  Returns : an  L<Bio::SeqFeature::FeaturePair>
 131  Args    : none
 132
 133
 134 =cut
 135
 136 sub next_result {
 137   my ($self) = @_;
 138   $self->_parse unless $self->_parsed;
 139   return $self->_next_result;
 140 }
 141
 142 sub _parse{
 143    my ($self) = @_;
 144    my (%hash,@fp);
 145    while (defined($_ = $self->_readline()) ) {
 146        chomp;
 147        my @array = split;
 148        push @{$hash{$array[-1]}}, \@array;
 149    }
 150    foreach my $key(keys %hash){
 151     my $sf1 = Bio::SeqFeature::Generic->new(-primary=>"conserved_element",
 152                                             -source_tag=>"promoterwise");
 153     $sf1->attach_seq($self->query1_seq) if $self->query1_seq;
 154     my $sf2 = Bio::SeqFeature::Generic->new(-primary=>"conserved_element",
 155                                             -source_tag=>"promoterwise");
 156     $sf2->attach_seq($self->query2_seq) if $self->query2_seq;
 157     foreach my $info(@{$hash{$key}}){
 158
 159       my ($score,$id1,$start_1,$end_1, $strand_1,$s1_len,
 160           $id2,$start_2,$end_2,$strand_2,$s2_len, $group);
 161       if( @{$info} == 12 ) {
 162           ($score,$id1,$start_1,$end_1, $strand_1,$s1_len,
 163            $id2,$start_2,$end_2,$strand_2,$s2_len, $group) = @{$info};
 164       } elsif( @{$info} == 10 ) {
 165           ($score,$id1,$start_1,$end_1, $strand_1,
 166            $id2,$start_2,$end_2,$s2_len, $group) = @{$info};
 167       } else {
 168           $self->throw("unknown promoterwise output, ", scalar @{$info},
 169                        " columns, expected 10 or 12\n");
 170       }
 171       if(!$sf1->strand && !$sf2->strand){
 172         $sf1->strand($strand_1);
 173         $sf2->strand($strand_2);
 174         $sf1->seq_id($id1);
 175         $sf2->seq_id($id2);
 176         $sf1->score($score);
 177         $sf2->score($score);
 178       }
 179
 180       my $sub1 = Bio::SeqFeature::Generic->new(-start=>$start_1,
 181                                               -seq_id=>$id1,
 182                                               -end  =>$end_1,
 183                                               -strand=>$strand_1,
 184                                               -primary=>"conserved_element",
 185                                               -source_tag=>"promoterwise",
 186                                               -score=>$score);
 187       $sub1->attach_seq($self->query1_seq) if $self->query1_seq;
 188
 189       my $sub2 = Bio::SeqFeature::Generic->new(-start=>$start_2,
 190                                               -seq_id=>$id2,
 191                                               -end  =>$end_2,
 192                                               -strand=>$strand_2,
 193                                               -primary=>"conserved_element",
 194                                               -source_tag=>"promoterwise",
 195                                               -score=>$score);
 196       $sub2->attach_seq($self->query2_seq) if $self->query2_seq;
 197       $sf1->add_SeqFeature($sub1,'EXPAND');
 198       $sf2->add_SeqFeature($sub2,'EXPAND');
 199     }
 200
 201     my $fp = Bio::SeqFeature::FeaturePair->new(-feature1=>$sf1,
 202                                                -feature2=>$sf2);
 203     push @fp, $fp;
 204   }
 205     $self->_feature_pairs(\@fp);
 206     $self->_parsed(1);
 207     return;
 208 }
 209
 210 sub _feature_pairs {
 211   my ($self,$fp) = @_;
 212   if($fp){
 213     $self->{'_feature_pairs'} = $fp;
 214   }
 215   return  $self->{'_feature_pairs'};
 216 }
 217
 218 sub _next_result {
 219   my ($self) = @_;
 220   return unless (exists($self->{'_feature_pairs'}) && @{$self->{'_feature_pairs'}});
 221   return shift(@{$self->{'_feature_pairs'}});
 222 }
 223 sub _parsed {
 224   my ($self,$flag) = @_;
 225   if($flag){
 226     $self->{'_flag'} = 1;
 227   }
 228   return $self->{'_flag'};
 229 }
 230
 231 sub query1_seq {
 232   my ($self,$val) = @_;
 233   if($val){
 234     $self->{'query1_seq'} = $val;
 235   }
 236   return $self->{'query1_seq'};
 237 }
 238 sub query2_seq {
 239   my ($self,$val) = @_;
 240   if($val){
 241     $self->{'query2_seq'} = $val;
 242   }
 243   return $self->{'query2_seq'};
 244 }
 245 1;