Sequence Collections and Alignments#

Note

Alpha Release of the New SequenceCollection API

We are pleased to announce an alpha release of our new SequenceCollection API. This version can be accessed by specifying the argument new_type=True in the load_unaligned_seqs() or make_unaligned_seqs() functions. The new API is designed for increased efficiency, offering access to the underlying data in multiple formats, including numpy arrays, strings, and bytes (via array(seqs), str(seqs) and bytes(seqs) respectively).

Please be aware that this alpha release has not been fully integrated with the library. Users are encouraged to explore its capabilities but should proceed with caution!

For loading collections of unaligned or aligned sequences see Loading an alignment from a file or url.

What’s the difference between `Alignment` and `ArrayAlignment`?#

The Alignment class can be annotated, meaning you can add annotations to an Alignment or it’s member sequences and you can slice the alignment via those objects. This capability is achieved, under the hood, by having the individual sequences represent gaps as a “span”, rather than explicitly as a “-” character in the sequence itself. This representation is also efficient for very long sequences.

The ArrayAlignment class cannot be annotated. The class represents its sequences as a numpy.ndarray instance. Thus, the gaps are included as a state in the array. This class is better at handling a lot of sequences and should typically be faster. This is the default class returned by the load_aligned_seqs and make_aligned_seqs functions. (See Loading an alignment from a file or url for details.)

You can change alignment types using the to_type() method.

Basic Collection objects#

Constructing a `SequenceCollection` or `Alignment` object from strings#

from cogent3 import make_aligned_seqs, make_unaligned_seqs

dna = {"seq1": "ATGACC", "seq2": "ATCGCC"}
seqs = make_aligned_seqs(data=dna, moltype="dna")
type(seqs)

cogent3.core.alignment.ArrayAlignment

seqs = make_unaligned_seqs(dna, moltype="dna")
type(seqs)

cogent3.core.alignment.SequenceCollection

Constructing a `ArrayAlignment` using `make_aligned_seqs`#

from cogent3 import make_aligned_seqs

dna = {"seq1": "ATGACC", "seq2": "ATCGCC"}
seqs = make_aligned_seqs(data=dna, moltype="dna", array_align=True)
print(type(seqs))
seqs

<class 'cogent3.core.alignment.ArrayAlignment'>

	0
seq1	ATGACC
seq2	..CG..

2 x 6 dna alignment

Converting a `SequenceCollection` to FASTA format#

from cogent3 import load_unaligned_seqs

seqs = load_unaligned_seqs("data/test.paml")
seqs

	0
NineBande	GCAAGGCGCCAACAGAGCAGATGGGCTGAAAGTAAGGAAACATGTAATGATAGGCAGACT
Mouse	GCAGTGAGCCAGCAGAGCAGATGGGCTGCAAGTAAAGGAACATGTAACGACAGGCAGGTT
Human	GCAAGGAGCCAACATAACAGATGGGCTGGAAGTAAGGAAACATGTAATGATAGGCGGACT
HowlerMon	GCAAGGAGCCAACATAACAGATGGGCTGAAAGTGAGGAAACATGTAATGATAGGCAGACT
DogFaced	GCAAGGAGCCAGCAGAACAGATGGGTTGAAACTAAGGAAACATGTAATGATAGGCAGACT

5 x {min=60, median=60, max=60} text sequence collection

Adding new sequences to an existing collection or alignment#

New sequences can be either appended or inserted using the add_seqs method. More than one sequence can be added at the same time. Note that add_seqs does not modify the existing collection/alignment, it creates a new one.

Appending the sequences#

add_seqs without additional parameters will append the sequences to the end of the collection/alignment.

from cogent3 import make_aligned_seqs

aln = make_aligned_seqs(
    [("seq1", "ATGAA------"), ("seq2", "ATG-AGTGATG"), ("seq3", "AT--AG-GATG")],
    moltype="dna",
)
aln

	0
seq2	ATG-AGTGATG
seq1	...A.------
seq3	..-...-....

3 x 11 dna alignment

new_seqs = make_aligned_seqs(
    [("seq0", "ATG-AGT-AGG"), ("seq4", "ATGCC------")], moltype="dna"
)
new_aln = aln.add_seqs(new_seqs)
new_aln

	0
seq2	ATG-AGTGATG
seq1	...A.------
seq3	..-...-....
seq0	.......-.G.
seq4	...CC------

5 x 11 dna alignment

Note

The order is not preserved if you use the to_fasta() method, which sorts sequences by name.

Inserting the sequences#

