File Coverage

Bio/SeqIO/fastq.pm

Criterion	Covered	Total	%
statement	125	144	86.8
branch	52	70	74.2
condition	33	39	84.6
subroutine	9	13	69.2
pod	8	8	100.0
total	227	274	82.8

# POD at __END__, let the code begin...

package Bio::SeqIO::fastq;
use strict;

use Bio::Seq::SeqFactory;

use base qw(Bio::SeqIO);

our %variant = (
        sanger     => {
            'offset'     => 33,
            'qual_start' => 0,
            'qual_end'   => 93
            },
        solexa     => {
            'offset'     => 64,
            'qual_start' => -5,
            'qual_end'   => 62
            },
        illumina   => {
            'offset'     => 64,
            'qual_start' => 0,
            'qual_end'   => 62
            },
    );

sub _initialize {
    my($self,@args) = @_;
    $self->SUPER::_initialize(@args);
    my ($variant, $validate, $header) = $self->_rearrange([qw(VARIANT
                                                              VALIDATE
                                                              QUALITY_HEADER)], @args);
    $variant ||= 'sanger';
    $self->variant($variant);
    $self->_init_tables($variant);
    $validate = defined $validate ? $validate : 1;
    $self->validate($validate);
    $header     && $self->quality_header($header);

    if( ! defined $self->sequence_factory ) {
        $self->sequence_factory(Bio::Seq::SeqFactory->new(
            -verbose => $self->verbose(),
            -type => 'Bio::Seq::Quality')
        );
    }
}

sub next_seq {
    my( $self ) = @_;
    while (defined(my $data = $self->next_dataset)) {
        # Are FASTQ sequences w/o any sequence valid?  Removing for now
        # -cjfields 6.22.09
        my $seq = $self->sequence_factory->create(%$data);
        return $seq;
    }
    return;
}

# pure perl version
sub next_dataset {
    my $self = shift;
    local $/ = "\n";
    my $data;
    my $mode = '-seq';
    # speed this up by directly accessing the filehandle and in-lining the
    # _readline stuff vs. making the repeated method calls. Tradeoff is speed
    # over repeated code.

    # we can probably normalize line endings using PerlIO::eol or
    # Encode::Newlines

    my $fh = $self->_fh;
    my $line = $self->{lastline} || <$fh>;

    FASTQ:
    while (defined $line) {
        $line =~ s/\015\012/\012/;
        $line =~ tr/\015/\n/;
        if ($mode eq '-seq' && $line =~ m{^@([^\n]+)$}xmso) {
            $data->{-descriptor} = $1;
            my ($id,$fulldesc);
            if ($data->{-descriptor} =~ /^\s*(\S+)\s*(.*)/) {
                ($id,$fulldesc) = ($1, $2);
            } else {
                $self->throw("Can't parse fastq header");
            }
            $data->{-id} = $id;
            $data->{-desc} = $fulldesc;
            $data->{-namespace} = $self->variant;
        } elsif ($mode eq '-seq' && $line =~ m{^\+([^\n]*)}xmso) {
            my $desc = $1;
            $self->throw("No description line parsed") unless $data->{-descriptor};
            if ($desc && $data->{-descriptor} ne $desc) {
                $self->throw("Quality descriptor [$desc] doesn't match seq ".
                    "descriptor ".$data->{-descriptor}.", line: $." );
            }
            $mode = '-raw_quality';
        } else {
            if ($mode eq '-raw_quality' && defined($data->{-raw_quality}) &&
                (length($data->{-raw_quality}) >= length($data->{-seq}))) {
                $self->{lastline} = $line;
                last FASTQ
            }
            chomp $line;
            if ($line =~ /^$/) {
                delete $self->{lastline};
                last FASTQ;
            }
            $data->{$mode} .= $line
        }
        $line = <$fh>;
        if (!defined $line) {
            delete $self->{lastline};
            last FASTQ;
        }
    }

    return unless $data;
    if (!$data->{-seq} || !defined($data->{-raw_quality})) {
        $self->throw("Missing sequence and/or quality data; line: $.");
    }

    # simple quality control tests
    if (length $data->{-seq} != length $data->{-raw_quality}) {
        $self->throw("Quality string [".$data->{-raw_quality}."] of length [".
            length($data->{-raw_quality})."]\ndoesn't match length of sequence ".
            $data->{-seq}."\n[".length($data->{-seq})."], line: $.");
    }

    $data->{-qual} = [map {
        if ($self->{_validate_qual} && !exists($self->{chr2qual}->{$_})) {
            $self->throw("Unknown symbol with ASCII value ".ord($_)." outside ".
                "of quality range")
            # TODO: fallback?
        }
        $self->variant eq 'solexa' ?
            $self->{sol2phred}->{$self->{chr2qual}->{$_}}:
            $self->{chr2qual}->{$_};
    } unpack("A1" x length($data->{-raw_quality}), $data->{-raw_quality})];
    return $data;
}

# This should be creating fastq output only.  Bio::SeqIO::fasta and
# Bio::SeqIO::qual should be used for that output

