The Perl Toolchain Summit needs more sponsors. If your company depends on Perl, please support this very important event.
# -*-Perl-*- Test Harness script for Bioperl
# $Id$

use strict;

BEGIN {
    use lib '.';
    use Bio::Root::Test;
    
    test_begin(-tests => 21);
	
	use_ok('Bio::SeqIO::phd');
}

my $DEBUG = test_debug();

print("Checking to see if Bio::Seq::Quality objects can be created from a file...\n") if ($DEBUG);
my $in_phd  = Bio::SeqIO->new('-file' => test_input_file('test.phd'),
			      '-format'  => 'phd',
			      '-verbose' => $DEBUG);
isa_ok($in_phd,'Bio::SeqIO::phd');



my $phd = $in_phd->next_seq();
is($phd->quality_levels,'99',"Did you get the 'QUALITY_LEVELS' comment?");
isa_ok($phd,"Bio::Seq::Quality");


if( $DEBUG ) {
    my $position = 6;
    print("I saw these in phredfile.phd:\n\n");
    print $_->tagname,": ",$_->display_text || 0," \n"
        for ($phd->annotation->get_Annotations('header'));

    print("What is the base at position $position (using subseq)?\n");
    print($phd->subseq($position,$position)."\n");
    print("What is the base at position $position (using baseat)?\n");
    print($phd->baseat($position)."\n");
    print("What is the quality at $position? (using subqual)\n");

my @qualsretr = @{$phd->subqual($position,$position)};
    print($qualsretr[0]."\n");
    print("What is the quality at $position? (using qualat)\n");
    print($phd->qualat($position)."\n");
    print("What is the trace at $position? (using trace_index_at)\n");
    print($phd->trace_index_at($position)."\n");
    print("What is the trace at $position? (using subtrace)\n");
    my @tracesretr = @{$phd->subtrace($position,$position)};
    print($tracesretr[0]."\n");
}

print("OK. Now testing write_phd...\n") if($DEBUG);

my $outfile = test_output_file();
my $out_phd = Bio::SeqIO->new(-file => ">$outfile",
			      '-format' => 'phd');
isa_ok($out_phd,"Bio::SeqIO::phd");

$out_phd->write_seq($phd);

ok( -s $outfile);

# Bug 2120

my @qual = q(9 9 12 12 8 8 9 8 8 8 9);
my @trace = q(113 121 130 145 153 169 177 203 210 218 234);

$in_phd  = Bio::SeqIO->new('-file' => test_input_file('bug2120.phd'),
			      '-format'  => 'phd',
			      '-verbose' => $DEBUG);

my $seq = $in_phd->next_seq;
is($seq->subseq(10,20),'gggggccttat','$seq->subseq()');
my @seq_qual =$seq->subqual_text(10,20);
is_deeply(\@seq_qual,\@qual,'$seq->subqual_tex()');
my @seq_trace = $seq->subtrace_text(10,20);
is_deeply(\@seq_trace,\@trace,'$seq->subqual_tex()');

if($DEBUG) {
    print "\nDefault header ... \n\n";
    use Bio::Seq::Quality;
    my $seq = Bio::Seq::Quality->new('-seq' => 'GAATTC');
    $out_phd->_fh(\*STDOUT);
    $out_phd->write_header($seq);
    print "Complete output\n\n";
    $out_phd->write_seq($seq);
}

print("Testing the header manipulation\n") if($DEBUG);
is($phd->chromat_file(),'ML4924R','$phd->chromat_file()');
$phd->chromat_file('ML4924R.esd');
is($phd->chromat_file(), 'ML4924R.esd','$phd->chromat_file()');
$phd->touch();

# Commented out 1/17/09.
# This isn't exactly a stable regression test as the comparison tests
# localtime() called from two different timepoints. They can differ if the calls
# occurred before and after a change in seconds, for example.

#my $localtime = localtime();
#is($phd->time, "$localtime", $phd->time.':'.$localtime);

