Site avoidance: restriction sites a genome avoids or favors

Ask which restriction sites a set of genomes carries far less often, or far more often, than chance predicts, with the false discovery rate controlled across every genome and every site at once.

What it doesCopy Anchor Link

Tick the records that make up the corpus (a set of bacterial genomes, phages, plasmids, viruses or any DNA you have), keep or drop the record open in the editor, choose the null that "chance" means, and run. Every exact recognition site in the NEB catalog is counted on both strands of every record and set against what the null predicts for that record. The result lists the sites each record avoids and the ones it enriches, most extreme first, each with its observed and expected count, the ratio between them and its q-value.

When you’d use thisCopy Anchor Link

• Reading a genome's history. Bacteria and phages that live with a restriction–modification system lose the sites it cuts. A site avoided across a lineage is a footprint of the systems it has met; E. coli's scarcity of CTAG-containing sites is the textbook case.
• Choosing an enzyme a design can rely on being rare. A site a chassis avoids is a safe cloning site in that chassis.
• Checking a finding is not composition. A genome depleted of CpG avoids every site that contains CG; the conditioned nulls say whether a site is rare beyond that.

What this computesCopy Anchor Link

For every record in the corpus and every exact recognition site of the catalog (degenerate sites are listed and left out):

• observed count, on both strands: a site and its reverse complement, once for a palindrome
• expected count under the chosen null, and the observed-to-expected ratio
• lower-tail and upper-tail p-values: is the site rarer, or commoner, than the null predicts
• q-values for both tails, adjusted across every (record, site) test of the run
• avoided and enriched calls: the tests whose q is at or below the false discovery rate you set
• site positions and strands, when asked for, so each occurrence can be found on the map

Three nulls, each a different meaning of "chance":

• base composition (order 0): independent bases at the record's own frequencies, with Poisson tails. The quickest look, and the one a CpG- or GC-skewed genome fools
• k-let-preserving shuffles: each record shuffled many times keeping every k-let count exact (uShuffle), with empirical tails (1 + shuffles at or beyond the observed count) / (shuffles + 1) and a 95% Clopper–Pearson interval on each tail's Monte Carlo error. At k = 2 every dinucleotide, CpG included, is held in the null
• maximal-order Markov model (R'MES): a site of length k is expected from the record's own (k−1)-mers given their (k−2)-mer cores, the site and its reverse complement scored as one family, with the R'MES score (observed − expected) / √variance, its asymptotic variance, and Gaussian tails on the score. The standard of published avoidance analyses

Two corrections for the false discovery rate:

• Benjamini–Hochberg, valid for independent or positively dependent tests
• Benjamini–Yekutieli, valid under any dependence (overlapping sites on one sequence included), with less power

Every result states, in words generated from the run: the null it stands on and what that null does not hold; the correction and why it is valid; and the caveats the data earned, such as avoided sites that carry CpG under base composition, tests whose expected count is too small for even zero occurrences to be significant, calls that sit at the shuffle null's smallest possible p, and tests the maximal-order null cannot score.

Worked exampleCopy Anchor Link

The eleven bundled genomes and genome segments (three E. coli segments, phage λ, SARS-CoV-2, SARS coronavirus, the Kryptoperidinium foliaceum chloroplast, two human chromosome segments and the human and orangutan mitochondria), all ticked, against the catalog's 101 exact sites of at least four bases: 1,111 tests.

• Base composition, Benjamini–Hochberg, q ≤ 0.05: 305 tests avoided and 181 enriched. AvrII CCTAGG is absent from the 250 kb E. coli K-12 segment where 61.9 are expected (q 5.5 × 10⁻²⁶), and SpeI ACTAGT from phage λ where 11.9 are expected (q 4.4 × 10⁻⁵). In the human segments every top hit carries CpG.
• Maximal-order Markov model, same settings: 150 avoided and 38 enriched. AvrII in E. coli and SpeI in λ are no longer called (λ SpeI: 0 against 0.31 expected given its own 5-mers): their rarity is the CTAG 4-mer's, which their 5-mers already carry. KasI/NarI GGCGCC is the strongest E. coli avoidance that survives, 6 of 118.2 expected (score −13.9).

LimitationsCopy Anchor Link

• Only exact recognition sites are swept; enzymes with degenerate sites are listed, not tested.
• The maximal-order null cannot score a site whose (k−2)-mer core never occurs in a record; that test has no p-value and is left out of the correction.
• The shuffle null's smallest p is 1 / (shuffles + 1). With few shuffles and many tests nothing can pass the correction; the result says so.
• A short record cannot show avoidance of a long site: when fewer than −ln α occurrences are expected, zero is not significant.
• Avoidance is evidence, not a mechanism: which restriction system caused it is a separate question.

Where the numbers come fromCopy Anchor Link

Counts and order-0 expectations are checked against seqinr and R's ppois; the shuffles against uShuffle's own C, shuffle for shuffle; the maximal-order expectations, variances and scores against R'MES 3.1.0, bit for bit; and the corrections against R's p.adjust. Each is checked on every site of the eleven genomes above.

AutomationCopy Anchor Link

The same sweep is a workflow node and a tool the assistant can call, with the same settings and defaults as the panel; the corpus is the list of record names.

RelatedCopy Anchor Link

In this group: Properties · BLAST Search · Search Databases · EMBL-EBI Services · Protein Domains · Protein Analysis · Protease Digest · Antibody Numbering · Mass spectrometry: peptide masses, b ladders, isotopes, spectrum scoring, q-values, protein inference · Codon Usage Analysis · Codon outliers: genes whose codon usage departs from their own genome · Tree incongruence: genes whose own alignment rejects the species tree · Reference stability: reference genes against genes of their own abundance · Annotation contradictions: UniProt features a predictor scores unlike confirmed ones · Hypothesis ledger: pre-register a hypothesis, attach runs, decide with e-BH · Reagent audit: does each annotated CDS code for what the map says, and are its cryptic promoters, splice sites and repeats worse than chance · Assay decay: does a published assay still match current genomes, and is the mismatch at the 3′ end · Gene incongruence: which loci of a multi-locus alignment have a different history · Target validation: test "gene G drives phenotype P" against gene-level evidence tables · Codon Tables · Enzyme Collections · Virtual Gel Bench · Phylogenetics · Dot Plot · Table · NumPy Array · Structure · PDF · Molecule · Image