API Reference

Alignment operation.

'M': non-specific match

'I': insertion into reference sequence

'D': deletion from reference sequence

'N': (typically long) deletion from the reference, e.g. due to RNA splicing

'S': sequence removed from the beginning or end of the query sequence but stored

'H': sequence removed from the beginning or end of the query sequence and not stored

'P': not currently supported, but present for SAM/BAM compatibility

'=': match operation with matching sequence positions

'X': match operation with mismatching sequence positions

'B': not currently supported, but present for SAM/BAM compatibility

'0': indicate the start of an alignment within the reference and query sequence

ismatchop(op::Operation)

Test if op is a match operation (i.e. op ∈ (OP_MATCH, OP_SEQ_MATCH, OP_SEQ_MISMATCH)).

isinsertop(op::Operation)

Test if op is a insertion operation (i.e. op ∈ (OP_INSERT, OP_SOFT_CLIP, OP_HARD_CLIP)).

isdeleteop(op::Operation)

Test if op is a deletion operation (i.e. op ∈ (OP_DELETE, OP_SKIP)).

Alignment operation with anchoring positions.

Alignment of two sequences.

Alignment(anchors::Vector{AlignmentAnchor}, check=true)

Create an alignment object from a sequence of alignment anchors.

Alignment(cigar::AbstractString, seqpos=1, refpos=1)

Make an alignment object from a CIGAR string.

seqpos and refpos specify the starting positions of two sequences.

seq2ref(aln::Alignment, i::Integer)::Tuple{Int,Operation}

Map a position i from sequence to reference.

ref2seq(aln::Alignment, i::Integer)::Tuple{Int,Operation}

Map a position i from reference to sequence.

cigar(aln::Alignment)

Make a CIGAR string encoding of aln.

This is not entirely lossless as it discards the alignments start positions.

Supertype of substitution matrix.

The required method:

• Base.getindex(submat, x, y): substitution score/cost from x to y

Substitution matrix.

Dichotomous substitution matrix.

EDNAFULL (or NUC4.4) substitution matrix

PAM30 substitution matrix

PAM70 substitution matrix

PAM250 substitution matrix

BLOSUM45 substitution matrix

BLOSUM50 substitution matrix

BLOSUM62 substitution matrix

BLOSUM80 substitution matrix

BLOSUM90 substitution matrix

Pairwise alignment

count(aln::PairwiseAlignment, target::Operation)

Count the number of positions where the target operation is applied.

count_matches(aln)

Count the number of matching positions.

count_mismatches(aln)

Count the number of mismatching positions.

count_insertions(aln)

Count the number of inserting positions.

count_deletions(aln)

Count the number of deleting positions.

count_aligned(aln)

Count the number of aligned positions.

Global-global alignment with end gap penalties.

Global-local alignment.

Global-global alignment without end gap penalties.

Local-local alignment.

Edit distance.

Hamming distance.

A special case of EditDistance with the costs of insertion and deletion are infinitely large.

Levenshtein distance.

A special case of EditDistance with the costs of mismatch, insertion, and deletion are 1.

Supertype of score model.

AffineGapScoreModel(submat, gap_open, gap_extend)
AffineGapScoreModel(submat, gap_open=, gap_extend=)
AffineGapScoreModel(match=, mismatch=, gap_open=, gap_extend=)

Affine gap scoring model.

This creates an affine gap scroing model object for alignment from a substitution matrix (submat), a gap opening score (gap_open), and a gap extending score (gap_extend). A consecutive gap of length k has a score of gap_open + gap_extend * k. Note that both of the gap scores should be non-positive. As a shorthand of creating a dichotomous substitution matrix, you can write as, for example, AffineGapScoreModel(match=5, mismatch=-3, gap_open=-2, gap_extend=-1).

Example

using BioSequences
using BioAlignments

# create an affine gap scoring model from a predefined substitution
# matrix and gap opening/extending scores.
affinegap = AffineGapScoreModel(BLOSUM62, gap_open=-10, gap_extend=-1)

