evospec

Tests of whether genome evolution was specific to the treatment it happened under.

Give it the mutations found in a set of evolved genomes, all descended from one ancestor, and the treatment each genome evolved under. It asks four questions:

• Dice similarity. Are genomes from the same treatment mutated in more of the same genes than genomes from different treatments? It computes the Dice similarity of the sets of mutated genes for every pair of genomes, averages it within and between treatments, and runs a randomization test that reshuffles which genome belongs to which treatment. This is done overall and, with --pairwise, for every pair of treatments.
• Signature genes. For each gene mutated more than once, are the genomes carrying those mutations concentrated in one treatment? This is a Fisher's exact test.
• Mutation counts. Did some treatments accumulate more mutations? This uses Mann-Whitney and Kruskal-Wallis tests, and can compare every treatment with a control treatment.
• Mutation types. How many SNPs, deletions, amplifications, insertions and IS insertions (by family) did each treatment accumulate, and how long were the deletions and amplifications?

It is the analysis from:

Deatherage, D.E., Kepner, J.L., Bennett, A.F., Lenski, R.E., and Barrick, J.E. (2017). Specificity of genome evolution in experimental populations of Escherichia coli evolved at different temperatures. Proc. Natl. Acad. Sci. U.S.A. 114:E1904–E1912. doi:10.1073/pnas.1616132114

Please cite it if you use this package.

evospec is the updated version of breseq-ext-specificity. That package held the script by Daniel Deatherage, Rohan Maddamsetti and Jeffrey Barrick that the paper used. The script is now a library with a command line on top of it, and the statistics can be called directly on data from anywhere. MutInt calls them this way.

Installing

pip install "evospec[cli] @ git+https://github.com/barricklab/evospec.git"

The cli extra adds Biopython, which reads the GenBank reference, and tqdm, which draws a progress bar. A program that assigns mutations to genes itself needs only numpy and SciPy: pip install "evospec @ git+https://github.com/barricklab/evospec.git".

The command line

There are two ways to lay out the input:

• -dt: the treatments are subdirectories of --samples. Each one holds the genomes as .gd files, or as breseq output directories containing evidence/annotated.gd.
• --header-treatments: every .gd sits directly in --samples, and its #=TREATMENT header line names its treatment.

evospec -s TEE-clone-curated --header-treatments -g REL606.gbk \
        -e TEE-clone-curated/REL1207.gd -x REL1206 -x REL1207 \
        -n 1000000 --pairwise -m genes.csv

• -e names mutations to leave out, usually the ancestor's.
• -x leaves a sample out entirely.
• -p assigns a mutation this many bases upstream of a gene to that gene.
• --dN_only counts only nonsynonymous SNPs as gene hits.
• -ct names a control treatment to compare every other one with.
• -m writes the genes-by-genomes matrix.

evospec --help lists everything.

The .gd files should be annotated (gdtools ANNOTATE), because a SNP's snp_type is what marks it synonymous. Without it, synonymous SNPs count as gene hits and the command warns you.

How a mutation is assigned to a gene. A mutation inside one gene is assigned to it, unless it is a synonymous SNP. With -p, a mutation within that many bases upstream of a gene's start is assigned to the nearest such gene. A mutation touching an annotated repeat region, or one that spans several genes, counts towards its genome's total but not towards any gene.

As a library

from evospec import Genome, Mutation, analyze

genomes = [
    Genome("REL2037", "32C", (Mutation("SNP", "nadR", "nonsynonymous"),
                              Mutation("DEL", None, None, 7554))),
    ...
]
result = analyze(genomes, permutations=1_000_000, seed=0, pairwise=True)

Inputs. A Mutation is (type, gene, snp_class, size, repeat_family). The gene is the one gene the mutation is assigned to, or None if it counts only towards the total.

Results. Everything comes back as plain dicts and lists, so it can be stored as JSON. dice_test, gene_fisher, count_tests and mutation_types run each analysis on its own.

Speed. The Dice matrix is computed once, and each randomization only regroups it. A million randomizations of thirty genomes take a few seconds.

Cancelling. progress(done, total) is called as randomizations complete. Raising from it stops the test.

Differences from breseq-ext-specificity

• Dice counts each gene once per genome. The old script counted a gene mutated twice in one genome twice, which could push a similarity above 1.
• A pair where either genome has no gene hits is left out of the averages. The old script scored such a pair −1 and averaged that in.
• The Bonferroni correction for pairwise Dice comparisons divides by the number of comparisons. The old script divided by the number of treatments minus one.
• Dice is recomputed without the significant genes. When any gene is individually significant, the Dice test runs a second time without those genes. This asks whether specificity remains beyond the genes that most obviously carry it. The old script's documentation described this, but it never ran it.
• Randomizations use a different random generator. They use numpy's generator with an explicit --seed, so a run can be repeated exactly. p-values agree with the old script's up to sampling noise.
• The old script cannot run on a current SciPy, and did not run as installed. It imported scipy.stats.binom_test, which was removed in SciPy 1.12, and its functions read variables its main() never made global.

Tests

pip install ".[cli]"
python -m unittest discover -s tests

License

MIT. See LICENSE.