File Coverage

blib/lib/Bio/JBrowse/Store/NCList/IntervalStore.pm

Criterion	Covered	Total	%
statement	19	125	15.2
branch	0	26	0.0
condition	0	10	0.0
subroutine	7	26	26.9
pod	6	12	50.0
total	32	199	16.0

package Bio::JBrowse::Store::NCList::IntervalStore;
BEGIN {
  $Bio::JBrowse::Store::NCList::IntervalStore::AUTHORITY = 'cpan:RBUELS';
}
{
  $Bio::JBrowse::Store::NCList::IntervalStore::VERSION = '0.1';
}
use strict;
use warnings;
use Carp;

use POSIX ();
use List::Util ();

use Bio::JBrowse::Store::NCList::LazyNCList ();


sub new {
    my ($class, $args) = @_;

    my $self = {
                store => $args->{store} || die,

                compress => $args->{compress},

                attrs => $args->{arrayRepr},
                urlTemplate => $args->{urlTemplate} || ("lf-{Chunk}"
                                                        . $args->{store}->ext),
                nclist => $args->{nclist},
                minStart => $args->{minStart},
                maxEnd => $args->{maxEnd},
                loadedChunks => {}
               };

    if( defined $args->{nclist} ) {
        # we're already loaded
        $self->{lazyNCList} =
          Bio::JBrowse::Store::NCList::LazyNCList->importExisting(
              $self->{attrs},
              $args->{count},
              $args->{minStart},
              $args->{maxEnd},
              sub { $self->_loadChunk( @_ ); },
              $args->{nclist} );
    }

    bless $self, $class;

    return $self;
}

sub _loadChunk {
    my ($self, $chunkId) = @_;
    my $chunk = $self->{loadedChunks}->{$chunkId};
    if (defined($chunk)) {
        return $chunk;
    } else {
        (my $path = $self->{urlTemplate}) =~ s/\{Chunk\}/$chunkId/g;
        $chunk = $self->{store}->get( $path );
        # TODO limit the number of chunks that we keep in memory
        $self->{loadedChunks}->{$chunkId} = $chunk;
        return $chunk;
    }
}


sub startLoad {
    my ( $self, $measure, $chunkBytes ) = @_;

    if (defined($self->{nclist})) {
        confess "loading into an already-loaded Bio::JBrowse::Store::NCList::IntervalStore";
    } else {
        my $lazyClass = $self->{attrs}->getClass( { start => 1, end => 2, chunk => 3 } );
        $self->{lazyClass} = $lazyClass->{index};

        # add a new class for "fake" features
        my $makeLazy = sub {
            my ($start, $end, $chunkId) = @_;
            my $f = { start => $start, end => $end, chunk => $chunkId };
            return [ $lazyClass->{index}, map { $f->{$_} } @{$lazyClass->{attributes}} ];
        };
        my $output = sub {
            my ($toStore, $chunkId) = @_;
            (my $path = $self->{urlTemplate}) =~ s/\{Chunk\}/$chunkId/ig;
            $self->{store}->put($path, $toStore);
        };
        $self->{lazyNCList} =
          Bio::JBrowse::Store::NCList::LazyNCList->new( $self->{attrs},
                          $self->{lazyClass},
                          $makeLazy,
                          sub { $self->_loadChunk( @_); },
                          $measure,
                          $output,
                          $chunkBytes);
    }
}


sub addSorted {
    my ($self, $feat) = @_;
    $self->{lazyNCList}->addSorted($feat);
}


sub finishLoad {
    my ($self) = @_;
    $self->{lazyNCList}->finish();
    $self->{nclist} = $self->lazyNCList->topLevelList();

    my $trackData = {
        featureCount => $self->count,
        intervals => $self->descriptor,
        histograms => $self->writeHistograms,
        formatVersion => 1
        };

    $self->store->put( 'trackData.'.( $self->{compress} ? 'jsonz' : 'json' ), $trackData);

}