# run global alignment between two amino acid sequenecs
pairalign(GlobalAlignment(), aa"IDGAAGQQL", aa"IDGATGQL", affinegap)

See also: SubstitutionMatrix, pairalign, CostModel

Supertype of cost model.

CostModel(submat, insertion, deletion)
CostModel(submat, insertion=, deletion=)
CostModel(match=, mismatch=, insertion=, deletion=)

Cost model.

This creates a cost model object for alignment from substitution matrix (submat), an insertion cost (insertion), and a deletion cost (deletion). Note that both of the insertion and deletion costs should be non-negative. As a shorthand of creating a dichotomous substitution matrix, you can write as, for example, CostModel(match=0, mismatch=1, insertion=2, deletion=2).

Example

using BioAlignments

# create a cost model from a substitution matrix and indel costs
cost = CostModel(ones(128, 128) - eye(128), insertion=.5, deletion=.5)

# run global alignment to minimize edit distance
pairalign(EditDistance(), "intension", "execution", cost)

See also: SubstitutionMatrix, pairalign, AffineGapScoreModel

Result of pairwise alignment

pairalign(type, seq, ref, model, [options...])

Run pairwise alignment between two sequences: seq and ref.

Available types are:

• GlobalAlignment()

• LocalAlignment()

• SemiGlobalAlignment()

• OverlapAlignment()

• EditDistance()

• LevenshteinDistance()

• HammingDistance()

GlobalAlignment, LocalAlignment, SemiGlobalAlignment, and OverlapAlignment are problem that maximizes alignment score between two sequences. Therefore, model should be an instance of AbstractScoreModel (e.g. AffineGapScoreModel).

EditDistance, LevenshteinDistance, and HammingDistance are problem that minimizes alignment cost between two sequences. As for EditDistance, model should be an instance of AbstractCostModel (e.g. CostModel). LevenshteinDistance and HammingDistance have predefined a cost model, so users cannot specify a cost model for these alignment types.

When you pass the score_only=true or distance_only=true option to pairalign, the result of pairwise alignment holds alignment score/distance only. This may enable some algorithms to run faster than calculating full alignment result. Other available options are documented for each alignemnt type.

Example

using BioSequences
using BioAlignments

# create affine gap scoring model
affinegap = AffineGapScoreModel(
    match=5,
    mismatch=-4,
    gap_open=-5,
    gap_extend=-3
)

# run global alignment between two DNA sequences
pairalign(GlobalAlignment(), dna"AGGTAG", dna"ATTG", affinegap)

# run local alignment between two DNA sequences
pairalign(LocalAlignment(), dna"AGGTAG", dna"ATTG", affinegap)

# you cannot specify a cost model in LevenshteinDistance
pairalign(LevenshteinDistance(), dna"AGGTAG", dna"ATTG")

See also: AffineGapScoreModel, CostModel

score(alignment_result)

Return score of alignment.

distance(alignment_result)

Retrun distance of alignment.

alignment(alignment_result)

Return alignment if any.

See also: hasalignment

hasalignment(alignment_result)

Check if alignment is stored or not.

seq2ref(aln::PairwiseAlignment, i::Integer)::Tuple{Int,Operation}

Map a position i from the first sequence to the second.

ref2seq(aln::PairwiseAlignment, i::Integer)::Tuple{Int,Operation}

Map a position i from the second sequence to the first.

SAM.Reader(input::IO)

Create a data reader of the SAM file format.

Arguments

• input: data source

header(reader::Reader)::Header

Get the header of reader.

SAM.Header()

Create an empty header.

find(header::Header, key::AbstractString)::Vector{MetaInfo}

Find metainfo objects satisfying SAM.tag(metainfo) == key.

Writer(output::IO, header::Header=Header())

Create a data writer of the SAM file format.

Arguments

• output: data sink

• header=Header(): SAM header object

MetaInfo(str::AbstractString)