if ($DEBUG){
    print "Testing the sequence ...\n";
    print ">",$phd->id," ",$phd->desc,"\n",$phd->seq,"\n";
    my $revcom = $phd->revcom;
    print ">revcom\n",$revcom->seq,"\n";
    print ">revcom_qual at 6\n",$revcom->qualat(6),"\n";
    print ">revcom_trace at 6 !!\n",$revcom->trace_index_at(6),"\n";
    my $trunc = $phd->trunc(10,20);
    print ">TRUNC 10,20\n",$trunc->seq,"\n>qual\n@{$trunc->qual}\n>trace\n@{$trunc->trace}\n";
}

# Multiple seqs in one file

$in_phd  = Bio::SeqIO->new('-file' => test_input_file('multi.phd'),
			      '-format'  => 'phd',
			      '-verbose' => $DEBUG);

@qual = qq(9 9 15 17 17 22 22 25 25 22 22);
@trace = qq(98 105 119 128 143 148 162 173 185 197 202);

$seq = $in_phd->next_seq;
is($seq->id, 'ML4924F');
is($seq->subseq(10,20),'tctcgagggta','$seq->subseq()');
@seq_qual =$seq->subqual_text(10,20);
is_deeply(\@seq_qual,\@qual,'$seq->subqual_tex()');
@seq_trace = $seq->subtrace_text(10,20);
is_deeply(\@seq_trace,\@trace,'$seq->subqual_tex()');

@qual = qq(11 9 6 6 9 19 20 32 34 34 39);
@trace = qq(98 104 122 128 140 147 159 167 178 190 200);

$seq = $in_phd->next_seq;
is($seq->id, 'ML4924R');
is($seq->subseq(10,20),'gcctgcaggta','$seq->subseq()');
@seq_qual =$seq->subqual_text(10,20);
is_deeply(\@seq_qual,\@qual,'$seq->subqual_tex()');
@seq_trace = $seq->subtrace_text(10,20);
is_deeply(\@seq_trace,\@trace,'$seq->subqual_tex()');

#if($DEBUG) {
#    print "\nDefault header ... \n\n";
#    use Bio::Seq::Quality;
#    my $seq = Bio::Seq::Quality->new('-seq' => 'GAATTC');
#    $out_phd->_fh(\*STDOUT);
#    $out_phd->write_header($seq);
#    print "Complete output\n\n";
#    $out_phd->write_seq($seq);
#}

##print("Testing the header manipulation\n") if($DEBUG);
#is($phd->chromat_file(),'ML4924R','$phd->chromat_file()');
#$phd->chromat_file('ML4924R.esd');
#is($phd->chromat_file(), 'ML4924R.esd','$phd->chromat_file()');
#$phd->touch();
#my $localtime = localtime();
#is($phd->time, "$localtime");
#if ($DEBUG){
#    print "Testing the sequence ...\n";
#    print ">",$phd->id," ",$phd->desc,"\n",$phd->seq,"\n";
#    my $revcom = $phd->revcom;
#    print ">revcom\n",$revcom->seq,"\n";
#    print ">revcom_qual at 6\n",$revcom->qualat(6),"\n";
#    print ">revcom_trace at 6 !!\n",$revcom->trace_index_at(6),"\n";
#    my $trunc = $phd->trunc(10,20);
#    print ">TRUNC 10,20\n",$trunc->seq,"\n>qual\n@{$trunc->qual}\n>trace\n@{$trunc->trace}\n";
#}
#

# Whole-read tags in the file
$in_phd  = Bio::SeqIO->new('-file' => test_input_file('multiseq_tags.phd'),
			      '-format'  => 'phd',
			      '-verbose' => $DEBUG);
isa_ok($in_phd,'Bio::SeqIO::phd');
my @seqs = ();
while (my $seq = $in_phd->next_seq){
  push @seqs, $seq;
}
is( scalar @seqs, 2 );