Sequences can be inserted into an alignment at the specified position using either the before_name or after_name arguments.

new_aln = aln.add_seqs(new_seqs, before_name="seq2")
new_aln

	0
seq2	ATG-AGTGATG
seq1	...A.------
seq0	.......-.G.
seq4	...CC------
seq3	..-...-....

5 x 11 dna alignment

new_aln = aln.add_seqs(new_seqs, after_name="seq2")
new_aln

	0
seq2	ATG-AGTGATG
seq1	...A.------
seq0	.......-.G.
seq4	...CC------
seq3	..-...-....

5 x 11 dna alignment

Inserting sequence(s) based on their alignment to a reference sequence#

Already aligned sequences can be added to an existing Alignment object and aligned at the same time using the add_from_ref_aln method. The alignment is performed based on their alignment to a reference sequence (which must be present in both alignments). The method assumes the first sequence in ref_aln.names[0] is the reference.

from cogent3 import make_aligned_seqs

aln = make_aligned_seqs(
    [("seq1", "ATGAA------"), ("seq2", "ATG-AGTGATG"), ("seq3", "AT--AG-GATG")],
    moltype="dna",
)
ref_aln = make_aligned_seqs(
    [("seq3", "ATAGGATG"), ("seq0", "ATG-AGCG"), ("seq4", "ATGCTGGG")],
    moltype="dna",
)
new_aln = aln.add_from_ref_aln(ref_aln)
new_aln

	0
seq2	ATG-AGTGATG
seq1	...A.------
seq3	..-...-....
seq0	..-.G--AGC.
seq4	..-.GC-TGG.

5 x 11 dna alignment

add_from_ref_aln has the same arguments as add_seqs so before_name and after_name can be used to insert the new sequences at the desired position.

Note

This method does not work with the ArrayAlignment class.

Removing all columns with gaps in a named sequence#

from cogent3 import make_aligned_seqs

aln = make_aligned_seqs(
    [("seq1", "ATGAA---TG-"), ("seq2", "ATG-AGTGATG"), ("seq3", "AT--AG-GATG")],
    moltype="dna",
)
new_aln = aln.get_degapped_relative_to("seq1")
new_aln

	0
seq1	ATGAATG
seq2	...-.AT
seq3	..--.AT

3 x 7 dna alignment

The elements of a collection or alignment#

Accessing individual sequences from a collection or alignment by name#

Using the get_seq method allows for extracting an unaligned sequence from a collection or alignment by name.

from cogent3 import make_aligned_seqs

aln = make_aligned_seqs(
    [("seq1", "ATGAA------"), ("seq2", "ATG-AGTGATG"), ("seq3", "AT--AG-GATG")],
    moltype="dna",
    array_align=False,
)
seq = aln.get_seq("seq1")
seq.name
type(seq)
seq.is_gapped()

False

Alternatively, if you want to extract the aligned (i.e., gapped) sequence from an alignment, you can use get_gapped_seq.

seq = aln.get_gapped_seq("seq1")
seq.is_gapped()
seq

	0
seq1	ATGAA------

DnaSequence, length=11

To see the names of the sequences in a sequence collection, use the names attribute.

aln.names

['seq1', 'seq2', 'seq3']

Slice the sequences from an alignment like a list#

The usual approach is to access a SequenceCollection or Alignment object as a dictionary, obtaining the individual sequences using the titles as “keys” (above). However, one can also iterate through the collection like a list.

from cogent3 import load_aligned_seqs, load_unaligned_seqs

fn = "data/long_testseqs.fasta"
seqs = load_unaligned_seqs(fn, moltype="dna")
my_seq = seqs.seqs[0]
my_seq[:24]

	0
Human	TGTGGCACAAATACTCATGCCAGC

DnaSequence, length=24

type(my_seq)

cogent3.core.sequence.DnaSequence

aln = load_aligned_seqs(fn, moltype="dna")
aln.seqs[0][:24]

	0
Human	TGTGGCACAAATACTCATGCCAGC

DnaSequence, length=24

Getting a subset of sequences from the alignment#

from cogent3 import load_aligned_seqs

aln = load_aligned_seqs("data/test.paml", moltype="dna")
aln.names

['NineBande', 'Mouse', 'Human', 'HowlerMon', 'DogFaced']

new = aln.take_seqs(["Human", "HowlerMon"])
new.names

['Human', 'HowlerMon']

Note

The Alignment class (which you get if you set array_align=False) is more memory efficient. The subset contain references to the original sequences, not copies.

Alignments#

