Bio/Tools/Blat.pm

   1 # $Id$
   2 #
   3 # BioPerl module for Bio::Tools::Blat
   4 #
   5 # Written by Balamurugan Kumarasamy
   6 #
   7 # You may distribute this module under the same terms as perl itself
   8 #
   9
  10 =head1 NAME
  11
  12 Bio::Tools::Blat - parser for Blat program
  13
  14 =head1 SYNOPSIS
  15
  16   use Bio::Tools::Blat;
  17   my $blat_parser = Bio::Tools::Blat->new(-fh =>$filehandle );
  18   while( my $blat_feat = $blat_parser->next_result ) {
  19         push @blat_feat, $blat_feat;
  20   }
  21
  22 =head1 DESCRIPTION
  23
  24  Parser for Blat program
  25
  26 =head1 FEEDBACK
  27
  28 =head2 Mailing Lists
  29
  30  User feedback is an integral part of the evolution of this and other
  31  Bioperl modules. Send your comments and suggestions preferably to
  32  the Bioperl mailing list.  Your participation is much appreciated.
  33
  34   bioperl-l@bioperl.org                  - General discussion
  35   http://bioperl.org/wiki/Mailing_lists  - About the mailing lists
  36
  37 =head2 Reporting Bugs
  38
  39  Report bugs to the Bioperl bug tracking system to help us keep track
  40  of the bugs and their resolution. Bug reports can be submitted the
  41  web:
  42
  43   http://bugzilla.open-bio.org/
  44
  45 =head1 AUTHOR - Balamurugan Kumarasamy
  46
  47  Email: bala@tll.org.sg
  48
  49 =head1 APPENDIX
  50
  51  The rest of the documentation details each of the object methods.
  52  Internal methods are usually preceded with a _
  53
  54 =cut
  55
  56 package Bio::Tools::Blat;
  57 use strict;
  58 use Bio::SeqFeature::Generic;
  59 use Bio::SeqFeature::FeaturePair;
  60 use Bio::SeqFeature::Generic;
  61 use base qw(Bio::Root::Root Bio::Root::IO);
  62
  63 =head2 new
  64
  65  Title   : new
  66  Usage   : my $obj = Bio::Tools::Blat->new(-fh=>$filehandle);
  67  Function: Builds a new Bio::Tools::Blat object
  68  Returns : Bio::Tools::Blat
  69  Args    : -filename
  70            -fh (filehandle)
  71
  72 =cut
  73
  74 sub new {
  75       my($class,@args) = @_;
  76
  77       my $self = $class->SUPER::new(@args);
  78       $self->_initialize_io(@args);
  79
  80       return $self;
  81 }
  82
  83 =head2 next_result
  84
  85  Title   : next_result
  86  Usage   : my $feat = $blat_parser->next_result
  87  Function: Get the next result set from parser data
  88  Returns : L<Bio::SeqFeature::Generic>
  89  Args    : none
  90
  91 =cut
  92
  93 sub next_result {
  94         my ($self) = @_;
  95         my $filehandle;
  96         my $line;
  97         my $id;
  98
  99         while ($_=$self->_readline()){
 100                 # first split on spaces:
 101                 $line = $_;
 102                 chomp $line;
 103
 104                 my ($matches, $mismatches, $rep_matches, $n_count, $q_num_insert,
 105                          $q_base_insert, $t_num_insert, $t_base_insert, $strand, $q_name,
 106                          $q_length, $q_start, $q_end, $t_name, $t_length,
 107                          $t_start, $t_end, $block_count, $block_sizes, $q_starts,
 108                          $t_starts
 109                         ) = split;
 110
 111                 my $superfeature = Bio::SeqFeature::Generic->new();
 112
 113                 # ignore any preceeding text
 114                 next unless ( $matches =~/^\d+$/ );
 115
 116                 # create as many features as blocks there are in each output line
 117                 my (%feat1, %feat2);
 118                 $feat1{name} = $t_name;
 119                 $feat2{name} = $q_name;
 120
 121                 $strand = $1 if ($strand =~/([+-])[+-]/);
 122
 123                 $feat2{strand} = 1;
 124                 $feat1{strand} = $strand;
 125
 126                 my $percent_id = sprintf "%.2f",
 127                 (100 * ($matches + $rep_matches)/( $matches + $mismatches + $rep_matches));
 128
 129                 unless ( $q_length ){
 130                         $self->warn("length of query is zero, something is wrong!");
 131                         next;
 132                 }
 133
 134                 my $score   = sprintf "%.2f",
 135                 (100 * ( $matches + $mismatches + $rep_matches ) / $q_length);
 136
 137                 # size of each block of alignment (inclusive)
 138                 my @block_sizes     = split ",",$block_sizes;
 139
 140                 # start position of each block (you must add 1 as psl output
 141                 # is off by one in the start coordinate)
 142                 my @q_start_positions = split ",",$q_starts;
 143                 my @t_start_positions = split ",",$t_starts;
 144
 145                 $superfeature->seq_id($q_name);
 146                 $superfeature->score( $score );
 147                 $superfeature->add_tag_value('percent_id',$percent_id);
 148
 149                 # each line of output represents one possible entire aligment
 150                 # of the query (feat1) and the target(feat2)
 151
 152                 for (my $i=0; $i<$block_count; $i++ ){
 153
 154                         my ($query_start,$query_end);
 155
 156                         if ( $strand eq '+' ){
 157                                 $query_start = $q_start_positions[$i] + 1;
 158                                 $query_end   = $query_start + $block_sizes[$i] - 1;
 159                         }else{
 160                                 $query_end   = $q_length  - $q_start_positions[$i];
 161                                 $query_start = $query_end - $block_sizes[$i] + 1;
 162                         }
 163
 164                         #$feat2 {start} = $q_start_positions[$i] + 1;
 165                         #$feat2 {end}   = $feat2{start} + $block_sizes[$i] - 1;
 166                         $feat2 {start} = $query_start;
 167                         $feat2 {end}   = $query_end;
 168                         if ( $query_end <  $query_start ){
 169                                 $self->warn("dodgy feature coordinates: end = $query_end, start = $query_start. Reversing...");
 170                                 $feat2 {end}   = $query_start;
 171                                 $feat2 {start} = $query_end;
 172                         }
 173
 174                         $feat1 {start} = $t_start_positions[$i] + 1;
 175                         $feat1 {end}   = $feat1{start} + $block_sizes[$i] - 1;
 176
 177                         # we put all the features with the same score and percent_id
 178                         $feat2 {score}   = $score;
 179                         $feat1 {score}   = $feat2 {score};
 180                         $feat2 {percent} = $percent_id;
 181                         $feat1 {percent} = $feat2 {percent};
 182
 183                         # other stuff:
 184                         $feat1 {db}         = undef;
 185                         $feat1 {db_version} = undef;
 186                         $feat1 {program}    = 'blat';
 187                         $feat1 {p_version}  = '1';
 188                         $feat1 {source}     = 'blat';
 189                         $feat1 {primary}    = 'similarity';
 190                         $feat2 {source}     = 'blat';
 191                         $feat2 {primary}    = 'similarity';
 192
 193                         my $feature_pair = $self->create_feature(\%feat1, \%feat2);
 194                         $superfeature->add_sub_SeqFeature( $feature_pair,'EXPAND');
 195                 }
 196                 return $superfeature;
 197         }
 198 }
 199
 200 =head2 create_feature
 201
 202  Title   : create_feature
 203  Usage   : my $feat=$blat_parser->create_feature($feature,$seqname)
 204  Function: creates a SeqFeature Generic object
 205  Returns : L<Bio::SeqFeature::Generic>
 206  Args    :
 207
 208
 209 =cut
 210
 211 sub create_feature {
 212     my ($self, $feat1,$feat2) = @_;
 213     my $feature1= Bio::SeqFeature::Generic->new(
 214                                                           -seq_id     =>$feat1->{name},
 215                                                           -start      =>$feat1->{start},
 216                        -end        =>$feat1->{end},
 217                        -strand     =>$feat1->{strand},
 218                        -score      =>$feat1->{score},
 219                        -source     =>$feat1->{source},
 220                        -primary    =>$feat1->{primary} );
 221
 222     my $feature2= Bio::SeqFeature::Generic->new(
 223                        -seq_id     =>$feat2->{name},
 224                                                           -start      =>$feat2->{start},
 225                        -end        =>$feat2->{end},
 226                        -strand     =>$feat2->{strand},
 227                        -score      =>$feat2->{score},
 228                        -source     =>$feat2->{source},
 229                        -primary    =>$feat2->{primary} );
 230
 231     my $featurepair = Bio::SeqFeature::FeaturePair->new;
 232     $featurepair->feature1 ($feature1);
 233     $featurepair->feature2 ($feature2);
 234
 235          $featurepair->add_tag_value('evalue',$feat2->{p});
 236          $featurepair->add_tag_value('percent_id',$feat2->{percent});
 237          $featurepair->add_tag_value("hid",$feat2->{primary});
 238     return  $featurepair;
 239 }
 240
 241 1;