t/phd.t

   1 # -*-Perl-*- Test Harness script for Bioperl
   2 # $Id$
   3
   4 use strict;
   5
   6 BEGIN {
   7     use lib 't/lib';
   8     use BioperlTest;
   9
  10     test_begin(-tests => 18);
  11
  12         use_ok('Bio::SeqIO');
  13 }
  14
  15 my $DEBUG = test_debug();
  16
  17 print("Checking to see if Bio::Seq::Quality objects can be created from a file...\n") if ($DEBUG);
  18 my $in_phd  = Bio::SeqIO->new('-file' => test_input_file('phredfile.phd'),
  19                               '-format'  => 'phd',
  20                               '-verbose' => $DEBUG);
  21 isa_ok($in_phd,'Bio::SeqIO::phd');
  22
  23
  24
  25 my $phd = $in_phd->next_seq();
  26 is($phd->quality_levels,'99',"Did you get the 'QUALITY_LEVELS' comment?");
  27 isa_ok($phd,"Bio::Seq::Quality");
  28
  29
  30 if( $DEBUG ) {
  31     my $position = 6;
  32     print("I saw these in phredfile.phd:\n\n");
  33     print $_->tagname,": ",$_->display_text || 0," \n"
  34         for ($phd->annotation->get_Annotations('header'));
  35
  36     print("What is the base at position $position (using subseq)?\n");
  37     print($phd->subseq($position,$position)."\n");
  38     print("What is the base at position $position (using baseat)?\n");
  39     print($phd->baseat($position)."\n");
  40     print("What is the quality at $position? (using subqual)\n");
  41
  42 my @qualsretr = @{$phd->subqual($position,$position)};
  43     print($qualsretr[0]."\n");
  44     print("What is the quality at $position? (using qualat)\n");
  45     print($phd->qualat($position)."\n");
  46     print("What is the trace at $position? (using trace_index_at)\n");
  47     print($phd->trace_index_at($position)."\n");
  48     print("What is the trace at $position? (using subtrace)\n");
  49     my @tracesretr = @{$phd->subtrace($position,$position)};
  50     print($tracesretr[0]."\n");
  51 }
  52
  53 print("OK. Now testing write_phd...\n") if($DEBUG);
  54
  55 my $outfile = test_output_file();
  56 my $out_phd = Bio::SeqIO->new(-file => ">$outfile",
  57                               '-format' => 'phd');
  58 isa_ok($out_phd,"Bio::SeqIO::phd");
  59
  60 $out_phd->write_seq($phd);
  61
  62 ok( -s $outfile);
  63
  64 # Bug 2120
  65
  66 my @qual = q(9 9 12 12 8 8 9 8 8 8 9);
  67 my @trace = q(113 121 130 145 153 169 177 203 210 218 234);
  68
  69 $in_phd  = Bio::SeqIO->new('-file' => test_input_file('bug2120.phd'),
  70                               '-format'  => 'phd',
  71                               '-verbose' => $DEBUG);
  72
  73 my $seq = $in_phd->next_seq;
  74 is($seq->subseq(10,20),'gggggccttat','$seq->subseq()');
  75 my @seq_qual =$seq->subqual_text(10,20);
  76 is_deeply(\@seq_qual,\@qual,'$seq->subqual_tex()');
  77 my @seq_trace = $seq->subtrace_text(10,20);
  78 is_deeply(\@seq_trace,\@trace,'$seq->subqual_tex()');
  79
  80 if($DEBUG) {
  81     print "\nDefault header ... \n\n";
  82     use Bio::Seq::Quality;
  83     my $seq = Bio::Seq::Quality->new('-seq' => 'GAATTC');
  84     $out_phd->_fh(\*STDOUT);
  85     $out_phd->write_header($seq);
  86     print "Complete output\n\n";
  87     $out_phd->write_seq($seq);
  88 }
  89
  90 print("Testing the header manipulation\n") if($DEBUG);
  91 is($phd->chromat_file(),'ML4924R','$phd->chromat_file()');
  92 $phd->chromat_file('ML4924R.esd');
  93 is($phd->chromat_file(), 'ML4924R.esd','$phd->chromat_file()');
  94 $phd->touch();
  95 my $localtime = localtime();
  96 is($phd->time, "$localtime");
  97 if ($DEBUG){
  98     print "Testing the sequence ...\n";
  99     print ">",$phd->id," ",$phd->desc,"\n",$phd->seq,"\n";
 100     my $revcom = $phd->revcom;
 101     print ">revcom\n",$revcom->seq,"\n";
 102     print ">revcom_qual at 6\n",$revcom->qualat(6),"\n";
 103     print ">revcom_trace at 6 !!\n",$revcom->trace_index_at(6),"\n";
 104     my $trunc = $phd->trunc(10,20);
 105     print ">TRUNC 10,20\n",$trunc->seq,"\n>qual\n@{$trunc->qual}\n>trace\n@{$trunc->trace}\n";
 106 }
 107
 108 # Multiple seqs in one file
 109
 110 $in_phd  = Bio::SeqIO->new('-file' => test_input_file('multi.phd'),
 111                               '-format'  => 'phd',
 112                               '-verbose' => $DEBUG);
 113
 114 @qual = qq(9 9 15 17 17 22 22 25 25 22 22);
 115 @trace = qq(98 105 119 128 143 148 162 173 185 197 202);
 116
 117 $seq = $in_phd->next_seq;
 118 is($seq->subseq(10,20),'tctcgagggta','$seq->subseq()');
 119 @seq_qual =$seq->subqual_text(10,20);
 120 is_deeply(\@seq_qual,\@qual,'$seq->subqual_tex()');
 121 @seq_trace = $seq->subtrace_text(10,20);
 122 is_deeply(\@seq_trace,\@trace,'$seq->subqual_tex()');
 123
 124 @qual = qq(11 9 6 6 9 19 20 32 34 34 39);
 125 @trace = qq(98 104 122 128 140 147 159 167 178 190 200);
 126
 127 $seq = $in_phd->next_seq;
 128 is($seq->subseq(10,20),'gcctgcaggta','$seq->subseq()');
 129 @seq_qual =$seq->subqual_text(10,20);
 130 is_deeply(\@seq_qual,\@qual,'$seq->subqual_tex()');
 131 @seq_trace = $seq->subtrace_text(10,20);
 132 is_deeply(\@seq_trace,\@trace,'$seq->subqual_tex()');
 133
 134 #if($DEBUG) {
 135 #    print "\nDefault header ... \n\n";
 136 #    use Bio::Seq::Quality;
 137 #    my $seq = Bio::Seq::Quality->new('-seq' => 'GAATTC');
 138 #    $out_phd->_fh(\*STDOUT);
 139 #    $out_phd->write_header($seq);
 140 #    print "Complete output\n\n";
 141 #    $out_phd->write_seq($seq);
 142 #}
 143
 144 ##print("Testing the header manipulation\n") if($DEBUG);
 145 #is($phd->chromat_file(),'ML4924R','$phd->chromat_file()');
 146 #$phd->chromat_file('ML4924R.esd');
 147 #is($phd->chromat_file(), 'ML4924R.esd','$phd->chromat_file()');
 148 #$phd->touch();
 149 #my $localtime = localtime();
 150 #is($phd->time, "$localtime");
 151 #if ($DEBUG){
 152 #    print "Testing the sequence ...\n";
 153 #    print ">",$phd->id," ",$phd->desc,"\n",$phd->seq,"\n";
 154 #    my $revcom = $phd->revcom;
 155 #    print ">revcom\n",$revcom->seq,"\n";
 156 #    print ">revcom_qual at 6\n",$revcom->qualat(6),"\n";
 157 #    print ">revcom_trace at 6 !!\n",$revcom->trace_index_at(6),"\n";
 158 #    my $trunc = $phd->trunc(10,20);
 159 #    print ">TRUNC 10,20\n",$trunc->seq,"\n>qual\n@{$trunc->qual}\n>trace\n@{$trunc->trace}\n";
 160 #}
 161 #