Creating an `Alignment` object from a `SequenceCollection`#

from cogent3 import load_unaligned_seqs
from cogent3.core.alignment import Alignment

seq = load_unaligned_seqs("data/test.paml", moltype="dna")
seq

	0
NineBande	GCAAGGCGCCAACAGAGCAGATGGGCTGAAAGTAAGGAAACATGTAATGATAGGCAGACT
Mouse	GCAGTGAGCCAGCAGAGCAGATGGGCTGCAAGTAAAGGAACATGTAACGACAGGCAGGTT
Human	GCAAGGAGCCAACATAACAGATGGGCTGGAAGTAAGGAAACATGTAATGATAGGCGGACT
HowlerMon	GCAAGGAGCCAACATAACAGATGGGCTGAAAGTGAGGAAACATGTAATGATAGGCAGACT
DogFaced	GCAAGGAGCCAGCAGAACAGATGGGTTGAAACTAAGGAAACATGTAATGATAGGCAGACT

5 x {min=60, median=60, max=60} dna sequence collection

aln = Alignment(seq)
aln

	0
DogFaced	GCAAGGAGCCAGCAGAACAGATGGGTTGAAACTAAGGAAACATGTAATGATAGGCAGACT
NineBande	......C....A....G........C.....G............................
Mouse	...GT...........G........C..C..G...A.G.........C..C......GT.
Human	...........A..T..........C..G..G.......................G....
HowlerMon	...........A..T..........C.....G.G..........................

5 x 60 dna alignment

Convert alignment to DNA, RNA or PROTEIN moltypes#

This is useful if you’ve loaded a sequence alignment without specifying the moltype and later need to convert it using the dedicated method

from cogent3 import make_aligned_seqs

data = [("a", "ACG---"), ("b", "CCTGGG")]
aln = make_aligned_seqs(data=data)
dna = aln.to_dna()
dna

	0
b	CCTGGG
a	A.G---

2 x 6 dna alignment

Or using the generic to_moltype() method

aln.to_moltype("dna")

	0
b	CCTGGG
a	A.G---

2 x 6 dna alignment

To RNA

from cogent3 import make_aligned_seqs

data = [("a", "ACG---"), ("b", "CCUGGG")]
aln = make_aligned_seqs(data=data)
rna = aln.to_rna()
rna

	0
b	CCUGGG
a	A.G---

2 x 6 rna alignment

To PROTEIN

from cogent3 import make_aligned_seqs

data = [("x", "TYV"), ("y", "TE-")]
aln = make_aligned_seqs(data=data)
prot = aln.to_protein()
prot

	0
x	TYV
y	.E-

2 x 3 protein alignment

Handling gaps#

Remove all gaps from an alignment in FASTA format#

This necessarily returns a SequenceCollection.

from cogent3 import load_aligned_seqs

aln = load_aligned_seqs("data/primate_cdx2_promoter.fasta")
degapped = aln.degap()
print(type(degapped))

<class 'cogent3.core.alignment.SequenceCollection'>

Writing sequences to file#

Both collection and alignment objects have a write method. The output format is inferred from the filename suffix,

from cogent3 import make_aligned_seqs

dna = {"seq1": "ATGACC", "seq2": "ATCGCC"}
aln = make_aligned_seqs(data=dna, moltype="dna")
aln.write("sample.fasta")

or by the format argument.

aln.write("sample", format="fasta")

from cogent3.util.io import remove_files

remove_files(["sample", "sample.fasta"], error_on_missing=False)

Converting an alignment to FASTA format#

from cogent3 import load_aligned_seqs
from cogent3.core.alignment import Alignment

seq = load_aligned_seqs("data/long_testseqs.fasta")
aln = Alignment(seq)
fasta_align = aln

Converting an alignment into Phylip format#

from cogent3 import load_aligned_seqs
from cogent3.core.alignment import Alignment

seq = load_aligned_seqs("data/test.paml")
aln = Alignment(seq)
got = aln.to_phylip()
print(got)

5  60
NineBande GCAAGGCGCCAACAGAGCAGATGGGCTGAAAGTAAGGAAACATGTAATGATAGGCAGACT
Mouse     GCAGTGAGCCAGCAGAGCAGATGGGCTGCAAGTAAAGGAACATGTAACGACAGGCAGGTT
Human     GCAAGGAGCCAACATAACAGATGGGCTGGAAGTAAGGAAACATGTAATGATAGGCGGACT
HowlerMon GCAAGGAGCCAACATAACAGATGGGCTGAAAGTGAGGAAACATGTAATGATAGGCAGACT
DogFaced  GCAAGGAGCCAGCAGAACAGATGGGTTGAAACTAAGGAAACATGTAATGATAGGCAGACT

