Bio::DB::Sam(3pm)
💖 Bio::DB::Sam 💖
hiii bestie, welcome to the Bio::DB::Sam(3pm) manual page 🫶 this is the User Contributed Perl Documentation situation, and it's literally the module that lets you read SAM/BAM database files from Perl. it's giving bioinformatics main character energy, no cap. grab a boba, we're going in 🌸
✨ NAME ✨
Bio::DB::Sam — read SAM/BAM database files 📖🧬
🌈 SYNOPSIS 🌈
okay so here's the whole vibe, both APIs at once, don't panic:
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);see? iconic. 💅
🌟 DESCRIPTION 🌟
so this module is your Perl lil gateway to the libbam library, for both indexed and unindexed SAM/BAM sequence alignment databases. it's got support for pulling info on individual alignments, read pairs, and alignment coverage across big ol' regions. it also does callback stuff for calling SNPs and doing other base-by-base functions. and the glow-up? most operations are compatible with the BioPerl Bio::SeqFeatureI interface, which means BAM files can be used as a backend to the GBrowse genome browser application (gmod.sourceforge.net). obsessed. 😍
💫 The high-level API 💫
the high-level API gives you a BioPerl-compatible interface to indexed BAM files. the BAM database just gets treated as a collection of Bio::SeqFeatureI features, and you can search for features by name, location, type, and combos of feature tags — like whether the alignment is part of a mate-pair. bestie, it's a whole vibe. 🦋
when you open a BAM database with the high-level API, you hand over the pathnames of two files: the FASTA file with the reference genome sequence, and the BAM file with the query sequences and their alignments. if either file needs indexing, the indexing just... happens automatically, no stress. then you can query the database for alignment features by combinations of name, position, type, and feature tag. ✨
the high-level API gives you 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 typeregionwill return a list of chromosomes or contigs in the reference sequence, rather than read alignments.chromosome— a synonym forregion.
features can be [fetched] en masse in a single call, retrieved in a memory-efficient streaming basis using an iterator, or interrogated using a filehandle that returns a series of TAM-format lines. pick your fighter 💀
SAM alignment flags can be retrieved using BioPerl's feature "tag" mechanism. for example, to interrogate the FIRST_MATE flag, you fetch 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. so versatile. 🎀
split alignments — like reads that cover introns — get dealt with in one of two ways. the default is to leave split alignments alone: you can detect them 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 gives you access to the samtools pileup API. this means you can 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. it's giving very demure, very mindful. 🧘
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 those tags will be used to reconstruct the reference sequence when necessary — in which case you can completely omit the -fasta argument. but heads up bestie: not all SAM/BAM files have MD tags, and the ones that do may not use them correctly, because this part of the SAM spec is new. 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'. periodt.
the main object classes you'll be dealing with in the high-level API are:
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. seeBio::DB::Bam::Alignmentfor 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 🌙
the low-level API closely mirrors the libbam library. it gives you the ability to open TAM and BAM files, read and write to them, build indexes, and perform searches across them. there's less overhead to using the API because there's very little Perl memory management — but ngl, the functions are less convenient to use. some operations, like writing BAM files, are only available through the low-level API. 😤
the classes you'll be interacting with in the low-level API are:
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.
👑 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.
🧚 Bio::DB::Sam Constructor and basic accessors 🧚
$sam = Bio::DB::Sam->new(%options)
the Bio::DB::Sam object combines a Fasta file of the reference sequences with a BAM file, so you can conveniently retrieve human-readable sequence IDs and reference sequences. the new() constructor accepts a -name=>value style list of options, like this:
-bam— path to the BAM file that contains the alignments (required). when using samtools 0.1.6 or higher, anhttp:orftp:URL is accepted.-fasta— path to the Fasta file that contains the reference sequences (optional). alternatively, you may pass any object that supports aseq()orfetch_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, bestie. If it's true, then alignments that are split across splices get broken out into a single alignment containing two sub-alignments. Default isfalse. ✨-split— literally the same thing as-split_splices. Same energy, different name. 💅-force_refseq— always use the reference sequence file to derive the reference sequence, even when the sequence could be derived from the MD tag instead. It's slower, yeah, but it's safer when you're working with BAM files that came from buggy aligners, or when the reference contains non-canonical (modified) bases. Protect yourself. 🫶-autoindex— creates a BAM index file if one doesn't exist yet, or if the current one has a modification date earlier than the BAM file. She keeps herself fresh. 💖
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'll be able to use the interface to fetch the reference sequence's bases. No -fasta? Then calls that return the reference sequence will just hand you sequences consisting entirely of "N". Which is giving "I don't know her." 💀
The -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 get concatenated in easily parseable form as a tag named "FLAGS". See get_all_tags() and get_tag_values() for more info, iykyk.
Any two-letter extension flags, such as H0 or H1, will always appear as separate tags regardless of the setting. They don't play by the rules, they just show up. 👑
-split_splices has the effect of breaking up alignments that contain an "N" operation into subparts for more convenient manipulation. So if you've got 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's some overhead to splitting up the spliced alignments, this option is false by default. It's an investment, not a freebie. 🎀
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. Not cute. Alarmingly not cute. 😭
It is not currently possible to refer to a remote FASTA file. You'll have to download those locally and index them before using. No remote FASTA, periodt.
$flag = $sam->expand_flags([$new_value])— get or set theexpand_flagsoption. You can do this after the object is created and it'll have an immediate effect on all alignments fetched from the BAM file. Instant glow-up. ✨$flag = $sam->split_splices([$new_value])— get or set thesplit_splicesoption. Also doable after object creation, and it'll affect all alignments fetched from the BAM file subsequently. So, moving forward only. 💅$header = $sam->header— returns theBio::DB::Bam::Headerobject associated with the BAM file. You can manipulate the header using the low-level API. Go off.$bam_path = $sam->bam_path— returns the path of the bam file used to create the sam object. This makes the sam object more portable, like it's ready to travel. 🧚$bam = $sam->bam— returns the low-levelBio::DB::Bamobject associated with the opened file.$fai = $sam->fai— returns theBio::DB::Sam::Faiobject associated with the Fasta file. You can then manipulate this object with the low-level API.
The index will be built automatically for you if it doesn't already exist. If index building is necessary, 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. And 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 {}. So like… walk carefully. 🙈
$bai = $sam->bam_index— returns theBio::DB::Bam::Indexobject associated with the BAM file.
The BAM file index will be built automatically for you if it doesn't already exist. On top of that, if the BAM file isn't 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 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::SAMobjects are not stable acrossfork()operations. If you fork, you must callclone()either in the parent or the child process before attempting to call any methods. Fork first, ask questions later, and clone. 🌸
🌈 Getting information about reference sequences 🌈
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 asseq_ids). This is the same as the identifier following the">"sign in the Fasta file (e.g."chr1"). 👑$num_targets = $sam->n_targets— returns 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 aseq_idused 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 aBio::PrimarySeqobject. And if no-fastapath was passed when the sam object was created, then you'll receive a series ofNnucleotides of the requested length. Just Ns. A whole vibe of Ns. 💀
🦋 Creating and querying segments 🦋
Bio::DB::Sam::Segment objects refer to 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. Iconic behavior, honestly.
$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'll be replaced with the beginning and end of the sequence respectively. Generous of her.
Alternatively, you may call segment() with the named -seq_id, -start and -end arguments.
All coordinates are 1-based. Don't come at me with 0-based energy. 🚫
$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 (always0). She's neutral. She's on her own axis.$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 aBio::PrimarySeqobject 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->seqThe advantage of working with a Bio::PrimarySeq object is that you can perform operations on it, including taking its reverse complement and subsequences. So it's redundant-looking but it's actually serving. 💅
@alignments = $segment->features(%args)— return alignments that overlap the segment in the associated BAM file. The optional%argslist lets you filter features by name, tag or other attributes. See the documentation of theBio::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. Team player behavior. 🫶
$segment->name,$segment->display_name,$segment->get_SeqFeatures,$segment->get_tag_values— these methods are provided forBio::SeqFeatureIcompatibility and don't do anything of interest. They're just there for the plot. 🎀
⭐ Retrieving alignments, mate pairs and coverage information ⭐
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(). She's the multitasker of the family.
These methods either return a list of features, an iterator across a list of features, or a filehandle opened on a pseudo-TAM file. Choose your own adventure. 🌈
@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:
okay bestie, here's the full option menu, no skips, no cuts 💅
-type— filter on features of a given type. You can hand it a single scalar typename, or a reference to an array of the feature types you want. The valid types are"match","read_pair","coverage"and"chromosome."(there is a whole explanation of what each one means right below, don't worry 🫶).-name— filter on reads that have the name you picked. Heads up tho: this can be slow unless you also give it the feature location.-seq_id,-start,-end— filter on features that align toseq_idbetweenstartandend.-startand-endmust be used together with-seq_id. If-startand/or-endare missing, they default to1and the end of the reference sequence, respectively.-flags— filter features that match a list of one or more flags. Format is explained below.-attributes,-tags— literally the same thing as-flags, just hanging around for compatibility with other APIs.-filter— filter on features with a coderef. That coderef gets exactly one argument, the feature, and should return true to keep the feature or false to yeet it.-iterator— instead of getting back a list of features, you get an iterator over the results. To get the actual results, call the iterator'snext_seq()method over and over until it returnsundef, which means no more matching features are left.-fh— instead of a list of features, you get a filehandle. Read from the filehandle to get each result in TAM format, one alignment per line read. This only works for features of type"match".
The high-level API brings in this whole concept of a feature "type" so it can hand you some very convenient functions ✨. You pick types with the optional -type argument. Right now, these are the types it supports:
match
The "match" type is the raw, unprocessed SAM alignment. It hands you single reads, mapped or unmapped, doesn't matter. Every match feature's primary_tag() method returns the string "match". The features you get are of type Bio::DB::Bam::AlignWrapper.
read_pair
The "paired_end" type makes the sam interface hunt down and merge mate pairs together. Fetching this type gives you a series of Bio::SeqFeatureI objects, each one as long as the total distance on the reference sequence that the mate pairs span. The top-level feature is type Bio::SeqFeature::Lite, and it holds two Bio::DB::Bam::AlignWrapper subparts inside it.
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 makes the sam interface calculate coverage across the region you designated. It only works properly if you also give it where you want the coverage graph; -seq_id is a mandatory argument for coverage calculation, and -start and -end are optional. The call returns one single Bio::SeqFeatureI object whose primary_tag() is "coverage". To recover the coverage data, call the object's coverage() method — you get an array in list context or an arrayref in scalar context, full of coverage counts across your 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 is at the base pair level. So for a 5000 bp wide region, coverage() gives you back an array or arrayref with exactly 5000 elements. But bestie, you also have the option of calculating coverage across bigger bins 🧚 — just append the number of intervals you care about onto the "coverage" typename. So for example, fetching "coverage:500" returns a feature whose coverage() method gives you the coverage across 500 intervals.
chromosome or region
The "chromosome" and "region" types are interchangeable, no cap. They ask the sam interface to build Bio::DB::Sam::Segments that represent the reference sequences. These two calls give you basically the same result:
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 as matches and other feature types. 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 you don't pass a -type argument at all, then features() defaults to finding features of type "match".
You can also call features() with a plain list of strings — positional arguments, not -type=>value arguments. That gets interpreted as a list of feature types to return:
my ($coverage) = $sam->features('coverage')For a description of all the methods available in the features that come back from this call, go peek at 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 you combine them with a location filter, they kick off an exhaustive search of the BAM database.
Name filters are case-insensitive, and they let you use shell-style * and ? wildcards. Flag filters made with the -flag, -attribute or -tag options have this syntax:
-flag => { FLAG_NAME_1 => ['list','of','possible','values'],
FLAG_NAME_2 => ['list','of','possible','values'],
...
}The value of -flag is a hash reference where the keys are flag names and the values are array references holding lists of acceptable values. The values in each list are OR'd with each other, and the flag names are AND'd with each other.
The -filter option gives you a completely generic filtering interface 🔥. Hand it a reference to a subroutine. It gets called once for every potential feature. Return true to keep the feature, false to discard it. Here's an example of finding every match whose alignment quality score is greater than 80:
@features = $sam->features(-filter=>sub {shift->qual > 80} );By default, features() returns a list of every matching feature. But you can request an iterator over the results list instead by passing -iterator=>1. That gives you an object with 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
}Same energy with -fh: pass it a true value and you get 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. Exactly the same as calling$sam->features(-name=>$name);$feature = $sam->get_feature_by_name($name)— convenience method. Exactly the same as($sam->features(-name=>$name))[0].@features = $sam->get_features_by_location($seqid,$start,$end)— convenience method. Exactly the same as calling$sam->features(-seq_id=>$seqid,-start=>$start,-end=>$end).@features = $sam->get_features_by_flag(%flags)— convenience method. Exactly the same as calling$sam->features(-flags=>\%flags). This method is also calledget_features_by_attribute()andget_features_by_tag(). Example:
@features = $sam->get_features_by_flag(H0=>1)$feature = $sam->get_feature_by_id($id)— the high-level API gives every feature a unique ID made out of its read name, position and strand, and returns it when you call the feature'sprimary_id()method. Give this method that ID, and it hands you back the feature 👑$iterator = $sam->get_seq_stream(%options)— convenience method. Same as calling$sam->features(%options,-iterator=>1).$fh = $sam->get_seq_fh(%options)— convenience method. Same as calling$sam->features(%options,-fh=>1).$fh = $sam->tam_fh— convenience method. Same as calling$sam->features(-fh=>1).@types = $sam->types— returns the list of feature types (like"read_pair") that the current version of the interface hands out.
The generic fetch() and pileup() methods 🌈
Last but absolutely not least, the high-level API supports two methods for rapidly traversing indexed BAM databases.
$sam->fetch($region,$callback)— this method, named after the nativebam_fetch()function in the C interface, traverses the region you indicate and invokes a callback code reference on every match 💫. Specify a region using the BAM syntax"seqid:start-end", or either of the alternative syntaxes"seqid:start..end"and"seqid:start,end". Ifstartandendare both absent, the entire reference sequence gets traversed. Ifendis absent, then the end of the reference sequence is assumed.
The callback gets called over and over 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. So bestie, heads up: your callback can totally get called for reads that align to the reference at positions starting before or ending after the region you asked for. That's not a bug, that's the vibe. 💅
$sam->pileup($region,$callback [,$keep_level])— okay so this one is named after the nativebam_lpileupfile()function in the C interfaces (iykyk), and what it does is traverse the region you gave it and build a whole "pileup" of every mapped read that covers it ✨. Then your user-provided callback function gets invoked on each position of the alignment, along with a data structure that hands you access to the individual aligned reads. It's giving very main character.
Just like with fetch(), you spell the region as a string in the format "seqid:start-end", "seqid:start..end" or "seqid:start,end".
The callback is a coderef that gets 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's the typical call signature:
sub {
my ($seqid,$pos,$pileup) = @_;
# do something
}So like, if you call pileup on the region "seq1:501-600", then the callback fires 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. So you may wish to ignore positions that are outside of the requested region, no cap. Also be aware (fr fr) that the reference sequence position uses 1-based coordinates, which is different from the low-level interface, which uses 0-based coordinates. Not me getting whiplash between the two 🫠.
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, so it 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, slay:
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 (aBio::DB::Bam::AlignWrapperobject).$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_levelis true. 👀$pileup->indel— length of the indel at this position:0for 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 thebam.hheader file.$pileup->is_tail— undocumented field in thebam.hheader file.
See "Examples" for a very simple SNP caller. 🧬
$sam->fast_pileup($region,$callback [,$keep_level])— literally identical topileup()except that the pileup object returns low-levelBio::DB::Bam::Alignmentobjects rather than the higher-levelBio::DB::Bam::AlignWrapperobjects. That makes it roughly 50% faster, but you lose the align objects'seq_id()andget_tag_values()methods 😔. As a compensation, the callback receives an additional argument corresponding to theBio::DB::Samobject. 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 to8000. 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 — nope, not a thing, don't ask bestie 💀.
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 toSTDOUT. You may provide a filehandle to redirect output to a file or pipe. Cute, right? 📊
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. Okay so it's giving ✨ internal integer energy ✨, but it's fine, we adapt.
🦋 Indexed Fasta Files 🦋
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 withBio::DB::Bam->open()this method is also calledopen().$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. 🧬✨
TAM Files
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 withgzipif desired.$header = $tam->header_read()— create and return aBio::DB::Bam::Headerobject from the information contained within@SQheader lines of the Sam file. If there are no@SQlines, then the header will not be useful, and you should callheader_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 aBio::DB::Bam::Headerobject from the information contained within the indexed Fasta file of the reference sequences. Note that you have to pass the path to the.faifile, and not the.fafile 🫶. The header object contains information on the reference sequence names and lengths.$bytes = $tam->read1($header,$alignment)— given aBio::DB::Bam::Headerobject, such as the one created byheader_read2(), and aBio::DB::Bam::Alignmentobject created byBio::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. 📖
BAM Files
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. Do not skip this step, no cap 😭.
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. Remote reads, we love to see it. 🌈
Example:
$bam = Bio::DB::Bam->open('http://some.site.com/nextgen/chr1_bowtie.bam');$header = $bam->header()— given an open BAM file, return aBio::DB::Bam::Headerobject containing information about the reference sequence(s). Note that you must invokeheader()at least once before callingread1(). Yes, once minimum, it's the law 💅.$status_code = $bam->header_write($header)— given aBio::DB::Bam::Headerobject 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.$posis one ofSEEK_SET,SEEK_CUR,SEEK_END. These constants can be obtained from theFcntlmodule by importing the":seek"group:
use Fcntl ':seek';$alignment = $bam->read1()— read one alignment from the BAM file and return it as aBio::DB::Bam::Alignmentobject. Note that you must invokeheader()at least once before callingread1(). (Yes, again. It's giving foreshadowing 🙈.)$bytes = $bam->write1($alignment)— given a BAM file that has been opened in write mode and aBio::DB::Bam::Alignmentobject, 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:$by_qname— bestie, if this is true we sort by read name instead of chromosomal location 💅$path— the path to the BAM file, obviously$prefix— the prefix for your shiny new sorted file. Like, pass"foo"and you get a BAM file named"foo.bam"✨$max_mem— the max core memory the sort is allowed to eat. If the sort needs more than that, intermediate sort files get written to disk instead. And if you don't pass it? The default is500M, no cap 🫶
🌸 BAM index methods 🌸
The Bio::DB::Bam::Index object is literally your key to BAM index (.bai) files, iykyk.
$status_code = Bio::DB::Bam->index_build('/path/to/file.bam')
Hand it the path to a .bam file and this function tries to build a ".bai" index for it. Heads up tho: the process the .bam file lives in has to be writable by the current process, and there must be enough disk space for the whole operation — otherwise the process gets terminated down in the C library layer. 😭 The result code is currently always zero, but in the future it may return a negative value to signal failure.
$index = Bio::DB::Bam->index('/path/to/file.bam',$reindex)
Tries to open the index for the given BAM file. If $reindex is true, and the index either doesn't exist or is out of date compared to the BAM file (it checks modification dates, sleuth behavior 💅), then it attempts to rebuild the index. It will throw an exception if the index doesn't exist or if the rebuild attempt didn't work out.
$index = Bio::DB::Bam->index_open('/path/to/file.bam')
Tries to open the index file for a BAM file and hands you back a Bio::DB::Bam::Index object. And bestie, listen: the filename path you pass is the .bam file, not the .bai file. Don't get it twisted 🙈
$index = Bio::DB::Bam->index_open_in_safewd('/path/to/file.bam' [,$mode])
When you're opening a remote BAM file, you might not want the index just... downloading itself into your current working directory. So this version of index_open copies the index into the directory named by the TMPDIR environment variable, or the system-defined /tmp directory if TMPDIR isn't there. Pro tip: you can change that environment variable right before the call to change how it behaves. 🌈
$code = $index->fetch($bam,$tid,$start,$end,$callback [,$callback_data])
This is the low-level equivalent of the $sam->fetch() function from the high-level API. You give it an open BAM file object, the numeric ID of the reference sequence, start and end ranges on the reference, and a coderef — and it will traverse the region and repeatedly invoke that coderef with every Bio::DB::Bam::Alignment object that overlaps the region. So many alignments, so little time 💖
Arguments:
$bam— theBio::DB::Bamobject that corresponds to the index object.$tid— the target ID of the reference sequence. You can get it by calling$header->parse_region()with an appropriate openedBio::DB::Bam::Headerobject.$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 gets called for each read overlapping the designated region.$callback_data— any arbitrary Perl data you wanna pass along to the$callback(optional).
The coderef's call signature should look like this:
my $callback = sub {
my ($alignment,$data) = @_;
...
}First argument is a Bio::DB::Bam::Alignment object. Second is the callback data (if any) that got passed to fetch().
Fetch() returns an integer code, but ngl, its meaning is not described in the SAM/BAM C library documentation. We just don't know bestie 💀
$index->pileup($bam,$tid,$start,$end,$callback [,$callback_data])
This is the low-level version of the pileup() method — it lets you invoke a coderef for every position in a BAM alignment. Obsessed honestly. Arguments are:
$bam— theBio::DB::Bamobject that corresponds to the index object.$tid— the target ID of the reference sequence. You can get it by calling$header->parse_region()with an appropriate openedBio::DB::Bam::Headerobject.$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 you wish to pass to the$callback(optional).
The callback gets invoked with four arguments: 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 looks like this:
$callback = sub {
my ($tid,$pos,$pileups,$callback_data) = @_;
for my $pileup (@$pileups) {
# do something
};Careful tho — that position argument is zero-based, not 1-based like it is in the high-level API. Don't say I didn't warn you 🎀
The Bio::DB::Bam::Pileup object was already described earlier, back in the description of the high-level pileup() method, so go peek there 👑
$coverage = $index->coverage($bam,$tid,$start,$end [,$bins [,maxcnt]])
Calculates coverage for the region on the target sequence given by $tid between positions $start and $end (zero-based coordinates, always zero-based with us). This method returns an array reference that by default is the size of the region. Each element is an integer telling you how many reads align over that position. But if you pass a binsize in $bins, the array will be $bins elements long, and each element will hold the average coverage over that region as a floating point number. Math but make it cute 🦋
By default the underlying Samtools library caps coverage counting at a fixed value of 8000. You can change that default by giving an optional numeric sixth value — which changes the cap for the duration of the call only — or by invoking Bio::DB::Sam->max_pileup_cnt($new_value), which changes the cap permanently. Unfortunately there is literally no way to say "bestie I want an unlimited cap". Tragic, I know 😭
🎀 BAM header methods 🎀
The Bio::DB::Bam::Header object holds all the info about the reference sequence(s) used to construct the corresponding TAM or BAM file. It's most frequently used to translate between numeric target IDs and human-readable seq_ids. Headers can be created either by 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's not much you can actually do with such objects at this point. 🙈
$header = Bio::DB::Bam::Header->new()
Returns a brand new, empty, header object. Blank slate energy ✨
$n_targets = $header->n_targets
Returns the number of reference sequences in the database.
$name_arrayref = $header->target_name
Returns a reference to an array of reference sequence names, which correspond to the high-level API's seq_ids.
To convert from a target ID to a seq_id, just index into this array, easy:
$seq_id = $header->target_name->[$tid];$length_arrayref = $header->target_len
Returns a reference to an array of reference sequence lengths. To get the length of the sequence corresponding to $tid, just index into the array that target_len() gave you:
$length = $header->target_len->[$tid];$text = $header->text
Reads the text portion of the BAM header.
$header->text("new value")
The text can be replaced by handing over the replacement string as an argument. Note tho: you should follow the header conventions when replacing the header text. No parsing, no error-checking, no nothing is performed. You're on your own bestie 💅
($tid,$start,$end) = $header->parse_region("seq_id:start-end")
Give it a string in the format "seqid:start-end" (human-readable seq_id, 1-based start and end coordinates) and it parses that string and returns the target ID plus start and end positions in 0-based coordinates. If you omit the range entirely, you get the start and end coordinates of the entire sequence. If you only omit the end position, then the end of the sequence is assumed. So forgiving, love that 🫶
$header->view1($alignment)
This one accepts a Bio::DB::Bam::Alignment object, converts it into a line of TAM output, and writes that output straight to STDOUT. And in the low-level API there's currently no way to send the output to a different filehandle or capture it as a string. So STDOUT it is, periodt 💀
💜 Bio::DB::Bam::Pileup methods 💜
An array of Bio::DB::Bam::Pileup objects gets passed to the pileup() callback for every position of a multi-read alignment. Each pileup object holds the info about the alignment of a single read at a single position. Main character energy, one read at a time 👑
$alignment = $pileup->alignment
Returns the Bio::DB::Bam::Alignment object at this level. That's your 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 they're exported to the Perl API just in case. Feral of them honestly, but we respect it 💅
🧬 The alignment objects 🧬
okay bestie, if you need the actual tea on the alignment objects — like the Bio::DB::Bam::Alignment and Bio::DB::Bam::AlignWrapper objects — the docs are living over in Bio::DB::Bam::Alignment. go read that one. it's giving source of truth. 💅
💥 EXAMPLES 💥
not me handing you an extremely stupid SNP caller, for illustrative purposes only. this one 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. periodt. ✨
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";🌈 GBrowse Compatibility 🌈
bestie, 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. it's giving range. 🦋
here is a typical configuration for a BAM database that contains information from a shotgun genomic sequencing project. some notes, okay?
- 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, and we are not doing all that. 💀 - 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 bporgreaterp, the track switches to a coverage graph. (yes,greaterp. iykyk.)
[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 = localto 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👀 SEE ALSO 👀
Bio::PerlBio::DB::Bam::AlignmentBio::DB::Bam::Constants
✍️ AUTHOR ✍️
Lincoln Stein <lincoln.stein@oicr.on.ca>. <lincoln.stein@bmail.com>. okay king. 👑
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.
and that's a wrap on Bio::DB::Sam(3pm) — from perl v5.42.3, dated 2026-08-19. mwah 💋 stay slaying, bestie. 🫶
📦 where this page comes from
- ubuntu stonking: libbio-samtools-perl 1.43-5build2 · GPL-1.0-or-later OR Artistic-2.0 OR Apache-2.0
- 💖 girlified by deepseek-ai/DeepSeek-V4.1-Flash on : a language model's rewrite of the original, under its terms