sub writeHistograms {
    my ( $self ) = @_;
    #this series of numbers is used in JBrowse for zoom level relationships
    my @multiples = (1, 2, 5, 10, 20, 50, 100, 200, 500,
                     1000, 2000, 5000, 10_000, 20_000, 50_000,
                     100_000, 200_000, 500_000, 1_000_000);
    my $histChunkSize = 10_000;

    my $attrs = $self->{attrs};
    my $getStart = $attrs->makeGetter("start");
    my $getEnd = $attrs->makeGetter("end");

    my $jsonStore = $self->store;
    my $refEnd = $self->lazyNCList->maxEnd || 0;
    my $featureCount = $self->count;

    # $histBinThresh is the approximate the number of bases per
    # histogram bin at the zoom level where FeatureTrack.js switches
    # to the histogram view by default
    my $histBinThresh = $featureCount ? ($refEnd * 2.5) / $featureCount : 999_999_999_999;
    my $histBinBases  = ( List::Util::first { $_ > $histBinThresh } @multiples ) || $multiples[-1];

    # initialize histogram arrays to all zeroes
    my @histograms;
    for (my $i = 0; $i < @multiples; $i++) {
        my $binBases = $histBinBases * $multiples[$i];
        $histograms[$i] = [(0) x POSIX::ceil($refEnd / $binBases)];
        # somewhat arbitrarily cut off the histograms at 100 bins
        last if $binBases * 100 > $refEnd;
    }

    my $processFeat = sub {
        my ($feature) = @_;
        my $curHist;
        my $start = List::Util::max(0, List::Util::min($getStart->($feature), $refEnd));
        my $end = List::Util::min($getEnd->($feature), $refEnd);
        return if ($end < 0);

        for (my $i = 0; $i <= $#multiples; $i++) {
            my $binBases = $histBinBases * $multiples[$i];
            $curHist = $histograms[$i];
            last unless defined($curHist);

            my $firstBin = int($start / $binBases);
            my $lastBin = int($end / $binBases);
            for (my $bin = $firstBin; $bin <= $lastBin; $bin++) {
                $curHist->[$bin] += 1;
            }
        }
    };

    $self->overlapCallback($self->lazyNCList->minStart,
                           $self->lazyNCList->maxEnd,
                           $processFeat);

    # find multiple of base hist bin size that's just over $histBinThresh
    my $i;
    for ($i = 1; $i <= $#multiples; $i++) {
        last if ($histBinBases * $multiples[$i]) > $histBinThresh;
    }

    my @histogramMeta;
    my @histStats;
    for (my $j = $i - 1; $j <= $#multiples; $j += 1) {
        my $curHist = $histograms[$j];
        last unless defined($curHist);
        my $histBases = $histBinBases * $multiples[$j];

        my $chunks = $self->chunkArray($curHist, $histChunkSize);
        for (my $k = 0; $k <= $#{$chunks}; $k++) {
            $jsonStore->put("hist-$histBases-$k" . $jsonStore->ext,
                            $chunks->[$k]);
        }
        push @histogramMeta,
            {
                basesPerBin => $histBases,
                arrayParams => {
                    length => $#{$curHist} + 1,
                    urlTemplate => "hist-$histBases-{Chunk}" . $jsonStore->ext,
                    chunkSize => $histChunkSize
                }
            };
        push @histStats,
            {
                'basesPerBin' => $histBases,
                'max'  => @$curHist ? List::Util::max( @$curHist ) : undef,
                'mean' => @$curHist ? ( List::Util::sum( @$curHist ) / @$curHist ) : undef,
            };
    }

    return { meta => \@histogramMeta,
             stats => \@histStats };
}


sub chunkArray {
    my ( $self, $bigArray, $chunkSize) = @_;

    my @result;
    for (my $start = 0; $start <= $#{$bigArray}; $start += $chunkSize) {
        my $lastIndex = $start + $chunkSize;
        $lastIndex = $#{$bigArray} if $lastIndex > $#{$bigArray};

        push @result, [@{$bigArray}[$start..$lastIndex]];
    }
    return \@result;
}



sub overlapCallback {
    my ($self, $start, $end, $cb) = @_;
    $self->lazyNCList->overlapCallback($start, $end, $cb);
}


sub lazyNCList   { shift->{lazyNCList}        }
sub count        { shift->{lazyNCList}->count }
sub hasIntervals { shift->count > 0           }
sub store        { shift->{store}             }


sub descriptor {
    my ( $self ) = @_;
    return {
            lazyClass => $self->{lazyClass},
            nclist => $self->{nclist},
            classes => $self->{attrs}->descriptor,
            urlTemplate => $self->{urlTemplate},
            count => $self->count,
            minStart => $self->lazyNCList->minStart,
            maxEnd => $self->lazyNCList->maxEnd
           };
}

1;

__END__