{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Prédiction de la charge ribosomique moyenne avec un réseau neuronal\n", "\n", "convolutif 1D\n", "\n", "CSI 4506 — Introduction à l’intelligence artificielle\n", "\n", "Marcel Turcotte \n", "Version : 31 août 2026, 18h44\n", "\n", "# Introduction\n", "\n", "Les réseaux neuronaux convolutifs sont le plus souvent présentés à\n", "l’aide d’images. Les images fournissent des exemples visuels\n", "convaincants, mais l’entraînement de classificateurs d’images peut\n", "exiger des ressources de calcul considérables. Dans ce notebook, nous\n", "étudions un problème unidimensionnel plus petit pour lequel la\n", "convolution est tout aussi naturelle.\n", "\n", "Notre tâche d’apprentissage automatique peut s’énoncer sans\n", "connaissances préalables en biologie :\n", "\n", "> Étant donné une chaîne de longueur 50 sur l’alphabet `{A, C, G, T}`,\n", "> prédire une valeur continue qui se situe généralement entre 0 et 13.\n", "\n", "Les chaînes sont des régions non traduites en 5′ (5′ UTR), et la\n", "variable cible est la charge ribosomique moyenne (*mean ribosome load*,\n", "MRL). Une 5′ UTR est un segment d’une molécule d’ARN qui précède la\n", "région codant une protéine. La MRL résume le nombre de ribosomes\n", "associés à l’ARN et fournit donc une mesure expérimentale liée à la\n", "traduction.\n", "\n", "L’expérience est adaptée de Sample et coll. et reprend des idées\n", "d’architecture employées dans des travaux ultérieurs, notamment ceux de\n", "Tang et coll. Notre stratégie de partitionnement est volontairement\n", "différente. Après avoir retiré les mesures de faible couverture, nous\n", "répartissons les données aléatoirement afin que l’ensemble de test\n", "provienne de la même population filtrée que l’ensemble d’entraînement.\n", "\n", "> **Pourquoi une convolution unidimensionnelle?**\n", ">\n", "> Un court noyau de convolution examine une sous-chaîne locale à chaque\n", "> position. Les mêmes poids appris sont réutilisés sur toute la\n", "> séquence. Un CNN 1D possède ainsi les deux mêmes biais inductifs\n", "> importants qu’un CNN appliqué aux images : la **connectivité locale**\n", "> et le **partage des paramètres**.\n", "\n", "## Objectifs d’apprentissage\n", "\n", "À la fin de ce notebook, vous devriez pouvoir :\n", "\n", "- représenter des chaînes de longueur fixe par des tenseurs numériques à\n", " l’aide d’un encodage à chaud (*one-hot encoding*);\n", "- expliquer comment une couche `Conv1D` traite une séquence;\n", "- distinguer les données d’entraînement, de validation et de test;\n", "- utiliser l’arrêt précoce pour contrôler l’entraînement d’un réseau\n", " neuronal;\n", "- interpréter les courbes d’apprentissage de la perte, de la corrélation\n", " de Pearson et de $R^2$;\n", "- expliquer comment des prédictions peuvent être fortement corrélées aux\n", " observations tout en demeurant mal calibrées.\n", "\n", "# Préparation\n", "\n", "Nous utilisons Keras avec TensorFlow comme moteur. NumPy et pandas\n", "servent au traitement des données, tandis que Matplotlib et seaborn\n", "permettent de les visualiser. Scikit-learn et SciPy fournissent des\n", "fonctions usuelles de prétraitement et d’évaluation." ], "id": "836b172a-96f9-47c5-a713-e277c506cc5a" }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [], "source": [ "from contextlib import redirect_stderr, redirect_stdout\n", "from io import StringIO\n", "from pathlib import Path\n", "\n", "import numpy as np\n", "import pandas as pd\n", "\n", "import matplotlib.pyplot as plt\n", "import seaborn as sns\n", "\n", "from scipy.stats import pearsonr, spearmanr\n", "from sklearn.metrics import mean_squared_error, r2_score\n", "from sklearn.preprocessing import StandardScaler\n", "\n", "import tensorflow as tf\n", "from tensorflow import keras\n", "from tensorflow.keras import callbacks, layers" ], "id": "imports" }, { "cell_type": "markdown", "metadata": {}, "source": [ "Tous les choix importants sont regroupés dans une seule cellule de\n", "configuration. Cette organisation facilite l’examen et la reproduction\n", "de l’expérience." ], "id": "03accff8-e05b-401c-9025-553f1d74106f" }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [ { "output_type": "stream", "name": "stdout", "text": [ "TensorFlow 2.21.0; dispositif d'entraînement : CPU" ] } ], "source": [ "SEED = 42\n", "SEQUENCE_LENGTH = 50\n", "\n", "N_RETAINED = 275_000\n", "N_TRAIN = 220_000\n", "N_VALIDATION = 27_500\n", "N_TEST = 27_500\n", "\n", "BATCH_SIZE = 256\n", "MAX_EPOCHS = 50\n", "PATIENCE = 5\n", "\n", "# Des limites communes rendent les graphiques de prédiction comparables.\n", "MRL_PLOT_LIMITS = (-0.5, 13.5)\n", "\n", "DATA_DIR = Path(\"data\")\n", "MODEL_DIR = Path(\"models\")\n", "DATA_DIR.mkdir(exist_ok=True)\n", "MODEL_DIR.mkdir(exist_ok=True)\n", "\n", "DATASETS = {\n", " \"GSM3130435_egfp_unmod_1.csv.gz\": (\n", " \"https://ftp.ncbi.nlm.nih.gov/geo/samples/GSM3130nnn/\"\n", " \"GSM3130435/suppl/GSM3130435_egfp_unmod_1.csv.gz\"\n", " ),\n", " \"GSE232927_processed_defined_end_tcell_r1.csv.gz\": (\n", " \"https://ftp.ncbi.nlm.nih.gov/geo/series/GSE232nnn/\"\n", " \"GSE232927/suppl/\"\n", " \"GSE232927_processed_defined_end_tcell_r1.csv.gz\"\n", " ),\n", " \"GSE232927_processed_defined_end_tcell_r2.csv.gz\": (\n", " \"https://ftp.ncbi.nlm.nih.gov/geo/series/GSE232nnn/\"\n", " \"GSE232927/suppl/\"\n", " \"GSE232927_processed_defined_end_tcell_r2.csv.gz\"\n", " ),\n", "}\n", "\n", "assert N_TRAIN + N_VALIDATION + N_TEST == N_RETAINED\n", "\n", "keras.utils.set_random_seed(SEED)\n", "sns.set_theme(style=\"ticks\", context=\"notebook\")\n", "\n", "accelerators = tf.config.list_physical_devices(\"GPU\")\n", "device_description = \"GPU\" if accelerators else \"CPU\"\n", "print(f\"TensorFlow {tf.__version__}; dispositif d'entraînement : {device_description}\")" ], "id": "configuration" }, { "cell_type": "markdown", "metadata": {}, "source": [ "La fonction auxiliaire ci-dessous télécharge un jeu de données seulement\n", "s’il n’est pas déjà disponible. `keras.utils.get_file` renvoie aussi le\n", "chemin local, que nous pouvons transmettre directement à pandas.\n", "L’affichage de la progression est masqué par défaut, car les caractères\n", "de contrôle du terminal peuvent produire une très longue sortie dans un\n", "notebook généré. Le paramètre `show_progress=True` rétablit la barre de\n", "progression interactive." ], "id": "6288ae7a-abef-4aca-babd-a224cf50de95" }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [], "source": [ "def download_dataset(filename, show_progress=False):\n", " \"\"\"Télécharge un jeu de données configuré et renvoie son chemin local.\"\"\"\n", " if filename not in DATASETS:\n", " raise KeyError(f\"Jeu de données inconnu : {filename}\")\n", "\n", " arguments = {\n", " \"fname\": filename,\n", " \"origin\": DATASETS[filename],\n", " \"cache_dir\": str(Path.cwd()),\n", " \"cache_subdir\": str(DATA_DIR),\n", " }\n", "\n", " if show_progress:\n", " path = keras.utils.get_file(**arguments)\n", " else:\n", " with redirect_stdout(StringIO()), redirect_stderr(StringIO()):\n", " path = keras.utils.get_file(**arguments)\n", "\n", " return Path(path)" ], "id": "download-helper" }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Chargement et validation du jeu de données principal\n", "\n", "Les fichiers sources utilisent des noms de colonnes légèrement\n", "différents. Nous les uniformisons à l’entrée de notre programme afin que\n", "toutes les fonctions suivantes puissent employer le même vocabulaire\n", "simple : `utr`, `mrl` et `read_count`." ], "id": "d59a4985-828b-4ff0-b220-496fcd982d5b" }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [], "source": [ "def load_utr_dataset(path, expected_length=SEQUENCE_LENGTH):\n", " \"\"\"Charge et valide un jeu de données UTR traité.\n", "\n", " Le DataFrame renvoyé contient toujours les colonnes ``utr``, ``mrl`` et\n", " ``read_count``.\n", " \"\"\"\n", " data = pd.read_csv(path)\n", "\n", " if \"utr\" not in data.columns and \"UTR\" in data.columns:\n", " data = data.rename(columns={\"UTR\": \"utr\"})\n", "\n", " if \"total_reads\" in data.columns:\n", " read_column = \"total_reads\"\n", " elif \"total\" in data.columns:\n", " read_column = \"total\"\n", " else:\n", " raise ValueError(\"Le jeu de données ne contient aucune colonne total_reads ou total.\")\n", "\n", " required = {\"utr\", \"rl\", read_column}\n", " missing = required.difference(data.columns)\n", " if missing:\n", " raise ValueError(f\"Colonnes requises manquantes : {sorted(missing)}\")\n", "\n", " data = (\n", " data.loc[:, [\"utr\", \"rl\", read_column]]\n", " .rename(columns={\"rl\": \"mrl\", read_column: \"read_count\"})\n", " .copy()\n", " )\n", "\n", " data[\"utr\"] = data[\"utr\"].astype(\"string\").str.upper()\n", " data[\"mrl\"] = pd.to_numeric(data[\"mrl\"], errors=\"coerce\")\n", " data[\"read_count\"] = pd.to_numeric(\n", " data[\"read_count\"], errors=\"coerce\"\n", " )\n", "\n", " if data[[\"utr\", \"mrl\", \"read_count\"]].isna().any().any():\n", " raise ValueError(\"Des valeurs manquantes ou non numériques ont été trouvées.\")\n", "\n", " valid_sequence = data[\"utr\"].str.fullmatch(\n", " rf\"[ACGTN]{{{expected_length}}}\"\n", " )\n", " if not valid_sequence.all():\n", " n_invalid = int((~valid_sequence).sum())\n", " raise ValueError(\n", " f\"{n_invalid} séquences ne sont pas des chaînes de longueur \"\n", " f\"{expected_length} sur A, C, G, T et N.\"\n", " )\n", "\n", " if (data[\"read_count\"] < 0).any():\n", " raise ValueError(\"Le nombre de lectures doit être positif ou nul.\")\n", "\n", " if data[\"utr\"].duplicated().any():\n", " raise ValueError(\"Le jeu de données contient des séquences UTR en double.\")\n", "\n", " return data.reset_index(drop=True)" ], "id": "load-helper" }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "output_type": "stream", "name": "stdout", "text": [ "326,033 séquences chargées." ] }, { "output_type": "display_data", "metadata": {}, "data": { "text/html": [ "\n", "" ] } } ], "source": [ "main_path = download_dataset(\"GSM3130435_egfp_unmod_1.csv.gz\")\n", "main_data = load_utr_dataset(main_path)\n", "\n", "print(f\"{len(main_data):,} séquences chargées.\")\n", "main_data.head()" ], "id": "cell-load-main-data" }, { "cell_type": "markdown", "metadata": {}, "source": [ "La cible et le nombre de lectures sont mesurés sur des échelles très\n", "différentes. Un résumé compact nous donne un premier aperçu des deux\n", "variables." ], "id": "9a95a353-fa6a-44c4-a4f8-9e496da0a552" }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [ { "output_type": "display_data", "metadata": {}, "data": { "text/html": [ "\n", "" ] } } ], "source": [ "main_data[[\"mrl\", \"read_count\"]].describe(\n", " percentiles=[0.01, 0.25, 0.50, 0.75, 0.99]\n", ").round(3)" ], "id": "cell-summarize-main-data" }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Filtrage des mesures selon le nombre de lectures\n", "\n", "Les mesures fondées sur un très petit nombre de lectures sont moins\n", "précises. Nous conservons donc exactement les 275 000 séquences ayant\n", "les nombres de lectures les plus élevés." ], "id": "5997921c-c9c5-4c2b-9c50-8b5204dc1a6e" }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [ { "output_type": "stream", "name": "stdout", "text": [ "Seuil du nombre de lectures : 294\n", "Séquences conservées : 275,000 (84.3%)" ] } ], "source": [ "def select_top_by_reads(data, n):\n", " \"\"\"Renvoie les n lignes les plus lues et le nombre de lectures à la limite.\"\"\"\n", " if not 0 < n <= len(data):\n", " raise ValueError(\"n doit être compris entre 1 et le nombre de lignes.\")\n", "\n", " ranked = data.sort_values(\n", " \"read_count\", ascending=False, kind=\"stable\"\n", " ).reset_index(drop=True)\n", " cutoff = float(ranked.loc[n - 1, \"read_count\"])\n", " return ranked.iloc[:n].copy(), cutoff\n", "\n", "\n", "retained_data, read_cutoff = select_top_by_reads(main_data, N_RETAINED)\n", "retained_fraction = N_RETAINED / len(main_data)\n", "\n", "print(f\"Seuil du nombre de lectures : {read_cutoff:,.0f}\")\n", "print(f\"Séquences conservées : {N_RETAINED:,} ({retained_fraction:.1%})\")" ], "id": "select-top-helper" }, { "cell_type": "markdown", "metadata": {}, "source": [ "La distribution du nombre de lectures est fortement asymétrique vers la\n", "droite. Nous traçons $\\log_{10}(\\text{nombre de lectures}+1)$ afin que\n", "les mesures de faible et de forte couverture demeurent visibles. La\n", "ligne rouge correspond au nombre de lectures de la 275 000e séquence du\n", "classement." ], "id": "b5e46d5a-8ac5-4abd-8a2d-1bed64680254" }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [ { "output_type": "display_data", "metadata": {}, "data": {} } ], "source": [ "log_read_counts = np.log10(main_data[\"read_count\"].to_numpy() + 1)\n", "log_cutoff = np.log10(read_cutoff + 1)\n", "\n", "fig, ax = plt.subplots(figsize=(9, 4.5))\n", "sns.histplot(log_read_counts, bins=80, color=\"steelblue\", ax=ax)\n", "ax.axvline(log_cutoff, color=\"red\", linewidth=2)\n", "ax.text(\n", " log_cutoff,\n", " ax.get_ylim()[1] * 0.92,\n", " f\" seuil = {read_cutoff:,.0f} lectures\",\n", " color=\"red\",\n", " va=\"top\",\n", ")\n", "ax.set(\n", " xlabel=r\"$\\log_{10}(\\mathrm{nombre\\ de\\ lectures}+1)$\",\n", " ylabel=\"Nombre de séquences\",\n", " title=\"Couverture avant le filtrage\",\n", ")\n", "sns.despine()\n", "plt.tight_layout()\n", "plt.show()" ], "id": "cell-plot-read-histogram" }, { "cell_type": "markdown", "metadata": {}, "source": [ "> **Note**\n", ">\n", "> Comme les nombres de lectures sont entiers, plusieurs séquences\n", "> peuvent être à égalité au seuil. Nous utilisons un classement stable\n", "> et conservons exactement 275 000 lignes. La ligne rouge représente\n", "> donc la valeur limite, et non une règle permettant de départager les\n", "> égalités.\n", "\n", "# Représentation des chaînes de nucléotides\n", "\n", "Un réseau neuronal manipule des nombres plutôt que des caractères. Nous\n", "représentons chaque nucléotide par quatre valeurs binaires :\n", "\n", "| Nucléotide | A | C | G | T |\n", "|:------------|----:|----:|----:|----:|\n", "| A | 1 | 0 | 0 | 0 |\n", "| C | 0 | 1 | 0 | 0 |\n", "| G | 0 | 0 | 1 | 0 |\n", "| T | 0 | 0 | 0 | 1 |\n", "| N (inconnu) | 0 | 0 | 0 | 0 |\n", "\n", "Une chaîne de longueur 50 devient donc une matrice de $50 \\times 4$." ], "id": "37951449-6ed5-416d-926c-33c1ca36331c" }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [], "source": [ "def one_hot_encode(sequences, sequence_length=SEQUENCE_LENGTH):\n", " \"\"\"Encode les chaînes de nucléotides dans un tableau (n, sequence_length, 4).\"\"\"\n", " sequences = pd.Series(sequences, dtype=\"string\").str.upper()\n", "\n", " valid = sequences.str.fullmatch(rf\"[ACGTN]{{{sequence_length}}}\")\n", " if not valid.all():\n", " raise ValueError(\"Chaque séquence doit respecter le format attendu.\")\n", "\n", " # Convertit tous les caractères ASCII à la fois et repère leur canal.\n", " characters = np.frombuffer(\n", " \"\".join(sequences).encode(\"ascii\"), dtype=np.uint8\n", " ).reshape(-1, sequence_length)\n", "\n", " channel_lookup = np.full(256, -1, dtype=np.int8)\n", " for channel, nucleotide in enumerate(\"ACGT\"):\n", " channel_lookup[ord(nucleotide)] = channel\n", "\n", " channels = channel_lookup[characters]\n", " encoded = np.zeros(\n", " (len(sequences), sequence_length, 4), dtype=np.float32\n", " )\n", "\n", " rows, positions = np.nonzero(channels >= 0)\n", " encoded[rows, positions, channels[rows, positions]] = 1.0\n", " return encoded" ], "id": "one-hot-helper" }, { "cell_type": "markdown", "metadata": {}, "source": [ "Nous visualisons le début d’une séquence avant d’encoder le jeu de\n", "données complet." ], "id": "e6eab2a7-e7a6-4d2f-85c9-288c47285d24" }, { "cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [ { "output_type": "display_data", "metadata": {}, "data": {} }, { "output_type": "stream", "name": "stdout", "text": [ "Longueur de la séquence complète : 50\n", "Forme de l'encodage : (50, 4)" ] } ], "source": [ "example_sequence = retained_data.loc[0, \"utr\"]\n", "example_encoding = one_hot_encode([example_sequence])[0]\n", "positions_to_show = 12\n", "\n", "fig, ax = plt.subplots(figsize=(9, 2.8))\n", "sns.heatmap(\n", " example_encoding[:positions_to_show].T,\n", " cmap=\"Blues\",\n", " vmin=0,\n", " vmax=1,\n", " cbar=False,\n", " linewidths=0.5,\n", " linecolor=\"lightgray\",\n", " xticklabels=np.arange(1, positions_to_show + 1),\n", " yticklabels=list(\"ACGT\"),\n", " ax=ax,\n", ")\n", "ax.set(\n", " xlabel=\"Position dans la séquence\",\n", " ylabel=\"Canal\",\n", " title=f\"Préfixe de la séquence : {example_sequence[:positions_to_show]}\",\n", ")\n", "plt.tight_layout()\n", "plt.show()\n", "\n", "print(\"Longueur de la séquence complète :\", len(example_sequence))\n", "print(\"Forme de l'encodage :\", example_encoding.shape)" ], "id": "cell-visualize-one-hot" }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Expérience 1 : partitionnement aléatoire\n", "\n", "## Création des ensembles d’entraînement, de validation et de test\n", "\n", "Nous mélangeons aléatoirement les séquences conservées, puis créons\n", "trois partitions disjointes :\n", "\n", "- 220 000 séquences pour ajuster le modèle;\n", "- 27 500 séquences pour l’arrêt précoce et la sélection du modèle;\n", "- 27 500 séquences pour l’évaluation finale.\n", "\n", "L’ensemble de test demeure intact jusqu’à la fin de l’entraînement." ], "id": "a5cc448e-4177-427c-86eb-a94dfd74dff6" }, { "cell_type": "code", "execution_count": 12, "metadata": {}, "outputs": [], "source": [ "def random_partition(data, seed=SEED):\n", " \"\"\"Crée les partitions aléatoires d'entraînement, de validation et de test.\"\"\"\n", " if len(data) != N_RETAINED:\n", " raise ValueError(f\"Exactement {N_RETAINED:,} lignes étaient attendues.\")\n", "\n", " shuffled = data.sample(frac=1, random_state=seed).reset_index(drop=True)\n", "\n", " train = shuffled.iloc[:N_TRAIN].copy()\n", " validation = shuffled.iloc[\n", " N_TRAIN:N_TRAIN + N_VALIDATION\n", " ].copy()\n", " test = shuffled.iloc[N_TRAIN + N_VALIDATION:].copy()\n", "\n", " if (len(train), len(validation), len(test)) != (\n", " N_TRAIN,\n", " N_VALIDATION,\n", " N_TEST,\n", " ):\n", " raise AssertionError(\"Tailles de partitions inattendues.\")\n", "\n", " train_sequences = set(train[\"utr\"])\n", " validation_sequences = set(validation[\"utr\"])\n", " test_sequences = set(test[\"utr\"])\n", "\n", " if not train_sequences.isdisjoint(validation_sequences):\n", " raise AssertionError(\"Les ensembles d'entraînement et de validation se chevauchent.\")\n", " if not train_sequences.isdisjoint(test_sequences):\n", " raise AssertionError(\"Les ensembles d'entraînement et de test se chevauchent.\")\n", " if not validation_sequences.isdisjoint(test_sequences):\n", " raise AssertionError(\"Les ensembles de validation et de test se chevauchent.\")\n", "\n", " return train, validation, test\n", "\n", "\n", "train_data, validation_data, test_data = random_partition(retained_data)" ], "id": "random-partition-helper" }, { "cell_type": "markdown", "metadata": {}, "source": [ "Les trois partitions devraient présenter des distributions semblables de\n", "la MRL et du nombre de lectures. Ce petit tableau permet de vérifier la\n", "randomisation." ], "id": "93eb3f9f-43c4-4551-9a4c-7c55d3a01536" }, { "cell_type": "code", "execution_count": 13, "metadata": {}, "outputs": [ { "output_type": "display_data", "metadata": {}, "data": { "text/html": [ "\n", "" ] } } ], "source": [ "def summarize_partitions(partitions):\n", " \"\"\"Résume la taille et la distribution cible des partitions nommées.\"\"\"\n", " rows = []\n", " for name, data in partitions.items():\n", " rows.append(\n", " {\n", " \"partition\": name,\n", " \"n\": len(data),\n", " \"MRL moyenne\": data[\"mrl\"].mean(),\n", " \"écart-type MRL\": data[\"mrl\"].std(),\n", " \"lectures médianes\": data[\"read_count\"].median(),\n", " }\n", " )\n", " return pd.DataFrame(rows).set_index(\"partition\")\n", "\n", "\n", "random_partition_summary = summarize_partitions(\n", " {\n", " \"Entraînement\": train_data,\n", " \"Validation\": validation_data,\n", " \"Test\": test_data,\n", " }\n", ")\n", "random_partition_summary.round(3)" ], "id": "cell-summarize-random-partitions" }, { "cell_type": "markdown", "metadata": {}, "source": [ "Nous standardisons les cibles d’entraînement pour qu’elles aient une\n", "moyenne de zéro et un écart-type de un. Nous appliquons ensuite la même\n", "transformation aux cibles de validation et de test. Ajuster cette\n", "transformation uniquement sur les données d’entraînement empêche\n", "l’information des ensembles de validation et de test de se propager dans\n", "la mise au point du modèle." ], "id": "55d5c966-d24f-4b4c-9a6a-530c4885ca51" }, { "cell_type": "code", "execution_count": 14, "metadata": {}, "outputs": [ { "output_type": "stream", "name": "stdout", "text": [ "Caractéristiques d'entraînement : (220000, 50, 4) float32\n", "Caractéristiques de validation : (27500, 50, 4) float32\n", "Caractéristiques de test : (27500, 50, 4) float32\n", "Cibles d'entraînement : (220000, 1) float32" ] } ], "source": [ "target_scaler = StandardScaler()\n", "\n", "y_train = target_scaler.fit_transform(\n", " train_data[[\"mrl\"]]\n", ").astype(np.float32)\n", "y_validation = target_scaler.transform(\n", " validation_data[[\"mrl\"]]\n", ").astype(np.float32)\n", "y_test = target_scaler.transform(\n", " test_data[[\"mrl\"]]\n", ").astype(np.float32)\n", "\n", "x_train = one_hot_encode(train_data[\"utr\"])\n", "x_validation = one_hot_encode(validation_data[\"utr\"])\n", "x_test = one_hot_encode(test_data[\"utr\"])\n", "\n", "print(\"Caractéristiques d'entraînement :\", x_train.shape, x_train.dtype)\n", "print(\"Caractéristiques de validation :\", x_validation.shape, x_validation.dtype)\n", "print(\"Caractéristiques de test :\", x_test.shape, x_test.dtype)\n", "print(\"Cibles d'entraînement :\", y_train.shape, y_train.dtype)" ], "id": "prepare-arrays" }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Construction du réseau neuronal convolutif\n", "\n", "La forme d’entrée est `(50, 4)` : 50 positions et quatre canaux de\n", "nucléotides. Keras utilise directement cette représentation où les\n", "canaux sont placés en dernier; aucune permutation du tenseur n’est donc\n", "nécessaire." ], "id": "889c432f-27b1-4f69-82a2-5a652dc4333b" }, { "cell_type": "code", "execution_count": 15, "metadata": {}, "outputs": [ { "output_type": "display_data", "metadata": {}, "data": { "text/html": [ "
Model: \"mrl_cnn\"\n",
              "
" ] } }, { "output_type": "display_data", "metadata": {}, "data": { "text/html": [ "
┏━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━┳━━━━━━━━━━━━━━━━━━━━━━━━┳━━━━━━━━━━━━━━━┓\n",
              "┃ Layer (type)                    ┃ Output Shape           ┃       Param # ┃\n",
              "┡━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━╇━━━━━━━━━━━━━━━━━━━━━━━━╇━━━━━━━━━━━━━━━┩\n",
              "│ convolution_1 (Conv1D)          │ (None, 50, 160)        │         5,280 │\n",
              "├─────────────────────────────────┼────────────────────────┼───────────────┤\n",
              "│ convolution_2 (Conv1D)          │ (None, 50, 160)        │       204,960 │\n",
              "├─────────────────────────────────┼────────────────────────┼───────────────┤\n",
              "│ batch_normalization_1           │ (None, 50, 160)        │           640 │\n",
              "│ (BatchNormalization)            │                        │               │\n",
              "├─────────────────────────────────┼────────────────────────┼───────────────┤\n",
              "│ dropout_1 (Dropout)             │ (None, 50, 160)        │             0 │\n",
              "├─────────────────────────────────┼────────────────────────┼───────────────┤\n",
              "│ convolution_3 (Conv1D)          │ (None, 50, 160)        │       204,960 │\n",
              "├─────────────────────────────────┼────────────────────────┼───────────────┤\n",
              "│ batch_normalization_2           │ (None, 50, 160)        │           640 │\n",
              "│ (BatchNormalization)            │                        │               │\n",
              "├─────────────────────────────────┼────────────────────────┼───────────────┤\n",
              "│ dropout_2 (Dropout)             │ (None, 50, 160)        │             0 │\n",
              "├─────────────────────────────────┼────────────────────────┼───────────────┤\n",
              "│ convolution_4 (Conv1D)          │ (None, 50, 80)         │       102,480 │\n",
              "├─────────────────────────────────┼────────────────────────┼───────────────┤\n",
              "│ batch_normalization_3           │ (None, 50, 80)         │           320 │\n",
              "│ (BatchNormalization)            │                        │               │\n",
              "├─────────────────────────────────┼────────────────────────┼───────────────┤\n",
              "│ dropout_3 (Dropout)             │ (None, 50, 80)         │             0 │\n",
              "├─────────────────────────────────┼────────────────────────┼───────────────┤\n",
              "│ flatten (Flatten)               │ (None, 4000)           │             0 │\n",
              "├─────────────────────────────────┼────────────────────────┼───────────────┤\n",
              "│ latent_representation (Dense)   │ (None, 80)             │       320,080 │\n",
              "├─────────────────────────────────┼────────────────────────┼───────────────┤\n",
              "│ scaled_mrl (Dense)              │ (None, 1)              │            81 │\n",
              "└─────────────────────────────────┴────────────────────────┴───────────────┘\n",
              "
" ] } }, { "output_type": "display_data", "metadata": {}, "data": { "text/html": [ "
 Total params: 839,441 (3.20 MB)\n",
              "
" ] } }, { "output_type": "display_data", "metadata": {}, "data": { "text/html": [ "
 Trainable params: 838,641 (3.20 MB)\n",
              "
" ] } }, { "output_type": "display_data", "metadata": {}, "data": { "text/html": [ "
 Non-trainable params: 800 (3.12 KB)\n",
              "
" ] } } ], "source": [ "def build_mrl_cnn(\n", " sequence_length=SEQUENCE_LENGTH,\n", " n_filters=160,\n", " kernel_size=8,\n", " latent_dimension=80,\n", "):\n", " \"\"\"Construit le CNN 1D servant à prédire la MRL standardisée.\"\"\"\n", " model = keras.Sequential(\n", " [\n", " keras.Input(shape=(sequence_length, 4), name=\"one_hot_utr\"),\n", " layers.Conv1D(\n", " n_filters,\n", " kernel_size,\n", " padding=\"same\",\n", " activation=\"relu\",\n", " name=\"convolution_1\",\n", " ),\n", " layers.Conv1D(\n", " n_filters,\n", " kernel_size,\n", " padding=\"same\",\n", " activation=\"relu\",\n", " name=\"convolution_2\",\n", " ),\n", " layers.BatchNormalization(name=\"batch_normalization_1\"),\n", " layers.Dropout(0.2, name=\"dropout_1\"),\n", " layers.Conv1D(\n", " n_filters,\n", " kernel_size,\n", " padding=\"same\",\n", " activation=\"relu\",\n", " name=\"convolution_3\",\n", " ),\n", " layers.BatchNormalization(name=\"batch_normalization_2\"),\n", " layers.Dropout(0.4, name=\"dropout_2\"),\n", " layers.Conv1D(\n", " n_filters // 2,\n", " kernel_size,\n", " padding=\"same\",\n", " activation=\"relu\",\n", " name=\"convolution_4\",\n", " ),\n", " layers.BatchNormalization(name=\"batch_normalization_3\"),\n", " layers.Dropout(0.2, name=\"dropout_3\"),\n", " layers.Flatten(name=\"flatten\"),\n", " layers.Dense(\n", " latent_dimension,\n", " activation=\"relu\",\n", " name=\"latent_representation\",\n", " ),\n", " layers.Dense(1, name=\"scaled_mrl\"),\n", " ],\n", " name=\"mrl_cnn\",\n", " )\n", " return model\n", "\n", "\n", "model = build_mrl_cnn()\n", "model.summary()" ], "id": "model-helper" }, { "cell_type": "markdown", "metadata": {}, "source": [ "Chaque convolution utilise un noyau de largeur huit et `padding=\"same\"`.\n", "La dimension spatiale demeure donc égale à 50 pendant que le réseau\n", "construit des représentations de plus en plus abstraites des motifs\n", "locaux de la séquence. La dernière unité linéaire produit une prédiction\n", "de MRL standardisée.\n", "\n", "## Entraînement du modèle\n", "\n", "Keras fournit une implémentation cumulative de $R^2$. Nous implémentons\n", "la corrélation de Pearson de la même façon, en accumulant les\n", "statistiques suffisantes pendant une époque complète. Calculer la\n", "corrélation séparément pour chaque mini-lot, puis faire la moyenne de\n", "ces valeurs ne donnerait pas la corrélation sur l’ensemble complet." ], "id": "81489792-2ffa-4a03-a794-656bab1494ed" }, { "cell_type": "code", "execution_count": 16, "metadata": {}, "outputs": [], "source": [ "@keras.utils.register_keras_serializable(package=\"CSI4106\")\n", "class PearsonCorrelation(keras.metrics.Metric):\n", " \"\"\"Corrélation de Pearson accumulée sur tous les exemples d'une époque.\"\"\"\n", "\n", " def __init__(self, name=\"pearson_r\", **kwargs):\n", " super().__init__(name=name, **kwargs)\n", " self.count = self.add_weight(name=\"count\", initializer=\"zeros\")\n", " self.sum_true = self.add_weight(name=\"sum_true\", initializer=\"zeros\")\n", " self.sum_pred = self.add_weight(name=\"sum_pred\", initializer=\"zeros\")\n", " self.sum_true_squared = self.add_weight(\n", " name=\"sum_true_squared\", initializer=\"zeros\"\n", " )\n", " self.sum_pred_squared = self.add_weight(\n", " name=\"sum_pred_squared\", initializer=\"zeros\"\n", " )\n", " self.sum_products = self.add_weight(\n", " name=\"sum_products\", initializer=\"zeros\"\n", " )\n", "\n", " def update_state(self, y_true, y_pred, sample_weight=None):\n", " y_true = tf.cast(tf.reshape(y_true, [-1]), self.dtype)\n", " y_pred = tf.cast(tf.reshape(y_pred, [-1]), self.dtype)\n", "\n", " self.count.assign_add(tf.cast(tf.size(y_true), self.dtype))\n", " self.sum_true.assign_add(tf.reduce_sum(y_true))\n", " self.sum_pred.assign_add(tf.reduce_sum(y_pred))\n", " self.sum_true_squared.assign_add(tf.reduce_sum(tf.square(y_true)))\n", " self.sum_pred_squared.assign_add(tf.reduce_sum(tf.square(y_pred)))\n", " self.sum_products.assign_add(tf.reduce_sum(y_true * y_pred))\n", "\n", " def result(self):\n", " covariance = (\n", " self.sum_products\n", " - self.sum_true * self.sum_pred / self.count\n", " )\n", " true_variation = (\n", " self.sum_true_squared\n", " - tf.square(self.sum_true) / self.count\n", " )\n", " pred_variation = (\n", " self.sum_pred_squared\n", " - tf.square(self.sum_pred) / self.count\n", " )\n", " denominator = tf.sqrt(true_variation * pred_variation)\n", " return tf.math.divide_no_nan(covariance, denominator)\n", "\n", " def reset_state(self):\n", " for variable in self.variables:\n", " variable.assign(0)" ], "id": "pearson-metric" }, { "cell_type": "markdown", "metadata": {}, "source": [ "L’erreur quadratique moyenne est à la fois notre objectif d’optimisation\n", "et la perte affichée pendant l’entraînement. La corrélation de Pearson\n", "et $R^2$ sont des métriques de suivi; elles ne modifient pas les mises à\n", "jour du gradient." ], "id": "57d31b18-336e-40ec-9b83-570d8fbcf11e" }, { "cell_type": "code", "execution_count": 17, "metadata": {}, "outputs": [], "source": [ "def compile_mrl_cnn(model):\n", " \"\"\"Configure l'optimiseur, la perte et les métriques calculées par époque.\"\"\"\n", " model.compile(\n", " optimizer=keras.optimizers.Adam(learning_rate=1e-3),\n", " loss=\"mean_squared_error\",\n", " metrics=[\n", " PearsonCorrelation(),\n", " keras.metrics.R2Score(name=\"r2\"),\n", " ],\n", " )\n", " return model\n", "\n", "\n", "model = compile_mrl_cnn(model)" ], "id": "compile-model" }, { "cell_type": "markdown", "metadata": {}, "source": [ "L’arrêt précoce surveille la perte de validation. Lorsque cette perte ne\n", "s’est pas améliorée pendant cinq époques, l’entraînement s’arrête et\n", "Keras restaure les poids de la meilleure époque. Un point de contrôle\n", "est également écrit sur le disque afin de pouvoir réutiliser le modèle\n", "dans des notebooks ultérieurs.\n", "\n", "> **Durée d’exécution prévue**\n", ">\n", "> Il s’agit de la cellule la plus exigeante en calcul du notebook. La\n", "> durée d’exécution dépend fortement de l’accélérateur disponible pour\n", "> TensorFlow. La cellule utilise les 220 000 exemples d’entraînement et\n", "> peut prendre plusieurs minutes." ], "id": "ab9b8d89-3d8c-4e7b-b4c8-75243a25d5b7" }, { "cell_type": "code", "execution_count": 18, "metadata": {}, "outputs": [ { "output_type": "stream", "name": "stdout", "text": [ "Epoch 1/50\n", "860/860 - 78s - 91ms/step - loss: 0.3765 - pearson_r: 0.7908 - r2: 0.6235 - val_loss: 0.2153 - val_pearson_r: 0.8923 - val_r2: 0.7849\n", "Epoch 2/50\n", "860/860 - 77s - 89ms/step - loss: 0.1669 - pearson_r: 0.9128 - r2: 0.8331 - val_loss: 0.1323 - val_pearson_r: 0.9366 - val_r2: 0.8678\n", "Epoch 3/50\n", "860/860 - 77s - 90ms/step - loss: 0.1300 - pearson_r: 0.9328 - r2: 0.8700 - val_loss: 0.1132 - val_pearson_r: 0.9443 - val_r2: 0.8869\n", "Epoch 4/50\n", "860/860 - 77s - 89ms/step - loss: 0.1187 - pearson_r: 0.9388 - r2: 0.8813 - val_loss: 0.1099 - val_pearson_r: 0.9471 - val_r2: 0.8902\n", "Epoch 5/50\n", "860/860 - 77s - 90ms/step - loss: 0.1127 - pearson_r: 0.9420 - r2: 0.8873 - val_loss: 0.1040 - val_pearson_r: 0.9485 - val_r2: 0.8961\n", "Epoch 6/50\n", "860/860 - 77s - 89ms/step - loss: 0.1089 - pearson_r: 0.9440 - r2: 0.8911 - val_loss: 0.1001 - val_pearson_r: 0.9504 - val_r2: 0.9000\n", "Epoch 7/50\n", "860/860 - 77s - 89ms/step - loss: 0.1057 - pearson_r: 0.9457 - r2: 0.8943 - val_loss: 0.0995 - val_pearson_r: 0.9506 - val_r2: 0.9006\n", "Epoch 8/50\n", "860/860 - 77s - 89ms/step - loss: 0.1036 - pearson_r: 0.9468 - r2: 0.8964 - val_loss: 0.1065 - val_pearson_r: 0.9500 - val_r2: 0.8936\n", "Epoch 9/50\n", "860/860 - 77s - 89ms/step - loss: 0.1011 - pearson_r: 0.9481 - r2: 0.8989 - val_loss: 0.0953 - val_pearson_r: 0.9524 - val_r2: 0.9047\n", "Epoch 10/50\n", "860/860 - 76s - 89ms/step - loss: 0.0992 - pearson_r: 0.9491 - r2: 0.9008 - val_loss: 0.0934 - val_pearson_r: 0.9527 - val_r2: 0.9067\n", "Epoch 11/50\n", "860/860 - 76s - 89ms/step - loss: 0.0973 - pearson_r: 0.9501 - r2: 0.9027 - val_loss: 0.0919 - val_pearson_r: 0.9533 - val_r2: 0.9082\n", "Epoch 12/50\n", "860/860 - 79s - 92ms/step - loss: 0.0955 - pearson_r: 0.9511 - r2: 0.9045 - val_loss: 0.0925 - val_pearson_r: 0.9536 - val_r2: 0.9075\n", "Epoch 13/50\n", "860/860 - 81s - 94ms/step - loss: 0.0933 - pearson_r: 0.9522 - r2: 0.9067 - val_loss: 0.0908 - val_pearson_r: 0.9537 - val_r2: 0.9093\n", "Epoch 14/50\n", "860/860 - 80s - 93ms/step - loss: 0.0918 - pearson_r: 0.9530 - r2: 0.9082 - val_loss: 0.0931 - val_pearson_r: 0.9534 - val_r2: 0.9070\n", "Epoch 15/50\n", "860/860 - 80s - 93ms/step - loss: 0.0902 - pearson_r: 0.9538 - r2: 0.9098 - val_loss: 0.0914 - val_pearson_r: 0.9537 - val_r2: 0.9086\n", "Epoch 16/50\n", "860/860 - 80s - 93ms/step - loss: 0.0883 - pearson_r: 0.9548 - r2: 0.9117 - val_loss: 0.0898 - val_pearson_r: 0.9544 - val_r2: 0.9102\n", "Epoch 17/50\n", "860/860 - 80s - 93ms/step - loss: 0.0870 - pearson_r: 0.9555 - r2: 0.9130 - val_loss: 0.0924 - val_pearson_r: 0.9536 - val_r2: 0.9076\n", "Epoch 18/50\n", "860/860 - 80s - 93ms/step - loss: 0.0849 - pearson_r: 0.9566 - r2: 0.9151 - val_loss: 0.0924 - val_pearson_r: 0.9530 - val_r2: 0.9077\n", "Epoch 19/50\n", "860/860 - 82s - 95ms/step - loss: 0.0836 - pearson_r: 0.9573 - r2: 0.9164 - val_loss: 0.0915 - val_pearson_r: 0.9532 - val_r2: 0.9086\n", "Epoch 20/50\n", "860/860 - 82s - 95ms/step - loss: 0.0821 - pearson_r: 0.9581 - r2: 0.9179 - val_loss: 0.0915 - val_pearson_r: 0.9533 - val_r2: 0.9086\n", "Epoch 21/50\n", "860/860 - 81s - 94ms/step - loss: 0.0806 - pearson_r: 0.9588 - r2: 0.9194 - val_loss: 0.0939 - val_pearson_r: 0.9531 - val_r2: 0.9062\n", "Epoch 21: early stopping\n", "Restoring model weights from the end of the best epoch: 16." ] } ], "source": [ "random_split_model_path = MODEL_DIR / \"mrl_cnn_random_split.keras\"\n", "\n", "training_callbacks = [\n", " callbacks.EarlyStopping(\n", " monitor=\"val_loss\",\n", " mode=\"min\",\n", " patience=PATIENCE,\n", " restore_best_weights=True,\n", " verbose=1,\n", " ),\n", " callbacks.ModelCheckpoint(\n", " filepath=str(random_split_model_path),\n", " monitor=\"val_loss\",\n", " mode=\"min\",\n", " save_best_only=True,\n", " verbose=0,\n", " ),\n", "]\n", "\n", "history = model.fit(\n", " x_train,\n", " y_train,\n", " validation_data=(x_validation, y_validation),\n", " epochs=MAX_EPOCHS,\n", " batch_size=BATCH_SIZE,\n", " shuffle=True,\n", " callbacks=training_callbacks,\n", " verbose=2,\n", ")" ], "id": "train-model" }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Visualisation de la dynamique d’apprentissage\n", "\n", "L’objet `History` de Keras enregistre chaque perte et chaque métrique\n", "nommée lors de la compilation du modèle. Nous affichons trois vues\n", "complémentaires de l’apprentissage." ], "id": "2e4d6fc7-0bc5-4b88-a439-6421a22dff3f" }, { "cell_type": "code", "execution_count": 19, "metadata": {}, "outputs": [], "source": [ "def plot_training_history(history):\n", " \"\"\"Trace la perte, r de Pearson et R carré pour l'entraînement et la validation.\"\"\"\n", " values = history.history\n", " epochs = np.arange(1, len(values[\"loss\"]) + 1)\n", " best_epoch = int(np.argmin(values[\"val_loss\"])) + 1\n", "\n", " panels = [\n", " (\"loss\", \"Erreur quadratique moyenne\", \"Perte\"),\n", " (\"pearson_r\", \"Corrélation de Pearson\", r\"Pearson $r$\"),\n", " (\"r2\", \"Coefficient de détermination\", r\"$R^2$\"),\n", " ]\n", "\n", " fig, axes = plt.subplots(1, 3, figsize=(15, 4.2))\n", " for ax, (metric, title, ylabel) in zip(axes, panels):\n", " ax.plot(epochs, values[metric], label=\"Entraînement\", linewidth=2)\n", " ax.plot(\n", " epochs,\n", " values[f\"val_{metric}\"],\n", " label=\"Validation\",\n", " linewidth=2,\n", " )\n", " ax.axvline(\n", " best_epoch,\n", " color=\"black\",\n", " linestyle=\":\",\n", " linewidth=1.5,\n", " label=\"Meilleure perte de validation\" if metric == \"loss\" else None,\n", " )\n", " ax.set(\n", " xlabel=\"Époque\",\n", " ylabel=ylabel,\n", " title=title,\n", " )\n", " ax.legend()\n", " sns.despine(ax=ax)\n", "\n", " fig.suptitle(\"Dynamique d'apprentissage\", fontsize=15)\n", " fig.tight_layout()\n", " return fig, axes" ], "id": "history-plot-helper" }, { "cell_type": "code", "execution_count": 20, "metadata": {}, "outputs": [ { "output_type": "display_data", "metadata": {}, "data": {} } ], "source": [ "plot_training_history(history)\n", "plt.show()" ], "id": "cell-plot-training-history" }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Évaluation de l’ensemble de test\n", "\n", "Le réseau prédit une MRL standardisée. Nous ramenons sa sortie à\n", "l’échelle expérimentale d’origine avant de calculer les métriques\n", "finales ou de tracer les résultats." ], "id": "01f2e729-618d-47e2-a789-802eaafa945e" }, { "cell_type": "code", "execution_count": 21, "metadata": {}, "outputs": [], "source": [ "def predict_mrl(model, encoded_sequences, scaler, batch_size=BATCH_SIZE):\n", " \"\"\"Prédit la MRL et renvoie les valeurs à l'échelle expérimentale d'origine.\"\"\"\n", " scaled_predictions = model.predict(\n", " encoded_sequences,\n", " batch_size=batch_size,\n", " verbose=0,\n", " )\n", " return scaler.inverse_transform(scaled_predictions).ravel()\n", "\n", "\n", "def regression_metrics(observed, predicted):\n", " \"\"\"Renvoie les principales métriques de régression d'un ensemble d'évaluation.\"\"\"\n", " observed = np.asarray(observed).ravel()\n", " predicted = np.asarray(predicted).ravel()\n", "\n", " return {\n", " \"n\": len(observed),\n", " \"Pearson r\": pearsonr(observed, predicted)[0],\n", " \"Spearman rho\": spearmanr(observed, predicted)[0],\n", " \"R2\": r2_score(observed, predicted),\n", " \"RMSE\": np.sqrt(mean_squared_error(observed, predicted)),\n", " }" ], "id": "evaluation-helpers" }, { "cell_type": "code", "execution_count": 22, "metadata": {}, "outputs": [ { "output_type": "display_data", "metadata": {}, "data": { "text/html": [ "\n", "" ] } } ], "source": [ "test_observed = test_data[\"mrl\"].to_numpy()\n", "test_predicted = predict_mrl(model, x_test, target_scaler)\n", "\n", "test_metrics = regression_metrics(test_observed, test_predicted)\n", "pd.DataFrame([test_metrics], index=[\"Ensemble de test aléatoire\"]).round(4)" ], "id": "cell-evaluate-test-set" }, { "cell_type": "markdown", "metadata": {}, "source": [ "Avec des dizaines de milliers de points, un nuage de points ordinaire\n", "souffre d’une forte superposition. Un diagramme de densité hexagonal\n", "montre où se concentrent les observations, tandis que les histogrammes\n", "marginaux présentent la distribution le long de chaque axe." ], "id": "b85a5eec-850e-467f-8a3a-c5329f136418" }, { "cell_type": "code", "execution_count": 23, "metadata": {}, "outputs": [], "source": [ "def plot_observed_vs_predicted(\n", " observed,\n", " predicted,\n", " title,\n", " limits=MRL_PLOT_LIMITS,\n", "):\n", " \"\"\"Trace la densité, les distributions marginales et les droites d'étalonnage.\"\"\"\n", " observed = np.asarray(observed).ravel()\n", " predicted = np.asarray(predicted).ravel()\n", " metrics = regression_metrics(observed, predicted)\n", "\n", " grid = sns.JointGrid(\n", " x=observed,\n", " y=predicted,\n", " height=7,\n", " ratio=5,\n", " space=0.05,\n", " xlim=limits,\n", " ylim=limits,\n", " )\n", "\n", " grid.ax_joint.hexbin(\n", " observed,\n", " predicted,\n", " gridsize=60,\n", " mincnt=1,\n", " bins=\"log\",\n", " cmap=\"viridis\",\n", " )\n", " sns.histplot(x=observed, bins=50, color=\"steelblue\", ax=grid.ax_marg_x)\n", " sns.histplot(y=predicted, bins=50, color=\"darkorange\", ax=grid.ax_marg_y)\n", "\n", " line_x = np.asarray(limits)\n", " grid.ax_joint.plot(\n", " line_x,\n", " line_x,\n", " color=\"red\",\n", " linestyle=\"--\",\n", " linewidth=2,\n", " label=\"Étalonnage parfait\",\n", " )\n", "\n", " slope, intercept = np.polyfit(observed, predicted, deg=1)\n", " grid.ax_joint.plot(\n", " line_x,\n", " intercept + slope * line_x,\n", " color=\"white\",\n", " linewidth=2,\n", " label=\"Relation ajustée\",\n", " )\n", "\n", " annotation = (\n", " f\"n = {metrics['n']:,}\\n\"\n", " f\"Pearson r = {metrics['Pearson r']:.3f}\\n\"\n", " f\"R² = {metrics['R2']:.3f}\\n\"\n", " f\"ajustement : y = {intercept:.2f} + {slope:.2f}x\"\n", " )\n", " grid.ax_joint.text(\n", " 0.03,\n", " 0.97,\n", " annotation,\n", " transform=grid.ax_joint.transAxes,\n", " ha=\"left\",\n", " va=\"top\",\n", " color=\"white\",\n", " bbox={\"boxstyle\": \"round\", \"facecolor\": \"black\", \"alpha\": 0.65},\n", " )\n", " grid.ax_joint.set(\n", " xlabel=\"MRL expérimentale\",\n", " ylabel=\"MRL prédite\",\n", " )\n", " grid.ax_joint.legend(loc=\"lower right\")\n", " grid.fig.suptitle(title, y=1.02, fontsize=15)\n", " return grid" ], "id": "prediction-plot-helper" }, { "cell_type": "code", "execution_count": 24, "metadata": {}, "outputs": [ { "output_type": "display_data", "metadata": {}, "data": {} } ], "source": [ "plot_observed_vs_predicted(\n", " test_observed,\n", " test_predicted,\n", " title=\"Ensemble de test aléatoire\",\n", ")\n", "plt.show()" ], "id": "cell-plot-test-predictions" }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Évaluation sur une expérience indépendante menée sur des lymphocytes T\n", "\n", "L’ensemble de test précédent provient de la même expérience et de la\n", "même population filtrée que les données d’entraînement. Nous posons\n", "maintenant une question plus difficile : le modèle se transfère-t-il à\n", "des mesures réalisées dans un contexte cellulaire différent?\n", "\n", "Deux réplicats biologiques sont disponibles. Une séquence mesurée avec\n", "une plus grande profondeur dans un réplicat devrait recevoir davantage\n", "d’influence de ce réplicat. Nous combinons donc les MRL des réplicats en\n", "utilisant leurs nombres de lectures comme poids :\n", "\n", "$$\n", "\\operatorname{MRL}_{\\mathrm{merged}}\n", "=\n", "\\frac{\n", " \\operatorname{MRL}_1 n_1 + \\operatorname{MRL}_2 n_2\n", "}{n_1+n_2},\n", "$$\n", "\n", "où $n_1$ et $n_2$ sont les nombres de lectures des réplicats." ], "id": "6c420c50-dff3-4e84-a4f4-4589a924de28" }, { "cell_type": "code", "execution_count": 25, "metadata": {}, "outputs": [], "source": [ "def merge_replicates(replicate_1, replicate_2):\n", " \"\"\"Joint deux réplicats et calcule la MRL pondérée par les lectures.\"\"\"\n", " first = replicate_1.rename(\n", " columns={\n", " \"mrl\": \"mrl_1\",\n", " \"read_count\": \"read_count_1\",\n", " }\n", " )\n", " second = replicate_2.rename(\n", " columns={\n", " \"mrl\": \"mrl_2\",\n", " \"read_count\": \"read_count_2\",\n", " }\n", " )\n", "\n", " merged = first.merge(second, how=\"outer\", on=\"utr\", validate=\"one_to_one\")\n", "\n", " for replicate in (1, 2):\n", " read_column = f\"read_count_{replicate}\"\n", " mrl_column = f\"mrl_{replicate}\"\n", "\n", " merged[read_column] = merged[read_column].fillna(0)\n", " inconsistent = (\n", " (merged[read_column] > 0) & merged[mrl_column].isna()\n", " )\n", " if inconsistent.any():\n", " raise ValueError(\n", " f\"Le réplicat {replicate} a des lectures positives, mais sa MRL est manquante.\"\n", " )\n", " merged[mrl_column] = merged[mrl_column].fillna(0)\n", "\n", " merged[\"read_count\"] = (\n", " merged[\"read_count_1\"] + merged[\"read_count_2\"]\n", " )\n", "\n", " positive_reads = merged[\"read_count\"] > 0\n", " merged = merged.loc[positive_reads].copy()\n", " merged[\"mrl\"] = (\n", " merged[\"mrl_1\"] * merged[\"read_count_1\"]\n", " + merged[\"mrl_2\"] * merged[\"read_count_2\"]\n", " ) / merged[\"read_count\"]\n", "\n", " return merged.loc[:, [\"utr\", \"mrl\", \"read_count\"]].reset_index(\n", " drop=True\n", " )" ], "id": "merge-replicates-helper" }, { "cell_type": "code", "execution_count": 26, "metadata": {}, "outputs": [ { "output_type": "stream", "name": "stdout", "text": [ "Réplicat 1 : 349,633 séquences\n", "Réplicat 2 : 352,724 séquences\n", "Union fusionnée : 387,794 séquences\n", "Séquences conservées : 100,000\n", "Seuil du nombre de lectures : 888" ] } ], "source": [ "tcell_r1_path = download_dataset(\n", " \"GSE232927_processed_defined_end_tcell_r1.csv.gz\"\n", ")\n", "tcell_r2_path = download_dataset(\n", " \"GSE232927_processed_defined_end_tcell_r2.csv.gz\"\n", ")\n", "\n", "tcell_r1 = load_utr_dataset(tcell_r1_path)\n", "tcell_r2 = load_utr_dataset(tcell_r2_path)\n", "\n", "tcell_merged = merge_replicates(tcell_r1, tcell_r2)\n", "tcell_data, tcell_read_cutoff = select_top_by_reads(tcell_merged, 100_000)\n", "\n", "print(f\"Réplicat 1 : {len(tcell_r1):,} séquences\")\n", "print(f\"Réplicat 2 : {len(tcell_r2):,} séquences\")\n", "print(f\"Union fusionnée : {len(tcell_merged):,} séquences\")\n", "print(f\"Séquences conservées : {len(tcell_data):,}\")\n", "print(f\"Seuil du nombre de lectures : {tcell_read_cutoff:,.0f}\")" ], "id": "load-and-merge-tcell" }, { "cell_type": "markdown", "metadata": {}, "source": [ "Les données indépendantes ne doivent pas influencer le prétraitement.\n", "Nous réutilisons l’encodeur à chaud et le normalisateur de la cible\n", "ajusté lors de l’expérience initiale." ], "id": "4ce80e0f-efe3-421f-b052-2e4f5ab4b2a3" }, { "cell_type": "code", "execution_count": 27, "metadata": {}, "outputs": [ { "output_type": "display_data", "metadata": {}, "data": { "text/html": [ "\n", "" ] } } ], "source": [ "x_tcell = one_hot_encode(tcell_data[\"utr\"])\n", "tcell_observed = tcell_data[\"mrl\"].to_numpy()\n", "tcell_predicted = predict_mrl(model, x_tcell, target_scaler)\n", "\n", "tcell_metrics = regression_metrics(tcell_observed, tcell_predicted)\n", "\n", "evaluation_summary = pd.DataFrame(\n", " [test_metrics, tcell_metrics],\n", " index=[\"Ensemble de test aléatoire\", \"Lymphocytes T indépendants\"],\n", ")\n", "evaluation_summary.round(4)" ], "id": "cell-evaluate-tcell" }, { "cell_type": "code", "execution_count": 28, "metadata": {}, "outputs": [ { "output_type": "display_data", "metadata": {}, "data": {} } ], "source": [ "plot_observed_vs_predicted(\n", " tcell_observed,\n", " tcell_predicted,\n", " title=\"Expérience indépendante sur des lymphocytes T\",\n", ")\n", "plt.show()" ], "id": "cell-plot-tcell-predictions" }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Corrélation, classement et étalonnage\n", "\n", "Les deux droites de référence, rouge et blanche, répondent à des\n", "questions différentes :\n", "\n", "- La droite d’identité rouge en tirets représente des prédictions\n", " numériques parfaites.\n", "- La droite blanche ajustée représente la relation linéaire réellement\n", " présente dans les données.\n", "\n", "La corrélation de Pearson mesure la force d’une relation linéaire. Elle\n", "demeure inchangée lorsque toutes les prédictions reçoivent le même\n", "décalage et la même mise à l’échelle positive. Le $R^2$ prédictif est\n", "plus strict :\n", "\n", "$$\n", "R^2\n", "=\n", "1-\n", "\\frac{\\sum_i(y_i-\\hat{y}_i)^2}\n", " {\\sum_i(y_i-\\bar{y})^2}.\n", "$$\n", "\n", "Il pénalise donc les décalages et les erreurs d’échelle. Les données\n", "indépendantes peuvent présenter une forte corrélation de Pearson, mais\n", "un $R^2$ plus faible lorsque la droite blanche ajustée diffère de la\n", "droite d’identité rouge.\n", "\n", "La corrélation de Spearman, incluse dans le tableau de résultats, mesure\n", "plus directement la concordance monotone des rangs. Cette propriété sera\n", "utile pour notre future application de conception : un algorithme\n", "génétique a surtout besoin que le modèle classe efficacement les\n", "séquences candidates prometteuses.\n", "\n", "> **En général, $R^2$ n’est pas égal à $r^2$**\n", ">\n", "> Certains articles présentent le carré de la corrélation de Pearson\n", "> sous le nom de « R carré ». Ici, `R2` désigne toujours le coefficient\n", "> de détermination prédictif calculé à partir des erreurs par rapport à\n", "> la droite d’identité.\n", "\n", "## Sauvegarde du prédicteur réutilisable\n", "\n", "Le meilleur modèle a été sauvegardé pendant l’entraînement. Nous\n", "sauvegardons aussi les deux nombres nécessaires pour ramener les\n", "prédictions à l’échelle originale de la MRL. Ces artefacts permettront à\n", "un notebook ultérieur sur les algorithmes génétiques d’évaluer de\n", "nouvelles chaînes sans charger les données expérimentales\n", "d’entraînement." ], "id": "6324a3f8-7494-4c51-86a5-f012a9b917c1" }, { "cell_type": "code", "execution_count": 29, "metadata": {}, "outputs": [ { "output_type": "stream", "name": "stdout", "text": [ "Modèle sauvegardé : models/mrl_cnn_random_split.keras\n", "Transformation de la cible sauvegardée : models/mrl_target_scaler.npz" ] } ], "source": [ "scaler_path = MODEL_DIR / \"mrl_target_scaler.npz\"\n", "np.savez(\n", " scaler_path,\n", " mean=target_scaler.mean_.astype(np.float32),\n", " scale=target_scaler.scale_.astype(np.float32),\n", ")\n", "\n", "print(\"Modèle sauvegardé :\", random_split_model_path)\n", "print(\"Transformation de la cible sauvegardée :\", scaler_path)" ], "id": "save-artifacts" }, { "cell_type": "markdown", "metadata": {}, "source": [ "L’interface de prédiction complète consiste maintenant à valider les\n", "chaînes, à les encoder à chaud, à appliquer le modèle, puis à annuler la\n", "standardisation de la cible." ], "id": "8f032131-3fa5-40a0-9e92-520acfc1eba8" }, { "cell_type": "code", "execution_count": 30, "metadata": {}, "outputs": [ { "output_type": "display_data", "metadata": {}, "data": { "text/html": [ "\n", "" ] } } ], "source": [ "example_sequences = pd.Series(\n", " [\n", " \"A\" * SEQUENCE_LENGTH,\n", " \"C\" * SEQUENCE_LENGTH,\n", " (\"ACGT\" * 13)[:SEQUENCE_LENGTH],\n", " ],\n", " name=\"utr\",\n", ")\n", "\n", "example_predictions = predict_mrl(\n", " model,\n", " one_hot_encode(example_sequences),\n", " target_scaler,\n", ")\n", "\n", "pd.DataFrame(\n", " {\n", " \"utr\": example_sequences,\n", " \"predicted_mrl\": example_predictions,\n", " }\n", ").round({\"predicted_mrl\": 3})" ], "id": "cell-predict-literal-sequences" }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Expérience 2 : partitionnement selon le nombre de lectures\n", "\n", "Notre expérience principale utilise un ensemble de test aléatoire, car\n", "elle cherche à déterminer dans quelle mesure le modèle se généralise à\n", "de nouvelles séquences de la même population filtrée. Sample et coll.\n", "posaient une question différente et réservaient à l’évaluation les\n", "séquences mesurées avec la plus grande profondeur.\n", "\n", "Pour isoler l’effet de ce choix expérimental, nous conservons exactement\n", "les mêmes tailles de partitions que précédemment :\n", "\n", "- les 27 500 séquences ayant le plus grand nombre de lectures forment\n", " l’ensemble de test;\n", "- 220 000 des séquences restantes forment l’ensemble d’entraînement;\n", "- les 27 500 dernières séquences forment l’ensemble de validation.\n", "\n", "Seule la règle de sélection de l’ensemble de test change." ], "id": "e9e61175-e6d1-4188-a193-c350096ea2fc" }, { "cell_type": "code", "execution_count": 31, "metadata": {}, "outputs": [ { "output_type": "display_data", "metadata": {}, "data": { "text/html": [ "\n", "" ] } } ], "source": [ "def read_count_partition(data, seed=SEED):\n", " \"\"\"Réserve les lignes les plus lues au test, puis répartit les autres.\"\"\"\n", " if len(data) != N_RETAINED:\n", " raise ValueError(f\"Exactement {N_RETAINED:,} lignes étaient attendues.\")\n", "\n", " ranked = data.sort_values(\n", " \"read_count\", ascending=False, kind=\"stable\"\n", " ).reset_index(drop=True)\n", "\n", " test = ranked.iloc[:N_TEST].copy()\n", " remaining = ranked.iloc[N_TEST:].sample(\n", " frac=1,\n", " random_state=seed,\n", " ).reset_index(drop=True)\n", "\n", " train = remaining.iloc[:N_TRAIN].copy()\n", " validation = remaining.iloc[N_TRAIN:].copy()\n", "\n", " if (len(train), len(validation), len(test)) != (\n", " N_TRAIN,\n", " N_VALIDATION,\n", " N_TEST,\n", " ):\n", " raise AssertionError(\"Tailles de partitions inattendues.\")\n", "\n", " return train, validation, test\n", "\n", "\n", "read_train_data, read_validation_data, read_test_data = (\n", " read_count_partition(retained_data)\n", ")\n", "\n", "read_partition_summary = summarize_partitions(\n", " {\n", " \"Entraînement\": read_train_data,\n", " \"Validation\": read_validation_data,\n", " \"Test (plus de lectures)\": read_test_data,\n", " }\n", ")\n", "read_partition_summary.round(3)" ], "id": "cell-read-count-partition-helper" }, { "cell_type": "markdown", "metadata": {}, "source": [ "Les distributions ci-dessous montrent que la sélection de l’ensemble de\n", "test selon le nombre de lectures modifie la couverture des mesures et,\n", "potentiellement, la distribution de la cible." ], "id": "4bfbd333-bd21-4976-b551-e83733c8cea0" }, { "cell_type": "code", "execution_count": 32, "metadata": {}, "outputs": [ { "output_type": "display_data", "metadata": {}, "data": {} } ], "source": [ "partition_plot_data = pd.concat(\n", " [\n", " read_train_data[[\"mrl\", \"read_count\"]].assign(\n", " partition=\"Entraînement\"\n", " ),\n", " read_validation_data[[\"mrl\", \"read_count\"]].assign(\n", " partition=\"Validation\"\n", " ),\n", " read_test_data[[\"mrl\", \"read_count\"]].assign(\n", " partition=\"Test (plus de lectures)\"\n", " ),\n", " ],\n", " ignore_index=True,\n", ")\n", "partition_plot_data[\"log_read_count\"] = np.log10(\n", " partition_plot_data[\"read_count\"] + 1\n", ")\n", "\n", "fig, axes = plt.subplots(1, 2, figsize=(13, 4.2))\n", "sns.histplot(\n", " data=partition_plot_data,\n", " x=\"log_read_count\",\n", " hue=\"partition\",\n", " bins=60,\n", " stat=\"density\",\n", " common_norm=False,\n", " element=\"step\",\n", " fill=False,\n", " ax=axes[0],\n", ")\n", "axes[0].set(\n", " xlabel=r\"$\\log_{10}(\\mathrm{nombre\\ de\\ lectures}+1)$\",\n", " ylabel=\"Densité\",\n", " title=\"Distribution du nombre de lectures\",\n", ")\n", "\n", "sns.histplot(\n", " data=partition_plot_data,\n", " x=\"mrl\",\n", " hue=\"partition\",\n", " bins=60,\n", " stat=\"density\",\n", " common_norm=False,\n", " element=\"step\",\n", " fill=False,\n", " ax=axes[1],\n", " legend=False,\n", ")\n", "axes[1].set(\n", " xlabel=\"MRL expérimentale\",\n", " ylabel=\"Densité\",\n", " title=\"Distribution de la cible\",\n", ")\n", "\n", "for ax in axes:\n", " sns.despine(ax=ax)\n", "fig.tight_layout()\n", "plt.show()" ], "id": "cell-compare-read-partitions" }, { "cell_type": "markdown", "metadata": {}, "source": [ "Nous réentraînons maintenant la même architecture à partir de zéro. La\n", "transformation de la cible est de nouveau ajustée uniquement à l’aide\n", "des données d’entraînement." ], "id": "9b002222-a53a-4c8f-a523-bcf2baeb008f" }, { "cell_type": "code", "execution_count": 33, "metadata": {}, "outputs": [ { "output_type": "stream", "name": "stdout", "text": [ "Caractéristiques d'entraînement : (220000, 50, 4) float32\n", "Caractéristiques de validation : (27500, 50, 4) float32\n", "Caractéristiques de test : (27500, 50, 4) float32" ] } ], "source": [ "read_target_scaler = StandardScaler()\n", "\n", "read_y_train = read_target_scaler.fit_transform(\n", " read_train_data[[\"mrl\"]]\n", ").astype(np.float32)\n", "read_y_validation = read_target_scaler.transform(\n", " read_validation_data[[\"mrl\"]]\n", ").astype(np.float32)\n", "\n", "read_x_train = one_hot_encode(read_train_data[\"utr\"])\n", "read_x_validation = one_hot_encode(read_validation_data[\"utr\"])\n", "read_x_test = one_hot_encode(read_test_data[\"utr\"])\n", "\n", "print(\"Caractéristiques d'entraînement :\", read_x_train.shape, read_x_train.dtype)\n", "print(\"Caractéristiques de validation :\", read_x_validation.shape, read_x_validation.dtype)\n", "print(\"Caractéristiques de test :\", read_x_test.shape, read_x_test.dtype)" ], "id": "prepare-read-count-arrays" }, { "cell_type": "markdown", "metadata": {}, "source": [ "La graine, l’architecture, l’optimiseur, la taille des lots, la règle\n", "d’arrêt précoce et le nombre maximal d’époques sont identiques à ceux de\n", "l’expérience 1. La stratégie de partitionnement des données constitue\n", "donc la seule différence expérimentale intentionnelle.\n", "\n", "> **Durée d’exécution prévue**\n", ">\n", "> Cette cellule effectue un second entraînement complet du modèle et\n", "> double donc approximativement la durée totale d’exécution du notebook." ], "id": "6f52f6fd-11d8-46d5-b983-e46baf21fa93" }, { "cell_type": "code", "execution_count": 34, "metadata": {}, "outputs": [ { "output_type": "stream", "name": "stdout", "text": [ "Epoch 1/50\n", "860/860 - 77s - 90ms/step - loss: 0.3927 - pearson_r: 0.7804 - r2: 0.6073 - val_loss: 0.2668 - val_pearson_r: 0.8870 - val_r2: 0.7331\n", "Epoch 2/50\n", "860/860 - 76s - 88ms/step - loss: 0.1746 - pearson_r: 0.9085 - r2: 0.8254 - val_loss: 0.1496 - val_pearson_r: 0.9342 - val_r2: 0.8503\n", "Epoch 3/50\n", "860/860 - 76s - 88ms/step - loss: 0.1392 - pearson_r: 0.9278 - r2: 0.8608 - val_loss: 0.1322 - val_pearson_r: 0.9383 - val_r2: 0.8677\n", "Epoch 4/50\n", "860/860 - 76s - 88ms/step - loss: 0.1279 - pearson_r: 0.9339 - r2: 0.8721 - val_loss: 0.1206 - val_pearson_r: 0.9416 - val_r2: 0.8793\n", "Epoch 5/50\n", "860/860 - 76s - 88ms/step - loss: 0.1217 - pearson_r: 0.9372 - r2: 0.8783 - val_loss: 0.1146 - val_pearson_r: 0.9433 - val_r2: 0.8854\n", "Epoch 6/50\n", "860/860 - 76s - 88ms/step - loss: 0.1175 - pearson_r: 0.9394 - r2: 0.8825 - val_loss: 0.1110 - val_pearson_r: 0.9446 - val_r2: 0.8889\n", "Epoch 7/50\n", "860/860 - 76s - 88ms/step - loss: 0.1136 - pearson_r: 0.9415 - r2: 0.8864 - val_loss: 0.1087 - val_pearson_r: 0.9462 - val_r2: 0.8913\n", "Epoch 8/50\n", "860/860 - 76s - 88ms/step - loss: 0.1110 - pearson_r: 0.9429 - r2: 0.8890 - val_loss: 0.1100 - val_pearson_r: 0.9458 - val_r2: 0.8899\n", "Epoch 9/50\n", "860/860 - 76s - 88ms/step - loss: 0.1087 - pearson_r: 0.9441 - r2: 0.8913 - val_loss: 0.1046 - val_pearson_r: 0.9471 - val_r2: 0.8954\n", "Epoch 10/50\n", "860/860 - 76s - 88ms/step - loss: 0.1067 - pearson_r: 0.9451 - r2: 0.8933 - val_loss: 0.1044 - val_pearson_r: 0.9478 - val_r2: 0.8955\n", "Epoch 11/50\n", "860/860 - 76s - 88ms/step - loss: 0.1041 - pearson_r: 0.9465 - r2: 0.8959 - val_loss: 0.1086 - val_pearson_r: 0.9475 - val_r2: 0.8914\n", "Epoch 12/50\n", "860/860 - 76s - 89ms/step - loss: 0.1023 - pearson_r: 0.9475 - r2: 0.8977 - val_loss: 0.1035 - val_pearson_r: 0.9482 - val_r2: 0.8965\n", "Epoch 13/50\n", "860/860 - 76s - 88ms/step - loss: 0.1006 - pearson_r: 0.9484 - r2: 0.8994 - val_loss: 0.1051 - val_pearson_r: 0.9476 - val_r2: 0.8949\n", "Epoch 14/50\n", "860/860 - 77s - 89ms/step - loss: 0.0984 - pearson_r: 0.9495 - r2: 0.9016 - val_loss: 0.1038 - val_pearson_r: 0.9483 - val_r2: 0.8962\n", "Epoch 15/50\n", "860/860 - 76s - 89ms/step - loss: 0.0964 - pearson_r: 0.9506 - r2: 0.9036 - val_loss: 0.1083 - val_pearson_r: 0.9468 - val_r2: 0.8916\n", "Epoch 16/50\n", "860/860 - 77s - 89ms/step - loss: 0.0946 - pearson_r: 0.9515 - r2: 0.9054 - val_loss: 0.1047 - val_pearson_r: 0.9482 - val_r2: 0.8953\n", "Epoch 17/50\n", "860/860 - 76s - 89ms/step - loss: 0.0931 - pearson_r: 0.9523 - r2: 0.9069 - val_loss: 0.1038 - val_pearson_r: 0.9475 - val_r2: 0.8961\n", "Epoch 17: early stopping\n", "Restoring model weights from the end of the best epoch: 12." ] } ], "source": [ "# Réinitialise l'état aléatoire avant de construire le second modèle.\n", "keras.utils.set_random_seed(SEED)\n", "read_count_model = compile_mrl_cnn(build_mrl_cnn())\n", "read_count_model_path = MODEL_DIR / \"mrl_cnn_read_count_split.keras\"\n", "\n", "read_history = read_count_model.fit(\n", " read_x_train,\n", " read_y_train,\n", " validation_data=(read_x_validation, read_y_validation),\n", " epochs=MAX_EPOCHS,\n", " batch_size=BATCH_SIZE,\n", " shuffle=True,\n", " callbacks=[\n", " callbacks.EarlyStopping(\n", " monitor=\"val_loss\",\n", " mode=\"min\",\n", " patience=PATIENCE,\n", " restore_best_weights=True,\n", " verbose=1,\n", " ),\n", " callbacks.ModelCheckpoint(\n", " filepath=str(read_count_model_path),\n", " monitor=\"val_loss\",\n", " mode=\"min\",\n", " save_best_only=True,\n", " ),\n", " ],\n", " verbose=2,\n", ")" ], "id": "train-read-count-model" }, { "cell_type": "markdown", "metadata": {}, "source": [ "Nous examinons d’abord la dynamique d’apprentissage et vérifions que\n", "l’arrêt précoce a sélectionné un modèle dont les performances de\n", "validation sont stables." ], "id": "79857753-0cce-4368-a05f-2d92cefaa999" }, { "cell_type": "code", "execution_count": 35, "metadata": {}, "outputs": [ { "output_type": "display_data", "metadata": {}, "data": {} } ], "source": [ "plot_training_history(read_history)\n", "plt.show()" ], "id": "cell-plot-read-count-training-history" }, { "cell_type": "markdown", "metadata": {}, "source": [ "L’évaluation finale utilise les 27 500 séquences ayant les nombres de\n", "lectures les plus élevés. Le tableau comparatif place ces résultats à\n", "côté de ceux de l’ensemble de test aléatoire." ], "id": "afeeb75f-45aa-412a-b39e-cdffc391ef9c" }, { "cell_type": "code", "execution_count": 36, "metadata": {}, "outputs": [ { "output_type": "display_data", "metadata": {}, "data": { "text/html": [ "\n", "" ] } } ], "source": [ "read_test_observed = read_test_data[\"mrl\"].to_numpy()\n", "read_test_predicted = predict_mrl(\n", " read_count_model,\n", " read_x_test,\n", " read_target_scaler,\n", ")\n", "\n", "split_comparison = pd.DataFrame(\n", " [\n", " test_metrics,\n", " regression_metrics(read_test_observed, read_test_predicted),\n", " ],\n", " index=[\"Sélection aléatoire du test\", \"Sélection du test selon les lectures\"],\n", ")\n", "split_comparison.round(4)" ], "id": "cell-evaluate-read-count-test-set" }, { "cell_type": "code", "execution_count": 37, "metadata": {}, "outputs": [ { "output_type": "display_data", "metadata": {}, "data": {} } ], "source": [ "plot_observed_vs_predicted(\n", " read_test_observed,\n", " read_test_predicted,\n", " title=\"Ensemble de test sélectionné selon le nombre de lectures\",\n", ")\n", "plt.show()" ], "id": "cell-plot-read-count-test-predictions" }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Évaluation sur l’expérience indépendante menée sur des lymphocytes T\n", "\n", "Nous appliquons maintenant le second modèle aux mêmes 100 000 séquences\n", "de lymphocytes T fusionnées qui ont servi à évaluer l’expérience 1. La\n", "réutilisation d’un ensemble d’évaluation externe identique permet de\n", "comparer directement l’effet de la stratégie de partitionnement de\n", "l’entraînement. Les prédictions du second modèle doivent être ramenées à\n", "l’échelle expérimentale de la MRL à l’aide de la transformation de la\n", "cible ajustée sur sa propre partition d’entraînement." ], "id": "d63817b4-baa1-44a3-9174-9f9ae8c33631" }, { "cell_type": "code", "execution_count": 38, "metadata": {}, "outputs": [ { "output_type": "display_data", "metadata": {}, "data": { "text/html": [ "\n", "" ] } } ], "source": [ "read_tcell_predicted = predict_mrl(\n", " read_count_model,\n", " x_tcell,\n", " read_target_scaler,\n", ")\n", "read_tcell_metrics = regression_metrics(\n", " tcell_observed,\n", " read_tcell_predicted,\n", ")\n", "\n", "cross_cell_comparison = pd.DataFrame(\n", " [\n", " tcell_metrics,\n", " read_tcell_metrics,\n", " ],\n", " index=[\n", " \"Expérience 1 : partition aléatoire\",\n", " \"Expérience 2 : partition selon les lectures\",\n", " ],\n", ")\n", "cross_cell_comparison.round(4)" ], "id": "cell-evaluate-read-count-model-on-tcell" }, { "cell_type": "code", "execution_count": 39, "metadata": {}, "outputs": [ { "output_type": "display_data", "metadata": {}, "data": {} } ], "source": [ "plot_observed_vs_predicted(\n", " tcell_observed,\n", " read_tcell_predicted,\n", " title=\"Lymphocytes T indépendants — modèle partitionné selon les lectures\",\n", ")\n", "plt.show()" ], "id": "cell-plot-read-count-model-tcell-predictions" }, { "cell_type": "markdown", "metadata": {}, "source": [ "Enfin, nous sauvegardons la transformation de la cible associée au\n", "second modèle. Le meilleur modèle lui-même a été sauvegardé par\n", "`ModelCheckpoint` pendant l’entraînement." ], "id": "a6f9efce-011c-456e-9fdb-18826ea54d41" }, { "cell_type": "code", "execution_count": 40, "metadata": {}, "outputs": [ { "output_type": "stream", "name": "stdout", "text": [ "Modèle sauvegardé : models/mrl_cnn_read_count_split.keras\n", "Transformation de la cible sauvegardée : models/mrl_target_scaler_read_count_split.npz" ] } ], "source": [ "read_scaler_path = MODEL_DIR / \"mrl_target_scaler_read_count_split.npz\"\n", "np.savez(\n", " read_scaler_path,\n", " mean=read_target_scaler.mean_.astype(np.float32),\n", " scale=read_target_scaler.scale_.astype(np.float32),\n", ")\n", "\n", "print(\"Modèle sauvegardé :\", read_count_model_path)\n", "print(\"Transformation de la cible sauvegardée :\", read_scaler_path)" ], "id": "save-read-count-artifacts" }, { "cell_type": "markdown", "metadata": {}, "source": [ "Les deux expériences estiment des formes différentes de généralisation.\n", "Leurs scores doivent donc être interprétés en parallèle avec les\n", "distributions de la cible affichées. En particulier, $R^2$ dépend de la\n", "variance des valeurs observées dans son ensemble de test; une variation\n", "de $R^2$ ne peut donc pas être automatiquement attribuée à une variation\n", "de la seule qualité du modèle.\n", "\n", "# Références\n", "\n", "- Sample, P. J., Wang, B., Reid, D. W., Presnyak, V., McFadyen, I. J.,\n", " Morris, D. R., & Seelig, G. (2019). Human 5′ UTR design and variant\n", " effect prediction from a massively parallel translation assay. *Nature\n", " Biotechnology*, *37*(7), 803–809.\n", " \n", "- Castillo-Hair, S., Fedak, S., Wang, B., Linder, J., Havens, K., Certo,\n", " M., & Seelig, G. (2024). Optimizing 5′ UTRs for mRNA-delivered gene\n", " editing using deep learning. *Nature Communications*, *15*(1), 5284.\n", " \n", "- Tang, X., Huo, M., Chen, Y., et al. (2024). A novel deep generative\n", " model for mRNA vaccine development: Designing 5′ UTRs with\n", " N1-methyl-pseudouridine modification. *Acta Pharmaceutica Sinica B*,\n", " *14*(4), 1814–1826. " ], "id": "9bece1ec-531b-4bd3-97bf-b8c646cfd5e5" } ], "nbformat": 4, "nbformat_minor": 5, "metadata": { "kernelspec": { "name": "python3", "display_name": "Python 3 (ipykernel)", "language": "python", "path": "/Users/turcotte/.virtualenvs/csi4106/share/jupyter/kernels/python3" } } }