code_perl/Maasha/Blast.pm

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