view release on metacpan or  search on metacpan

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

    my $index      = $self->bam_index;
    my $coverage   = $index->coverage($self->{bam},
				      $id,$s,$e,
				      $bins);

    return Bio::SeqFeature::Coverage->new(
	-display_name => "$seqid coverage",
	-seq_id       => $seqid,
	-start        => $start,
	-end          => $end,
	-strand       => 0,
	-type         => "coverage:$bins",
	-class        => "coverage:$bins",
	-attributes   => { coverage => [$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 bam_index {
    my $self = shift;
    return $self->{bai} ||= Bio::DB::Bam->index($self->{bam_path},$self->autoindex);
}

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::Bam->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;
}

# 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::Sam::Fai;

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::Coverage;

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

sub coverage {
    my $self       = shift;
    my ($coverage) = $self->get_tag_values('coverage');
    return wantarray ? @$coverage : $coverage;