Converting an alignment to a list of strings#

from cogent3 import load_aligned_seqs
from cogent3.core.alignment import Alignment

seq = load_aligned_seqs("data/test.paml")
aln = Alignment(seq)
string_list = aln.to_dict().values()

Slicing an alignment#

By rows (sequences)#

An Alignment can be sliced

from cogent3 import load_aligned_seqs

fn = "data/long_testseqs.fasta"
aln = load_aligned_seqs(fn, moltype="dna")
aln[:24]

	0
DogFaced	TGTGGCACAAATACTCATGCCAAC
Human	......................G.
HowlerMon	......................G.
Mouse	.........G..G.........G.
NineBande	........................

5 x 24 dna alignment

but a SequenceCollection cannot be sliced

from cogent3 import load_unaligned_seqs

fn = "data/long_testseqs.fasta"
seqs = load_unaligned_seqs(fn)
seqs[:24]

---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
Cell In[34], line 5
      3 fn = "data/long_testseqs.fasta"
      4 seqs = load_unaligned_seqs(fn)
----> 5 seqs[:24]

TypeError: 'SequenceCollection' object is not subscriptable

Getting a single column from an alignment#

from cogent3 import load_aligned_seqs

seq = load_aligned_seqs("data/test.paml")
column_four = aln[3]

Getting a region of contiguous columns#

from cogent3 import load_aligned_seqs

aln = load_aligned_seqs("data/long_testseqs.fasta")
region = aln[50:70]

Iterating over alignment positions#

from cogent3 import load_aligned_seqs

aln = load_aligned_seqs("data/primate_cdx2_promoter.fasta")
col = aln[113:115].iter_positions()
type(col)
list(col)

[[ByteSequence(A), ByteSequence(A), ByteSequence(A)],
 [ByteSequence(T), ByteSequence(-), ByteSequence(-)]]

Getting codon 3rd positions from `Alignment`#

We’ll do this by specifying the position indices of interest, creating a sequence Feature and using that to extract the positions.

from cogent3 import make_aligned_seqs

aln = make_aligned_seqs(
    data={"seq1": "ATGATGATG---", "seq2": "ATGATGATGATG"}, array_align=False
)
list(range(len(aln))[2::3])
indices = [(i, i + 1) for i in range(len(aln))[2::3]]
indices

[(2, 3), (5, 6), (8, 9), (11, 12)]

pos3 = aln.add_feature(biotype="pos3", name="pos3", spans=indices)
pos3 = pos3.get_slice()
pos3

	0
seq2	GGGG
seq1	...-

2 x 4 text alignment

Getting codon 3rd positions from `ArrayAlignment`#

We can use more conventional slice notation in this instance. Note, because Python counts from 0, the 3rd position starts at index 2.

from cogent3 import make_aligned_seqs

aln = make_aligned_seqs(
    data={"seq1": "ATGATGATG---", "seq2": "ATGATGATGATG"}, array_align=True
)
pos3 = aln[2::3]
pos3

	0
seq2	GGGG
seq1	...-

2 x 4 bytes alignment

Filtering positions#

Trim terminal stop codons#

For evolutionary analyses that use codon models we need to exclude terminating stop codons. For the case where the sequences are all of length divisible by 3.

from cogent3 import make_aligned_seqs

aln = make_aligned_seqs(
    data={"seq1": "ACGTAA---", "seq2": "ACGACA---", "seq3": "ACGCAATGA"},
    moltype="dna",
)
new = aln.trim_stop_codons()
new

	0
seq2	ACGACA---
seq1	...---...
seq3	...CA....

3 x 9 dna alignment

To detect if the alignment contains sequences not divisible by 3, use the strict argument. This argument covers both allowing partial terminating codons / not divisible by 3.

aln = make_aligned_seqs(
    data={
        "seq1": "ACGTAA---",
        "seq2": "ACGAC----",  # terminal codon incomplete
        "seq3": "ACGCAATGA",
    },
    moltype="dna",
)
new = aln.trim_stop_codons(strict=True)

Eliminating columns with non-nucleotide characters#

We sometimes want to eliminate ambiguous or gap data from our alignments. We show how to exclude alignment columns by the characters they contain. In the first instance we do this just for single nucleotide columns, then for trinucleotides (equivalent for handling codons). Both are done using the no_degenerates method.

