natsel_zhang – a branch-site test

Note

These docs now use the new_type core objects via the following setting.

import os

# using new types without requiring an explicit argument
os.environ["COGENT3_NEW_TYPE"] = "1"

This is the hypothesis test presented in Zhang et al. It evaluates the hypothesis that a set of sites have undergone positive natural selection on a pre-specified set of lineages.

For this model class, there are groups of branches for which all positions are evolving neutrally but some proportion of those neutrally evolving sites change to adaptively evolving on so-called foreground edges. For the current example, we’ll define the Chimpanzee and Human branches as foreground and everything else as background. The following table defines the parameter scopes.

<IPython.core.display.HTML object>

Note

Our implementation is not as parametrically succinct as that of Zhang et al, we have 1 additional bin probability.

from cogent3 import get_app

loader = get_app("load_aligned", format="fasta", moltype="dna")
aln = loader("data/primate_brca1.fasta")

zhang_test = get_app("natsel_zhang",
    "GNC",
    tree="data/primate_brca1.tree",
    optimise_motif_probs=False,
    tip1="Human",
    tip2="Chimpanzee",
)

result = zhang_test(aln)
result

Statistics

LR	df	pvalue
4.9553	3	0.1751

hypothesis	key	lnL	nfp	DLC	unique_Q
null	'GNC-null'	-6708.3129	24	True	True
alt	'GNC-alt'	-6705.8352	27	True	True

result.alt.lf

GNC-alt

log-likelihood = -6705.8352

number of free parameters = 27

Global params

A>C	A>G	A>T	C>A	C>G	C>T	G>A	G>C	G>T	T>A	T>C
0.86	3.53	0.97	1.66	2.19	6.26	8.01	1.24	0.79	1.27	2.96

Bin params

bin	bprobs
0	0.05
1	0.01
2a	0.04
2b	0.89

Edge params

edge	parent	length
Galago	root	0.54
HowlerMon	root	0.14
Rhesus	edge.3	0.06
Orangutan	edge.2	0.02
Gorilla	edge.1	0.01
Human	edge.0	0.02
Chimpanzee	edge.0	0.01
edge.0	edge.1	0.00
edge.1	edge.2	0.01
edge.2	edge.3	0.04
edge.3	root	0.02

Edge bin params

edge	bin	omega
Galago	0	0.00
Galago	1	1.00
Galago	2a	0.00
Galago	2b	1.00
HowlerMon	0	0.00
HowlerMon	1	1.00
HowlerMon	2a	0.00
HowlerMon	2b	1.00
Rhesus	0	0.00
Rhesus	1	1.00
Rhesus	2a	0.00
Rhesus	2b	1.00
Orangutan	0	0.00
Orangutan	1	1.00
Orangutan	2a	0.00
Orangutan	2b	1.00
Gorilla	0	0.00
Gorilla	1	1.00
Gorilla	2a	0.00
Gorilla	2b	1.00
Human	0	0.00
Human	1	1.00
Human	2a	3.76
Human	2b	3.76
Chimpanzee	0	0.00
Chimpanzee	1	1.00
Chimpanzee	2a	3.76
Chimpanzee	2b	3.76
edge.0	0	0.00
edge.0	1	1.00
edge.0	2a	0.00
edge.0	2b	1.00
edge.1	0	0.00
edge.1	1	1.00
edge.1	2a	0.00
edge.1	2b	1.00
edge.2	0	0.00
edge.2	1	1.00
edge.2	2a	0.00
edge.2	2b	1.00
edge.3	0	0.00
edge.3	1	1.00
edge.3	2a	0.00
edge.3	2b	1.00

Motif params

AAA	AAC	AAG	AAT	ACA	ACC	ACG	ACT	AGA	AGC	AGG	AGT	ATA
0.06	0.02	0.03	0.06	0.02	0.00	0.00	0.03	0.02	0.03	0.01	0.04	0.02

continuation

ATC	ATG	ATT	CAA	CAC	CAG	CAT	CCA	CCC	CCG	CCT	CGA	CGC
0.01	0.01	0.02	0.02	0.01	0.02	0.02	0.02	0.01	0.00	0.03	0.00	0.00

continuation

CGG	CGT	CTA	CTC	CTG	CTT	GAA	GAC	GAG	GAT	GCA	GCC	GCG
0.00	0.00	0.01	0.01	0.01	0.01	0.08	0.01	0.03	0.03	0.02	0.01	0.00

continuation

GCT	GGA	GGC	GGG	GGT	GTA	GTC	GTG	GTT	TAC	TAT	TCA	TCC
0.01	0.02	0.01	0.01	0.01	0.01	0.01	0.01	0.02	0.00	0.01	0.02	0.01

continuation

TCG	TCT	TGC	TGG	TGT	TTA	TTC	TTG	TTT
0.00	0.03	0.00	0.00	0.02	0.02	0.01	0.01	0.02

Getting the posterior probabilities of site-class membership

bprobs = result.alt.lf.get_bin_probs()
bprobs[:, :20]

	0	1	2	3	4	5	6	7	8	9	10	11	12	13	14	15	16	17	18	19
0	0.0723	0.0407	0.0000	0.0638	0.0558	0.0763	0.0409	0.0579	0.0494	0.0391	0.0373	0.0763	0.0457	0.0000	0.0759	0.2494	0.0391	0.0338	0.0558	0.0590
1	0.0117	0.0124	0.0135	0.0119	0.0121	0.0116	0.0124	0.0120	0.0122	0.0125	0.0125	0.0116	0.0123	0.0135	0.0116	0.0072	0.0125	0.0126	0.0121	0.0120
2a	0.0603	0.0349	0.0000	0.0536	0.0471	0.0635	0.0351	0.0488	0.0420	0.0337	0.0321	0.0635	0.0390	0.0000	0.0632	0.2146	0.0337	0.0292	0.0471	0.0498
2b	0.8556	0.9120	0.9865	0.8707	0.8850	0.8486	0.9116	0.8812	0.8964	0.9148	0.9180	0.8486	0.9030	0.9865	0.8492	0.5287	0.9148	0.9244	0.8850	0.8792

Getting all the statistics in tabular form

tab = get_app("tabulate_stats")
stats = tab(result.alt)
stats

5x tabular_result('global params': Table, 'bin params': Table, 'edge params': Table, 'edge bin params': Table, 'motif params': Table)

stats["edge bin params"][:10]  # truncating the table

edge bin params

edge	bin	omega
Galago	0	0.00
Galago	1	1.00
Galago	2a	0.00
Galago	2b	1.00
HowlerMon	0	0.00
HowlerMon	1	1.00
HowlerMon	2a	0.00
HowlerMon	2b	1.00
Rhesus	0	0.00
Rhesus	1	1.00

10 rows x 3 columns