Edit on GitHub

starplast.codons

Codon usage bias, computed from the coding sequences and nothing else.

Three numbers per gene, all sequence-derived, which is what the slot asks for:

  • ENC (Wright's effective number of codons) -- 20 when a gene uses one codon per amino acid, 61 when it uses all of them evenly. Reference-free: it is computed from the gene's own codon homozygosity, so it needs no notion of which genes are highly expressed and cannot smuggle expression into a sequence column.
  • GC3 -- GC content at synonymous third positions. The compositional half of the same story, and the thing ENC is usually plotted against.
  • CAI -- codon adaptation index against ribosomal proteins as the reference set. Ribosomal proteins are the standard choice because they are translated heavily in every organism, and using them rather than this map's own expression columns is what keeps the column from being a restatement of expr_tachy. That choice is the one judgement call here and it is declared in the registry as well as stated here.

Why this is not circular, and where it nearly was

CAI's reference set is normally "the most highly expressed genes", which for this map would mean building a sequence column out of an expression column and then being surprised when they correlate. The ribosomal set is chosen by product annotation instead. The correlation with expression that remains -- and there is one, which is the whole reason CAI is interesting -- is then a finding rather than an artefact of construction.

TABLE = 'toxodb_cds.tsv.gz'
BASES = 'TCAG'
AMINO = 'FFLLSSSSYY**CC*WLLLLPPPPHHQQRRRRIIIMTTTTNNKKSSRRVVVVAAAADDEEGGGG'
CODE = {'TTT': 'F', 'TTC': 'F', 'TTA': 'L', 'TTG': 'L', 'TCT': 'S', 'TCC': 'S', 'TCA': 'S', 'TCG': 'S', 'TAT': 'Y', 'TAC': 'Y', 'TAA': '*', 'TAG': '*', 'TGT': 'C', 'TGC': 'C', 'TGA': '*', 'TGG': 'W', 'CTT': 'L', 'CTC': 'L', 'CTA': 'L', 'CTG': 'L', 'CCT': 'P', 'CCC': 'P', 'CCA': 'P', 'CCG': 'P', 'CAT': 'H', 'CAC': 'H', 'CAA': 'Q', 'CAG': 'Q', 'CGT': 'R', 'CGC': 'R', 'CGA': 'R', 'CGG': 'R', 'ATT': 'I', 'ATC': 'I', 'ATA': 'I', 'ATG': 'M', 'ACT': 'T', 'ACC': 'T', 'ACA': 'T', 'ACG': 'T', 'AAT': 'N', 'AAC': 'N', 'AAA': 'K', 'AAG': 'K', 'AGT': 'S', 'AGC': 'S', 'AGA': 'R', 'AGG': 'R', 'GTT': 'V', 'GTC': 'V', 'GTA': 'V', 'GTG': 'V', 'GCT': 'A', 'GCC': 'A', 'GCA': 'A', 'GCG': 'A', 'GAT': 'D', 'GAC': 'D', 'GAA': 'E', 'GAG': 'E', 'GGT': 'G', 'GGC': 'G', 'GGA': 'G', 'GGG': 'G'}
SYNONYMS = defaultdict(<class 'list'>, {'F': ['TTT', 'TTC'], 'L': ['TTA', 'TTG', 'CTT', 'CTC', 'CTA', 'CTG'], 'S': ['TCT', 'TCC', 'TCA', 'TCG', 'AGT', 'AGC'], 'Y': ['TAT', 'TAC'], 'C': ['TGT', 'TGC'], 'W': ['TGG'], 'P': ['CCT', 'CCC', 'CCA', 'CCG'], 'H': ['CAT', 'CAC'], 'Q': ['CAA', 'CAG'], 'R': ['CGT', 'CGC', 'CGA', 'CGG', 'AGA', 'AGG'], 'I': ['ATT', 'ATC', 'ATA'], 'M': ['ATG'], 'T': ['ACT', 'ACC', 'ACA', 'ACG'], 'N': ['AAT', 'AAC'], 'K': ['AAA', 'AAG'], 'V': ['GTT', 'GTC', 'GTA', 'GTG'], 'A': ['GCT', 'GCC', 'GCA', 'GCG'], 'D': ['GAT', 'GAC'], 'E': ['GAA', 'GAG'], 'G': ['GGT', 'GGC', 'GGA', 'GGG']})
DEGENERATE = {'F': ['TTT', 'TTC'], 'L': ['TTA', 'TTG', 'CTT', 'CTC', 'CTA', 'CTG'], 'S': ['TCT', 'TCC', 'TCA', 'TCG', 'AGT', 'AGC'], 'Y': ['TAT', 'TAC'], 'C': ['TGT', 'TGC'], 'P': ['CCT', 'CCC', 'CCA', 'CCG'], 'H': ['CAT', 'CAC'], 'Q': ['CAA', 'CAG'], 'R': ['CGT', 'CGC', 'CGA', 'CGG', 'AGA', 'AGG'], 'I': ['ATT', 'ATC', 'ATA'], 'T': ['ACT', 'ACC', 'ACA', 'ACG'], 'N': ['AAT', 'AAC'], 'K': ['AAA', 'AAG'], 'V': ['GTT', 'GTC', 'GTA', 'GTG'], 'A': ['GCT', 'GCC', 'GCA', 'GCG'], 'D': ['GAT', 'GAC'], 'E': ['GAA', 'GAG'], 'G': ['GGT', 'GGC', 'GGA', 'GGG']}
def translate(cds: str) -> str:

A coding sequence to its protein, in frame 1, with the terminator trimmed.

Here because a phosphosite table can key on nothing this project's identity index carries -- the CDPK1 supplement reports numeric ids from a 2017 annotation -- while giving the 15-residue window around each site. A window is an identifier if it occurs in exactly one protein, and checking that requires the proteome the CDS table already ships. Anything that is not a clean triplet of ACGT becomes X, so a match is never made across a gap that was silently closed.

def codon_counts(cds: str) -> dict:

Codons in frame, skipping anything that is not a clean triplet of ACGT.

def gc3(counts: dict) -> float:

GC at the third position of degenerate codons. NaN when there are none to look at.

def enc(counts: dict) -> float:

Wright's effective number of codons, from per-family homozygosity.

Families with too few observations to estimate homozygosity are dropped and their class average is used in their place -- Wright's own prescription. A gene with no usable family at all gets NaN rather than 61, because "we could not tell" is not "unbiased".

def reference_weights(counts_by_gene: dict) -> dict:

Relative adaptiveness per codon, from a pooled reference set of genes.

def cai(counts: dict, weights: dict) -> float:

Geometric mean of relative adaptiveness over a gene's degenerate codons.

def codon_usage( base: str, resolve=None, reference=None, log=<built-in function print>, table: str | None = None, nodes: str | None = None) -> pandas.DataFrame:

ENC, GC3 and CAI per gene.

reference is the gene ids of the CAI reference set. Left None it is taken from the product descriptions in the shipped node table -- the ribosomal proteins -- and if that table is not present the CAI column is left out rather than computed against an arbitrary set.

table and nodes name the CDS file and the node table to read, so the Plasmodium arm can use this construction instead of writing its own. Sharing it is the point: ENC, GC3 and a CAI against ribosomal proteins are definitions, and two arms computing them differently would make a difference between the arms unreadable.