from cogent3 import make_aligned_seqs

aln = make_aligned_seqs(
    data=[
        ("seq1", "ATGAAGGTG---"),
        ("seq2", "ATGAAGGTGATG"),
        ("seq3", "ATGAAGGNGATG"),
    ],
    moltype="dna",
)

We apply to nucleotides,

nucs = aln.no_degenerates()
nucs

	0
seq1	ATGAAGGG
seq2	........
seq3	........

3 x 8 dna alignment

Applying the same filter to trinucleotides (specified by setting motif_length=3).

trinucs = aln.no_degenerates(motif_length=3)
trinucs

	0
seq1	ATGAAG
seq2	......
seq3	......

3 x 6 dna alignment

Getting all variable positions from an alignment#

from cogent3 import load_aligned_seqs

aln = load_aligned_seqs("data/long_testseqs.fasta")
pos = aln.variable_positions()
just_variable_aln = aln.take_positions(pos)
just_variable_aln[:10]

	0
DogFaced	AAACAAAATA
Human	..G.....CT
HowlerMon	..G...G.CT
Mouse	GGG.CC.GCT
NineBande	...T....CT

5 x 10 bytes alignment

Getting all constant positions from an alignment#

from cogent3 import load_aligned_seqs

aln = load_aligned_seqs("data/long_testseqs.fasta")
pos = aln.variable_positions()
just_constant_aln = aln.take_positions(pos, negate=True)
just_constant_aln[:10]

	0
DogFaced	TGTGGCACAA
Human	..........
HowlerMon	..........
Mouse	..........
NineBande	..........

5 x 10 bytes alignment

Getting all variable codons from an alignment#

This is done using the filtered method using the motif_length argument. We demonstrate this first for the ArrayAlignment.

from cogent3 import load_aligned_seqs

aln = load_aligned_seqs("data/long_testseqs.fasta")
variable_codons = aln.filtered(
    lambda x: len(set(map(tuple, x))) > 1, motif_length=3
)
just_variable_aln[:9]

	0
DogFaced	AAACAAAAT
Human	..G.....C
HowlerMon	..G...G.C
Mouse	GGG.CC.GC
NineBande	...T....C

5 x 9 bytes alignment

Then for the standard Alignment by first converting the ArrayAlignment.

aln = aln.to_type(array_align=False)
variable_codons = aln.filtered(lambda x: len(set("".join(x))) > 1, motif_length=3)
just_variable_aln[:9]

	0
DogFaced	AAACAAAAT
Human	..G.....C
HowlerMon	..G...G.C
Mouse	GGG.CC.GC
NineBande	...T....C

5 x 9 bytes alignment

Filtering sequences#

Extracting sequences by sequence identifier into a new alignment object#

You can use take_seqs to extract some sequences by sequence identifier from an alignment to a new alignment object:

from cogent3 import load_aligned_seqs

aln = load_aligned_seqs("data/long_testseqs.fasta")
aln.take_seqs(["Human", "Mouse"])

	0
Human	TGTGGCACAAATACTCATGCCAGCTCATTACAGCATGAGAACAGCAGTTTATTACTCACT
Mouse	.........G..G.....................C.....C............G....T.

2 x 2532 (truncated to 2 x 60) bytes alignment

Alternatively, you can extract only the sequences which are not specified by passing negate=True:

aln.take_seqs(["Human", "Mouse"], negate=True)

	0
DogFaced	TGTGGCACAAATACTCATGCCAACTCATTACAGCATGAGAACAGCAGTTTATTATACACT
HowlerMon	......................G...........................G...CT....
NineBande	.........................T............................CT....

3 x 2532 (truncated to 3 x 60) bytes alignment

Extracting sequences using an arbitrary function into a new alignment object#

You can use take_seqs_if to extract sequences into a new alignment object based on whether an arbitrary function applied to the sequence evaluates to True. For example, to extract sequences which don’t contain any N bases you could do the following:

from cogent3 import make_aligned_seqs

aln = make_aligned_seqs(
    data=[
        ("seq1", "ATGAAGGTG---"),
        ("seq2", "ATGAAGGTGATG"),
        ("seq3", "ATGAAGGNGATG"),
    ],
    moltype="dna",
)

def no_N_chars(s):
    return "N" not in s

aln.take_seqs_if(no_N_chars)

	0
seq2	ATGAAGGTGATG
seq1	.........---

2 x 12 dna alignment

You can additionally get the sequences where the provided function evaluates to False:

aln.take_seqs_if(no_N_chars, negate=True)

	0
