t/SearchIO/Tiling.t

   1 #-*-perl-*-
   2 #$Id$
   3 use strict;
   4 use warnings;
   5     use vars qw($EXHAUSTIVE $VERBOSE);
   6 BEGIN {
   7     use Bio::Root::Test;
   8     $EXHAUSTIVE = $ENV{BIOPERL_TILING_EXHAUSTIVE_TESTS};
   9     $VERBOSE    = $ENV{BIOPERL_TILING_VERBOSE_TESTS};
  10     test_begin(-tests => ($EXHAUSTIVE ? 6519 : 1141) );
  11 }
  12
  13 use_ok('Bio::Search::Tiling::MapTiling');
  14 use_ok('Bio::Search::Tiling::MapTileUtils');
  15 use_ok('Bio::SearchIO');
  16 use_ok('Bio::Search::Hit::BlastHit');
  17 use_ok('File::Spec');
  18
  19 my ($blio, $result, $hit, $tiling, $hsp);
  20
  21 my @normal_formats = qw( blast  wublast
  22                          blastn wublastn
  23                          blastp wublastp
  24                          multiblast
  25                          megablast
  26                          rpsblast
  27                          psiblast );
  28 my @xltd_formats  = qw(  blastx wublastx
  29                          tblastn wutblastn
  30                          tblastx wutblastx );
  31
  32
  33 # an exhaustive listing of search reports in
  34 # t/data
  35
  36 my %test_files = (
  37     'blast' => [qw(
  38                ecolitst.bls
  39                frac_problems.blast
  40                frac_problems2.blast
  41                frac_problems3.blast
  42                bl2seq.out
  43                )],
  44     'multiblast' => [qw(
  45                multi_blast.bls
  46                )],
  47     'blastn' => [qw(
  48                a_thaliana.blastn
  49                bl2seq.blastn
  50                new_blastn.txt
  51                hsinsulin.blastcl3.blastn
  52                )],
  53     'wublastn' =>[qw(
  54                brassica_ATH.WUBLASTN
  55                echofilter.wublastn
  56                )],
  57     'blastp' => [qw(
  58                blastp2215.blast
  59                no_hsps.blastp
  60                catalase-webblast.BLASTP
  61                )],
  62     'wublastp' => [qw(
  63                dcr1_sp.WUBLASTP
  64                ecolitst.wublastp
  65                contig-by-hand.wublastp
  66                ecolitst.noseqs.wublastp
  67                )],
  68     'blastx' => [qw(
  69                bl2seq.blastx.out
  70                )],
  71     'wublastx' => [qw(
  72                dnaEbsub_ecoli.wublastx
  73                )],
  74     'wublast' => [qw(
  75                tricky.wublast
  76                )],
  77     'tblastn' => [qw(
  78                tblastn.out
  79                1ZZ19XR301R-Alignment.tblastn
  80                )],
  81     'wutblastn' => [qw(
  82                dnaEbsub_ecoli.wutblastn
  83                )],
  84     'tblastx' => [qw(
  85                bl2seq.tblastx.out
  86                HUMBETGLOA.tblastx
  87                )],
  88     'wutblastx' => [qw(
  89                dnaEbsub_ecoli.wutblastx
  90                )],
  91     'megablast' => [qw(
  92                503384.MEGABLAST.2
  93                )],
  94     'rpsblast' => [qw(
  95                ecoli_domains.rpsblast
  96                )],
  97     'psiblast' => [qw(
  98                psiblastreport.out
  99                )],
 100     'bug2942'  => [qw(
 101                bug2942.blastx
 102                )]
 103     );
 104
 105 # a subset of search reports for
 106 # run-o-the-mill regression tests
 107
 108 my %example_files = (
 109     'blast' => [qw(
 110                ecolitst.bls
 111                )],
 112     'blastn' => [qw(
 113                a_thaliana.blastn
 114                )],
 115     'wublastn' =>[qw(
 116                brassica_ATH.WUBLASTN
 117                )],
 118     'blastp' => [qw(
 119                no_hsps.blastp
 120                catalase-webblast.BLASTP
 121                )],
 122     'wublastp' => [qw(
 123                dcr1_sp.WUBLASTP
 124                )],
 125     'blastx' => [qw(
 126                bl2seq.blastx.out
 127                )],
 128     'wublastx' => [qw(
 129                dnaEbsub_ecoli.wublastx
 130                )],
 131     'wublast' => [qw(
 132                tricky.wublast
 133                )],
 134     'tblastn' => [qw(
 135                tblastn.out
 136                )],
 137     'wutblastn' => [qw(
 138                dnaEbsub_ecoli.wutblastn
 139                )],
 140     'tblastx' => [qw(
 141                HUMBETGLOA.tblastx
 142                )],
 143     'wutblastx' => [qw(
 144                dnaEbsub_ecoli.wutblastx
 145                )],
 146     'megablast' => [qw(
 147                503384.MEGABLAST.2
 148                )]
 149     );
 150
 151 ok( $blio = new Bio::SearchIO(
 152         -file=>test_input_file('dcr1_sp.WUBLASTP'),
 153         -format=>'blast'), 'parse data file');
 154
 155 $result = $blio->next_result;
 156 while ( $_ = $result->next_hit ) {
 157     last if $_->name =~ /ASPTN/;
 158 }
 159 ok($hit = $_, 'got test hit');
 160 ok($tiling = Bio::Search::Tiling::MapTiling->new($hit), 'create tiling');
 161
 162
 163 # TilingI compliance
 164
 165 isa_ok($tiling, 'Bio::Search::Tiling::TilingI');
 166 foreach ( qw( next_tiling rewind_tilings identities conserved length ) ) {
 167     ok( $tiling->$_, "implements '$_'" );
 168 }
 169
 170 # regression test on original calculations
 171
 172 my @orig_id_results = ( 387,388,388,381,382,389 );
 173 my @orig_cn_results = ( 622,619,628,608,611,613 );
 174 my @id_results = (
 175     $tiling->identities('query', 'exact'),
 176     $tiling->identities('query', 'est'),
 177     $tiling->identities('query', 'max'),
 178     $tiling->identities('subject', 'exact'),
 179     $tiling->identities('subject', 'est'),
 180     $tiling->identities('subject', 'max')
 181     );
 182 my @cn_results = (
 183     $tiling->conserved('query', 'exact'),
 184     $tiling->conserved('query', 'est'),
 185     $tiling->conserved('query', 'max'),
 186     $tiling->conserved('subject', 'exact'),
 187     $tiling->conserved('subject', 'est'),
 188     $tiling->conserved('subject', 'max')
 189     );
 190 map { $_ = int($_) } @id_results, @cn_results;
 191
 192 is_deeply(\@id_results, \@orig_id_results, 'identities regression test');
 193 is_deeply(\@cn_results, \@orig_cn_results, 'conserved regression test');
 194
 195 # tiling iterator regression tests
 196
 197 my ($qn, $sn)=(0,0);
 198 while ($tiling->next_tiling('query')) {$qn++};
 199 while ($tiling->next_tiling('subject')) {$sn++};
 200 is ($qn, 8, 'tiling iterator regression test(1)');
 201 is ($sn, 128, 'tiling iterator regression test(2)');
 202 $tiling->rewind('subject');
 203 while ($tiling->next_tiling('subject')) {$sn++};
 204 is ($sn, 256, 'tiling iterator regression test(3, rewind)');
 205
 206 diag("Old blast.t tiling tests") if $VERBOSE;
 207
 208 ok($blio = Bio::SearchIO->new(
 209     '-format' => 'blast',
 210     '-file'   => test_input_file('ecolitst.wublastp')
 211    ), "ecolitst.wublastp");
 212 $result = $blio->next_result;
 213 $result->next_hit;
 214 $hit = $result->next_hit;
 215 $tiling = Bio::Search::Tiling::MapTiling->new($hit);
 216 # Test HSP contig data returned by SearchUtils::tile_hsps()
 217 # Second hit has two hsps that overlap.
 218
 219 # compare with the contig made by hand for these two contigs
 220 # in t/data/contig-by-hand.wublastp
 221 # (in this made-up file, the hsps from ecolitst.wublastp
 222 #  were aligned and contiged, and Length, Identities, Positives
 223 #  were counted, by a human (maj) )
 224
 225 my $hand_hit = Bio::SearchIO->new(
 226     -format=>'blast',
 227     -file=>test_input_file('contig-by-hand.wublastp')
 228     )->next_result->next_hit;
 229 my $hand_hsp = $hand_hit->next_hsp;
 230 my @hand_qrng = $hand_hsp->range('query');
 231 my @hand_srng = $hand_hsp->range('hit');
 232 my @hand_matches = $hand_hit->matches;
 233
 234 is(($tiling->range('query'))[0], $hand_qrng[0]);
 235 is(($tiling->range('query'))[1], $hand_qrng[1]);
 236 is(sprintf("%d",$tiling->identities('query')), $hand_matches[0]);
 237 is(sprintf("%d",$tiling->conserved('query')), $hand_matches[1]);
 238 is(($tiling->range('hit'))[0], $hand_srng[0]);
 239 is(($tiling->range('hit'))[1], $hand_srng[1]);
 240 is(sprintf("%d",$tiling->identities('hit')), $hand_matches[0]);
 241 is(sprintf("%d",$tiling->conserved('hit')), $hand_matches[1]);
 242
 243 ok( $blio = Bio::SearchIO->new(
 244         '-format' => 'blast',
 245         '-file'   => test_input_file('dnaEbsub_ecoli.wublastx')
 246     ), "dnaEbsub_ecoli.wublastx");
 247
 248 $hit = $blio->next_result->next_hit;
 249 $tiling = Bio::Search::Tiling::MapTiling->new($hit);
 250 is(sprintf("%.3f",$tiling->frac_identical(-type=>'query',-denom=>'aligned',-context=>'p2')), '0.364');
 251 is(sprintf("%.3f",$tiling->frac_identical(-type=>'hit',-denom=>'aligned',-context=>'all')), '0.366');
 252 is(sprintf("%.3f",$tiling->frac_conserved(-type=>'query',-denom=>'aligned',-context=>'p2')), '0.537');
 253 is(sprintf("%.3f",$tiling->frac_conserved(-type=>'hit',-denom=>'aligned',-context=>'all')), '0.540');
 254 is(sprintf("%.2f",$tiling->frac_aligned_query(-context=>'p2')), '0.62');
 255 is(sprintf("%.2f",$tiling->frac_aligned_hit(-context=>'all')), '0.71');
 256
 257 ok( $blio = Bio::SearchIO->new(
 258         '-format' => 'blast',
 259         '-file'   => test_input_file('tricky.wublast')
 260     ), "tricky.wublast");
 261
 262 $hit = $blio->next_result->next_hit;
 263 $tiling = Bio::Search::Tiling::MapTiling->new($hit);
 264 cmp_ok sprintf("%.3f",$tiling->frac_identical(-denom => 'aligned')), '>', 0.2, 'tricky.wublast(1)';
 265 cmp_ok sprintf("%.3f",$tiling->frac_conserved(-denom => 'aligned')), '<=', 1, 'tricky.wublast(2)';
 266 is(sprintf("%.2f",$tiling->frac_aligned_query), '0.92', 'tricky.wublast(3)');
 267 is(sprintf("%.2f",$tiling->frac_aligned_hit), '0.91','tricky.wublast(4)');
 268
 269 diag("New tiling tests") if $VERBOSE;
 270
 271 # select test file set based on the environment variable
 272 # BIOPERL_TILING_EXHAUSTIVE_TESTS
 273
 274 my $files = ($EXHAUSTIVE ? \%test_files : \%example_files);
 275
 276 foreach my $alg (@normal_formats, @xltd_formats) {
 277     diag("*******$alg files*******") if ($files->{$alg} && $VERBOSE);
 278     foreach my $tf (@{$files->{$alg}}) {
 279         ok( $blio = Bio::SearchIO->new( -format=>'blast',
 280                                         -file=>test_input_file($tf)
 281             ), "$tf" );
 282         $result = $blio->next_result;
 283         my $hit_count = 0;
 284         # compare the per-aligned-base identity avg over hsps
 285         # with frac_identical (bzw, conserved)
 286
 287       HIT:
 288         while ( $hit = $result->next_hit ) {
 289             ++$hit_count;
 290             # quiet the "No HSPs" warning with -verbose => -1
 291             ok( $tiling = Bio::Search::Tiling::MapTiling->new(-hit=>$hit,-verbose=>-1), "tile $tf hit $hit_count #hsps ".scalar $tiling->hsps );
 292             my @hsps = $tiling->hsps;
 293
 294             unless (@hsps) {
 295                 diag( "--no hsps for $tf hit $hit_count") if $VERBOSE;
 296                 next HIT;
 297             }
 298             my ($dpct, $est, $fast,$exact, $max);
 299             my $tol = 0.10; # % difference accepted as approx. equal
 300
 301             ## loop through contexts:
 302             for my $type (qw( query hit )) {
 303                 for my $context ($tiling->contexts($type)) {
 304                     diag(" --- $type $context ---") if $VERBOSE;
 305                     if (scalar($tiling->contexts($type, $context)) == 1) {
 306                         # equality
 307                         ($dpct, $est, $fast) = $tiling->cmp_frac($type,'identical','aligned', 'est', 'fast', $context);
 308                         is( $est,$fast, substr($type,0,1)." id: est ($est) = fast ($fast)");
 309                         ($dpct, $est, $fast) = $tiling->cmp_frac($type,'conserved','aligned', 'est', 'fast', $context);
 310                         is( $est,$fast, substr($type,0,1)." cn: est ($est) = fast ($fast)");
 311                     }
 312                     else {
 313                         # comparisons
 314                         ($dpct, $est, $fast) = $tiling->cmp_frac($type,'identical','aligned', 'est', 'fast', $context);
 315 #                       cmp_ok( $dpct, "<", $tol,
 316 #                               substr($type,0,1)." id: est ($est) ~ fast ($fast)");
 317                         ($dpct, $exact, $max) = $tiling->cmp_frac($type,'identical','aligned', 'exact', 'max', $context);
 318                         cmp_ok( abs($exact-$est)/$exact, "<" , $tol,
 319                                 substr($type,0,1)." id: exact ($exact) ~ est ($est)");
 320                         cmp_ok( $exact, "<=" , $max,
 321                                 substr($type,0,1)." id: exact ($exact) <= max ($max)");
 322
 323                         ($dpct, $est, $fast) = $tiling->cmp_frac($type,'conserved','aligned', 'est', 'fast', $context);
 324 #                       cmp_ok( $dpct, "<", $tol,
 325 #                               substr($type,0,1)." cn: est ($est) ~ fast ($fast)");
 326                         ($dpct, $exact, $max) = $tiling->cmp_frac($type,'conserved','aligned', 'exact', 'max', $context);
 327                         cmp_ok(  abs($exact-$est)/$exact, "<" , $tol,
 328                                  substr($type,0,1)." cn: exact ($exact) ~ est ($est)");
 329                         cmp_ok( $exact, "<=" , $max,
 330                                 substr($type,0,1)." cn: exact ($exact) <= max ($max)");
 331                     }
 332                 }
 333             }
 334         }
 335     }
 336 }
 337
 338 # bug 2942
 339
 340 my %expected_ranges = ( 'm0' => [7, 11037], #query
 341                         'm1' => [1770, 10865], #query
 342                         'm2' => [2462, 14599], #query
 343                         'all' => [231, 3563] #subject
 344     );
 345 $blio = Bio::SearchIO->new( -file=>test_input_file( $test_files{'bug2942'}->[0] ),
 346                             -format => 'blast' );
 347 $hit = $blio->next_result->next_hit;
 348 $tiling = Bio::Search::Tiling::MapTiling->new($hit);
 349 for ( 'm0', 'm1', 'm2' ) {
 350     is_deeply( [$tiling->range('query',$_)], $expected_ranges{$_}, "bug2942: query $_: range correct");
 351 }
 352 is_deeply( [$tiling->range('subject', 'all')], $expected_ranges{'all'}, "bug2942: subject all : range correct" );
 353
 354 # test get_tiled_alns
 355
 356 $blio = Bio::SearchIO->new( -file=>test_input_file( 'dcr1_sp.WUBLASTP' ) );
 357 $result = $blio->next_result;
 358 while ($hit = $result->next_hit) {
 359     last if $hit->name =~ /ASPTN/;
 360 }
 361
 362 $tiling = Bio::Search::Tiling::MapTiling->new($hit);
 363
 364 ok my @alns = $tiling->get_tiled_alns, "get_tiled_alns";
 365 is scalar @alns, 6, "got all alns";
 366
 367 for my $aln ( @alns ) {
 368     my (@aint, @qint, @sint);
 369     my $qs = $aln->get_seq_by_id('query');
 370     my $ss = $aln->get_seq_by_id('subject');
 371     ok my @qfeats = $qs->get_SeqFeatures;
 372     foreach (@qfeats) {
 373         push @aint, [$_->start, $_->end];
 374         push @qint, [($_->get_tag_values('query_start'))[0],
 375                      ($_->get_tag_values('query_end'))[0] ];
 376     }
 377     is( eval(join('+', map {$$_[1]-$$_[0]+1} @aint)),
 378         eval(join('+', map {$$_[1]-$$_[0]+1} @qint)),
 379         "aln and qfeat lengths correspond" );
 380     is( $qs->length - $qs->num_gaps('-'), eval(join('+', map {$$_[1]-$$_[0]+1} @qint)), "q length correct");
 381     ok my @hfeats = $ss->get_SeqFeatures;
 382     @aint = ();
 383     ok ( @qfeats == @hfeats, "features on q and s correspond");
 384     foreach (@hfeats) {
 385         push @aint, [$_->start, $_->end];
 386         push @sint, [($_->get_tag_values('subject_start'))[0],
 387                      ($_->get_tag_values('subject_end'))[0] ];
 388     }
 389     is( eval(join('+', map {$$_[1]-$$_[0]+1} @aint)),
 390         eval(join('+', map {$$_[1]-$$_[0]+1} @sint)),
 391         "aln and hfeat lengths correspond" );
 392     is( $ss->length - $ss->num_gaps('-'), eval(join('+', map {$$_[1]-$$_[0]+1} @sint)), "s length correct");
 393
 394 }
 395 1;
 396
 397
 398
 399 package Bio::Search::Tiling::MapTiling;
 400
 401 sub cmp_frac {
 402     my ($tiling, $type, $method, $denom, @actions) = @_;
 403     my ($a, $b);
 404     my $context = ($actions[2] ? $actions[2] : 'all');
 405     $a = $tiling->frac(-type=>$type,
 406                        -method=>$method,
 407                        -denom=>$denom,
 408                        -action=>$actions[0],
 409                        -context=>$context);
 410     $b = $tiling->frac(-type=>$type,
 411                        -method=>$method,
 412                        -denom=>$denom,
 413                        -action=>$actions[1],
 414                        -context=>$context);
 415     return ( abs($a-$b)/$a, f(5,$a), f(5,$b) );
 416 }
 417
 418 sub f { my ($d,$val) = @_; sprintf("%.${d}f",$val) }
 419
 420
 421
 422
 423 1;
 424