Bio::DB::HTS(3pm)
okay so bestie, you've stumbled onto Bio::DB::HTS(3pm) — a User Contributed Perl Documentation page — and honestly? it's giving main character energy. Buckle up, we're doing six parts, this is part one, and no I will not be normal about it. 💅
📜 LICENSE 📜
Copyright [2015-2018] EMBL-European Bioinformatics Institute
Licensed under the Apache License, Version 2.0 (the "License"); you may not use this file except in compliance with the License. You may obtain a copy of the License at
http://www.apache.org/licenses/LICENSE-2.0Unless required by applicable law or agreed to in writing, software distributed under the License is distributed on an "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. See the License for the specific language governing permissions and limitations under the License.
translating: it's the Apache 2.0 license, it's free, use it, don't sue, keep the notice. iconic behaviour. 🫶
✨ NAME ✨
Bio::DB::HTS — Read files using HTSlib including BAM/CRAM, Tabix and BCF database files. yes girl, all of them.
🌸 SYNOPSIS 🌸
use Bio::DB::HTS;
# high level API
# Note that the high level API does not reset the CRAM file pointer to the start
# of the file as the method to do so is (at time or writing) not easily accessible.
# Therefore a new HTS object may be needed to repeat a query.
my $hts = Bio::DB::HTS->new(-bam =>"data/ex1.bam",
-fasta=>"data/ex1.fa",
);
my @targets = $hts->seq_ids;
my @alignments = $hts->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 = $hts->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 $hfile = Bio::DB::HTSfile->open('/path/to/alignment_file');
my $header = $hfile->header_read;
my $target_count = $header->n_targets;
my $target_names = $header->target_name;
while (my $align = $hfile->read1($header))
{
my $seqid = $target_names->[$align->tid];
my $start = $align->pos+1;
my $end = $align->calend;
my $cigar = $align->cigar_str;
}
Bio::DB::HTSfile->index_build($bamfile);
my $index = Bio::DB::HTSfile->index_load($hfile);
my $index = Bio::DB::HTSfile->index_open_in_safewd($hfile);
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($hfile,$header->parse_region('seq2'),$callback);heads up bestie ⚠️ the high level API does not reset the CRAM file pointer to the start of the file, because the method that does that is (at the time of writing) not easily accessible. So if you wanna repeat a query, you might need to make a whole new HTS object. it's giving commitment issues but we move. 💀
🦋 DESCRIPTION 🦋
This module hands you a Perl interface to the HTSlib library, for both indexed and unindexed SAM/BAM/CRAM sequence alignment databases. It lets you pull info on individual alignments, read pairs, and alignment coverage across big regions. There's also callback functionality for calling SNPs and doing other base-by-base little tasks. she does it all, no notes.
💖 The high-level API 💖
The high-level API is your BioPerl-compatible bestie for indexed BAM and CRAM files. It treats the alignment file database like a whole collection of Bio::SeqFeatureI features, and you can search it by name, location, type, and combos of feature tags — like whether an alignment is part of a mate-pair.
When you open an alignment 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 an index, that indexing has to be built. After that you can query the database for alignment features by any combination of name, position, type, and feature tag. slay.
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; one single composite feature containing two subfeatures, one for each mate's alignment in the pair. they're a package deal fr fr."coverage"— a feature spanning a region of interest that carries numeric info on the coverage of reads across that region."region"— a way to get info about the reference sequence. Searching for"region"features gives you back a list of chromosomes or contigs in the reference sequence instead of read alignments."chromosome"— a synonym for"region".
(yes, the list is longer than "four". we don't gatekeep here. 💅)
You can grab features en masse in one single call, pull them in a memory-efficient streaming way with an iterator, or interrogate them through a filehandle that hands back a series of SAM-format lines. choose your fighter.
SAM alignment flags are retrievable using BioPerl's feature "tag" mechanism. So 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 so it can retrieve all flags as a compact human-readable string, and return the CIGAR alignment in a bunch of different formats. versatile queen.
Split alignments — like reads that cover introns — get handled one of two ways. The default is to leave them alone: you can spot them by one or more "N" operations in the CIGAR string. Optionally, you can choose to have the API split those alignments across two or more subfeatures, and the CIGAR strings of these split alignments get 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 gets invoked on every column of the alignment, for the purpose of calculating coverage, quality score metrics, or SNP calling. per-base realness.
Access to the reference sequence — when you create the Bio::DB::HTS object you can pass the path to a FASTA file holding the reference sequence. Alternatively you can pass an object that knows how to retrieve DNA sequences across a range via the seq() or fetch_seq() methods, as described under new(). 👑
If the SAM/BAM file has MD tags, those tags will be used to reconstruct the reference sequence when necessary — and in that case you can completely omit the -fasta argument. but!! not all SAM/BAM files have MD tags, and the ones that do might not use them correctly, because this part of the SAM spec is still new. You may wish to populate these tags using samtools' "calmd" command. (it's giving "trust but verify".)
If the -fasta argument is omitted and there are no MD tags present, then the reference sequence comes back as 'N'.
The main object classes you'll be dealing with in the high-level API are:
Bio::DB::HTS— a collection of alignments and reference sequences.Bio::DB::HTS::Alignment— the alignment between a query and the reference.Bio::DB::HTS::Query— an object corresponding to the query sequence in which both (+) and (-) strand alignments are shown in the reference (+) strand.Bio::DB::HTS::Target— an interface to the query sequence in which (-) strand alignments are shown in reverse complement.
You may encounter other classes too, including:
Bio::DB::HTS::Segment— corresponds to a region on the reference sequence.Bio::DB::HTS::Constants— defines CIGAR symbol constants and flags.Bio::DB::HTS::AlignWrapper— an alignment helper object that adds split alignment functionality. SeeBio::DB::HTS::Alignmentfor the documentation on using it.Bio::DB::HTS::ReadIterator— an iterator that mediates the one-feature-at-a-time retrieval mechanism.Bio::DB::HTS::FetchIterator— another iterator for feature-at-a-time retrieval.
💜 The low-level API 💜
The low-level API closely mirrors the HTSlib library itself. It gives you the ability to open and read SAM, BAM and CRAM files, build indexes, and perform searches across them. raw, unfiltered, close to the metal.
The classes you'll be interacting with in the low-level API are:
Bio::DB::HTS— methods that read and write SAM, BAM and CRAM files.Bio::DB::HTS::Header— methods for manipulating the BAM file header.Bio::DB::HTS::Alignment— methods for manipulating alignment data.Bio::DB::HTS::Pileup— methods for manipulating the pileup data structure.Bio::DB::HTS::Fai— methods for creating and reading from indexed Fasta files.
💅 METHODS 💅
We cover the high-level API first, bestie. The high-level API code lives in the files Bio/DB/HTS.pm and Bio/DB/HTS/*.pm. iykyk.
🎀 Bio::DB::HTS Constructor and basic accessors 🎀
$sam = Bio::DB::HTS->new(%options)
The Bio::DB::HTS object pairs a Fasta file of reference sequences with an SAM/BAM/CRAM alignment file, so you can conveniently retrieve human-readable sequence IDs and reference sequences. The new() constructor takes a -name=>value style list of options, like this:
-bam— Path to the SAM/BAM/CRAM alignment file that contains the alignments (required). Ahttp:orftp:URL is accepted. yes, a URL, she's flexible like that. 🦋-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).
and that's the vibe for part one, the rest is coming right up — don't touch that dial 🫶
-expand_flags— a boolean. If it's true, the standard alignment flags get broken out into their own individual tags, like'M_UNMAPPED'(default isfalse).-split_splices— a boolean. If it's true, then alignments that are split across splices get broken out into a single alignment that holds two sub-alignments (default isfalse).-split— literally just the same thing as-split_splices.-force_refseq— always use the reference sequence file to derive the reference sequence, even when the sequence could be derived from theMDtag. Yeah it's slower, but it's safer when you're working with BAM files that came out of buggy aligners, or when the reference has non-canonical (modified) bases.-autoindex— create an alignment index file if there isn't one already, or if the current one has a modification date earlier than the alignment file.
An example of a typical new() constructor invocation is:
$hts = Bio::DB::HTS->new(-fasta => '/home/projects/genomes/hu17.fa',
-bam => '/home/projects/alignments/ej88.bam',
-expand_flags => 1,
-split_splices => 1);If the -fasta argument is there, then you'll be able to use the interface to fetch the reference sequence's bases. If it's not? Then any call that returns the reference sequence will hand you sequences made up entirely of "N". 🙃
The -expand_flags option, when true, turns each of the standard SAM flags into a separately retrievable tag in the Bio::SeqFeatureI interface. Otherwise, the standard flags get concatenated in an easily parseable form as a tag named "FLAGS". Go peek at get_all_tags() and get_tag_values() for more info.
Any two-letter extension flags, like H0 or H1, will always show up as separate tags no matter what the setting is. No cap.
-split_splices breaks up alignments that contain an "N" operation into subparts so they're more convenient to manipulate. Say you've got both paired reads and spliced alignments in the BAM file — the following code shows the subpart relationships:
$pair = $hts->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. 💅
Remote access to alignment files living on an HTTP or FTP server is possible. Just swap the path to the BAM file for the appropriate URL. But heads up bestie — incorrect URLs may lead to a core dump. 💀
It is not currently possible to refer to a remote FASTA file. Those have to be downloaded locally and indexed before you can use them.
💖 $flag = $hts->expand_flags([$new_value])
Get or set the expand_flags option. You can do this after object creation and it'll have an immediate effect on all alignments fetched from the alignment file. ✨
🦋 $flag = $hts->split_splices([$new_value])
Get or set the split_splices option. You can do this after object creation too, and it'll affect all alignments fetched from the alignment file subsequently.
🎀 $header = $hts->header
Return the Bio::DB::HTS::Header object associated with the alignment file. You can manipulate the header using the low-level API. iykyk.
🌸 $hts_path = $hts->hts_path
Return the path of the alignment file that was used to create the hts object. This makes the object more portable. 🫶
⭐ $hts_file = $hts->$hts_file
Returns the low-level Bio::DB::HTSfile object associated with the opened file.
🧚 $fai = $hts->fai
Returns the Bio::DB::HTS::Fai object associated with the Fasta file. You can then manipulate this object with the low-level API.
The index can 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 doesn't have write permission, then the call will fail. 😭
🔥 $hts_idx = $hts->hts_index
Return the Bio::DB::HTS::Index object associated with the alignment file.
The index is not automatically built. Just so you know.
👑 $hts->clone
Bio::DB::HTS objects are not stable across fork() operations. So if you fork, you must call clone() either in the parent or the child process before you try to call any methods. Do not skip this, fr fr.
🌈 Getting information about reference sequences
The Bio::DB::HTS object gives you the following methods for getting info about the reference sequence(s) contained in the associated Fasta file. 📖
💜 @seq_ids = $hts->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 = $hts->n_targets
Return the number of reference sequences.
🦋 $length = $hts->length('seqid')
Returns the length of the reference sequence named "seqid".
🎀 $seq_id = $hts->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 = $hts->target_len($tid)
Translates a numeric target ID (TID) from the low-level API into a sequence length.
⭐ $dna = $hts->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. And if no -fasta path was passed when the sam object was created, then you'll receive a series of N nucleotides of the requested length. 🙈
🧚 Creating and querying segments
Bio::DB::HTS::Segment objects refer to regions on the reference sequence. You can use them to retrieve the sequence of the reference, plus alignments that overlap with the region. Obsessed.
$segment = $hts->segment($seqid,$start,$end);
$segment = $hts->segment(-seq_id=>'chr1',-start=>5000,-end=>6000);Segments are created using the Bio:DB::HTS->segment() method. You can call it with 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.
Alternatively, you may call segment() with named -seq_id, -start and -end arguments.
All coordinates are 1-based. Don't forget, bestie.
💅 $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 grab the actual DNA string in this redundant-looking way:
$dna = $segment->seq->seqThe whole 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 alignment file. The optional %args list lets you filter features by name, tag or other attributes. See the docs for the Bio::DB::HTS->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 SAM 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 let the segment 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 honestly don't do anything of interest. 🙈
🌟 Retrieving alignments, mate pairs and coverage information
The features() method is the all-purpose tool for retrieving alignment information from the SAM/BAM/CRAM alignment file database. On top of that, the methods get_features_by_name(), get_features_by_location() and friends give you 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-SAM file.
@features = $hts->features(%options)
$iterator = $hts->features(-iterator=>1,%more_options)
$filehandle = $hts->features(-fh=>1,%more_options)
@features = $hts->features('type1','type2'...)This is the all-purpose interface for fetching alignments and other types of features from the database. The arguments are a -name=>value option list picked from the following list of options: 👇
-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— filters on reads that have the name you name 💅. Heads up bestie: this is a slow operation unless you also give it a feature location to hang out with. Solo? It's gonna go rummaging through everything.-seq_id,-start,-end— filters on features that align toseq_idbetweenstartandend.-startand-endmust be used together with-seq_id, no exceptions. If-startand/or-endghost you, they default to1and the end of the reference sequence respectively. Simple as.-flags— filters features that match a list of one or more flags. Format is down below, keep scrolling 👇.-attributes,-tags— literally the same thing as-flags, just here for compatibility with other APIs. Same energy, different name 🫶.-filter— filter on features using a coderef. That coderef gets one single argument: the feature. Return true to keep it, false to yeet it into the void. That's it, that's the whole job.-iterator— instead of handing you a whole list of features, this gives you an iterator across the results. To actually get your results, call the iterator'snext_seq()method over and over until it returnsundef, which is the "no more matching features, we're done here" signal.-fh— instead of a list of features, you get a filehandle. Read from it to pull each result in TAM format, one alignment per line. Warning though: this only works for features of type"match."No match, no filehandle, periodt.
The high-level API also invented the concept of a feature "type" so it could hand out a bunch of convenience functions like candy 🍓. You set types with the optional -type argument. Currently supported types are:
match
The "match" type is the unprocessed SAM alignment, raw and unfiltered. It retrieves single reads — mapped or unmapped, it does not care. Every match feature's primary_tag() method returns the string "match." The features coming back from this call are Bio::DB::HTS::AlignWrapper objects.
read_pair
The "paired_end" type makes the sam interface go find mate pairs and merge them together like the ultimate friendship reunion 💞. Fetching this type yields a series of Bio::SeqFeatureI objects, each one as long as the total distance on the reference sequence spanned by the mate pairs. The top-level feature is a Bio::SeqFeature::Lite; it contains two Bio::DB::HTS::AlignWrapper subparts.
Call get_SeqFeatures() to get the two individual reads. Like so:
my @pairs = $hts->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 picked. It only behaves properly if you also tell it where the coverage graph goes: -seq_id is a mandatory argument for coverage calculation, and -start and -end are optional. The call returns a 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 the region of interest:
my ($coverage) = $hts->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 sits at the base pair level. So for a region 5000 bp wide, coverage() returns an array or arrayref with exactly 5000 elements. But ✨ plot twist ✨ you can also calculate coverage across bigger bins. Just append the number of intervals you want onto the "coverage" typename. For example, fetching "coverage:500" returns a feature whose coverage() method gives you coverage across 500 intervals. Big brain.
chromosome or region
The "chromosome" or "region" types are interchangeable, no drama between them. They ask the sam interface to construct Bio::DB::HTS::Segment representing the reference sequences. These two calls give you similar results:
my $segment = $hts->segment('seq2',1=>500);
my ($seg) = $hts->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 exactly as expected:my @chromosomes = $hts->features (-type=>'chromosome');But this does not (as of 18 June 2009):my @chromosomes_and_matches = $hts->features(-type=>['match','chromosome']);Not me reading a bug note from 2009 😭 it's still there bestie.
If you don't hand over a -type argument, 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) = $hts->features('coverage')For a description of the methods available in the features returned from this call, please see Bio::SeqfeatureI and Bio::DB::HTS::Alignment. Go read those, they're the receipts 🧾.
You can filter "match" and "read_pair" features by name, location and/or flags. Fair warning: the name and flag filters are not very efficient. Unless you combine them with a location filter, they will initiate an exhaustive search of the BAM database. That's a full sweep of the whole thing, no shortcut, no cap.
Name filters are case-insensitive and let you use shell-style "*" and "?" wildcards. Flag filters created 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 list of values are OR'd with each other, and the flag names are AND'd with each other. That's the whole logic, slay 🔥.
The -filter option gives you a completely generic filtering interface. Hand it a reference to a subroutine. It gets called once for each potential feature. Return true to keep the feature, false to discard it. Here's how you'd find all matches whose alignment quality scores are greater than 80:
@features = $hts->features(-filter=>sub {shift->qual > 80} );By default, features() returns a list of all matching features. You can instead request an iterator across the results list by passing -iterator=>1. That gives you an object with a single method, next_seq():
my $high_qual = $hts->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 SAM file. This only works with the "match" feature type:
my $high_qual = $hts->features(-filter => sub {shift->qual > 80},
-fh => 1 );
while (my $tam_line = <$high_qual>) {
chomp($tam_line);
# do something with it
}@features = $hts->get_features_by_name($name)
Convenience method, doing the absolute most by doing the least. It's the same as calling $hts->features(-name=>$name);
$feature = $hts->get_feature_by_name($name)
Convenience method. Same as ($hts->features(-name=>$name))[0]. Just the first one, no entourage.
@features = $hts->get_features_by_location($seqid,$start,$end)
Convenience method. Same as calling $hts->features(-seq_id=>$seqid,-start=>$start,-end=>$end).
@features = $hts->get_features_by_flag(%flags)
Convenience method. Same as calling $hts->features(-flags=>\%flags). This method is also called get_features_by_attribute() and get_features_by_tag() — three names, one icon 💅. Example:
@features = $hts->get_features_by_flag(H0=>1)$feature = $hts->get_feature_by_id($id)
The high-level API assigns every feature a unique ID made out of its read name, position and strand, and hands it back when you call the feature's primary_id() method. Give that ID to this method and it returns the feature. Like a little lost-and-found for alignments 🎀.
$iterator = $hts->get_seq_stream(%options)
Convenience method. Same as calling $hts->features(%options,-iterator=>1).
$fh = $hts->get_seq_fh(%options)
Convenience method. Same as calling $hts->features(%options,-fh=>1).
$fh = $hts->tam_fh
Convenience method. Same as calling $hts->features(-fh=>1).
@types = $hts->types
This method returns the list of feature types (e.g. "read_pair") returned by the current version of the interface. Know your options bestie 🌈.
🌸 The generic fetch() and pileup() methods 🌸
Lastly, the high-level API supports two methods for rapidly traversing indexed BAM databases. Speed, speed, speed 🚀.
$hts->fetch($region,$callback)
This method traverses the indicated region and invokes a callback code reference on each match. Specify a region with the 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 gets called repeatedly with a Bio::DB::HTS::AlignWrapper on the argument list.
Example:
$hts->fetch('seq1:600-700',
sub {
my $a = shift;
print $a->display_name,' ',$a->cigar_str,"\n";
});Heads up: fetch() operates on reads that overlap the indicated region. So your callback may get called for reads that align to the reference at positions that start before or end after the region you asked for. Not a bug, just how overlap works — don't come for me 🙈.
$hts->pileup($region,$callback [,$keep_level])
This method, 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. Then your user-provided callback function gets invoked on each position of the alignment, along with a data structure that gives you access to the individual aligned reads. Stack 'em up, bestie 🧚.
okay so literally same rules as fetch(): you hand it a region string like "seqid:start-end", "seqid:start..end" or "seqid:start,end". 💖
The callback is a coderef that gets called with three things: the seq_id of the reference sequence, the current position on the reference (in 1-based coordinates, bestie!!), and a reference to an array of Bio::DB::HTS::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", the callback fires for all reads that overlap the indicated region. The first invocation 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 uses 0-based coordinates. 🙈 (why do they do this to us)
The size of the $pileup array reference tells you the read coverage at that position. Here's 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;
}
$hts->pileup('seq1:501-600',$callback);
print "coverage = ",$depth/$positions;Each Bio::DB::HTS::Pileup object describes the position of a read in the alignment. Briefly, Bio::DB::HTS::Pileup has the following methods:
$pileup->alignment— The alignment at this level (aBio::DB::HTS::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, so it may not be relevant post htslib move. 👀$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— True if this is the first base in the query sequence.$pileup->is_tail— True if this is the last base in the query sequence.
See "Examples" for a very simple SNP caller. 🧚
🌸 $hts->fast_pileup($region,$callback [,$keep_level]) 🌸
This is literally pileup()'s twin except that the pileup object hands back low-level Bio::DB::HTS::Alignment objects instead of the higher-level Bio::DB::HTS::AlignWrapper objects. That 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::HTS object. You can use this to create AlignWrapper objects on an as needed basis:
my $callback = sub {
my($seqid,$pos,$pileup,$hts) = @_;
for my $p (@$pileup) {
my $alignment = $p->alignment;
my $wrapper = Bio::DB::HTS::AlignWrapper->new($alignment,$hts);
my $has_mate = $wrapper->get_tag_values('PAIRED');
}
};✨ Bio::DB::HTS->max_pileup_cnt([$new_cnt]) ✨
✨ $hts->max_pileup_cnt([$new_cnt]) ✨
The HTSlib 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. slay.
🎀 $hts->coverage2BedGraph([$fh]) 🎀
This special-purpose method will compute a four-column BED graph of the coverage across the entire alignment 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 alignment file and is a small 0-based integer index. The Bio::DB::HTS::Header object provides methods for converting from seq_ids to tids. 🦋
🌈 Indexed Fasta Files 🌈
These methods relate to the indexed Fasta (".fai") files.
$fai = Bio::DB::HTS::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::HTS->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.
💜 Alignment Files 💜
These methods provide interfaces to alignment files in SAM/BAM/CRAM format.
$hts_file = Bio::DB::HTSfile->open('/path/to/file.bam' [,$mode])— Open the alignment file at the indicated path. Mode, if present, must be one of the file stream open flags ("r","w","wb","wc","a","r+", etc.). If absent, mode defaults to"r". [write formats:w= SAM,wb= BAM,wc= CRAM]
⚠️ Note thatBio::DB::HTSobjects are not stable acrossfork()operations. If you fork, and intend to use the object in both parent and child, you must reopen theBio::DB::HTSin either the child or the parent (but not both) before attempting to call any of the object's methods. fr fr do not skip this one. 😭
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:
$hfile = Bio::DB::HTSfile->open('http://some.site.com/nextgen/chr1_bowtie.bam');$header = $hfile->header_read()— Given an open alignment file, return aBio::DB::HTS::Headerobject containing information about the reference sequence(s). Note that you must invokeheader_read()at least once before callingread1(). 👑$status_code = $hfile->header_write($header, [$reference])— Given aBio::DB::HTSfile::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. 💀 If$hfileis CRAM formated a second argument$reference, which is the path to the reference Fasta file, must be passed. The result code is (currently) always zero.$alignment = $hfile->read1($header)— Read one alignment from the alignment file and return it as aBio::DB::HTS::Alignmentobject. The$headerparameter is returned by invokingheader().$bytes = $hfile->write1($header, $alignment)— Given a BAM file that has been opened in write mode and aBio::DB::HTS::Alignmentobject, write the alignment to the BAM file and return the number of bytes successfully written.
🫶 Index methods 🫶
The Bio::DB::HTS::Index object provides access to index (.bai|.csi, .crai) files.
$status_code = Bio::DB::HTS->index_build('/path/to/file.?am')— Given the path to an alignment file, this function attempts to build an index. The process in which the alignment 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.
The index file built will depend on the alignment file type specified. For CRAM this will be a .crai file, for BAM .bai.
$index = Bio::DB::HTS->index('/path/to/file.?am',$reindex)— Attempt to open the index for the indicated alignment file. If$reindexis true, and the index either does not exist or is out of date with respect to the alignment 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::HTS->index_load('/path/to/file.?am')— Attempt to open the index file for an alignment file, returning aBio::DB::HTS::Indexobject. The filename path to use is the alignment file, not the index file (i.e..bamor.cram, not.bai|.csior.crai). iykyk ✨$index = Bio::DB::HTS->index_open_in_safewd('/path/to/file.?am' [,$mode])— When opening a remote alignmentfile, you may not wish for the index to be downloaded to the current working directory. This version ofindex_opencopies the index into the directory indicated by theTMPDIRenvironment variable or the system-defined/tmpdirectory if not present. You may change the environment variable just before the call to change its behavior. 🎀$code = $index->fetch($hfile,$tid,$start,$end,$callback [,$callback_data])— This is the low-level equivalent of the$hts->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 eachBio::DB::HTS::Alignmentobject that overlaps the region.
Arguments:
Argument Description
-------- -----------
$hts_file The Bio::DB::HTSfile 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::HTS::Header object.$start,$end— where the range you want starts and ends on the reference sequence that$tidpoints at, in 0-based coordinates. Same deal as$tid— you can grab both straight out of$header->parse_region(). 📍$callback— a coderef that gets pinged for every single read that overlaps the region you chose. 💖$callback_data— literally any arbitrary Perl data you wanna hand over to the$callback(optional, bestie). ✨
The coderef's call signature should look exactly like this:
my $callback = sub {
my ($alignment,$data) = @_;
...
}First argument? A Bio::DB::HTS::Alignment object, periodt. Second? The callback data (if any) you passed into fetch(). 🫶
Fetch() returns an integer code, but ngl — its meaning is not described in the SAM/BAM C library documentation. So like, it is what it is. 🙈
$index->pileup($htsfile,$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. Arguments are:
$hts_file— theBio::DB::HTSfileobject that corresponds to the index object.$tid— the target ID of the reference sequence. You get this by calling$header->parse_region()with an appropriately openedBio::DB::HTS::Headerobject.$start,$end— the start and end positions of the desired range on the reference sequence given by$tid, in 0-based coordinates. Just like$tid, these can be snagged from$header->parse_region().$callback— a coderef that gets 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::HTS::Pileup objects, and the callback data, if any. A typical call signature will look like this:
$callback = sub {
my ($tid,$pos,$pileups,$callback_data) = @_;
for my $pileup (@$pileups) {
# do something
};Heads up bestie: the position argument is zero-based rather than 1-based, the way it is in the high-level API. Don't get caught slippin' 😭
The Bio::DB::HTS::Pileup object was described earlier in the description of the high-level pileup() method.
$coverage = $index->coverage($hfile,$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 returns an array reference equal to the size of the region (by default). Each element of the array is an integer indicating how many reads are aligning over that position. If you provide the 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 can 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::HTS->max_pileup_cnt($new_value), which changes the cap permanently. Unfortunately there is no way of specifying that you want an unlimited cap. 💀
✨ HTS header methods ✨
The Bio::DB::HTS::Header object holds all the info about the reference sequence(s) used to build the corresponding alignment file. It's most frequently used to translate between numeric target IDs and human-readable seq_ids. You can create headers by reading from a BAM file using Bio::DB::HTS->header(). You can also create header objects from scratch, although ngl there isn't much you can do with those objects at this point. 🦋
$header = Bio::DB::HTS::Header->new()
Return a new, empty, header object. Fresh out the box. ✨
$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 $header->text("new value")
Read the text portion of the 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 accepts a Bio::DB::HTS::Alignment object, converts it to a line of TAM output, and writes 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. F in the chat. 😭
🌸 Bio::DB::HTS::Pileup methods 🌸
An array of Bio::DB::HTS::Pileup objects gets passed to the pileup() callback for each position of a multi-read alignment. Each pileup object contains info about the alignment of a single read at a single position.
$alignment = $pileup->alignment
Return the Bio::DB::HTS::Alignment object at this level. This gives you 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
True if this is the first base in the query sequence.
$flag = $pileup->is_tail
True if this is the last base in the query sequence.
🎀 The alignment objects 🎀
Please see Bio::DB::HTS::Alignment for documentation of the Bio::DB::HTS::Alignment and Bio::DB::HTS::AlignWrapper objects. Go read it, it's giving main character energy. 👑
🫶 DEPENDENCIES 🫶
Module::Build, Carp, Bio::Perl (>=1.006001), Test::More
✨ EXPORT ✨
None
👑 AUTHORS 👑
Rishi Nag <rishi@ebi.ac.uk>, original author.
Alessandro Vullo "<avullo at cpan.org>", the current developer and maintainer.
🌈 CONTRIBUTORS 🌈
Andy Yates, Keiran Raine, John Marshall, Zhicheng Liu, Can Wood, Dietmar Rieder, Chris Fields, David Jones, James Gilbert, Alex Hodgkins (Congenica Ltd.), Rob Aganrab
💀 KNOWN BUGS 💀
- SAM file reading and iterating over alignments does not work with older htslib versions (
<1.5) - The
padded_alignment()function with CRAM files may produce invalid output: unequal lenght of the strings that specify the pairwise alignment
Please report any bugs or feature requests to "bug-bio-db-hts at rt.cpan.org", or through the web interface at <http://rt.cpan.org/NoAuth/ReportBug.html?Queue=Bio-DB-HTS>. I will be notified, and then you'll automatically be notified of progress on your bug as I make changes. 🫶
🧚 TESTING AND CONTRIBUTING 🧚
You can obtain the most recent development version of this module via the GitHub repository at https://github.com/Ensembl/Bio-DB-HTS. Please feel free to submit bug reports, patches etc. We love a contribution era. 💖
⭐ SUPPORT ⭐
You can find documentation for this module with the perldoc command.
perldoc Bio::DB::HTSYou can also look for information at:
- RT: CPAN's request tracker (report bugs here)
<http://rt.cpan.org/NoAuth/Bugs.html?Dist=Bio-DB-HTS>
- AnnoCPAN: Annotated CPAN documentation
<http://annocpan.org/dist/Bio-DB-HTS>
- CPAN Ratings
<http://cpanratings.perl.org/d/Bio-DB-HTS>
- Search CPAN
<http://search.cpan.org/dist/Bio-DB-HTS/>
🔥 EXAMPLES 🔥
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. It's not it, but it is iconic. 💅
my @SNPs; # this will be list of SNPs
my $snp_caller = sub {
my ($seqid,$pos,$p) = @_;
my $refbase = $hts->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";
}
};
$hts->pileup('seq1',$snp_caller);
print "Found SNPs: @SNPs\n";🦋 GBrowse Compatibility 🦋
okay so real talk, bestie — the Bio::DB::HTS interface can totally moonlight as a backend for GBrowse (gmod.sourceforge.net/gbrowse). like, that's a whole thing it can do. and GBrowse is out here being versatile: it can calculate and display coverage graphs across big ol' regions, alignment cartoons across medium regions, and detailed base-pair level alignments across tiny regions. it's giving range. 💅
Here's a typical configuration for a BAM database stuffed with info from a shotgun genomic sequencing project. some notes before you dive in:
- 🚨 it is important to set
search options = none— this is not optional, this is the vibe check. otherwise GBrowse will try to scan through the entire BAM database to match read names, and that operation is time-consuming. we do not stan. no ma'am. - the callback to
bgcolorrenders pairs whose mates are unmapped in red. - the callback to
balloon hovermakes a balloon pop up with the read name when the user hovers over each paired read. otherwise the default behavior would be to give you info about the pair as a whole. so it's that personal touch. 🎈 - when the user zooms out to
1001bp or greater, the track switches to a coverage graph. slay.
[bamtest:database]
db_adaptor = Bio::DB::HTSfile
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✨ and to show alignment data correctly when the user is zoomed in, you should also hand over a pointer to the FASTA file that holds the reference genome. in that case, modify the db_args line so it reads like this:
db_args = -bam /var/www/gbrowse2/databases/bamtest/ex1.bam
-fasta /var/www/gbrowse2/databases/bamtest/ex1.fa🌈 SEE ALSO 🌈
go read these when you're done here, they're the friend group:
Bio::PerlBio::DB::HTS::AlignmentBio::DB::HTS::Constants
and that's the whole page, bestie!! 🫶 you made it to the end, which honestly is main character energy. this was Bio::DB::HTS(3pm) from perl v5.40.0, dated 2024-10-31. go forth and index some genomes. periodt. 💖✨🎀
📦 where this page comes from
- debian trixie: libbio-db-hts-perl 3.01-5 · Apache-2.0
- ubuntu questing: libbio-db-hts-perl 3.01-5
- 💖 girlified by deepseek-ai/DeepSeek-V4.1-Flash on : a language model's rewrite of the original, under its terms