seq3	ATGAAGGNGATG

1 x 12 dna alignment

Computing alignment statistics#

Getting motif counts#

We state the motif length we want and whether to allow gap or ambiguous characters. The latter only has meaning for IPUAC character sets (the DNA, RNA or PROTEIN moltypes). We illustrate this for the DNA moltype with motif lengths of 1 and 3.

from cogent3 import make_aligned_seqs

aln = make_aligned_seqs(
    data=[
        ("seq1", "ATGAAGGTG---"),
        ("seq2", "ATGAAGGTGATG"),
        ("seq3", "ATGAAGGNGATG"),
    ],
    moltype="dna",
)
counts = aln.counts()
counts

A	C	G	T
11	0	14	7

counts = aln.counts(motif_length=3)
counts

AAA	AAC	AAG	AAT	ACA	ACC	ACG	ACT	AGA	AGC	AGG	AGT	ATA	ATC	ATG	ATT	CAA	CAC	CAG	CAT	CCA	CCC	CCG	CCT	CGA	CGC	CGG	CGT	CTA	CTC	CTG	CTT	GAA	GAC	GAG	GAT	GCA	GCC	GCG	GCT	GGA	GGC	GGG	GGT	GTA	GTC	GTG	GTT	TAA	TAC	TAG	TAT	TCA	TCC	TCG	TCT	TGA	TGC	TGG	TGT	TTA	TTC	TTG	TTT
0	0	3	0	0	0	0	0	0	0	0	0	0	0	5	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	2	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0

counts = aln.counts(include_ambiguity=True)
counts

A	C	G	N	T
11	0	14	1	7

Note

Only the observed motifs are returned, rather than all defined by the alphabet.

Computing motif probabilities from an alignment#

The method get_motif_probs of Alignment objects returns the probabilities for all motifs of a given length. For individual nucleotides:

from cogent3 import load_aligned_seqs

aln = load_aligned_seqs("data/primate_cdx2_promoter.fasta", moltype="dna")
motif_probs = aln.get_motif_probs()
motif_probs

{'T': 0.25520602569782896,
 'C': 0.25808595480726626,
 'A': 0.24390784226849802,
 'G': 0.24280017722640673}

For dinucleotides or longer, we need to pass in an Alphabet with the appropriate word length. Here is an example with trinucleotides:

from cogent3 import DNA, load_aligned_seqs

trinuc_alphabet = DNA.alphabet.get_word_alphabet(3)
aln = load_aligned_seqs("data/primate_cdx2_promoter.fasta", moltype="dna")
motif_probs = aln.get_motif_probs(alphabet=trinuc_alphabet)
for m in sorted(motif_probs, key=lambda x: motif_probs[x], reverse=True):
    print("%s  %.3f" % (m, motif_probs[m]))

CAG  0.037
CCT  0.034
CGC  0.030
TCT  0.029
AAC  0.028
GCC  0.028
AAA  0.027
AAG  0.027
GGG  0.026
TAA  0.025
CCC  0.025
TTC  0.024
CTT  0.024
GAG  0.024
TGT  0.023
CTG  0.023
CTC  0.021
GGA  0.021
TTT  0.020
CAA  0.019
GGT  0.019
TCC  0.019
AGT  0.019
GAA  0.019
ATT  0.018
TAG  0.017
TTG  0.017
AAT  0.017
TCA  0.015
TGG  0.015
TGC  0.015
GTT  0.014
GGC  0.014
CAT  0.013
GTA  0.013
TGA  0.013
ATA  0.012
AGA  0.012
TTA  0.011
CCA  0.011
CAC  0.011
CCG  0.011
GCT  0.011
GCA  0.011
GAT  0.011
TAC  0.010
AGC  0.009
GAC  0.009
CGG  0.008
ATC  0.008
AGG  0.008
GTC  0.008
ACT  0.007
ACG  0.007
CGT  0.007
GCG  0.007
GTG  0.006
TCG  0.005
TAT  0.005
ATG  0.005
CGA  0.005
ACC  0.005
CTA  0.003
ACA  0.002

The same holds for other arbitrary alphabets, as long as they match the alignment MolType.

Some calculations in cogent3 require all non-zero values in the motif probabilities, in which case we use a pseudo-count. We illustrate that here with a simple example where T is missing. Without the pseudo-count, the frequency of T is 0.0, with the pseudo-count defined as 1e-6 then the frequency of T will be slightly less than 1e-6.