Create a SAM metainfo from str.

Examples

julia> SAM.MetaInfo("@CO	some comment")
BioAlignments.SAM.MetaInfo:
    tag: CO
  value: some comment

julia> SAM.MetaInfo("@SQ	SN:chr1	LN:12345")
BioAlignments.SAM.MetaInfo:
    tag: SQ
  value: SN=chr1 LN=12345

MetaInfo(tag::AbstractString, value)

Create a SAM metainfo with tag and value.

tag is a two-byte ASCII string. If tag is "CO", value must be a string; otherwise, value is an iterable object with key and value pairs.

Examples

julia> SAM.MetaInfo("CO", "some comment")
BioAlignments.SAM.MetaInfo:
    tag: CO
  value: some comment

julia> string(ans)
"@CO	some comment"

julia> SAM.MetaInfo("SQ", ["SN" => "chr1", "LN" => 12345])
BioAlignments.SAM.MetaInfo:
    tag: SQ
  value: SN=chr1 LN=12345

julia> string(ans)
"@SQ	SN:chr1	LN:12345"

iscomment(metainfo::MetaInfo)::Bool

Test if metainfo is a comment (i.e. its tag is "CO").

tag(metainfo::MetaInfo)::String

Get the tag of metainfo.

value(metainfo::MetaInfo)::String

Get the value of metainfo as a string.

keyvalues(metainfo::MetaInfo)::Vector{Pair{String,String}}

Get the values of metainfo as string pairs.

SAM.Record()

Create an unfilled SAM record.

SAM.Record(data::Vector{UInt8})

Create a SAM record from data. This function verifies the format and indexes fields for accessors. Note that the ownership of data is transferred to a new record object.

SAM.Record(str::AbstractString)

Create a SAM record from str. This function verifies the format and indexes fields for accessors.

flag(record::Record)::UInt16

Get the bitwise flag of record.

ismapped(record::Record)::Bool

Test if record is mapped.

isprimary(record::Record)::Bool

Test if record is a primary line of the read.

This is equivalent to flag(record) & 0x900 == 0.

refname(record::Record)::String

Get the reference sequence name of record.

position(record::Record)::Int

Get the 1-based leftmost mapping position of record.

rightposition(record::Record)::Int

Get the 1-based rightmost mapping position of record.

isnextmapped(record::Record)::Bool

Test if the mate/next read of record is mapped.

nextrefname(record::Record)::String

Get the reference name of the mate/next read of record.

nextposition(record::Record)::Int

Get the position of the mate/next read of record.

mappingquality(record::Record)::UInt8

Get the mapping quality of record.

cigar(record::Record)::String

Get the CIGAR string of record.

alignment(record::Record)::BioAlignments.Alignment

Get the alignment of record.

alignlength(record::Record)::Int

Get the alignment length of record.

tempname(record::Record)::String

Get the query template name of record.

templength(record::Record)::Int

Get the template length of record.

sequence(record::Record)::BioSequences.DNASequence

Get the segment sequence of record.

sequence(::Type{String}, record::Record)::String

Get the segment sequence of record as String.

seqlength(record::Record)::Int

Get the sequence length of record.

quality(record::Record)::Vector{UInt8}

Get the Phred-scaled base quality of record.

quality(::Type{String}, record::Record)::String

Get the ASCII-encoded base quality of record.

auxdata(record::Record)::Dict{String,Any}

Get the auxiliary data (optional fields) of record.

0x0001: the read is paired in sequencing, no matter whether it is mapped in a pair

0x0002: the read is mapped in a proper pair

0x0004: the read itself is unmapped; conflictive with SAM.FLAG_PROPER_PAIR

0x0008: the mate is unmapped

0x0010: the read is mapped to the reverse strand

0x0020: the mate is mapped to the reverse strand

0x0040: this is read1

0x0080: this is read2

0x0100: not primary alignment

0x0200: QC failure

0x0400: optical or PCR duplicate

0x0800: supplementary alignment