sub write_seq {
    my ($self,@seq) = @_;
    my $var = $self->variant;
    foreach my $seq (@seq) {
        unless ($seq->isa("Bio::Seq::Quality")){
            $self->warn("You can't write FASTQ without supplying a Bio::Seq::".
                "Quality object! ".ref($seq)."\n");
            next;
        }
        my $str = $seq->seq || '';
        my @qual = @{$seq->qual};

        # this should be the origin of the sequence (illumina, solexa, sanger)
        my $ns= $seq->namespace;

        my $top = $seq->display_id();
        if (my $desc = $seq->desc()) {
            $desc =~ s/\n//g;
            $top .= " $desc";
        }
        my $qual = '';
        my $qual_map =
            ($ns eq 'solexa' && $var eq 'solexa') ? $self->{phred_fp2chr} :
            ($var eq 'solexa')                    ? $self->{phred_int2chr} :
            $self->{qual2chr};

        my %bad_qual;
        for my $q (@qual) {
            $q = sprintf("%.0f", $q) if ($var ne 'solexa' && $ns eq 'solexa');
            if (exists $qual_map->{$q}) {
                $qual .= $qual_map->{$q};
                next;
            } else {
                # fuzzy mapping, for edited qual scores
                my $qr = sprintf("%.0f",$q);
                my $bounds = sprintf("%.1f-%.1f",$qr-0.5, $qr+0.5);
                if (exists $self->{fuzzy_qual2chr}->{$bounds}) {
                    $qual .= $self->{fuzzy_qual2chr}->{$bounds};
                    next;
                } else {
                    my $rep = ($q <= $self->{qual_start}) ?
                        $qual_map->{$self->{qual_start}} : $qual_map->{$self->{qual_end}};
                    $qual .= $rep;
                    $bad_qual{$q}++;
                }
            }
        }
        if ($self->{_validate_qual} && %bad_qual) {
            $self->warn("Data loss for $var: following values not found\n".
                        join(',',sort {$a <=> $b} keys %bad_qual))
        }
        $self->_print("\@",$top,"\n",$str,"\n") or return;
        $self->_print("+",($self->{_quality_header} ? $top : ''),"\n",$qual,"\n") or return;
    }
    return 1;
}

sub write_fastq {
    my ($self,@seq) = @_;
    return $self->write_seq(@seq);
}

sub write_fasta {
    my ($self,@seq) = @_;
    if (!exists($self->{fasta_proxy})) {
        $self->{fasta_proxy} = Bio::SeqIO->new(-format => 'fasta', -fh => $self->_fh);
    }
    return $self->{fasta_proxy}->write_seq(@seq);
}

sub write_qual {
    my ($self,@seq) = @_;
    if (!exists($self->{qual_proxy})) {
        $self->{qual_proxy} = Bio::SeqIO->new(-format => 'qual', -fh => $self->_fh);
    }
    return $self->{qual_proxy}->write_seq(@seq);
}

# variant() method inherited from Bio::Root::IO

sub _init_tables {
    my ($self, $var) = @_;
    # cache encode/decode values for quicker accession
    ($self->{qual_start}, $self->{qual_end}, $self->{qual_offset}) =
        @{ $variant{$var} }{qw(qual_start qual_end offset)};
    if ($var eq 'solexa') {
        for my $q ($self->{qual_start} .. $self->{qual_end}) {
            my $char = chr($q + $self->{qual_offset});
            $self->{chr2qual}->{$char} = $q;
            $self->{qual2chr}->{$q} = $char;
            my $s2p = 10 * log(1 + 10 ** ($q / 10.0)) / log(10);

            # solexa <=> solexa mapping speedup (retain floating pt precision)
            $self->{phred_fp2chr}->{$s2p} = $char;
            $self->{sol2phred}->{$q} = $s2p;

            # this is for mapping values fuzzily (fallback)
            $self->{fuzzy_qual2chr}->{sprintf("%.1f-%.1f",$q - 0.5, $q + 0.5)} = $char;

            next if $q < 0; # skip loop; PHRED scores greater than 0
            my $p2s = sprintf("%.0f",($q <= 1) ? -5 : 10 * log(-1 + 10 ** ($q / 10.0)) / log(10));
            # sanger/illumina PHRED <=> Solexa char mapping speedup
            $self->{phred_int2chr}->{$q} = chr($p2s + $self->{qual_offset});
        }
    } else {
        for my $c ($self->{qual_start}..$self->{qual_end}) {
            # PHRED mapping
            my $char = chr($c + $self->{qual_offset});
            $self->{chr2qual}->{$char} = $c;
            $self->{qual2chr}->{$c} = $char;
            # this is for mapping values not found with above
            $self->{fuzzy_qual2chr}->{sprintf("%.1f-%.1f",$c - 0.5, $c + 0.5)} = $char;
        }
    }
}

sub validate {
    my ($self, $val) = @_;
    if (defined $val) {
        $self->{_validate_qual} = $val;
    }
    return $self->{_validate_qual};
}

sub quality_header{
    my ($self, $val) = @_;
    if (defined $val) {
        $self->{_quality_header} = $val;
    }
    return $self->{_quality_header} || 0;
}

1;

__END__