aln = make_aligned_seqs(data=[("a", "AACAAC"), ("b", "AAGAAG")], moltype="dna")
motif_probs = aln.get_motif_probs()
assert motif_probs["T"] == 0.0
motif_probs = aln.get_motif_probs(pseudocount=1e-6)
assert 0 < motif_probs["T"] <= 1e-6

It is important to notice that motif probabilities are computed by treating sequences as non-overlapping tuples. Below is a very simple pair of identical sequences where there are clearly 2 ‘AA’ dinucleotides per sequence but only the first one is ‘in-frame’ (frame width = 2).

We then create a dinucleotide Alphabet object and use this to get dinucleotide probabilities. These frequencies are determined by breaking each aligned sequence up into non-overlapping dinucleotides and then doing a count. The expected value for the ‘AA’ dinucleotide in this case will be 2/8 = 0.25.

seqs = [("a", "AACGTAAG"), ("b", "AACGTAAG")]
aln = make_aligned_seqs(data=seqs, moltype="dna")
dinuc_alphabet = DNA.alphabet.get_word_alphabet(2)
motif_probs = aln.get_motif_probs(alphabet=dinuc_alphabet)
assert motif_probs["AA"] == 0.25

What about counting the total incidence of dinucleotides including those not in-frame? A naive application of the Python string object’s count method will not work as desired either because it “returns the number of non-overlapping occurrences”.

seqs = [("my_seq", "AAAGTAAG")]
aln = make_aligned_seqs(data=seqs, moltype="dna")
my_seq = aln.get_seq("my_seq")
my_seq.count("AA")
"AAA".count("AA")
"AAAA".count("AA")

To count all occurrences of a given dinucleotide in a DNA sequence, one could use a standard Python approach such as list comprehension:

from cogent3 import make_seq

seq = make_seq(moltype="dna", seq="AAAGTAAG")
seq
di_nucs = [seq[i : i + 2] for i in range(len(seq) - 1)]
sum([nn == "AA" for nn in di_nucs])

Working with alignment gaps#

Filtering extracted columns for the gap character#

from cogent3 import load_aligned_seqs

aln = load_aligned_seqs("data/primate_cdx2_promoter.fasta")
col = aln[113:115].iter_positions()
c1, c2 = list(col)
c1, c2
list(filter(lambda x: x == "-", c1))
list(filter(lambda x: x == "-", c2))

[ByteSequence(-), ByteSequence(-)]

Calculating the gap fraction#

from cogent3 import load_aligned_seqs

aln = load_aligned_seqs("data/primate_cdx2_promoter.fasta")
for column in aln[113:150].iter_positions():
    ungapped = list(filter(lambda x: x == "-", column))
    gap_fraction = len(ungapped) * 1.0 / len(column)
    print(gap_fraction)

0
6666666666666666
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0

Extracting maps of aligned to unaligned positions (i.e., gap maps)#

It’s often important to know how an alignment position relates to a position in one or more of the sequences in the alignment. The gap_maps method of the individual sequences is useful for this. To get a map of sequence to alignment positions for a specific sequence in your alignment, do the following:

from cogent3 import make_aligned_seqs

aln = make_aligned_seqs(
    data=[
        ("seq1", "ATGAAGG-TG--"),
        ("seq2", "ATG-AGGTGATG"),
        ("seq3", "ATGAAG--GATG"),
    ],
    moltype="dna",
)
seq_to_aln_map = aln.get_gapped_seq("seq1").gap_maps()[0]

It’s now possible to look up positions in the seq1, and find out what they map to in the alignment:

seq_to_aln_map[3]
seq_to_aln_map[8]

This tells us that in position 3 in seq1 corresponds to position 3 in aln, and that position 8 in seq1 corresponds to position 9 in aln.

Notice that we grabbed the first result from the call to gap_maps. This is the sequence position to alignment position map. The second value returned is the alignment position to sequence position map, so if you want to find out what sequence positions the alignment positions correspond to (opposed to what alignment positions the sequence positions correspond to) for a given sequence, you would take the following steps:

aln_to_seq_map = aln.get_gapped_seq("seq1").gap_maps()[1]
aln_to_seq_map[3]
aln_to_seq_map[8]

If an alignment position is a gap, and therefore has no corresponding sequence position, you’ll get a KeyError.

seq_pos = aln_to_seq_map[7]

---------------------------------------------------------------------------
KeyError                                  Traceback (most recent call last)
Cell In[68], line 1
----> 1 seq_pos = aln_to_seq_map[7]

KeyError: 7

Note

The first position in alignments and sequences is always numbered position 0.

