Bio::D...ignment(3pm)

📦
version 3.01 · 6 versions
, in debian trixie, ubuntu questing

💖 Bio::DB::HTS::Alignment 💖

okay bestie, welcome to the Bio::DB::HTS::Alignment(3pm) page — this is the manual page for the HTS alignment object, and it's a whole lot, but we're gonna get through it together 🫶

📜 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.0

Unless 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.

fr fr, read the license, it's giving legally binding 💅

✨ NAME ✨

Bio::DB::HTS::Alignment -- the HTS alignment object. that's it, that's the name, she's iconic 👑

🎀 SYNOPSIS 🎀

here's the vibe, copy it, run it, live it:

use Bio::DB::HTS;

my $sam = Bio::DB::HTS->new(-fasta=>"data/ex1.fa",
                            -bam  =>"data/ex1.bam");

my @alignments = $sam->get_features_by_location(-seq_id => 'seq2',
                                                -start  => 500,
                                                -end    => 800);
for my $a (@alignments) {
   my $seqid  = $a->seq_id;
   my $start  = $a->start;
   my $end    = $a->end;
   my $strand = $a->strand;
   my $ref_dna= $a->dna;

   my $query_start  = $a->query->start;
   my $query_end    = $a->query->end;
   my $query_strand = $a->query->strand;
   my $query_dna    = $a->query->dna;

   my $cigar     = $a->cigar_str;
   my @scores    = $a->qscore;     # per-base quality scores
   my $score     = $a->qstring;    # TAM-style quality string
   my $match_qual= $a->qual;       # quality of the match

   my $paired = $a->get_tag_values('PAIRED');
}

like, the comments in there are literal and load-bearing, don't delete them 😌

🦋 DESCRIPTION 🦋

so. Bio::DB::HTS::Alignment and Bio::DB::HTS::AlignWrapper are a duo — together they represent an alignment between a sequence read (the "query") and a reference sequence (the "target"). 💫

  • Bio::DB::HTS::Alignment adheres strictly to the C-level BAM library's definition of a bam1_t* and is used in the Bio::DB::HTS low-level API. no cap, strictly.
  • The latter (Bio::DB::HTS::AlignWrapper) adds convenience methods that make it similar to a BioPerl Bio::SeqFeatureI object.

and this manual page describes both of them 🌈

👑 High-level Bio::DB::HTS::Alignment methods 👑

these methods are provided by Bio::DB::HTS::Alignment, and are intended to be compatible with the Bio::SeqFeatureI interfaces.

⚠️ heads up bestie: these objects are not compatible with Bio::Align::AlignI, because the BAM API is fundamentally incompatible with the BioPerl API for alignments. the first deals with the alignment of a single read against the reference sequence, while the second deals with a multiple alignment. different universes, different rules 💀

