Genetic distance calculation

Fast pairwise distance estimation

For a limited number of evolutionary models a fast implementation is available.

from cogent3 import available_distances

available_distances()

Specify a pairwise genetic distance calculator using 'Abbreviation' (case insensitive).
Abbreviation	Suitable for moltype
paralinear	dna, rna, protein
logdet	dna, rna, protein
jc69	dna, rna
tn93	dna, rna
hamming	dna, rna, protein, text, bytes
pdist	dna, rna, protein, text, bytes

6 rows x 2 columns

Computing genetic distances using the `Alignment` object

Abbreviations listed from available_distances() can be used as values for the distance_matrix(calc=<abbreviation>).

from cogent3 import load_aligned_seqs

aln = load_aligned_seqs("data/primate_brca1.fasta", moltype="dna")
dists = aln.distance_matrix(calc="tn93")
dists

names	Galago	HowlerMon	Rhesus	Orangutan	Gorilla	Human	Chimpanzee
Galago	0.0000	0.2155	0.1959	0.1943	0.1922	0.1963	0.1920
HowlerMon	0.2155	0.0000	0.0736	0.0719	0.0700	0.0736	0.0704
Rhesus	0.1959	0.0736	0.0000	0.0411	0.0393	0.0423	0.0396
Orangutan	0.1943	0.0719	0.0411	0.0000	0.0137	0.0173	0.0140
Gorilla	0.1922	0.0700	0.0393	0.0137	0.0000	0.0086	0.0054
Human	0.1963	0.0736	0.0423	0.0173	0.0086	0.0000	0.0089
Chimpanzee	0.1920	0.0704	0.0396	0.0140	0.0054	0.0089	0.0000

Including ambiguity codes

By default, a site where either sequence has an ambiguity code, such as R (A or G), is excluded. Here the sequences differ only at the last site, so all the distances are 0.

from cogent3 import make_aligned_seqs

aln = make_aligned_seqs(
    {"s1": "ACGTACGTAC", "s2": "ACGTACGTAR", "s3": "ACGTACGTAY"},
    moltype="dna",
)
aln.distance_matrix(calc="tn93")

names	s1	s2	s3
s1	0.0000	0.0000	0.0000
s2	0.0000	0.0000	0.0000
s3	0.0000	0.0000	0.0000

Setting include_ambiguity=True includes those sites, for hamming, pdist, jc69, tn93, logdet and paralinear. Each code is weighted by how often the bases it represents occur in the alignment. In this alignment C occurs 7 times and T 6 times, so Y (C or T) counts as 7/13 C and 6/13 T. Against the C in s1 it is 6/13 of a difference, giving a proportion of (6/13) / 10 = 0.046 differences. R against C is a whole difference whichever base R stands for.

aln.distance_matrix(calc="pdist", include_ambiguity=True)

names	s1	s2	s3
s1	0.0000	0.1000	0.0462
s2	0.1000	0.0000	0.1000
s3	0.0462	0.1000	0.0000

The same code in both sequences counts as a match. Gaps, and codes that represent every base such as N, are still excluded.

aln.distance_matrix(calc="tn93", include_ambiguity=True)

names	s1	s2	s3
s1	0.0000	0.1086	0.0515
s2	0.1086	0.0000	0.1086
s3	0.0515	0.1086	0.0000

Note

include_ambiguity is available through distance_matrix(), not the distance calculator object described next.

Using the distance calculator directly

from cogent3 import get_distance_calculator, load_aligned_seqs

aln = load_aligned_seqs("data/primate_brca1.fasta", moltype="dna")
dist_calc = get_distance_calculator("tn93", alignment=aln)
dist_calc

<cogent3.evolve.fast_distance.TN93Pair at 0x74f7a4483a10>

dist_calc.run(show_progress=False)
dists = dist_calc.get_pairwise_distances()
dists

names	Chimpanzee	Galago	Gorilla	HowlerMon	Human	Orangutan	Rhesus
Chimpanzee	0.0000	0.1921	0.0054	0.0704	0.0089	0.0140	0.0396
Galago	0.1921	0.0000	0.1923	0.2157	0.1965	0.1944	0.1962
Gorilla	0.0054	0.1923	0.0000	0.0700	0.0086	0.0137	0.0393
HowlerMon	0.0704	0.2157	0.0700	0.0000	0.0736	0.0719	0.0736
Human	0.0089	0.1965	0.0086	0.0736	0.0000	0.0173	0.0423
Orangutan	0.0140	0.1944	0.0137	0.0719	0.0173	0.0000	0.0411
Rhesus	0.0396	0.1962	0.0393	0.0736	0.0423	0.0411	0.0000

