DA-Example-Tutorial.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "b2d9f710-01cb-4972-b9f9-273e95285687",
   "metadata": {},
   "source": [
    "## Tutorial of DeepAdapter\n",
    "### A self-adaptive and versatile tool for eliminating multiple undesirable variations from transcriptome\n",
    "In this notebook, you will learn how to re-train DeepAdapter with the example dataset.\n",
    "\n",
    "## 1. Installation and requirements\n",
    "### 1.1. Installation\n",
    "To run locally, please open a terminal window and download the code with:\n",
    "```sh\n",
    "$ # create a new conda environment\n",
    "$ conda create -n DA python=3.9\n",
    "$ # activate environment\n",
    "$ conda activate DA\n",
    "$ # Install dependencies\n",
    "$ pip install deepadapter\n",
    "$ # Launch jupyter notebook\n",
    "$ jupyter notebook\n",
    "```\n",
    "### 1.2. Download datasets\n",
    "Please download the open datasets in [Zenodo](https://zenodo.org/records/10494751).\n",
    "These datasets are collected from literatures to demonstrate multiple unwanted variations, including:\n",
    "* batch datasets: LINCS-DToxS ([van Hasselt et al. Nature Communications, 2020](https://www.nature.com/articles/s41467-020-18396-7)) and Quartet project ([Yu, Y. et al. Nature Biotechnology, 2023](https://www.nature.com/articles/s41587-023-01867-9)).\n",
    "* platform datasets: profiles from microarray ([Iorio, F. et al. Cell, 2016](https://www.cell.com/cell/pdf/S0092-8674(16)30746-2.pdf)) and RNA-seq ([Ghandi, M. et al. Nature, 2019](https://www.nature.com/articles/s41586-019-1186-3)).\n",
    "* purity datasets: profiles from cancer cell lines ([Ghandi, M. et al. Nature, 2019](https://www.nature.com/articles/s41586-019-1186-3)) and tissues ([Weinstein, J.N. et al. Nature genetics, 2013](https://www.nature.com/articles/ng.2764)).\n",
    "\n",
    "After downloading, place the datasets in the `data/` directory located in the same hierarchy as this tutorial.\n",
    "* batch datasets: `data/batch_data/`\n",
    "* platform datasets: `data/platform_data/`\n",
    "* purity datasets: `data/purity_data/`\n",
    "  \n",
    "**Putting datasets in the right directory is important for loading the example datasets successfully.**\n",
    "\n",
    "To execute a \"cell\", please press Shift+Enter"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e98c179c-882b-4e1e-aba7-0594a6647bc5",
   "metadata": {},
   "source": [
    "## 2. Load the datasets and preprocess\n",
    "### 2.1. load the modules"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fbc07082-e02d-431f-9fba-b5b78c85fe24",
   "metadata": {},
   "outputs": [],
   "source": [
    "%load_ext autoreload\n",
    "%autoreload 2\n",
    "\n",
    "import os, sys\n",
    "import pandas as pd\n",
    "import numpy as np\n",
    "\n",
    "from deepadapter.utils import data_utils as DT\n",
    "from deepadapter.utils import utils as UT\n",
    "from deepadapter.utils import triplet as TRP\n",
    "from deepadapter.models.trainer import Trainer\n",
    "from deepadapter.models.data_loader import TransData, DataLoader\n",
    "from deepadapter.models.dl_utils import AE, FBatch\n",
    "from deepadapter.params import dl_params as DLPARAM\n",
    "from deepadapter.run import run as RUN"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "56913165-41eb-4ce8-a8dc-496d7dc71f07",
   "metadata": {},
   "source": [
    "### 2.2. Load the demonstrated datasets\n",
    "We ultilize Batch-LINCS for demonstration. To load datasets of platform and purity variations, please download them in Zenodo (https://zenodo.org/records/10494751).\n",
    "  * In the tutorial, we have **data** for gene expression, **batches** for unwanted variations, and **donors** for biological signals.\n",
    "  * In training your own DeepAdapter, please refer to `DeepAdapter-YourOwnData-Tutorial.ipynb`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "45c4e81d-f681-43bb-b299-e7953ff7e2ad",
   "metadata": {},
   "outputs": [],
   "source": [
    "loadTransData = DT.LoadTransData()\n",
    "data, batches, wells, donors, infos, test_infos = loadTransData.load_lincs_lds1593()\n",
    "ids = np.arange(len(data))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7845814c-c265-4bf1-b662-a685d52933a9",
   "metadata": {},
   "source": [
    "### 2.3. Preprocess the transcriptomic data\n",
    "The gene expression profiles are preprocessed by sample normalization, gene ranking, and log normalization. Let $S_i = \\sum_l x_{i l}$ denote the sum over all genes. In sample normalization, we divide $S_i$ for every sample and multiply a constant 10000 ([Xiaokang Yu et al. Nature communications, 2023](https://www.nature.com/articles/s41467-023-36635-5)):\n",
    "$$x_{i l} = \\frac{x_{i l}}{S_i} 10^4.$$\n",
    "Then, we sort genes by their expression levels and perform the log transformation $x_{i l} = \\log {(x_{i l} + 1)}$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d387d965-c88f-4aed-95cd-b6dbf2d0017b",
   "metadata": {},
   "outputs": [],
   "source": [
    "prepTransData = DT.PrepTransData()\n",
    "raw_df = prepTransData.sample_norm(data)\n",
    "raw_df, sorted_cols = prepTransData.sort_genes_sgl_df(raw_df)\n",
    "input_arr = prepTransData.sample_log(raw_df)\n",
    "bat2label, label2bat, unwanted_labels, unwanted_onehot = prepTransData.label2onehot(batches)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4ed7bdaf-271f-4754-af3d-fe7bb4d93f25",
   "metadata": {},
   "source": [
    "## 3. Train DeepAdapter\n",
    "### 3.1. Adjust DeepAdapter's parameters\n",
    "The parameters for DeepAdapter are as follows (**Note: you can revise parameter directly in `net_args`, e.g., `net_args.batch_size = 512`.**):\n",
    "* **epochs**: the total training epochs of DeepAdapter, default = $150000$\n",
    "* **ae_epochs**: the warmup epochs of autoencoder in DeepAdapter, default = $400$\n",
    "* **batch_epochs**: the warmup epochs of discriminator in DeepAdapter, default = $50$\n",
    "* **batch_size**: the batch size of dataloader, default = $256$\n",
    "* **hidden_dim**: the hidden units of autoencoder in DeepAdapter, default = $256$\n",
    "* **z_dim**: the latent units of autoencoder in DeepAdapter, default = $128$\n",
    "* **drop**: the dropout rate of DeepAdapter, default = $0.3$\n",
    "* **lr_lower_ae**: the lower learning rate of autoencoder in DeepAdapter, default = $1e-5$\n",
    "* **lr_upper_ae**: the upper learning rate of autoencoder in DeepAdapter, default = $5e-4$\n",
    "* **lr_lower_batch**: the lower learning rate of discriminator in DeepAdapter, default = $1e-5$\n",
    "* **lr_upper_batch**: the upper learning rate of discriminator in DeepAdapter, default = $5e-4$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "64a84343-cb36-4a4d-93c6-7637eef175db",
   "metadata": {},
   "outputs": [],
   "source": [
    "net_args = DLPARAM.load_dl_params()\n",
    "\n",
    "## demonstration abou how to modify parameteres\n",
    "print(net_args.batch_size)\n",
    "net_args.batch_size = 512\n",
    "print(net_args.batch_size)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "61d92d21-32e7-438c-94a2-606045a2ac21",
   "metadata": {},
   "source": [
    "### 3.2. Split dataset\n",
    "* For the tutorial, we extract the biosamples across all batches as the test set; then split the rest into training and validation set randomly.</br>\n",
    "That means the training data seen by DeepAdapter doesn't disperse across all unwanted variations while the testing data does.</br>\n",
    "Acutally, this split method could increase the training difficulty.\n",
    "* For your own dataset, you can split the dataset randomly using the function `DT.data_split_random`.</br>\n",
    "Please refer to `DeepAdapter-YourOwnData-Tutorial.ipynb`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fdeb151d-5a9f-414a-a25b-4dca43a7b7cf",
   "metadata": {},
   "outputs": [],
   "source": [
    "train_data, train_labels, train_labels_hot, \\\n",
    "    val_data, val_labels, val_labels_hot, \\\n",
    "    test_data, test_labels, test_labels_hot, \\\n",
    "    train_ids, val_ids, test_ids, \\\n",
    "    tot_train_val_idxs, tot_train_idxs, tot_val_idxs, tot_test_idxs = DT.data_split_lds1593(input_arr, unwanted_labels, unwanted_onehot, ids, infos, test_infos)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "32181bab-39f6-49c5-8b11-65a00562fa18",
   "metadata": {},
   "outputs": [],
   "source": [
    "train_bios, val_bios, test_bios = donors[tot_train_idxs], donors[tot_val_idxs], donors[tot_test_idxs]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c6b4302d-0e87-40a3-9793-e10fc3e66bfb",
   "metadata": {},
   "outputs": [],
   "source": [
    "bio_label2bat = {t:t for t in set(train_bios)}"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "15aaab98-836f-4f74-9b4a-e2af4c65fd52",
   "metadata": {},
   "source": [
    "### 3.3. Train DeepAdapter\n",
    "Two options are provided for training DeepAdapter. If you want to learn the training process, please train it step by step. If you want to skip these initializations, please use the one-line code :)!\n",
    "* To train DeepAdapter step by step, you need to initialize models, dataloaders, trainer, and the mutual nearest neighbors.\n",
    "* To train DeepAdapter in one-line code, just utilize `deepAdapter.run.train()`."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2ce7a975-f0d9-48ed-8968-87e40d85642e",
   "metadata": {},
   "source": [
    "#### 3.3.1. Train it step by step"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6ba951b1-797d-409e-ab41-e5a0e048471e",
   "metadata": {},
   "outputs": [],
   "source": [
    "db_name = \"LDS1593\"\n",
    "out_dir = os.path.join(\"DA_model/batch\", \"stepByStep_{}/\".format(db_name))\n",
    "os.makedirs(out_dir, exist_ok = True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f8d2801b-cce2-4d53-921f-4875ec4140cc",
   "metadata": {},
   "outputs": [],
   "source": [
    "## initialize models\n",
    "in_dim = input_arr.shape[1]\n",
    "num_unw_vars = len(bat2label)\n",
    "ae = AE(in_dim, net_args.hidden_dim, num_unw_vars, net_args.z_dim, net_args.drop).cuda()\n",
    "fbatch = FBatch(net_args.hidden_dim, num_unw_vars, net_args.z_dim, net_args.drop).cuda()\n",
    "\n",
    "## initialize dataloaders\n",
    "train_trans = TransData(train_data, train_labels, train_bios, train_ids, train_labels_hot)\n",
    "train_loader = DataLoader(train_trans, batch_size = net_args.batch_size, collate_fn = train_trans.collate_fn, shuffle = True, drop_last = False)\n",
    "val_trans = TransData(val_data, val_labels, val_bios, val_ids, val_labels_hot)\n",
    "val_loader = DataLoader(val_trans, batch_size = net_args.batch_size, collate_fn = val_trans.collate_fn, shuffle = False, drop_last = False)\n",
    "test_trans = TransData(test_data, test_labels, test_bios, test_ids, test_labels_hot)\n",
    "test_loader = DataLoader(test_trans, batch_size = net_args.batch_size, collate_fn = test_trans.collate_fn, shuffle = False, drop_last = False)\n",
    "\n",
    "## initialize trainer\n",
    "trainer = Trainer(train_loader, val_loader, test_loader, ae, fbatch, bio_label2bat, label2bat, net_args, out_dir)\n",
    "\n",
    "## initialize mutual nearest neighbors\n",
    "train_mutuals = TRP.find_MNN_cosine_kSources(train_data, train_labels, train_ids)\n",
    "val_mutuals = TRP.find_MNN_cosine_kSources(val_data, val_labels, val_ids)\n",
    "\n",
    "## begin training!\n",
    "trainer.fit(train_mutuals, val_mutuals)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "eeb97409-d395-4289-86db-36e6fd5a341d",
   "metadata": {},
   "source": [
    "#### 3.3.2. Train it in one-line code\n",
    "Parameters for one-line code training:\n",
    "* **train_list**: the list of training transcriptomic profiles, unwanted variations, biological signals, data ids, and onehot representations of unwanted variations.\n",
    "* **val_list**: the list of validation transcriptomic profiles, unwanted variations, biological signals, data ids, and onehot representations of unwanted variations.\n",
    "* **test_list**: the list of testing transcriptomic profiles, unwanted variations, biological signals, data ids, and onehot representations of unwanted variations.\n",
    "* **label2unw**: the dictionary which maps unwanted labels (e.g., 0, 1 ...) to unwanted variations (e.g., batch1, batch2 ...)\n",
    "* **label2wnt**: the dictionary which maps biological labels (e.g., 0, 1 ...) to biological annotations (e.g., donor1, donor2 ...)\n",
    "* **net_args**: the parameters to construct DeepAdapter\n",
    "* **out_dir**: the out directory for saved models and logged losses."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ca6f06f3-34bd-4124-8fce-4ac24f7ad1d9",
   "metadata": {},
   "outputs": [],
   "source": [
    "db_name = \"LDS1593\"\n",
    "out_dir = os.path.join(\"DA_model/batch\", \"oneLineCode_{}/\".format(db_name))\n",
    "os.makedirs(out_dir, exist_ok = True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4aa8eed2-8ef2-4105-b352-ebdbfb6763d4",
   "metadata": {},
   "outputs": [],
   "source": [
    "train_list = [train_data, train_labels, train_bios, train_ids, train_labels_hot]\n",
    "val_list = [val_data, val_labels, val_bios, val_ids, val_labels_hot]\n",
    "test_list = [test_data, test_labels, test_bios, test_ids, test_labels_hot]\n",
    "\n",
    "trainer = RUN.train(\n",
    "    train_list = train_list, \n",
    "    val_list = val_list, \n",
    "    test_list = test_list, \n",
    "    label2unw = label2bat, \n",
    "    label2wnt = bio_label2bat, \n",
    "    net_args = net_args, \n",
    "    out_dir = out_dir)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1660e05c-41e7-4673-843c-08ab3c35a9c5",
   "metadata": {},
   "source": [
    "## 4. Align the data\n",
    "### 4.1. Load trained model & quantatitive evaluation\n",
    "* Step 1: load the best-trained model\n",
    "* Step 2: utilize `trainer.evaluate()`\n",
    "\n",
    "In `trainer.evaluate()`, we perform decomposition analysis of aligned data and perform the quantatitive analysis including alignment score, ASW, NMI, and ARI calcuation. The quantatitive results are recorded in `record_path`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5b4e13a5-9df6-458d-84e1-1343ed016454",
   "metadata": {},
   "outputs": [],
   "source": [
    "trainer.load_trained_ae(os.path.join(out_dir, \"ae.tar\"))\n",
    "\n",
    "record_path = os.path.join(out_dir, \"test_res.csv\")\n",
    "test_data, test_aligned_data, test_wnt_infs, test_unw_infs = trainer.evaluate(record_path, db_name, test_loader)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "dbe16151-a72a-4b70-8285-ea408df3f3ce",
   "metadata": {},
   "source": [
    "Additionally, you can perform any other analysis you like with the aligned data `aligned_data`!"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "47764a01-b524-4073-9a33-036a3dce6e9d",
   "metadata": {},
   "source": [
    "### 4.2. Save the aligned data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e4c40232-ca88-440f-a4f4-5bd74d11a597",
   "metadata": {},
   "outputs": [],
   "source": [
    "all_trans = TransData(np.vstack((train_data, val_data, test_data)), \n",
    "                      np.array(list(train_labels) + list(val_labels) + list(test_labels)),\n",
    "                      np.array(list(train_bios) + list(val_bios) + list(test_bios)),\n",
    "                      np.array(list(train_ids) + list(val_ids) + list(test_ids)), \n",
    "                      np.vstack((train_labels_hot, val_labels_hot, test_labels_hot)))\n",
    "all_loader = DataLoader(all_trans, batch_size = net_args.batch_size, collate_fn = all_trans.collate_fn, shuffle = False, drop_last = False)\n",
    "record_path = os.path.join(out_dir, \"res.csv\")\n",
    "data, aligned_data, wnt_infs, unw_infs = trainer.evaluate(record_path, db_name, all_loader)\n",
    "\n",
    "save_path = os.path.join(out_dir, \"DA_data.csv\")\n",
    "df = pd.DataFrame(data, columns = sorted_cols)\n",
    "df[\"ID\"] = np.array(list(train_ids) + list(val_ids) + list(test_ids))\n",
    "df[\"wantInfo\"] = wnt_infs\n",
    "df[\"unwantInfo\"] = unw_infs\n",
    "df.to_csv(save_path, index = False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ee68202d-ec0e-4eae-b249-c0e5ee665b71",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "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.9.19"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}