code_perl/Maasha/Blast.pm

   1 package Maasha::Blast;
   2
   3 # Copyright (C) 2007-2009 Martin A. Hansen.
   4
   5 # This program is free software; you can redistribute it and/or
   6 # modify it under the terms of the GNU General Public License
   7 # as published by the Free Software Foundation; either version 2
   8 # of the License, or (at your option) any later version.
   9
  10 # This program is distributed in the hope that it will be useful,
  11 # but WITHOUT ANY WARRANTY; without even the implied warranty of
  12 # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  13 # GNU General Public License for more details.
  14
  15 # You should have received a copy of the GNU General Public License
  16 # along with this program; if not, write to the Free Software
  17 # Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
  18
  19 # http://www.gnu.org/copyleft/gpl.html
  20
  21
  22 # >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> DESCRIPTION <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
  23
  24
  25 # Routines to run BLAST and parse results.
  26
  27
  28 # >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
  29
  30
  31 use warnings;
  32 use strict;
  33 use Storable qw( dclone );
  34 use Data::Dumper;
  35 use vars qw( @ISA @EXPORT );
  36
  37 @ISA = qw( Exporter );
  38
  39
  40 # >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
  41
  42
  43 # The blast report structured like this:
  44 #
  45 #    for the first entry:
  46 #
  47 #    1  -  blast program name
  48 #    2  -  blast reference
  49 #    3  -  query sequence name and length
  50 #    4  -  subject database
  51 #    5  -  sequences producing significant alignments
  52 #    6  -  one or more HSP for each subject sequence
  53 #    7  -  blast statistics
  54 #
  55 #    for subsequent entries:
  56 #
  57 #    3  -  query sequence name and length
  58 #    5  -  sequences producing significant alignments
  59 #    6  -  one or more HSP for each subject sequence
  60 #    7  -  blast statistics
  61 #
  62 # ________________________
  63 #
  64 # info
  65 # query
  66 # parems
  67 # subject
  68 #     hit1
  69 #         hsp1
  70 #         hsp2
  71 #         hsp3
  72 #     hit2
  73 #         hsp1
  74 #     hit3
  75 #         hsp1
  76 #         hsp2
  77 # stats
  78 # ________________________
  79
  80
  81 # >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> SUBROUTINES <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
  82
  83
  84 sub xml_parse_blast
  85 {
  86     # Martin A. Hansen, March 2007.
  87
  88     # determines if the results is from ncbi blast or blastcl3
  89     # and parses the results in accordance.
  90
  91     my ( $fh,
  92        ) = @_;
  93
  94     # returns list
  95
  96     my ( $results, $line,$doctype );
  97
  98     while ( $line = <$fh> )
  99     {
 100         chomp $line;
 101
 102         if ( $line =~ /^<!DOCTYPE/ )
 103         {
 104             $doctype = $line;
 105             last;
 106         }
 107     }
 108
 109     if ( $doctype eq '<!DOCTYPE BlastOutput PUBLIC "-//NCBI//NCBI BlastOutput/EN" "NCBI_BlastOutput.dtd">' )
 110     {
 111         print STDERR qq(Parsing blastcl3 results ...\n);
 112         $results = xml_parse_blast_blastcl3( $fh );
 113     }
 114     elsif ( $doctype eq '<!DOCTYPE BlastOutput PUBLIC "-//NCBI//NCBI BlastOutput/EN" "http://www.ncbi.nlm.nih.gov/dtd/NCBI_BlastOutput.dtd">' )
 115     {
 116         print STDERR qq(Parsing NCBI blast results ...\n);
 117         $results = xml_parse_blast_ncbi( $fh );
 118     }
 119     else
 120     {
 121         die qq(ERROR: Could not determine doctype\n);
 122     }
 123
 124     return wantarray ? @{ $results } : $results;
 125 }
 126
 127
 128 sub xml_parse_blast_ncbi
 129 {
 130     # Martin A. Hansen, February 2007.
 131
 132     my ( $fh,
 133        ) = @_;
 134
 135     my ( $blast_record, $line, @blast_query, @blast_subject, $query, $subject, @results );
 136
 137     while ( $blast_record = xml_get_blast_record( $fh ) and scalar @{ $blast_record } > 0 )
 138     {
 139         foreach $line ( @{ $blast_record } )
 140         {
 141             if ( $line =~ /<Iteration_query-ID>|<Iteration_query-def>|<Iteration_query-len>/ )
 142             {
 143                 push @blast_query, $line;
 144             }
 145             elsif ( @blast_query )
 146             {
 147                 push @blast_subject, $line;
 148
 149                 if ( $line =~ /<\/Iteration_hits>/ )
 150                 {
 151                     $query   = xml_parse_blast_query( \@blast_query );
 152                     $subject = xml_parse_blast_subject( \@blast_subject );
 153
 154                     push @results, {
 155                         "QUERY"   => $query,
 156                         "SUBJECT" => $subject,
 157                     };
 158
 159                     undef @blast_query;
 160                     undef @blast_subject;
 161                 }
 162             }
 163         }
 164     }
 165
 166     return wantarray ? @results : \@results;
 167 }
 168
 169
 170 sub xml_parse_blast_blastcl3
 171 {
 172     # Martin A. Hansen, February 2007.
 173
 174     my ( $fh,
 175        ) = @_;
 176
 177     my ( $blast_record, $line, @blast_query, @blast_subject, $query, $subject, @results );
 178
 179     while ( $blast_record = xml_get_blast_record( $fh ) and scalar @{ $blast_record } > 0 )
 180     {
 181         foreach $line ( @{ $blast_record } )
 182         {
 183             if ( $line =~ /<BlastOutput_query-ID>|<BlastOutput_query-def>|<BlastOutput_query-len>/ )
 184             {
 185                 push @blast_query, $line;
 186             }
 187             elsif ( @blast_query )
 188             {
 189                 push @blast_subject, $line;
 190
 191                 if ( $line =~ /<\/Iteration_hits>/ )
 192                 {
 193                     $query   = xml_parse_blast_query( \@blast_query );
 194                     $subject = xml_parse_blast_subject( \@blast_subject );
 195
 196                     push @results, {
 197                         "QUERY"   => $query,
 198                         "SUBJECT" => $subject,
 199                     };
 200
 201                     undef @blast_query;
 202                     undef @blast_subject;
 203                 }
 204             }
 205         }
 206     }
 207
 208     return wantarray ? @results : \@results;
 209 }
 210
 211
 212 sub xml_get_blast_record
 213 {
 214     # Martin A. Hansen, March 2007.
 215
 216     my ( $fh,   # file handle to BLAST file in XML format
 217        ) = @_;
 218
 219     # returns list of lines
 220
 221     my ( $line, @blast_record );
 222
 223     while ( $line = <$fh> )
 224     {
 225         chomp $line;
 226
 227         push @blast_record, $line;
 228
 229         last if $line =~ /<\/BlastOutput>/;
 230     }
 231
 232     return wantarray ? @blast_record : \@blast_record;
 233 }
 234
 235
 236 sub xml_parse_blast_query
 237 {
 238     my ( $lines,
 239        ) = @_;
 240
 241     my ( $line, %hash );
 242
 243     foreach $line ( @{ $lines } )
 244     {
 245         if ( $line =~ /<Iteration_query-ID>([^<]+)/ ) {
 246             $hash{ "Q_ID" } = $1;
 247         } elsif ( $line =~ /<Iteration_query-def>([^<]+)/ ) {
 248             $hash{ "Q_DEF" } = $1;
 249         } elsif ( $line =~ /<Iteration_query-len>([^<]+)/ ) {
 250             $hash{ "Q_LEN" } = $1;
 251         }
 252     }
 253
 254     return wantarray ? %hash : \%hash;
 255 }
 256
 257
 258 sub xml_parse_blast_subject
 259 {
 260     # Martin A. Hansen, March 2007.
 261
 262     my ( $lines,   #
 263        ) = @_;
 264
 265     # returns
 266
 267     my ( $line, @blast_hit, @blast_hsps, $hit, $hsps, @hits );
 268
 269     foreach $line ( @{ $lines } )
 270     {
 271         if ( $line =~ /<Hit_id>|<Hit_def>|<Hit_accession>|<Hit_len>|<Hit_num>/ )
 272         {
 273             push @blast_hit, $line;
 274         }
 275         elsif ( @blast_hit )
 276         {
 277             push @blast_hsps, $line;
 278
 279             if ( $line =~ /<\/Hit_hsps>/ )
 280             {
 281                 $hit  = xml_parse_blast_hit( \@blast_hit );
 282                 $hsps = xml_parse_blast_hsps( \@blast_hsps );
 283
 284                 $hit->{ "HSPS" } = $hsps;
 285
 286                 push @hits, $hit;
 287
 288                 undef @blast_hit;
 289                 undef @blast_hsps;
 290             }
 291         }
 292     }
 293
 294     return wantarray ? @hits : \@hits;
 295 }
 296
 297
 298 sub xml_parse_blast_hit
 299 {
 300     my ( $lines
 301        ) = @_;
 302
 303     my ( $line, %hash );
 304
 305     foreach $line ( @{ $lines } )
 306     {
 307         if ( $line =~ /<Hit_num>([^<]+)/ ) {
 308             $hash{ "S_NUM" } = $1;
 309         } elsif ( $line =~ /<Hit_id>([^<]+)/ ) {
 310             $hash{ "S_ID" } = $1;
 311         } elsif ( $line =~ /<Hit_def>([^<]+)/ ) {
 312             $hash{ "S_DEF" } = $1;
 313         } elsif ( $line =~ /<Hit_accession>([^<]+)/ ) {
 314             $hash{ "S_ACC" } = $1;
 315         } elsif ( $line =~ /<Hit_len>([^<]+)/ ) {
 316             $hash{ "S_LEN" } = $1;
 317         }
 318     }
 319
 320     return wantarray ? %hash : \%hash;
 321 }
 322
 323
 324 sub xml_parse_blast_hsps
 325 {
 326     # Martin A. Hansen, March 2007.
 327
 328     my ( $blast_hits,   #
 329        ) = @_;
 330
 331     # returns
 332
 333     my ( $line, %hash, @hsps );
 334
 335     foreach $line ( @{ $blast_hits } )
 336     {
 337         if ( $line =~ /<Hsp_num>([^<]+)/ ) {
 338             $hash{ "NUM" } = $1;
 339         } elsif ( $line =~ /<Hsp_evalue>([^<]+)/ ) {
 340             $hash{ "E_VAL" } = $1;
 341         } elsif ( $line =~ /<Hsp_query-from>([^<]+)/ ) {
 342             $hash{ "Q_BEG" } = $1 - 1;
 343         } elsif ( $line =~ /<Hsp_query-to>([^<]+)/ ) {
 344             $hash{ "Q_END" } = $1 - 1;
 345         } elsif ( $line =~ /<Hsp_hit-from>([^<]+)/ ) {
 346             $hash{ "S_BEG" } = $1 - 1;
 347         } elsif ( $line =~ /<Hsp_hit-to>([^<]+)/ ) {
 348             $hash{ "S_END" } = $1 - 1;
 349         } elsif ( $line =~ /<Hsp_query-frame>([^<]+)/ ) {
 350             $hash{ "Q_FRAME" } = $1;
 351         } elsif ( $line =~ /<Hsp_hit-frame>([^<]+)/ ) {
 352             $hash{ "S_FRAME" } = $1;
 353         } elsif ( $line =~ /<Hsp_qseq>([^<]+)/ ) {
 354             $hash{ "Q_ALIGN" } = $1;
 355         } elsif ( $line =~ /<Hsp_hseq>([^<]+)/ ) {
 356             $hash{ "S_ALIGN" } = $1;
 357         } elsif ( $line =~ /<Hsp_midline>([^<]+)/ ) {
 358             $hash{ "MIDLINE" } = $1;
 359         } elsif ( $line =~ /<\/Hsp>/ ) {
 360             push @hsps, dclone \%hash;
 361         }
 362     }
 363
 364     return wantarray ? @hsps : \@hsps;
 365 }
 366
 367
 368 # >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
 369
 370
 371 __END__