The distance calculation object can provide more information, for instance, the standard errors.

dist_calc.stderr

Standard Error of Pairwise Distances
Seq1 \ Seq2	Galago	HowlerMon	Rhesus	Orangutan	Gorilla	Human	Chimpanzee
Galago	0	0.010274827058395896	0.009616307832648562	0.009535646532276787	0.009491382249540176	0.009615033091864917	0.009469268026590141
HowlerMon	0.010274827058395896	0	0.005411811712554772	0.0053348584951611175	0.005265612474694246	0.005406760238748984	0.005273572620183854
Rhesus	0.009616307832648562	0.005411811712554772	0	0.0039408549417865755	0.003852798161903095	0.004005045920100125	0.0038665597157698894
Orangutan	0.009535646532276787	0.0053348584951611175	0.0039408549417865755	0	0.0022291124743011375	0.0025151838791803655	0.0022606571679022955
Gorilla	0.009491382249540176	0.005265612474694246	0.003852798161903095	0.0022291124743011375	0	0.0017596919902326876	0.0013848543487237903
Human	0.009615033091864917	0.005406760238748984	0.004005045920100125	0.0025151838791803655	0.0017596919902326876	0	0.0017949285088691988
Chimpanzee	0.009469268026590141	0.005273572620183854	0.0038665597157698894	0.0022606571679022955	0.0013848543487237903	0.0017949285088691988	0

7 rows x 8 columns

Likelihood based pairwise distance estimation

The standard cogent3 likelihood function can also be used to estimate distances. Because these require numerical optimisation they can be significantly slower than the fast estimation approach above.

The following will use the F81 nucleotide substitution model and perform numerical optimisation.

from cogent3 import get_model, load_aligned_seqs
from cogent3.evolve import distance

aln = load_aligned_seqs("data/primate_brca1.fasta", moltype="dna")
d = distance.EstimateDistances(aln, submodel=get_model("F81"))
d.run(show_progress=False)
dists = d.get_pairwise_distances()
dists

names	Chimpanzee	Galago	Gorilla	HowlerMon	Human	Orangutan	Rhesus
Chimpanzee	0.0000	0.1892	0.0054	0.0697	0.0089	0.0140	0.0395
Galago	0.1892	0.0000	0.1891	0.2112	0.1934	0.1915	0.1930
Gorilla	0.0054	0.1891	0.0000	0.0693	0.0086	0.0136	0.0391
HowlerMon	0.0697	0.2112	0.0693	0.0000	0.0729	0.0713	0.0729
Human	0.0089	0.1934	0.0086	0.0729	0.0000	0.0173	0.0421
Orangutan	0.0140	0.1915	0.0136	0.0713	0.0173	0.0000	0.0410
Rhesus	0.0395	0.1930	0.0391	0.0729	0.0421	0.0410	0.0000

Get the names of sequences with max pairwise distance

Given a DistanceMatrix object, finding the sequences that have the maximum pairwise distance is achieved through the max_pair method.

from cogent3 import load_aligned_seqs

aln = load_aligned_seqs("data/primate_brca1.fasta", moltype="dna")
dists = aln.distance_matrix(calc="tn93")
dists.max_pair()

('Galago', 'HowlerMon')

To find the maximum distance, index the DistanceMatrix with the result of max_pair.

dists[dists.max_pair()]

np.float64(0.21548371016979218)

Get the names of sequences with min pairwise distance

Given a DistanceMatrix object, finding the sequences that have the minimum pairwise distance is achieved through the min_pair method.

Note

As the distance between a sequence and itself is zero, and this is not informative, min_pair will return the smallest distance not on the diagonal.

from cogent3 import load_aligned_seqs

aln = load_aligned_seqs("data/primate_brca1.fasta", moltype="dna")
dists = aln.distance_matrix(calc="tn93")
dists.min_pair()

('Gorilla', 'Chimpanzee')

To find the minimum distance, index the DistanceMatrix with the result of min_pair.

dists[dists.min_pair()]

np.float64(0.005354191176593304)

Fast pairwise distance estimation

Computing genetic distances using the Alignment object

Including ambiguity codes

Using the distance calculator directly

Likelihood based pairwise distance estimation

Get the names of sequences with max pairwise distance

Get the names of sequences with min pairwise distance

Computing genetic distances using the `Alignment` object