Filtering alignments based on gaps#

Note

An alternate, computationally faster, approach to removing gaps is to use the filtered method as discussed in Filtering positions.

The omit_gap_runs method can be applied to remove long stretches of gaps in an alignment. In the following example, we remove sequences that have more than two adjacent gaps anywhere in the aligned sequence.

aln = make_aligned_seqs(
    data=[
        ("seq1", "ATGAA---TG-"),
        ("seq2", "ATG-AGTGATG"),
        ("seq3", "AT--AG-GATG"),
    ],
    moltype="dna",
)
aln.omit_gap_runs(2)

	0
seq2	ATG-AGTGATG
seq3	..-...-....

2 x 11 dna alignment

If instead, we just wanted to remove positions from the alignment which are gaps in more than a certain percentage of the sequences, we could use the omit_gap_pos function. For example:

aln = make_aligned_seqs(
    data=[
        ("seq1", "ATGAA---TG-"),
        ("seq2", "ATG-AGTGATG"),
        ("seq3", "AT--AG-GATG"),
    ],
    moltype="dna",
)
aln.omit_gap_pos(0.40)

	0
seq2	ATGAGGATG
seq1	....--TG-
seq3	..-......

3 x 9 dna alignment

If you wanted to remove sequences which contain more than a certain percent gap characters, you could use the omit_gap_seqs method. This is commonly applied to filter partial sequences from an alignment.

aln = make_aligned_seqs(
    data=[
        ("seq1", "ATGAA------"),
        ("seq2", "ATG-AGTGATG"),
        ("seq3", "AT--AG-GATG"),
    ],
    moltype="dna",
)
filtered_aln = aln.omit_gap_seqs(0.50)
filtered_aln

	0
seq2	ATG-AGTGATG
seq3	..-...-....

2 x 11 dna alignment

Note that following this call to omit_gap_seqs, the 4th column of filtered_aln is 100% gaps. This is generally not desirable, so a call to omit_gap_seqs is frequently followed with a call to omit_gap_pos with no parameters – this defaults to removing positions which are all gaps:

filtered_aln.omit_gap_pos()

	0
seq2	ATGAGTGATG
seq3	..-..-....

2 x 10 dna alignment

Sequence Collections and Alignments#

What’s the difference between Alignment and ArrayAlignment?#

Basic Collection objects#

Constructing a SequenceCollection or Alignment object from strings#

Constructing a ArrayAlignment using make_aligned_seqs#

Converting a SequenceCollection to FASTA format#

Adding new sequences to an existing collection or alignment#

Appending the sequences#

Inserting the sequences#

Inserting sequence(s) based on their alignment to a reference sequence#

Removing all columns with gaps in a named sequence#

The elements of a collection or alignment#

Accessing individual sequences from a collection or alignment by name#

Slice the sequences from an alignment like a list#

Getting a subset of sequences from the alignment#

Alignments#

Creating an Alignment object from a SequenceCollection#

Convert alignment to DNA, RNA or PROTEIN moltypes#

Handling gaps#

Remove all gaps from an alignment in FASTA format#

Writing sequences to file#

Converting an alignment to FASTA format#

Converting an alignment into Phylip format#

Converting an alignment to a list of strings#

Slicing an alignment#

By rows (sequences)#

Getting a single column from an alignment#

Getting a region of contiguous columns#

Iterating over alignment positions#

Getting codon 3rd positions from Alignment#

Getting codon 3rd positions from ArrayAlignment#

Filtering positions#

Trim terminal stop codons#

Eliminating columns with non-nucleotide characters#

Getting all variable positions from an alignment#

Getting all constant positions from an alignment#

Getting all variable codons from an alignment#

Filtering sequences#

Extracting sequences by sequence identifier into a new alignment object#

Extracting sequences using an arbitrary function into a new alignment object#

Computing alignment statistics#

Getting motif counts#

Computing motif probabilities from an alignment#

Working with alignment gaps#

Filtering extracted columns for the gap character#

Calculating the gap fraction#

Extracting maps of aligned to unaligned positions (i.e., gap maps)#

Filtering alignments based on gaps#

What’s the difference between `Alignment` and `ArrayAlignment`?#

Constructing a `SequenceCollection` or `Alignment` object from strings#

Constructing a `ArrayAlignment` using `make_aligned_seqs`#

Converting a `SequenceCollection` to FASTA format#

Creating an `Alignment` object from a `SequenceCollection`#

Getting codon 3rd positions from `Alignment`#

Getting codon 3rd positions from `ArrayAlignment`#