also, quick note: the high-level API returns Bio::DB::HTS::AlignWrapper objects, except in the case of the callback to the fast_pileup() method. in that one case only, the object returned by calling $pileup->b() is a Bio::DB::HTS::Alignment object — for performance reasons. she's doing it for speed 🏃♀️💨

  • $seq_id = $align->seq_id — return the seq_id of the reference (target) sequence. this method is only available in the Bio::DB::HTS::AlignWrapper extension.
  • $start = $align->start — return the start of the alignment in 1-based reference sequence coordinates.
  • $end = $align->end — return the end of the alignment in 1-based reference sequence coordinates.
  • $len = $align->length — return the length of the alignment on the reference sequence.
  • $mseqid = $align->mate_seq_id — return the seq_id of the mate's reference (target) sequence. only available in the Bio::DB::AlignWrapper extension.
  • $mstart = $align->mate_start — for paired reads, return the start of the mate's alignment in 1-based reference sequence coordinates.
  • $mend = $align->mate_end — for paired reads, return the end position of the mate's alignment in 1-based reference sequence coordinates.
  • $mlen = $align->mate_len — for mate-pairs, retrieve the length of the mate's alignment on the reference sequence.
  • $strand = $align->strand — return the strand of the alignment as -1 for reversed, +1 for forward.
  • $mstrand = $align->mstrand — if the read has a mate pair, return the strand of the mate in the format -1 or +1.
  • $ref_dna = $align->dna — returns the reference sequence's DNA across the aligned region. if an MD tag is present in the alignment, it will be used preferentially to reconstruct the reference sequence. otherwise the reference DNA access object passed to Bio::DB::HTS->new() will be used.
  • $ref_dna = $align->seq — the reference sequence's DNA as a Bio::PrimarySeqI object (useful for passing to BioPerl functions and for calculating subsequences and reverse complements).
  • $query = $align->query — this method returns a Bio::DB::Alignment::Query object that can be used to retrieve information about the query sequence. the next few entries show how to use this object 👇
  • $read_name = $align->query->name — the name of the read.
  • $q_start = $align->query->start — this returns the start position of the query (read) sequence in 1-based coordinates. it acts via a transient Bio::DB::HTS::Query object that is provided for Bio::Graphics compatibility (see Bio::Graphics).
  • $q_end = $align->query->end — this returns the end position of the query sequence in 1-based coordinates.
  • $q_len = $align->query->length — return the length of the alignment on the read.
  • $scores = $align->query->score — return an array reference containing the unpacked quality scores for each base of the query sequence. the length of this array reference will be equal to the length of the read.
  • $read_dna = $align->query->dna — the read's DNA string.
  • $read_seq = $align->query->seq — the read's DNA as a Bio::PrimarySeqI object.
  • $target = $align->target; — the target() method is similar to query(), except that it follows Bio::AlignIO conventions for how to represent minus strand alignments. the object returned has start(), end(), qscore(), dna() and seq() methods, but for minus strand alignments the sequence will be represented as it appears on the reverse strand, rather than on the forward strand. the glow-up here: you get the read as it came off the machine, before being reverse complemented for use in the SAM file. periodt 💅
  • $query = $align->hit — the hit() method is identical to target() and returns information about the read. it is present for compatibility with some of the Bio::Graphics glyphs, which use hit() to represent the non-reference sequence in aligned sequences.
  • $primary_id = $align->primary_id — this method synthesizes a unique ID for the alignment which can be passed to $sam->get_feature_by_id() to retrieve the alignment at a later date.
  • @tags = $align->get_all_tags — return all tag names known to this alignment. this includes SAM flags such as M_UNMAPPED, as well as auxiliary flags such as H0. the behavior of this method depends on the value of -expand_flags when the SAM object was created. if false (the default), then the standard SAM flags will be concatenated together into a single string and stored in a tag named 'FLAGS'. the format of this tag value is the list of one or more flag constants separated by the "|" character, as in: "PAIRED|MAP_PAIR|REVERSED|SECOND_MATE". if -expand_flags was true, then each flag becomes its own named tag, such as "MAP_PAIR".
  • @values = $align->get_tag_values($tag) — given a tag name, such as 'PAIRED' or 'H0', return its value(s). -expand_flags must be true in order to use the standard SAM flag constants as tags. otherwise, they can be fetched by asking for the "FLAGS" tag, or by using the low-level methods described below.
  • $is_true = $align->has_tag($tag) — return true if the alignment has the indicated tag.
  • $string = $align->cigar_str — return the CIGAR string for this alignment in conventional human readable format (e.g. "M34D1M1").
  • $arrayref = $align->cigar_array — return a reference to an array representing the CIGAR string. this is an array of arrays, in which each subarray consists of a CIGAR operation and a count. example:
 [ ['M',34], ['D',1], ['M',1] ]
  • ($ref,$matches,$query) = $align->padded_alignment — return three strings that show the alignment between the reference sequence (the target) and the query. it will look like this:
 $ref     AGTGCCTTTGTTCA-----ACCCCCTTGCAACAACC
 $matches ||||||||||||||     |||||||||||||||||
 $query   AGTGCCTTTGTTCACATAGACCCCCTTGCAACAACC
  • $str = $align->aux — returns the text version of the SAM tags, e.g. "XT:A:M NM:i:2 SM:i:37 AM:i:37 XM:i:1 XO:i:1 XG:i:1 MD:Z:6^C0A47"
  • $str = $align->tam_line — returns the TAM (text) representation of the alignment (available in the high-level "AlignWrapper" interface only).
  • $tag = $align->primary_tag — this is provided for Bio::SeqFeatureI compatibility. return the string "match".
  • $tag = $align->source_tag — this is provided for Bio::SeqFeatureI compatibility. return the string "sam/bam".
  • @parts = $align->get_SeqFeatures — return subfeatures of this alignment. if you have fetched a "read_pair" feature, this will be the two mate pair objects (both of type Bio::DB::HTS::AlignWrapper). if you have -split_splices set to true in the Bio::DB::HTS database, calling get_SeqFeatures() will return the components of split alignments. see "Bio::DB::HTS Constructor and basic accessors" in Bio::DB::HTS for an example of how to use this. 🧚

💫 Low-level Bio::DB::HTS::Alignment methods 💫