BAM.Reader(input::IO; index=nothing)

Create a data reader of the BAM file format.

Arguments

• input: data source

• index=nothing: filepath to a random access index (currently bai is supported)

header(reader::Reader; fillSQ::Bool=false)::SAM.Header

Get the header of reader.

If fillSQ is true, this function fills missing "SQ" metainfo in the header.

BAM.Writer(output::BGZFStream, header::SAM.Header)

Create a data writer of the BAM file format.

Arguments

• output: data sink

• header: SAM header object

BAM.Record()

Create an unfilled BAM record.

flag(record::Record)::UInt16

Get the bitwise flag of record.

ismapped(record::Record)::Bool

Test if record is mapped.

isprimary(record::Record)::Bool

Test if record is a primary line of the read.

This is equivalent to flag(record) & 0x900 == 0.

ispositivestrand(record::Record)::Bool

Test if record is aligned to the positive strand.

This is equivalent to flag(record) & 0x10 == 0.

refid(record::Record)::Int

Get the reference sequence ID of record.

The ID is 1-based (i.e. the first sequence is 1) and is 0 for a record without a mapping position.

See also: BAM.rname

refname(record::Record)::String

Get the reference sequence name of record.

See also: BAM.refid

position(record::Record)::Int

Get the 1-based leftmost mapping position of record.

rightposition(record::Record)::Int

Get the 1-based rightmost mapping position of record.

isnextmapped(record::Record)::Bool

Test if the mate/next read of record is mapped.

nextrefid(record::Record)::Int

Get the next/mate reference sequence ID of record.

nextrefname(record::Record)::String

Get the reference name of the mate/next read of record.

nextposition(record::Record)::Int

Get the 1-based leftmost mapping position of the next/mate read of record.

mappingquality(record::Record)::UInt8

Get the mapping quality of record.

cigar(record::Record)::String

Get the CIGAR string of record.

Note that in the BAM specification, the field called cigar typically stores the cigar string of the record. However, this is not always true, sometimes the true cigar is very long, and due to some constraints of the BAM format, the actual cigar string is stored in an extra tag: CG:B,I, and the cigar field stores a pseudo-cigar string.

Calling this method with checkCG set to true (default) this method will always yield the true cigar string, because this is probably what you want the vast majority of the time.

If you have a record that stores the true cigar in a CG:B,I tag, but you still want to access the pseudo-cigar that is stored in the cigar field of the BAM record, then you can set checkCG to false.

See also BAM.cigar_rle.

cigar_rle(record::Record, checkCG::Bool = true)::Tuple{Vector{BioAlignments.Operation},Vector{Int}}

Get a run-length encoded tuple (ops, lens) of the CIGAR string in record.

Note that in the BAM specification, the field called cigar typically stores the cigar string of the record. However, this is not always true, sometimes the true cigar is very long, and due to some constraints of the BAM format, the actual cigar string is stored in an extra tag: CG:B,I, and the cigar field stores a pseudo-cigar string.

Calling this method with checkCG set to true (default) this method will always yield the true cigar string, because this is probably what you want the vast majority of the time.

If you have a record that stores the true cigar in a CG:B,I tag, but you still want to access the pseudo-cigar that is stored in the cigar field of the BAM record, then you can set checkCG to false.

See also BAM.cigar.

alignment(record::Record)::BioAlignments.Alignment

Get the alignment of record.

alignlength(record::Record)::Int

Get the alignment length of record.

tempname(record::Record)::String

Get the query template name of record.

templength(record::Record)::Int

Get the template length of record.

sequence(record::Record)::BioSequences.DNASequence

Get the segment sequence of record.

seqlength(record::Record)::Int

Get the sequence length of record.

quality(record::Record)::Vector{UInt8}

Get the base quality of record.

auxdata(record::Record)::BAM.AuxData

Get the auxiliary data of record.

BAI(filename::AbstractString)

Load a BAI index from filename.

BAI(input::IO)

Load a BAI index from input.