Sequence Alignments (skbio.alignment)#

This module provides functionality for computing and manipulating sequence alignments. DNA, RNA, and protein sequences can be aligned, as well as sequences with or without custom alphabets.

Alignment structures#

Tabular format of aligned sequences.

TabularMSA(sequences[, metadata, ...])

Store a multiple sequence alignment in tabular (row/column) form.

Compact format of alignment paths.

AlignPath(lengths, states, *[, ranges, ...])

Store an alignment path between sequences.

PairAlignPath(lengths, states, *[, ranges, ...])

Store a pairwise alignment path between two sequences.

Pairwise alignment#

A versatile, efficient and generalizable pairwise alignment algorithm.

pair_align(seq1, seq2, /[, mode, sub_score, ...])

Perform pairwise alignment of two sequences.

Convenience wrappers with preset scoring schemes.

pair_align_nucl(seq1, seq2, /, **kwargs)

Align two nucleotide sequences.

pair_align_prot(seq1, seq2, /, **kwargs)

Align two protein sequences.

Alignment statistics#

align_score(alignment[, sub_score, ...])

Calculate the alignment score of two or more aligned sequences.

Deprecated functionality#

Pure Python algorithms (slow; educational-purposes only)

global_pairwise_align_nucleotide(seq1, seq2)

Globally align nucleotide seqs or alignments with Needleman-Wunsch.

global_pairwise_align_protein(seq1, seq2[, ...])

Globally align pair of protein seqs or alignments with Needleman-Wunsch.

global_pairwise_align(seq1, seq2, ...[, ...])

Globally align a pair of seqs or alignments with Needleman-Wunsch.

local_pairwise_align_nucleotide(seq1, seq2)

Locally align exactly two nucleotide seqs with Smith-Waterman.

local_pairwise_align_protein(seq1, seq2[, ...])

Locally align exactly two protein seqs with Smith-Waterman.

local_pairwise_align(seq1, seq2, ...)

Locally align exactly two seqs with Smith-Waterman.

Tutorial#

Alignment data structure#

Load two DNA sequences that have been previously aligned into a TabularMSA object, using sequence IDs as the MSA’s index:

>>> from skbio import TabularMSA, DNA
>>> seqs = [DNA("ACC--G-GGTA..", metadata={'id': "seq1"}),
...         DNA("TCC--G-GGCA..", metadata={'id': "seq2"})]
>>> msa = TabularMSA(seqs, minter='id')
>>> msa
TabularMSA[DNA]
----------------------
Stats:
    sequence count: 2
    position count: 13
----------------------
ACC--G-GGTA..
TCC--G-GGCA..
>>> msa.index
Index(['seq1', 'seq2'], dtype='object')

Alignment Algorithm Examples#

scikit-bio provides pure-Python implementations of Smith-Waterman and Needleman-Wunsch alignment. These are much slower than implementations in languages such as C, but serve as useful educational examples as they’re simpler to understand for programmers not experienced with C and are easy to experiment with. These are _far_ too slow to use in production applications when run on long sequences or run many times.

Functions are provided for local and global alignment of protein and nucleotide sequences. The global* and local* functions differ in the underlying algorithm that is applied (global* uses Needleman- Wunsch while local* uses Smith-Waterman), and *protein and *nucleotide differ in their default scoring of matches, mismatches, and gaps.

Here we locally align a pair of protein sequences using gap open penalty of 11 and a gap extend penalty of 1 (in other words, it is much more costly to open a new gap than extend an existing one).

>>> from skbio import Protein
>>> from skbio.alignment import local_pairwise_align_protein
>>> s1 = Protein("HEAGAWGHEE")
>>> s2 = Protein("PAWHEAE")
>>> alignment, score, start_end_positions = local_pairwise_align_protein(
...     s1, s2, 11, 1)

This returns an skbio.TabularMSA object, the alignment score, and start/end positions of each aligned sequence:

>>> alignment
TabularMSA[Protein]
---------------------
Stats:
    sequence count: 2
    position count: 5
---------------------
AWGHE
AW-HE
>>> score
25.0
>>> start_end_positions
[(4, 8), (1, 4)]

Similarly, we can perform global alignment of nucleotide sequences:

>>> from skbio import DNA
>>> from skbio.alignment import global_pairwise_align_nucleotide
>>> s1 = DNA("GCGTGCCTAAGGTATGCAAG")
>>> s2 = DNA("ACGTGCCTAGGTACGCAAG")
>>> alignment, score, start_end_positions = global_pairwise_align_nucleotide(
...     s1, s2)
>>> alignment
TabularMSA[DNA]
----------------------
Stats:
    sequence count: 2
    position count: 20
----------------------
GCGTGCCTAAGGTATGCAAG
ACGTGCCTA-GGTACGCAAG