File Coverage

blib/lib/BioX/Seq/Fetch.pm

Criterion	Covered	Total	%
statement	146	149	97.9
branch	51	64	81.2
condition	16	20	80.0
subroutine	12	12	100.0
pod	5	5	100.0
total	230	250	92.4

package BioX::Seq::Fetch;

use 5.012;
use strict;
use warnings;

use BioX::Seq;
use File::stat;

use constant MAGIC_GZIP => pack('C3', (0x1f, 0x8b, 0x08));

sub new {

    my ($class,$fn,%args) = @_;

    die "Must define input filename for random access"
        if (! defined $fn);

    my $self = bless {fn => $fn}, $class;

    open my $fh, '<', $fn
        or die "Error opening $fn for reading: $!\n";

    # read magic bytes and reset filehandle
    my $old_layers = join '', map {":$_"} PerlIO::get_layers($fh);
    binmode($fh);
    read( $fh, my $magic_start, 3 );
    read( $fh, my $magic_end,   1 );
    binmode($fh,$old_layers);
    seek($fh,0,0);

    # detect compression (must be BGZIP)
    if ($magic_start eq MAGIC_GZIP) {
        die "Gzip files must be compressed with bgzip for non-sequential access"
            if (ord($magic_end) != 0x04);
        # TODO: find way to test missing library during unit testing
        require Compress::BGZF::Reader
            or die "Compress::BGZF::Reader not installed, no gzip support";
        close $fh;
        $fh = Compress::BGZF::Reader->new_filehandle($fn);
    }

    $self->{with_description} = 1 if ($args{with_description});

    $self->{fh} = $fh;

    $self->_load_faidx;

    return $self;

}

sub write_index {

    my ($self, $fn_idx) = @_;

    $fn_idx //= $self->{fn} . '.fai';

    if (-e $fn_idx) {
        warn "Index file exists, won't overwrite\n";
        return 0;
    }

    open my $idx, '>:raw', $fn_idx
        or die "Error opening $fn_idx for writing: $!\n";

    my $fh = $self->{fh};

    my $curr_id;
    my $curr_len;
    my $curr_bases;
    my $curr_bytes;
    my $curr_offset;
    my $bl_mismatch = 0;
    my $ll_mismatch = 0;
    while (my $line =  <$fh>) {
        if ($line =~ /^>(\S+)/) {
            if (defined $curr_id) {
                say {$idx} join "\t",
                    $curr_id,
                    $curr_len,
                    $curr_offset,
                    $curr_bases,
                    $curr_bytes,
                ;
            }
            $curr_id = $1;
            $curr_offset = tell $fh;
            $curr_len = 0;
            $curr_bases = undef;
            $curr_bytes = undef;
            $bl_mismatch = 0;
            $ll_mismatch = 0;
        }
        elsif ($line =~ /^[A-Za-z\-]+(\r?\n?)$/) {
            if ($bl_mismatch) {
                close $idx;
                unlink $fn_idx;
                die "Base length mismatch\n";
            }
            if ($ll_mismatch) {
                close $idx;
                unlink $fn_idx;
                die "Line length mismatch\n";
            }
            my $ll  = length $line;
            my $bl  = $ll - length $1;
            $curr_len += $bl;
            $curr_bases //= $bl;
            $curr_bytes //= $ll;
            $bl_mismatch = 1 if ($bl != $curr_bases);
            $ll_mismatch = 1 if ($ll != $curr_bytes);
        }
        elsif ($line =~ /\S/) {
            close $idx;
            unlink $fn_idx;
            die "Unexpected content: $line";
        }

    }

    # write remaining index
    # uncoverable branch false
    if (defined $curr_id) {
        say {$idx} join "\t",
            $curr_id,
            $curr_len,
            $curr_offset,
            $curr_bases,
            $curr_bytes,
        ;
    }

    close $idx;

    #reset filehandle
    seek $self->{fh}, 0, 0;

    return 1;

}


sub _load_faidx {

    my ($self) = @_;

    my $fn_idx  = $self->{fn} . '.fai';
    $self->write_index if (! -e $fn_idx);

    # make sure FASTA file is not newer than its index; otherwise index file
    # may be out of date and cause silent corruption downstream
    my $mtime_fas = stat($self->{fn})->mtime;
    my $mtime_idx = stat($fn_idx)->mtime;
    if ($mtime_fas > $mtime_idx) {
        die "Index file exists but is older than FASTA file and probably out"
            . " of date. Refresh or remove the existing index before"
            . " proceeding.\n"
    }

    my @ids;

    open my $in, '<', $fn_idx
        or die "Error opening index file: $!\n";
    while (my $line = <$in>) {
        chomp $line;
        my ($name, $len, $offset, $bases_per_line, $bytes_per_line)
            = split "\t", $line;
        die "ERROR: $fn_idx contains duplicate entries"
            if ( defined $self->{idx}->{$name} );
        my $eol = $bytes_per_line - $bases_per_line;
        $self->{idx}->{$name} = [$len, $offset, $bases_per_line, $eol];
        push @ids, $name;
    }
    close $in;
    $self->{ids} = \@ids;

    return;

}

sub ids {

    my ($self) = @_;
    return @{ $self->{ids} };

}

sub length {

    my ($self, $id) = @_;
    return $self->{idx}->{$id}->[0];

}

sub fetch_seq {

    my ($self, $id, $start_bp, $end_bp) = @_;

    return undef if (! defined $self->{idx}->{$id});
    my ($len, $off, $bpl, $eol) = @{ $self->{idx}->{$id} };

    my $fh = $self->{fh};

    $start_bp = $start_bp // 1;
    $end_bp   = $end_bp   // $len;
    --$start_bp; #make 0-based
    die "Coordinates out of bounds" if ($start_bp < 0 || $end_bp > $len);
    my $l = $end_bp - $start_bp;

    seek $fh, $off + $start_bp + int(($start_bp)/$bpl)*$eol, 0;
    read($fh, my $seq, $l + int(($l+$start_bp%$bpl)/$bpl)*$eol);
    $seq =~ s/\s//g;

    ++$start_bp;

    my $desc;
    if ($start_bp > 1 || $end_bp < $len) {
        $desc = "($start_bp-$end_bp)";
    }

    # backtrack to find defline, if asked
    if ($self->{with_description}) {

        my $string = '';
        my $p = $off;
        my $chr;

        BACK:
        while ($p > 0) {
            --$p;
            seek $fh, $p, 0;
            read($fh, $chr, 1)
                or die "Error during backtrack read: $@\n";
            # stop at record separator
            if ($chr eq '>') {
                last BACK if ($p == 0);
                my $prev;
                seek $fh, $p-1, 0;
                read($fh, $prev, 1);
                last BACK if ($prev =~ /[\r\n]/);
            }

            $string = $chr . $string;
        }

        die "Backtracking inexplicably failed\n"
            if (! CORE::length $string);
        if ($string =~ /^(\S+)(?:\s+(.*))?$/) {
            my $nid  = $1;
            $desc = join ' ', grep {defined $_} $2, $desc;
            die "ID mismatch ($id vs $nid)!\n"
                if ($nid ne $id);
        }

    }

    return BioX::Seq->new($seq, $id, $desc);

}


1;


__END__