{ "cells": [ { "cell_type": "markdown", "id": "4ac6391a", "metadata": {}, "source": [ "# Tutorial 2 - Differential expression analysis\n", "\n", "This tutorial shows how to perform differential expression analysis (DEA) between a set of \"treated\" samples and their corresponding in silico \"control\" samples, found using the bulkdgd model (see :doc:`Tutorial 1 ` for how representations are found). For each sample, the model's prediction acts as the \"normal\"/control counterpart against which the observed gene counts are compared.\n", "\n", "As a concrete example, we use 10 real breast cancer (TCGA-BRCA) samples: their representations, predicted means, and predicted r-values were obtained by running `bulkdgd_find_representations` with the default (pre-trained) bulkdgd model on these samples ahead of time." ] }, { "cell_type": "code", "execution_count": 1, "id": "50aeea55", "metadata": { "execution": { "iopub.execute_input": "2026-07-05T18:44:08.142177Z", "iopub.status.busy": "2026-07-05T18:44:08.141753Z", "iopub.status.idle": "2026-07-05T18:44:13.850612Z", "shell.execute_reply": "2026-07-05T18:44:13.849118Z" } }, "outputs": [], "source": [ "# Import the 'logging' module\n", "import logging as log\n", "# Import the 'analysis.dea' module\n", "from bulkdgd.analysis import dea\n", "# Import the 'ioutil' package\n", "from bulkdgd import ioutil\n", "\n", "# Set the logging options so that every message of level INFO or above\n", "# is emitted.\n", "log.basicConfig(level = \"INFO\")" ] }, { "cell_type": "markdown", "id": "446513a2", "metadata": {}, "source": [ "## Load the samples\n", "\n", "Load the preprocessed samples (the observed gene counts) produced in Tutorial 1, keeping only the first ten for this example." ] }, { "cell_type": "code", "execution_count": 2, "id": "4753948d", "metadata": { "execution": { "iopub.execute_input": "2026-07-05T18:44:13.854893Z", "iopub.status.busy": "2026-07-05T18:44:13.854374Z", "iopub.status.idle": "2026-07-05T18:44:14.058428Z", "shell.execute_reply": "2026-07-05T18:44:14.057014Z" } }, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "INFO:bulkdgd.ioutil.samplesio:Since 'keep_samples_names = True', the samples will be identified using the values contained in the first column of the data frame.\n" ] }, { "name": "stderr", "output_type": "stream", "text": [ "INFO:bulkdgd.ioutil.samplesio:14895 column(s) containing gene expression data was (were) found in the input data frame.\n" ] } ], "source": [ "# Load the preprocessed samples into a data frame.\n", "df_samples = \\\n", " ioutil.load_samples(# The CSV file where the samples are stored\n", " csv_file = \"samples_preprocessed.csv\",\n", " # The field separator used in the CSV file\n", " sep = \",\",\n", " # Whether to keep the original samples' names/\n", " # indexes (if True, they are assumed to be in\n", " # the first column of the data frame) \n", " keep_samples_names = True,\n", " # Whether to split the input data frame into\n", " # two data frames, one containing only gene\n", " # expression data and the other containing\n", " # additional information about the samples\n", " split = False)\n", "\n", "# Get only the first ten rows.\n", "df_samples = df_samples.iloc[:10,:]" ] }, { "cell_type": "markdown", "id": "f503c85e", "metadata": {}, "source": [ "## Load the predicted scaled means\n", "\n", "Load the decoder's predicted scaled means for each sample found in Tutorial 1 - these represent the in silico control counterpart of each treated sample." ] }, { "cell_type": "code", "execution_count": 3, "id": "73dff7d4", "metadata": { "execution": { "iopub.execute_input": "2026-07-05T18:44:14.062633Z", "iopub.status.busy": "2026-07-05T18:44:14.062423Z", "iopub.status.idle": "2026-07-05T18:44:14.286162Z", "shell.execute_reply": "2026-07-05T18:44:14.284973Z" } }, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "INFO:bulkdgd.ioutil.decoutio:14895 column(s) containing the decoder's outputs was (were) found in the input data frame.\n" ] } ], "source": [ "# Load the predicted scaled means into a data frame.\n", "df_pred_means = \\\n", " ioutil.load_decoder_outputs(# The CSV file where the predicted\n", " # scaled means are stored\n", " csv_file = \"pred_means.csv\",\n", " # The field separator used in the CSV\n", " # file\n", " sep = \",\",\n", " # Whether to split the input data frame\n", " # into two data frames, one containing\n", " # only the predicted scaled means and\n", " # the other containing additional\n", " # information about the original samples\n", " split = False)\n", "\n", "# Get only the first ten rows.\n", "df_pred_means = df_pred_means.iloc[:10,:]" ] }, { "cell_type": "markdown", "id": "0dc119ca", "metadata": {}, "source": [ "## Load the predicted r-values\n", "\n", "Load the decoder's predicted r-values for each sample, needed because the genes' counts are modeled with negative binomial distributions." ] }, { "cell_type": "code", "execution_count": 4, "id": "13cebf0c", "metadata": { "execution": { "iopub.execute_input": "2026-07-05T18:44:14.289314Z", "iopub.status.busy": "2026-07-05T18:44:14.288683Z", "iopub.status.idle": "2026-07-05T18:44:14.814208Z", "shell.execute_reply": "2026-07-05T18:44:14.812968Z" } }, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "INFO:bulkdgd.ioutil.decoutio:14895 column(s) containing the decoder's outputs was (were) found in the input data frame.\n" ] } ], "source": [ "# Load the predicted r-values into a data frame.\n", "df_pred_r_values = \\\n", " ioutil.load_decoder_outputs(# The CSV file where the predicted\n", " # r-values are stored\n", " csv_file = \"pred_r_values.csv\",\n", " # The field separator used in the CSV\n", " # file\n", " sep = \",\",\n", " # Whether to split the input data frame\n", " # into two data frames, one containing\n", " # only the predicted r-values and\n", " # the other containing additional\n", " # information about the original samples\n", " split = False)\n", "\n", "# Get only the first ten rows.\n", "df_pred_r_values = df_pred_r_values.iloc[:10,:]" ] }, { "cell_type": "markdown", "id": "48c722e6", "metadata": {}, "source": [ "## Differential expression analysis\n", "\n", "For each sample, compare its observed gene counts against the model's in silico control prediction, computing p-values, q-values (adjusted p-values), and log2-fold changes for every gene, and save the results to one CSV file per sample." ] }, { "cell_type": "code", "execution_count": 5, "id": "1b7a4fdc", "metadata": { "execution": { "iopub.execute_input": "2026-07-05T18:44:14.818195Z", "iopub.status.busy": "2026-07-05T18:44:14.817704Z", "iopub.status.idle": "2026-07-05T18:45:35.764322Z", "shell.execute_reply": "2026-07-05T18:45:35.762915Z" } }, "outputs": [], "source": [ "# For each sample\n", "for sample in df_samples.index:\n", "\n", " # Perform differential expression analysis.\n", " dea_results, _ = \\\n", " dea.get_statistics(# The observed gene counts for the current\n", " # sample\n", " obs_counts = df_samples.loc[sample,:],\n", " # The predicted scaled means for the current\n", " # sample\n", " pred_means = df_pred_means.loc[sample,:],\n", " # The r-values for the current sample\n", " r_values = df_pred_r_values.loc[sample,:],\n", " # Which statistics should be computed and\n", " # included in the results\n", " statistics = \\\n", " [\"p_values\", \"q_values\",\n", " \"log2_fold_changes\"],\n", " # The resolution for the p-values calculation\n", " # (the higher, the more accurate the\n", " # calculation; set to 'None' for an exact\n", " # calculation)\n", " resolution = 1e4,\n", " # The family-wise error rate for the\n", " # calculation of the q-values\n", " alpha = 0.05,\n", " # The method used to calculate the q-values\n", " method = \"fdr_bh\")\n", "\n", " # Save the results.\n", " dea_results.to_csv(# The CSV file where to save the results\n", " # for the current sample\n", " f\"dea_sample_{sample}.csv\",\n", " # The field separator to use in the output\n", " # CSV file\n", " sep = \",\",\n", " # Whether to keep the rows' names\n", " index = True,\n", " # Whether to keep the columns' names\n", " header = True)" ] } ], "metadata": { "kernelspec": { "display_name": "bulkdgd2-env", "language": "python", "name": "bulkdgd2-env" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.13.14" } }, "nbformat": 4, "nbformat_minor": 5 }