Bio::DB::Sam -- Read SAM/BAM database files
use Bio::DB::Sam;
# high level API
my $sam = Bio::DB::Sam->new(-bam =>"data/ex1.bam",
-fasta=>"data/ex1.fa",
);
my @targets = $sam->seq_ids;
my @alignments = $sam->get_features_by_location(-seq_id => 'seq2',
-start => 500,
-end => 800);
for my $a (@alignments) {
# where does the alignment start in the reference sequence
my $seqid = $a->seq_id;
my $start = $a->start;
my $end = $a->end;
my $strand = $a->strand;
my $cigar = $a->cigar_str;
my $paired = $a->get_tag_values('PAIRED');
# where does the alignment start in the query sequence
my $query_start = $a->query->start;
my $query_end = $a->query->end;
my $ref_dna = $a->dna; # reference sequence bases
my $query_dna = $a->query->dna; # query sequence bases
my @scores = $a->qscore; # per-base quality scores
my $match_qual= $a->qual; # quality of the match
}
my @pairs = $sam->get_features_by_location(-type => 'read_pair',
-seq_id => 'seq2',
-start => 500,
-end => 800);
for my $pair (@pairs) {
my $length = $pair->length; # insert length
my ($first_mate,$second_mate) = $pair->get_SeqFeatures;
my $f_start = $first_mate->start;
my $s_start = $second_mate->start;
}
# low level API
my $bam = Bio::DB::Bam->open('/path/to/bamfile');
my $header = $bam->header;
my $target_count = $header->n_targets;
my $target_names = $header->target_name;
while (my $align = $bam->read1) {
my $seqid = $target_names->[$align->tid];
my $start = $align->pos+1;
my $end = $align->calend;
my $cigar = $align->cigar_str;
}
my $index = Bio::DB::Bam->index_open('/path/to/bamfile');
my $index = Bio::DB::Bam->index_open_in_safewd('/path/to/bamfile');
my $callback = sub {
my $alignment = shift;
my $start = $alignment->start;
my $end = $alignment->end;
my $seqid = $target_names->[$alignment->tid];
print $alignment->qname," aligns to $seqid:$start..$end\n";
}
my $header = $index->header;
$index->fetch($bam,$header->parse_region('seq2'),$callback);
This module provides a Perl interface to the libbam library for indexed and
unindexed SAM/BAM sequence alignment databases. It provides support for
retrieving information on individual alignments, read pairs, and alignment
coverage information across large regions. It also provides callback
functionality for calling SNPs and performing other base-by-base functions.
Most operations are compatible with the BioPerl Bio::SeqFeatureI interface,
allowing BAM files to be used as a backend to the GBrowse genome browser
application (gmod.sourceforge.net).
The high-level API provides a BioPerl-compatible interface to indexed BAM files.
The BAM database is treated as a collection of Bio::SeqFeatureI features, and
can be searched for features by name, location, type and combinations of
feature tags such as whether the alignment is part of a mate-pair.
When opening a BAM database using the high-level API, you provide the pathnames
of two files: the FASTA file that contains the reference genome sequence, and
the BAM file that contains the query sequences and their alignments. If either
of the two files needs to be indexed, the indexing will happen automatically.
You can then query the database for alignment features by combinations of
name, position, type, and feature tag.
The high-level API provides access to up to four feature "types":
* "match": The "raw" unpaired alignment between a read and the
reference sequence.
* "read_pair": Paired alignments; a single composite
feature that contains two subfeatures for the alignments of each
of the mates in a mate pair.
* "coverage": A feature that spans a region of interest that contains
numeric information on the coverage of reads across the region.
* "region": A way of retrieving information about the reference
sequence. Searching for features of type "region" will return a
list of chromosomes or contigs in the reference sequence, rather
than read alignments.
* "chromosome": A synonym for "region".
Features can be en masse in a single call, retrieved in a
memory-efficient streaming basis using an iterator, or interrogated using a
filehandle that return a series of TAM-format lines.
SAM alignment flags can be retrieved using BioPerl's feature
"tag" mechanism. For example, to interrogate the FIRST_MATE flag,
one fetches the "FIRST_MATE" tag:
warn "aye aye captain!" if $alignment->get_tag_values('FIRST_MATE');
The Bio::SeqFeatureI interface has been extended to retrieve all flags as a
compact human-readable string, and to return the CIGAR alignment in a variety
of formats.
Split alignments, such as reads that cover introns, are dealt with in one
of two ways. The default is to leave split alignments alone: they can be
detected by one or more "N" operations in the CIGAR string.
Optionally, you can choose to have the API split these alignments across two
or more subfeatures; the CIGAR strings of these split alignments will be
adjusted accordingly.
Interface to the pileup routines The API provides you with access to the
samtools "pileup" API. This gives you the ability to write a
callback that will be invoked on every column of the alignment for the purpose
of calculating coverage, quality score metrics, or SNP calling.
Access to the reference sequence When you create the Bio::DB::Sam object,
you can pass the path to a FASTA file containing the reference sequence.
Alternatively, you may pass an object that knows how to retrieve DNA sequences
across a range via the
seq() of
fetch_seq() methods, as
described under
new().
If the SAM/BAM file has MD tags, then these tags will be used to reconstruct the
reference sequence when necessary, in which case you can completely omit the
-fasta argument. Note that not all SAM/BAM files have MD tags, and those that
do may not use them correctly due to the newness of this part of the SAM spec.
You may wish to populate these tags using samtools' "calmd" command.
If the -fasta argument is omitted and no MD tags are present, then the reference
sequence will be returned as 'N'.
The
main object classes that you will be dealing with in the high-level
API are as follows:
* Bio::DB::Sam -- A collection of alignments and reference sequences.
* Bio::DB::Bam::Alignment -- The alignment between a query and the reference.
* Bio::DB::Bam::Query -- An object corresponding to the query sequence in
which both (+) and (-) strand alignments are
shown in the reference (+) strand.
* Bio::DB::Bam::Target -- An interface to the query sequence in which
(-) strand alignments are shown in reverse
complement
You may encounter other classes as well. These include:
* Bio::DB::Sam::Segment -- This corresponds to a region on the reference
sequence.
* Bio::DB::Sam::Constants -- This defines CIGAR symbol constants and flags.
* Bio::DB::Bam::AlignWrapper -- An alignment helper object that adds split
alignment functionality. See Bio::DB::Bam::Alignment
for the documentation on using it.
* Bio::DB::Bam::ReadIterator -- An iterator that mediates the one-feature-at-a-time
retrieval mechanism.
* Bio::DB::Bam::FetchIterator -- Another iterator for feature-at-a-time retrieval.
The low-level API closely mirrors that of the libbam library. It provides the
ability to open TAM and BAM files, read and write to them, build indexes, and
perform searches across them. There is less overhead to using the API because
there is very little Perl memory management, but the functions are less
convenient to use. Some operations, such as writing BAM files, are only
available through the low-level API.
The classes you will be interacting with in the low-level API are as follows:
* Bio::DB::Tam -- Methods that read and write TAM (text SAM) files.
* Bio::DB::Bam -- Methods that read and write BAM (binary SAM) files.
* Bio::DB::Bam::Header -- Methods for manipulating the BAM file header.
* Bio::DB::Bam::Index -- Methods for retrieving data from indexed BAM files.
* Bio::DB::Bam::Alignment -- Methods for manipulating alignment data.
* Bio::DB::Bam::Pileup -- Methods for manipulating the pileup data structure.
* Bio::DB::Sam::Fai -- Methods for creating and reading from indexed Fasta
files.
=head1 METHODS
We cover the high-level API first. The high-level API code can be found in the
files Bio/DB/Sam.pm, Bio/DB/Sam/*.pm, and Bio/DB/Bam/*.pm.
- $sam = Bio::DB::Sam->new(%options)
- The Bio::DB::Sam object combines a Fasta file of the
reference sequences with a BAM file to allow for convenient retrieval of
human-readable sequence IDs and reference sequences. The new()
constructor accepts a -name=>value style list of options as follows:
Option Description
------ -------------
-bam Path to the BAM file that contains the
alignments (required). When using samtools 0.1.6
or higher, an http: or ftp: URL is accepted.
-fasta Path to the Fasta file that contains
the reference sequences (optional). Alternatively,
you may pass any object that supports a seq()
or fetch_seq() method and takes the three arguments
($seq_id,$start,$end).
-expand_flags A boolean value. If true then the standard
alignment flags will be broken out as
individual tags such as 'M_UNMAPPED' (default
false).
-split_splices A boolean value. If true, then alignments that
are split across splices will be broken out
into a single alignment containing two sub-
alignments (default false).
-split The same as -split_splices.
-force_refseq Always use the reference sequence file to derive the
reference sequence, even when the sequence can be
derived from the MD tag. This is slower, but safer
when working with BAM files derived from buggy aligners
or when the reference contains non-canonical (modified)
bases.
-autoindex Create a BAM index file if one does not exist
or the current one has a modification date
earlier than the BAM file.
An example of a typical new() constructor invocation is:
$sam = Bio::DB::Sam->new(-fasta => '/home/projects/genomes/hu17.fa',
-bam => '/home/projects/alignments/ej88.bam',
-expand_flags => 1,
-split_splices => 1);
If the -fasta argument is present, then you will be able to use the
interface to fetch the reference sequence's bases. Otherwise, calls that
return the reference sequence will return sequences consisting entirely of
"N".
-expand_flags option, if true, has the effect of turning each of the
standard SAM flags into a separately retrievable tag in the
Bio::SeqFeatureI interface. Otherwise, the standard flags will be
concatenated in easily parseable form as a tag named "FLAGS".
See get_all_tags() and get_tag_values() for more
information.
Any two-letter extension flags, such as H0 or H1, will always appear as
separate tags regardless of the setting.
-split_splices has the effect of breaking up alignments that contain
an "N" operation into subparts for more convenient manipulation.
For example, if you have both paired reads and spliced alignments in the
BAM file, the following code shows the subpart relationships:
$pair = $sam->get_feature_by_name('E113:01:01:23');
@mates = $pair->get_SeqFeatures;
@mate1_parts = $mates[0]->get_SeqFeatures;
@mate2_parts = $mates[1]->get_SeqFeatures;
Because there is some overhead to splitting up the spliced alignments, this
option is false by default.
Remote access to BAM files located on an HTTP or FTP server is
possible when using the Samtools library version 0.1.6 or higher. Simply
replace the path to the BAM file with the appropriate URL. Note that
incorrect URLs may lead to a core dump.
It is not currently possible to refer to a remote FASTA file. These will
have to be downloaded locally and indexed before using.
- $flag = $sam->expand_flags([$new_value])
- Get or set the expand_flags option. This can be done after
object creation and will have an immediate effect on all alignments
fetched from the BAM file.
- $flag = $sam->split_splices([$new_value])
- Get or set the split_splices option. This can be done after
object creation and will affect all alignments fetched from the BAM file
subsequently.
- $header = $sam->header
- Return the Bio::DB::Bam::Header object associated with the
BAM file. You can manipulate the header using the low-level API.
- $bam_path = $sam->bam_path
- Return the path of the bam file used to create the sam
object. This makes the sam object more portable.
- $bam = $sam->bam
- Returns the low-level Bio::DB::Bam object associated with
the opened file.
- $fai = $sam->fai
- Returns the Bio::DB::Sam::Fai object associated with the
Fasta file. You can then manipuate this object with the low-level API.
The index will be built automatically for you if it does not already
exist. If index building is necessarily, the process will need
write privileges to the same directory in which the Fasta file
resides.> If the process does not have write permission, then the call
will fail. Unfortunately, the BAM library does not do great error recovery
for this condition, and you may experience a core dump. This is not
trappable via an eval {}.
- $bai = $sam->bam_index
- Return the Bio::DB::Bam::Index object associated with the
BAM file.
The BAM file index will be built automatically for you if it does
not already exist. In addition, if the BAM file is not already
sorted by chromosome and coordinate, it will be sorted automatically, an
operation that consumes significant time and disk space. The current
process must have write permission to the directory in which the BAM file
resides in order for this to work.> In case of a permissions problem,
the Perl library will catch the error and die. You can trap it with an
eval {}.
- $sam->clone
- Bio::DB::SAM objects are not stable across fork()
operations. If you fork, you must call clone() either in the parent
or the child process before attempting to call any methods.
The Bio::DB::Sam object provides the following methods for getting information
about the reference sequence(s) contained in the associated Fasta file.
- @seq_ids = $sam->seq_ids
- Returns an unsorted list of the IDs of the reference
sequences (known elsewhere in this document as seq_ids). This is the same
as the identifier following the ">" sign in the Fasta file
(e.g. "chr1").
- $num_targets = $sam->n_targets
- Return the number of reference sequences.
- $length = $sam->length('seqid')
- Returns the length of the reference sequence named
"seqid".
- $seq_id = $sam->target_name($tid)
- Translates a numeric target ID (TID) returned by the
low-level API into a seq_id used by the high-level API.
- $length = $sam->target_len($tid)
- Translates a numeric target ID (TID) from the low-level API
to a sequence length.
- $dna = $sam->seq($seqid,$start,$end)
- Returns the DNA across the region from start to end on
reference seqid. Note that this is a string, not a Bio::PrimarySeq object.
If no -fasta path was passed when the sam object was created, then you
will receive a series of N nucleotides of the requested length.
Bio::DB::Sam::Segment objects refer regions on the reference sequence. They can
be used to retrieve the sequence of the reference, as well as alignments that
overlap with the region.
- $segment = $sam->segment($seqid,$start,$end);
- $segment =
$sam->segment(-seq_id=>'chr1',-start=>5000,-end=>6000);
- Segments are created using the
Bio:DB::Sam->segment() method. It can be called using one to
three positional arguments corresponding to the seq_id of the reference
sequence, and optionally the start and end positions of a subregion on the
sequence. If the start and/or end are undefined, they will be replaced
with the beginning and end of the sequence respectively.
Alternatively, you may call segment() with named -seq_id, -start and
-end arguments.
All coordinates are 1-based.
- $seqid = $segment->seq_id
- Return the segment's sequence ID.
- $start = $segment->start
- Return the segment's start position.
- $end = $segment->end
- Return the segment's end position.
- $strand = $segment->strand
- Return the strand of the segment (always 0).
- $length = $segment->length
- Return the length of the segment.
- $dna = $segment->dna
- Return the DNA string for the reference sequence under this
segment.
- $seq = $segment->seq
- Return a Bio::PrimarySeq object corresponding to the
sequence of the reference under this segment. You can get the actual DNA
string in this redundant-looking way:
$dna = $segment->seq->seq
The advantage of working with a Bio::PrimarySeq object is that you can
perform operations on it, including taking its reverse complement and
subsequences.
- @alignments = $segment->features(%args)
- Return alignments that overlap the segment in the
associated BAM file. The optional %args list allows you to filter features
by name, tag or other attributes. See the documentation of the
Bio::DB::Sam-> features() method for the full list of options.
Here are some typical examples:
# get all the overlapping alignments
@all_alignments = $segment->features;
# get an iterator across the alignments
my $iterator = $segment->features(-iterator=>1);
while (my $align = $iterator->next_seq) { do something }
# get a TAM filehandle across the alignments
my $fh = $segment->features(-fh=>1);
while (<$fh>) { print }
# get only the alignments with unmapped mates
my @unmapped = $segment->features(-flags=>{M_UNMAPPED=>1});
# get coverage across this region
my ($coverage) = $segment->features('coverage');
my @data_points = $coverage->coverage;
# grep through features using a coderef
my @reverse_alignments = $segment->features(
-filter => sub {
my $a = shift;
return $a->strand < 0;
});
- $tag = $segment->primary_tag
- $tag = $segment->source_tag
- Return the strings "region" and
"sam/bam" respectively. These methods allow the segment to be
passed to BioPerl methods that expect Bio::SeqFeatureI objects.
- $segment->name, $segment->display_name,
$segment->get_SeqFeatures, $segment->get_tag_values
- These methods are provided for Bio::SeqFeatureI
compatibility and don't do anything of interest.
The
features() method is an all-purpose tool for retrieving alignment
information from the SAM/BAM database. In addition, the methods
get_features_by_name(),
get_features_by_location() and others
provide convenient shortcuts to
features().
These methods either return a list of features, an iterator across a list of
features, or a filehandle opened on a pseudo-TAM file.
- @features = $sam->features(%options)
- $iterator =
$sam->features(-iterator=>1,%more_options)
- $filehandle =
$sam->features(-fh=>1,%more_options)
- @features = $sam->features('type1','type2'...)
- This is the all-purpose interface for fetching alignments
and other types of features from the database. Arguments are a
-name=>value option list selected from the following list of options:
Option Description
------ -------------
-type Filter on features of a given type. You may provide
either a scalar typename, or a reference to an
array of desired feature types. Valid types are
"match", "read_pair", "coverage" and "chromosome."
See below for a full explanation of feature types.
-name Filter on reads with the designated name. Note that
this can be a slow operation unless accompanied by
the feature location as well.
-seq_id Filter on features that align to seq_id between start
-start and end. -start and -end must be used in conjunction
-end with -seq_id. If -start and/or -end are absent, they
will default to 1 and the end of the reference
sequence, respectively.
-flags Filter features that match a list of one or more
flags. See below for the format.
-attributes The same as -flags, for compatibility with other
-tags APIs.
-filter Filter on features with a coderef. The coderef will
receive a single argument consisting of the feature
and should return true to keep the feature, or false
to discard it.
-iterator Instead of returning a list of features, return an
iterator across the results. To retrieve the results,
call the iterator's next_seq() method repeatedly
until it returns undef to indicate that no more
matching features remain.
-fh Instead of returning a list of features, return a
filehandle. Read from the filehandle to retrieve
each of the results in TAM format, one alignment
per line read. This only works for features of type
"match."
The high-level API introduces the concept of a feature
"type" in order to provide several convenience functions.
You specify types by using the optional -type argument. The
following types are currently supported:
match. The "match" type corresponds to the unprocessed SAM
alignment. It will retrieve single reads, either mapped or unmapped. Each
match feature's primary_tag() method will return the string
"match." The features returned by this call are of type
Bio::DB::Bam::AlignWrapper.
read_pair. The "paired_end" type causes the sam interface
to find and merge together mate pairs. Fetching this type of feature will
yield a series of Bio::SeqFeatureI objects, each as long as the total
distance on the reference sequence spanned by the mate pairs. The
top-level feature is of type Bio::SeqFeature::Lite; it contains two
Bio::DB::Bam::AlignWrapper subparts.
Call get_SeqFeatures() to get the two individual reads. Example:
my @pairs = $sam->features(-type=>'read_pair');
my $p = $pairs[0];
my $i_length = $p->length;
my @ends = $p->get_SeqFeatures;
my $left = $ends[0]->start;
my $right = $ends[1]->end;
coverage. The "coverage" type causes the sam interface to
calculate coverage across the designated region. It only works properly if
accompanied by the desired location of the coverage graph; -seq_id is a
mandatory argument for coverage calculation, and -start and -end are
optional. The call will return a single Bio::SeqFeatureI object whose
primary_tag() is "coverage." To recover the coverage
data, call the object's coverage() method to obtain an array (list
context) or arrayref (scalar context) of coverage counts across the region
of interest:
my ($coverage) = $sam->features(-type=>'coverage',-seq_id=>'seq1');
my @data = $coverage->coverage;
my $total;
for (@data) { $total += $_ }
my $average_coverage = $total/@data;
By default the coverage graph will be at the base pair level. So for a
region 5000 bp wide, coverage() will return an array or arrayref
with exactly 5000 elements. However, you also have the option of
calculating the coverage across larger bins. Simply append the number of
intervals you are interested to the "coverage" typename. For
example, fetching "coverage:500" will return a feature whose
coverage() method will return the coverage across 500 intervals.
chromosome or region. The "chromosome" or
"region" type are interchangeable. They ask the sam interface to
construct Bio::DB::Sam::Segment representing the reference sequences.
These two calls give similar results:
my $segment = $sam->segment('seq2',1=>500);
my ($seg) = $sam->features(-type=>'chromosome',
-seq_id=>'seq2',-start=>1,-end=>500);
Due to an unresolved bug, you cannot fetch chromosome features in the same
call with matches and other feature types call. Specifically, this works
as expected:
my @chromosomes = $sam->features (-type=>'chromosome');
But this doesn't (as of 18 June 2009):
my @chromosomes_and_matches = $sam->features(-type=>['match','chromosome']);
If no -type argument is provided, then features() defaults to finding
features of type "match."
You may call features() with a plain list of strings (positional
arguments, not -type=>value arguments). This will be interpreted as a
list of feature types to return:
my ($coverage) = $sam->features('coverage')
For a description of the methods available in the features returned from
this call, please see Bio::SeqfeatureI and Bio::DB::Bam::Alignment.
You can filter "match" and "read_pair" features
by name, location and/or flags. The name and flag filters are not very
efficient. Unless they are combined with a location filter, they will
initiate an exhaustive search of the BAM database.
Name filters are case-insensitive, and allow you to use shell-style
"*" and "?" wildcards. Flag filters created with the
-flag, -attribute or -tag options have the following
syntax:
-flag => { FLAG_NAME_1 => ['list','of','possible','values'],
FLAG_NAME_2 => ['list','of','possible','values'],
...
}
The value of -flag is a hash reference in which the keys are flag
names and the values are array references containing lists of acceptable
values. The list of values are OR'd with each other, and the flag names
are AND'd with each other.
The -filter option provides a completely generic filtering interface.
Provide a reference to a subroutine. It will be called once for each
potential feature. Return true to keep the feature, or false to discard
it. Here is an example of how to find all matches whose alignment quality
scores are greater than 80.
@features = $sam->features(-filter=>sub {shift->qual > 80} );
By default, features() returns a list of all matching features. You
may instead request an iterator across the results list by passing
-iterator=>1. This will give you an object that has a single method,
next_seq():
my $high_qual = $sam->features(-filter => sub {shift->qual > 80},
-iterator=> 1 );
while (my $feature = $high_qual->next_seq) {
# do something with the alignment
}
Similarly, by passing a true value to the argument -fh, you can
obtain a filehandle to a virtual TAM file. This only works with the
"match" feature type:
my $high_qual = $sam->features(-filter => sub {shift->qual > 80},
-fh => 1 );
while (my $tam_line = <$high_qual>) {
chomp($tam_line);
# do something with it
}
- @features = $sam->get_features_by_name($name)
- Convenience method. The same as calling
$sam->features(-name=>$name);
- $feature = $sam->get_feature_by_name($name)
- Convenience method. The same as
($sam->features(-name=>$name))[0].
- @features =
$sam->get_features_by_location($seqid,$start,$end)
- Convenience method. The same as calling
$sam->features(-seq_id=>$seqid,-start=>$start,-end=>$end).
- @features = $sam->get_features_by_flag(%flags)
- Convenience method. The same as calling
$sam->features(-flags=>\%flags). This method is also called
get_features_by_attribute() and get_features_by_tag().
Example:
@features = $sam->get_features_by_flag(H0=>1)
- $feature = $sam->get_feature_by_id($id)
- The high-level API assigns each feature a unique ID
composed of its read name, position and strand and returns it when you
call the feature's primary_id() method. Given that ID, this method
returns the feature.
- $iterator = $sam->get_seq_stream(%options)
- Convenience method. This is the same as calling
$sam->features(%options,-iterator=>1).
- $fh = $sam->get_seq_fh(%options)
- Convenience method. This is the same as calling
$sam->features(%options,-fh=>1).
- $fh = $sam->tam_fh
- Convenience method. It is the same as calling
$sam->features(-fh=>1).
- @types = $sam->types
- This method returns the list of feature types (e.g.
"read_pair") returned by the current version of the
interface.
Lastly, the high-level API supports two methods for rapidly traversing indexed
BAM databases.
- $sam->fetch($region,$callback)
- This method, which is named after the native
bam_fetch() function in the C interface, traverses the indicated
region and invokes a callback code reference on each match. Specify a
region using the BAM syntax "seqid:start-end", or either of the
alternative syntaxes "seqid:start..end" and
"seqid:start,end". If start and end are absent, then the entire
reference sequence is traversed. If end is absent, then the end of the
reference sequence is assumed.
The callback will be called repeatedly with a Bio::DB::Bam::AlignWrapper on
the argument list.
Example:
$sam->fetch('seq1:600-700',
sub {
my $a = shift;
print $a->display_name,' ',$a->cigar_str,"\n";
});
Note that the fetch() operation works on reads that overlap
the indicated region. Therefore the callback may be called for reads that
align to the reference at positions that start before or end after the
indicated region.
- $sam->pileup($region,$callback [,$keep_level])
- This method, which is named after the native
bam_lpileupfile() function in the C interfaces, traverses the
indicated region and generates a "pileup" of all the mapped
reads that cover it. The user-provided callback function is then invoked
on each position of the alignment along with a data structure that
provides access to the individual aligned reads.
As with fetch(), the region is specified as a string in the format
"seqid:start-end", "seqid:start..end" or
"seqid:start,end".
The callback is a coderef that will be invoked with three arguments: the
seq_id of the reference sequence, the current position on the reference
(in 1-based coordinates!), and a reference to an array of
Bio::DB::Bam::Pileup objects. Here is the typical call signature:
sub {
my ($seqid,$pos,$pileup) = @_;
# do something
}
For example, if you call pileup on the region "seq1:501-600", then
the callback will be invoked for all reads that overlap the indicated
region. The first invocation of the callback will typically have a $pos
argument somewhat to the left of the desired region and the last call will
be somewhat to the right. You may wish to ignore positions that are
outside of the requested region. Also be aware that the reference sequence
position uses 1-based coordinates, which is different from the low-level
interface, which use 0-based coordinates.
The optional $keep_level argument, if true, asks the BAM library to keep
track of the level of the read in the multiple alignment, an operation
that generates some overhead. This is mostly useful for text alignment
viewers, and so is off by default.
The size of the $pileup array reference indicates the read coverage at that
position. Here is a simple average coverage calculator:
my $depth = 0;
my $positions = 0;
my $callback = sub {
my ($seqid,$pos,$pileup) = @_;
next unless $pos >= 501 && $pos <= 600;
$positions++;
$depth += @$pileup;
}
$sam->pileup('seq1:501-600',$callback);
print "coverage = ",$depth/$positions;
Each Bio::DB::Bam::Pileup object describes the position of a read in the
alignment. Briefly, Bio::DB::Bam::Pileup has the following methods:
$pileup->alignment The alignment at this level (a
Bio::DB::Bam::AlignWrapper object).
$pileup->qpos The position of the read base at the pileup site,
in 0-based coordinates.
$pileup->pos The position of the read base at the pileup site,
in 1-based coordinates;
$pileup->level The level of the read in the multiple alignment
view. Note that this field is only valid when
$keep_level is true.
$pileup->indel Length of the indel at this position: 0 for no indel, positive
for an insertion (relative to the reference), negative for a
deletion (relative to the reference.)
$pileup->is_del True if the base on the padded read is a deletion.
$pileup->is_refskip True if the base on the padded read is a gap relative to the reference (denoted as < or > in the pileup)
$pileup->is_head Undocumented field in the bam.h header file.
$pileup->is_tail Undocumented field in the bam.h header file.
See "Examples" for a very simple SNP caller.
- $sam->fast_pileup($region,$callback [,$keep_level])
- This is identical to pileup() except that the pileup
object returns low-level Bio::DB::Bam::Alignment objects rather than the
higher-level Bio::DB::Bam::AlignWrapper objects. This makes it roughly 50%
faster, but you lose the align objects' seq_id() and
get_tag_values() methods. As a compensation, the callback receives
an additional argument corresponding to the Bio::DB::Sam object. You can
use this to create AlignWrapper objects on an as needed basis:
my $callback = sub {
my($seqid,$pos,$pileup,$sam) = @_;
for my $p (@$pileup) {
my $alignment = $p->alignment;
my $wrapper = Bio::DB::Bam::AlignWrapper->new($alignment,$sam);
my $has_mate = $wrapper->get_tag_values('PAIRED');
}
};
- Bio::DB::Sam->max_pileup_cnt([$new_cnt])
- $sam->max_pileup_cnt([$new_cnt])
- The Samtools library caps pileups at a set level,
defaulting to 8000. The callback will not be invoked on a single position
more than the level set by the cap, even if there are more reads. Called
with no arguments, this method returns the current cap value. Called with
a numeric argument, it changes the cap. There is currently no way to
specify an unlimited cap.
This method can be called as an instance method or a class method.
- $sam->coverage2BedGraph([$fh])
- This special-purpose method will compute a four-column BED
graph of the coverage across the entire SAM/BAM file and print it to
STDOUT. You may provide a filehandle to redirect output to a file or
pipe.
The next sections correspond to the low-level API, which let you create and
manipulate Perl objects that correspond directly to data structures in the C
interface. A major difference between the high and low level APIs is that in
the high-level API, the reference sequence is identified using a
human-readable seq_id. However, in the low-level API, the reference is
identified using a numeric target ID ("tid"). The target ID is
established during the creation of the BAM file and is a small 0-based integer
index. The Bio::DB::Bam::Header object provides methods for converting from
seq_ids to tids.
These methods relate to the BAM library's indexed Fasta (".fai")
files.
- $fai = Bio::DB::Sam::Fai->load('/path/to/file.fa')
- Load an indexed Fasta file and return the object
corresponding to it. If the index does not exist, it will be created
automatically. Note that you pass the path to the Fasta file, not the
index.
For consistency with Bio::DB::Bam-> open() this method is also
called open().
- $dna_string =
$fai->fetch("seqid:start-end")
- Given a sequence ID contained in the Fasta file and
optionally a subrange in the form "start-end", finds the
indicated subsequence and returns it as a string.
These methods provide interfaces to the "TAM" text version of SAM
files; they often have a .sam extension.
- $tam = Bio::DB::Tam->open('/path/to/file.sam')
- Given the path to a SAM file, opens it for reading. The
file can be compressed with gzip if desired.
- $header = $tam->header_read()
- Create and return a Bio::DB::Bam::Header object from the
information contained within @SQ header lines of the Sam file. If there
are no @SQ lines, then the header will not be useful, and you should call
header_read2() to generate the missing information from the
appropriate indexed Fasta file. Here is some code to illustrate the
suggested logic:
my $header = $tam->header_read;
unless ($header->n_targets > 0) {
$header = $tam->header_read2('/path/to/file.fa.fai');
}
- $header =
$tam->header_read2('/path/to/file.fa.fai')
- Create and return a Bio::DB::Bam::Header object from the
information contained within the indexed Fasta file of the reference
sequences. Note that you have to pass the path to the .fai file, and not
the .fa file. The header object contains information on the reference
sequence names and lengths.
- $bytes = $tam->read1($header,$alignment)
- Given a Bio::DB::Bam::Header object, such as the one
created by header_read2(), and a Bio::DB::Bam::Alignment object
created by Bio::DB::Bam::Alignment-> new(), reads one line of
alignment information into the alignment object from the TAM file and
returns a status code. The result code will be the number of bytes
read.
These methods provide interfaces to the "BAM" binary version of SAM.
They usually have a .bam extension.
- $bam = Bio::DB::Bam->open('/path/to/file.bam'
[,$mode])
- Open up the BAM file at the indicated path. Mode, if
present, must be one of the file stream open flags ("r",
"w", "a", "r+", etc.). If absent, mode
defaults to "r".
Note that Bio::DB::Bam objects are not stable across fork()
operations. If you fork, and intend to use the object in both parent and
child, you must reopen the Bio::DB::Bam in either the child or the parent
(but not both) before attempting to call any of the object's methods.
The path may be an http: or ftp: URL, in which case a copy of the index file
will be downloaded to the current working directory (see below) and all
accesses will be performed on the remote BAM file.
Example:
$bam = Bio::DB::Bam->open('http://some.site.com/nextgen/chr1_bowtie.bam');
- $header = $bam->header()
- Given an open BAM file, return a Bio::DB::Bam::Header
object containing information about the reference sequence(s). Note that
you must invoke header() at least once before calling
read1().
- $status_code = $bam->header_write($header)
- Given a Bio::DB::Bam::Header object and a BAM file opened
in write mode, write the header to the file. If the write fails the
process will be terminated at the C layer. The result code is (currently)
always zero.
- $integer = $bam->tell()
- Return the current position of the BAM file read/write
pointer.
- $bam->seek($integer,$pos)
- Set the current position of the BAM file read/write
pointer. $pos is one of SEEK_SET, SEEK_CUR, SEEK_END. These constants can
be obtained from the Fcntl module by importing the ":seek"
group:
use Fcntl ':seek';
- $alignment = $bam->read1()
- Read one alignment from the BAM file and return it as a
Bio::DB::Bam::Alignment object. Note that you must invoke header()
at least once before calling read1().
- $bytes = $bam->write1($alignment)
- Given a BAM file that has been opened in write mode and a
Bio::DB::Bam::Alignment object, write the alignment to the BAM file and
return the number of bytes successfully written.
- Bio::DB::Bam->sort_core($by_qname,$path,$prefix,$max_mem)
- Attempt to sort a BAM file by chromosomal location or name
and create a new sorted BAM file. Arguments are as follows:
Argument Description
-------- -----------
$by_qname If true, sort by read name rather than chromosomal
location.
$path Path to the BAM file
$prefix Prefix to use for the new sorted file. For example,
passing "foo" will result in a BAM file named
"foo.bam".
$max_mem Maximum core memory to use for the sort. If the sort
requires more than this amount of memory, intermediate
sort files will be written to disk. The default, if not
provided is 500M.
The Bio::DB::Bam::Index object provides access to BAM index (.bai) files.
- $status_code =
Bio::DB::Bam->index_build('/path/to/file.bam')
- Given the path to a .bam file, this function attempts to
build a ".bai" index. The process in which the .bam file exists
must be writable by the current process and there must be sufficient disk
space for the operation or the process will be terminated in the C library
layer. The result code is currently always zero, but in the future may
return a negative value to indicate failure.
- $index =
Bio::DB::Bam->index('/path/to/file.bam',$reindex)
- Attempt to open the index for the indicated BAM file. If
$reindex is true, and the index either does not exist or is out of date
with respect to the BAM file (by checking modification dates), then
attempt to rebuild the index. Will throw an exception if the index does
not exist or if attempting to rebuild the index was unsuccessful.
- $index =
Bio::DB::Bam->index_open('/path/to/file.bam')
- Attempt to open the index file for a BAM file, returning a
Bio::DB::Bam::Index object. The filename path to use is the .bam file, not
the .bai file.
- $index =
Bio::DB::Bam->index_open_in_safewd('/path/to/file.bam' [,$mode])
- When opening a remote BAM file, you may not wish for the
index to be downloaded to the current working directory. This version of
index_open copies the index into the directory indicated by the TMPDIR
environment variable or the system-defined /tmp directory if not present.
You may change the environment variable just before the call to change its
behavior.
- $code = $index->fetch($bam,$tid,$start,$end,$callback
[,$callback_data])
- This is the low-level equivalent of the
$sam->fetch() function described for the high-level API. Given a
open BAM file object, the numeric ID of the reference sequence, start and
end ranges on the reference, and a coderef, this function will traverse
the region and repeatedly invoke the coderef with each
Bio::DB::Bam::Alignment object that overlaps the region.
Arguments:
Argument Description
-------- -----------
$bam The Bio::DB::Bam object that corresponds to the
index object.
$tid The target ID of the reference sequence. This can
be obtained by calling $header->parse_region() with
an appropriate opened Bio::DB::Bam::Header object.
$start The start and end positions of the desired range on
the reference sequence given by $tid, in 0-based
$end coordinates. Like the $tid, these can be obtained from
$header->parse_region().
$callback A coderef that will be called for each read overlapping
the designated region.
$callback_data Any arbitrary Perl data that you wish to pass to the
$callback (optional).
The coderef's call signature should look like this:
my $callback = sub {
my ($alignment,$data) = @_;
...
}
The first argument is a Bio::DB::Bam::Alignment object. The second is the
callback data (if any) passed to fetch().
Fetch() returns an integer code, but its meaning is not described in
the SAM/BAM C library documentation.
- $index->pileup($bam,$tid,$start,$end,$callback
[,$callback_data])
- This is the low-level version of the pileup()
method, which allows you to invoke a coderef for every position in a BAM
alignment. Arguments are:
Argument Description
-------- -----------
$bam The Bio::DB::Bam object that corresponds to the
index object.
$tid The target ID of the reference sequence. This can
be obtained by calling $header->parse_region() with
an appropriate opened Bio::DB::Bam::Header object.
$start The start and end positions of the desired range on
the reference sequence given by $tid, in 0-based
$end coordinates. Like the $tid, these can be obtained from
$header->parse_region().
$callback A coderef that will be called for each position of the
alignment across the designated region.
$callback_data Any arbitrary Perl data that you wish to pass to the
$callback (optional).
The callback will be invoked with four arguments corresponding to the
numeric sequence ID of the reference sequence, the zero-based
position on the alignment, an arrayref of Bio::DB::Bam::Pileup objects,
and the callback data, if any. A typical call signature will be this:
$callback = sub {
my ($tid,$pos,$pileups,$callback_data) = @_;
for my $pileup (@$pileups) {
# do something
};
Note that the position argument is zero-based rather than 1-based, as it is
in the high-level API.
The Bio::DB::Bam::Pileup object was described earlier in the description of
the high-level pileup() method.
- $coverage = $index->coverage($bam,$tid,$start,$end
[,$bins [,maxcnt]])
- Calculate coverage for the region on the target sequence
given by $tid between positions $start and $end (zero-based coordinates).
This method will return an array reference equal to the size of the region
(by default). Each element of the array will be an integer indicating the
number of reads aligning over that position. If you provide an option
binsize in $bins, the array will be $bins elements in length, and each
element will contain the average coverage over that region as a floating
point number.
By default, the underlying Samtools library caps coverage counting at a
fixed value of 8000. You may change this default by providing an optional
numeric sixth value, which changes the cap for the duration of the call,
or by invoking Bio::DB::Sam->max_pileup_cnt($new_value), which changes
the cap permanently. Unfortunately there is no way of specifying that you
want an unlimited cap.
The Bio::DB::Bam::Header object contains information regarding the reference
sequence(s) used to construct the corresponding TAM or BAM file. It is most
frequently used to translate between numeric target IDs and human-readable
seq_ids. Headers can be created either from reading from a .fai file with the
Bio::DB::Tam->
header_read2() method, or by reading from a BAM file
using Bio::DB::Bam->
header(). You can also create header objects
from scratch, although there is not much that you can do with such objects at
this point.
- $header = Bio::DB::Bam::Header->new()
- Return a new, empty, header object.
- $n_targets = $header->n_targets
- Return the number of reference sequences in the
database.
- $name_arrayref = $header->target_name
- Return a reference to an array of reference sequence names,
corresponding to the high-level API's seq_ids.
To convert from a target ID to a seq_id, simply index into this array:
$seq_id = $header->target_name->[$tid];
- $length_arrayref = $header->target_len
- Return a reference to an array of reference sequence
lengths. To get the length of the sequence corresponding to $tid, just
index into the array returned by target_len():
$length = $header->target_len->[$tid];
- $text = $header->text =item $header->text("new
value")
- Read the text portion of the BAM header. The text can be
replaced by providing the replacement string as an argument. Note that you
should follow the header conventions when replacing the header text. No
parsing or other error-checking is performed.
- ($tid,$start,$end) =
$header->parse_region("seq_id:start-end")
- Given a string in the format "seqid:start-end"
(using a human-readable seq_id and 1-based start and end coordinates),
parse the string and return the target ID and start and end positions in
0-based coordinates. If the range is omitted, then the start and end
coordinates of the entire sequence is returned. If only the end position
is omitted, then the end of the sequence is assumed.
- $header->view1($alignment)
- This method will accept a Bio::DB::Bam::Alignment object,
convert it to a line of TAM output, and write the output to STDOUT. In the
low-level API there is currently no way to send the output to a different
filehandle or capture it as a string.
An array of Bio::DB::Bam::Pileup object is passed to the
pileup()
callback for each position of a multi-read alignment. Each pileup object
contains information about the alignment of a single read at a single
position.
- $alignment = $pileup->alignment
- Return the Bio::DB::Bam::Alignment object at this level.
This provides you with access to the aligning read.
- $alignment = $pileup->b
- An alias for alignment(), provided for compatibility
with the C API.
- $pos = $pileup->qpos
- The position of the aligning base in the read in zero-based
coordinates.
- $pos = $pileup->pos
- The position of the aligning base in 1-based
coordinates.
- $level = $pileup->level
- The "level" of the read in the BAM-generated text
display of the alignment.
- $indel = $pileup->indel
- Length of the indel at this position: 0 for no indel,
positive for an insertion (relative to the reference), negative for a
deletion (relative to the reference sequence.)
- $flag = $pileup->is_del
- True if the base on the padded read is a deletion.
- $flag = $pileup->is_refskip
- True if the base on the padded read is a gap relative to
the reference (denoted as < or > in the pileup)
- $flag = $pileup->is_head
- $flag = $pileup->is_del
- These fields are undocumented in the BAM documentation, but
are exported to the Perl API just in case.
Please see Bio::DB::Bam::Alignment for documentation of the
Bio::DB::Bam::Alignment and Bio::DB::Bam::AlignWrapper objects.
For illustrative purposes only, here is an extremely stupid SNP caller that
tallies up bases that are q>20 and calls a SNP if there are at least 4
non-N/non-indel bases at the position and at least 25% of them are a
non-reference base.
my @SNPs; # this will be list of SNPs
my $snp_caller = sub {
my ($seqid,$pos,$p) = @_;
my $refbase = $sam->segment($seqid,$pos,$pos)->dna;
my ($total,$different);
for my $pileup (@$p) {
my $b = $pileup->alignment;
next if $pileup->indel or $pileup->is_refskip; # don't deal with these ;-)
my $qbase = substr($b->qseq,$pileup->qpos,1);
next if $qbase =~ /[nN]/;
my $qscore = $b->qscore->[$pileup->qpos];
next unless $qscore > 25;
$total++;
$different++ if $refbase ne $qbase;
}
if ($total >= 4 && $different/$total >= 0.25) {
push @SNPs,"$seqid:$pos";
}
};
$sam->pileup('seq1',$snp_caller);
print "Found SNPs: @SNPs\n";
The Bio::DB::Sam interface can be used as a backend to GBrowse
(gmod.sourceforge.net/gbrowse). GBrowse can calculate and display coverage
graphs across large regions, alignment cartoons across intermediate size
regions, and detailed base-pair level alignments across small regions.
Here is a typical configuration for a BAM database that contains information
from a shotgun genomic sequencing project. Some notes:
* It is important to set "search options = none" in order to avoid
GBrowse trying to scan through the BAM database to match read
names. This is a time-consuming operation.
* The callback to "bgcolor" renders pairs whose mates are unmapped in
red.
* The callback to "balloon hover" causes a balloon to pop up with the
read name when the user hovers over each paired read. Otherwise the
default behavior would be to provide information about the pair as
a whole.
* When the user zooms out to 1001 bp or greaterp, the track switches
to a coverage graph.
[bamtest:database]
db_adaptor = Bio::DB::Sam
db_args = -bam /var/www/gbrowse2/databases/bamtest/ex1.bam
search options= default
[Pair]
feature = read_pair
glyph = segments
database = bamtest
draw_target = 1
show_mismatch = 1
bgcolor = sub {
my $f = shift;
return $f->get_tag_values('M_UNMAPPED') ? 'red' : 'green';
}
fgcolor = green
height = 3
label = sub {shift->display_name}
label density = 50
bump = fast
connector = dashed
balloon hover = sub {
my $f = shift;
return '' unless $f->type eq 'match';
return 'Read: '.$f->display_name.' : '.$f->flag_str;
}
key = Read Pairs
[Pair:1000]
feature = coverage:1001
glyph = wiggle_xyplot
height = 50
min_score = 0
autoscale = local
To show alignment data correctly when the user is zoomed in, you should also
provide a pointer to the FASTA file containing the reference genome. In this
case, modify the db_args line to read:
db_args = -bam /var/www/gbrowse2/databases/bamtest/ex1.bam
-fasta /var/www/gbrowse2/databases/bamtest/ex1.fa
Bio::Perl, Bio::DB::Bam::Alignment, Bio::DB::Bam::Constants
Lincoln Stein <
[email protected]>. <
[email protected]>
Copyright (c) 2009-2015 Ontario Institute for Cancer Research.
This package and its accompanying libraries are free software; you can
redistribute it and/or modify it under the terms of the Artistic License 2.0,
the Apache 2.0 License, or the GNU General Public License (version 1 or
higher). Refer to LICENSE for the full license text.