Analyse d’expression génique (RNA-seq)

Le RNA-seq donne, après alignement et comptage, une matrice de comptages : combien de reads se sont alignés sur chaque gène, pour chaque échantillon. C’est la donnée primaire de toute étude transcriptomique.

Voici un mini-jeu de données pédagogique : 10 gènes oncologiques, mesurés dans 3 contrôles (ctrl_1 à ctrl_3) et 3 traités (treat_1 à treat_3). Sauvegardez-le dans comptages.tsv :

gene	ctrl_1	ctrl_2	ctrl_3	treat_1	treat_2	treat_3
ENSG_BRCA1	820	790	850	1240	1180	1310
ENSG_TP53	450	480	470	920	880	950
ENSG_GAPDH	2100	2150	2080	2080	2120	2110
ENSG_MYC	320	310	340	130	120	140
ENSG_KRAS	180	170	190	60	70	80
ENSG_BAX	500	520	490	790	810	780
ENSG_BCL2	600	620	590	240	260	220
ENSG_CDK2	380	360	400	380	390	410
ENSG_RB1	290	310	280	295	300	305
ENSG_VEGFA	150	160	140	450	470	440

Premier regard sur les données : pour chaque gène, quel est son niveau d’expression moyen ? Est-il très variable d’un échantillon à l’autre ? Cette première inspection révèle déjà les gènes très exprimés (housekeeping comme GAPDH) et les gènes faiblement détectés (à filtrer plus tard).

On ne peut pas comparer directement les comptages bruts entre deux échantillons. Pourquoi ? Parce que l’échantillon A peut avoir été séquencé deux fois plus profondément que B — son gène GAPDH aura mécaniquement deux fois plus de reads, sans changement biologique réel.

Solution la plus simple : normaliser par la taille de librairie (total de reads par échantillon), exprimée en CPM (Counts Per Million) :

CPM(gene, échantillon) = comptage_gene / total_échantillon * 1_000_000

(Méthode rudimentaire mais didactique. En vrai, on utilise TPM ou DESeq2/TMM pour des analyses sérieuses.)

Un gène détecté à 2-3 reads dans seulement un échantillon n’apporte qu’un bruit statistique. Avant le calcul du fold-change, on élimine les gènes faibles.

Règle classique (edgeR et consorts) : un gène est gardé s’il atteint au moins N CPM dans au moins M échantillons.

Le fold-change mesure de combien un gène est surexprimé (ou sous-exprimé) entre deux conditions :

FC = moyenne(traité) / moyenne(contrôle)

On prend ensuite le log2 : log2(2) = 1 (2× plus), log2(4) = 2 (4× plus), log2(0.5) = −1 (2× moins). Cette transformation rend la distribution symétrique autour de 0.

Astuce numérique : on ajoute un pseudo-comptage (+1) pour éviter log2(0).

Un DEG (Differentially Expressed Gene) est, dans notre version simplifiée, un gène dont |log2FC| ≥ 1 — soit au moins 2× de différence entre les deux conditions.

Une vraie analyse ajouterait un test statistique (t-test, Wald via DESeq2) puis une correction multiple (Benjamini-Hochberg FDR). On reste sur le critère du fold-change pour rester en Python pur.

Une heatmap est l’outil de visualisation classique d’une analyse d’expression : chaque ligne = un gène, chaque colonne = un échantillon, intensité de couleur = niveau d’expression. On peut faire une version terminal-only avec des blocs Unicode █ — utile sur un serveur sans interface graphique.

Avant l’affichage, on z-scale par ligne : chaque gène est centré-réduit, ce qui rend visible le motif de variation indépendamment du niveau absolu d’expression.

On a sept briques. Reste à les assembler en un outil en ligne de commande qui prend une matrice de comptages, deux listes d’échantillons (contrôle et traité), et qui produit un rapport complet.