view release on metacpan or  search on metacpan

lib/Bio/DB/HTS.pm  view on Meta::CPAN

    return
      Bio::SeqFeature::HTSCoverage->new( -display_name => "$seqid coverage",
                                      -seq_id       => $seqid,
                                      -start        => $start,
                                      -end          => $end,
                                      -strand       => 0,
                                      -type         => "coverage:$bins",
                                      -class        => "coverage:$bins",
                                      -attributes => { coverage => [$coverage] }
      );
} ## end sub _coverage

sub _segment_search {
    my $self = shift;
    my $name = shift;

    my $targets = $self->_cache_targets;
    return $self->segment($name) if $targets->{$name};

    if ( my $regexp = $self->_glob_match($name) ) {
        my @results = grep { /^$regexp$/i } keys %$targets;
        return map { $self->segment($_) } @results;
    }

    return;
}

sub hts_index {
    my $self = shift;
    if ( defined $self->{hts_idx} ) {
        return $self->{hts_idx};
    }
    $self->{hts_idx} = Bio::DB::HTSfile->index($self);
    return $self->{hts_idx};
}

sub _features_fh {
    my $self = shift;
    my ( $seqid, $start, $end, $filter ) = @_;

    my $result = open my $fh, "-|";
    if ( !$result ) {    # in child
        $self->_filter_features( $seqid, $start, $end, $filter, 'do_fh' )
          ;              # will print TAM to stdout
        exit 0;
    }
    return $fh;

}

sub tam_fh {
    my $self = shift;
    return $self->features( -fh => 1 );
}

sub max_pileup_cnt {
    my $self = shift;
    return Bio::DB::HTSfile->max_pileup_cnt(@_);
}

# return a fragment of code that will be placed in the eval "" filter
# to eliminate alignments that don't match by name
sub _filter_by_name {
    my $self = shift;
    my $name = shift;

    my $frag = "my \$name=\$a->qname; defined \$name or return; ";

    if ( my $regexp = $self->_glob_match($name) ) {
        $frag .= "return unless \$name =~ /^$regexp\$/i;\n";
    }
    else {
        $frag .= "return unless lc \$name eq '$name';\n";
    }
}

# return a fragment of code that will be placed in the eval "" filter
# to eliminate alignments that don't match by attribute
sub _filter_by_attribute {
    my $self       = shift;
    my $attributes = shift;
    my $result;
    for my $tag ( keys %$attributes ) {
        $result .= "my \$value = lc \$a->get_tag_values('$tag');\n";
        $result .= "return unless defined \$value;\n";
        my @comps =
          ref $attributes->{$tag} eq 'ARRAY' ? @{ $attributes->{$tag} } :
          $attributes->{$tag};
        my @matches;
        for my $c (@comps) {
            if ( $c =~ /^[+-]?[\deE.]+$/ ) {    # numeric-looking argument
                push @matches, "CORE::length \$value && \$value == $c";
            }
            elsif ( my $regexp = $self->_glob_match($c) ) {
                push @matches, "\$value =~ /^$regexp\$/i";
            }
            else {
                push @matches, "\$value eq lc '$c'";
            }
        }
        $result .= "return unless " . join( ' OR ', @matches ) . ";\n";
    }
    return $result;
} ## end sub _filter_by_attribute

# turn a glob expression into a regexp
sub _glob_match {
    my $self = shift;
    my $term = shift;
    return unless $term =~ /(?:^|[^\\])[*?]/;
    $term =~ s/(^|[^\\])([+\[\]^{}\$|\(\).])/$1\\$2/g;
    $term =~ s/(^|[^\\])\*/$1.*/g;
    $term =~ s/(^|[^\\])\?/$1./g;
    return $term;
}

package Bio::DB::HTS::Fai;

$Bio::DB::HTS::Fai::VERSION = '3.01';

sub open { shift->load(@_) }

sub seq {
    my $self = shift;
    my ( $seqid, $start, $end ) = @_;
    my $region = $seqid;
    $region .= ":$start" if defined $start;
    $region .= "-$end"   if defined $end;
    return $self->fetch($region);
}

package Bio::SeqFeature::HTSCoverage;

use base 'Bio::SeqFeature::Lite';

$Bio::SeqFeature::HTSCoverage::VERSION = '3.01';