okay bestie, these are the down-in-the-weeds ones 🫶 They live on objects of type Bio::DB::HTS::Alignment and on Bio::DB::HTS::AlignWrapper, and they lowkey just mirror the native C API. No cap, they're basically the C API wearing a Perl hoodie 💅

  • $align = Bio::DB::HTS::Alignment->new — makes a brand-new, empty alignment object. You'll usually need this when you're iterating through an HTS file with Bio::DB::HTS->read1(). ✨
  • $tid = $align->tid( [$new_tid] ) — hands you the target ID of the alignment. You can optionally change the tid by passing it in as an argument — but heads up, that's currently the only field you're allowed to change, because the feature was implemented as a proof of principle 🙈
  • $read_name = $align->qname — the name of the read. iconic.
  • $pos = $align->pos — the 0-based leftmost coordinate of the aligned sequence on the reference sequence. 📍
  • $end = $align->calend — the 0-based rightmost coordinate of the aligned sequence on the reference sequence, after alignment gaps get taken into account.
  • $len = $align->cigar2qlen — the length of the query sequence, calculated from the CIGAR string. 🎀
  • $quality = $align->qual — the quality score for the alignment as a whole.
  • $flag = $align->flag — the bitwise flag field (go check the SAM documentation, iykyk).
  • $mtid = $align->mtid — for paired reads, the target ID of the mate's alignment.
  • $mpos = $align->mpos — for paired reads, the 0-based leftmost coordinate of the mate's alignment on the reference sequence.
  • $n_cigar = $align->n_cigar — how many CIGAR operations are in this alignment.
  • $length = $align->l_qseq — the length of the query sequence (that's the read, babe).
  • $dna = $align->qseq — the actual DNA sequence of the query. And just like in the SAM file, reads aligned to the minus strand of the reference come back in reverse complemented form. 🧬
  • $score_str = $align->_qscore — a packed binary string holding the quality scores for every base of the read. It's the same length as the DNA. Unpack it with unpack('C*',$score_str), or just use the high-level qscore() method instead. 💖
  • $score_arry = $align->qscore
  • @score_arry = $align->qscore — in scalar context you get an array reference with the unpacked quality scores for each base of the query sequence; in list context you get an actual list of the scores. And this array is in the same orientation as the reference sequence. 👑
  • $score_str = $align->qstring — the quality string, in the exact same format used in the SAM (TAM) file.
  • $length = $align->isize — the calculated insert size for mapped paired reads.
  • $length = $align->l_aux — the length of the alignment "auxiliary" data.
  • $value = $align->aux_get("tag") — give it an auxiliary tag, like "H0", and it gives you back its value. 🍓
  • @keys = $align->aux_keys — the list of auxiliary tags this alignment knows about.
  • $data = $align->data — a packed string with the alignment data (sequence, quality scores and cigar string all bundled up).
  • $length = $align->data_len — the current length of the alignment data.
  • $length = $align->m_data — the maximum length of the alignment data.
  • $is_paired = $align->paired — true if the aligned read is part of a mate/read pair, regardless of whether the mate actually mapped. 💅
  • $is_proper = $align->proper_pair — true if the aligned read is part of a mate/read pair and both partners mapped to the reference sequence.
  • $is_unmapped = $align->unmapped — true if the read failed to align. 😭
  • $mate_is_unmapped = $align->munmapped — true if the read's mate failed to align.
  • $reversed = $align->reversed — true if the aligned read was reverse complemented before aligning.
  • $mate_reversed = $align->mreversed — true if the aligned read's mate was reverse complemented before aligning.
  • $isize = $align->isize — for mate-pairs, gives you the computed insert size.
  • $arrayref = $align->cigar — this one returns the CIGAR data in its native BAM format 🌈 You get an arrayref where each operation and count are packed together into an 8-bit structure. To decode each element you gotta use these operations:
use Bio::DB::HTS::Constants;
my $c   = $align->cigar;
my $op  = $c->[0] & BAM_CIGAR_MASK;
my $len = $c->[0] >> BAM_CIGAR_SHIFT;

🌸 AUTHOR 🌸

Rishi Nag <rishi@ebi.ac.uk<gt>

🦋 SEE ALSO 🦋

  • Bio::Perl
  • Bio::DB::HTS
  • Bio::DB::HTS::Constants

and that's the whole thing bestie, we ate and left zero crumbs 💅✨ this was perl v5.40.0 energy, dated 2024-10-31, and yeah the page name got squished into Bio::D...ignment(3pm) but we all know what it is 🎀 go forth and align, mwah 🫶

📦 where this page comes from