tutorials/Prepare_input_from_user_defined_trajectory.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Get started\n",
    "\n",
    "In this tutorial, we will go through how to prepare for the inputs required by TraSig, using user defined trajectory. Regardless of the pseudotime trajectory tool you used, to prepare the inputs for TraSig, you need:\n",
    "\n",
    "    * Expression matrix (row: cells; column: genes)\n",
    "    * Cluster / branch / edge label assigned to each cell\n",
    "    * Progression on the branch assigned to each cell, a value in [0, 1]\n",
    "    * Cell type label assigned to each cell (for visualization purpose only)\n",
    "    * Sampling time of cells (optional; if unknown, then all set to 0)\n",
    "\n",
    "To prepare for the expression matrix starting from the counts, after filtering no-quality cells / genes. You may normalize the counts using TPM (we use a scaling factor of 1e4) and then transform the values using log (x + 1).\n",
    "\n",
    "After obtaining these, follow the steps below to prepare the inputs for TraSig. We will use ti_slingshot (Slingshot in dynverse) on the dataset \"oligodendrocyte-differentiation-clusters_marques.rds\" as an example. \n",
    "\n",
    "Note that the emphasis of this tutorial is to **illustarte the input file formats**. You may use your own prefered pseudotime trajectory tools (not necessarily those implemented in dynverse). If you would like to learn the details on how to prepare inputs from dynverse trajectory outputs, refer to [Prepare_input_from_dynverse_ti_methods](Prepare_input_from_dynverse_ti_methods.ipynb) in the tutorials. \n",
    "\n",
    "**Table of Content**\n",
    "1. [Load expression and trajectory results](#1)\n",
    "2. [Prepare pseudotime results per cell for TraSig](#2)\n",
    "3. [Subsetting expression data (keep only relevant genes)](#3)\n",
    "4. [Save correspondence from sampling time to paths](#4)\n",
    "\n",
    "**Extra Package Requirements**\n",
    "* h5py >= 3.1.0 (required to load dynverse trajectory results)\n",
    "* rpy2 >= 3.3.6 (required to load dynverse datasets)\n",
    "* matplotlib-base >= 3.3.4 (required for plotting)\n",
    "* scikit-learn >= 0.23.2 (required for evaluating trajectory results)\n",
    "* scipy >= 1.5.4 (required to prepare sampling time input)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os, sys\n",
    "import argparse\n",
    "import time\n",
    "from os.path import exists\n",
    "import collections\n",
    "from typing import Iterable\n",
    "import pickle\n",
    "from collections import Counter\n",
    "import requests\n",
    "\n",
    "import pandas as pd\n",
    "import numpy as np\n",
    "from matplotlib import pyplot as plt\n",
    "import h5py \n",
    "import rpy2.robjects as robjects  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# example data set \n",
    "project = \"oligodendrocyte-differentiation-clusters_marques\"\n",
    "\n",
    "# set the path to the inputs for the trajectory inference (e.g. expression)\n",
    "input_path = \"../trajectory/input\"\n",
    "\n",
    "# set the path to save the outputs of this script (place to save inputs for TraSig)\n",
    "output_path = \"../example/input\"\n",
    "\n",
    "# set the path to the trajectory output\n",
    "trajectory_filename = f\"../trajectory/output/output.h5\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "# set the names for output files \n",
    "preprocess = \"None\"\n",
    "model_name = \"ti_slingshot\"\n",
    "others = \"None\"\n",
    "\n",
    "if preprocess != \"None\":\n",
    "    _preprocess = f\"_{preprocess}\"\n",
    "else:\n",
    "    _preprocess = \"\"\n",
    "    \n",
    "if others == \"None\":\n",
    "    condition = \"\"\n",
    "        \n",
    "suffix = f\"{_preprocess}_{model_name}{condition}\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 1. Load expression and trajectory results\n",
    "\n",
    "* You may need to customize your code in this part to load the expression matrix and the trajectory results from the specific pseudotime inference tool you used. \n",
    "\n",
    "<a id=1></a>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "filepath = f\"{input_path}/{project}.rds\"\n",
    "\n",
    "if os.path.exists(filepath):\n",
    "    pass \n",
    "else:\n",
    "    url = f\"https://zenodo.org/record/1443566/files/real/silver/{project}.rds?download=1\"\n",
    "    r = requests.get(url)  \n",
    "\n",
    "    with open(filepath, 'wb') as f:\n",
    "        f.write(r.content)\n",
    "        \n",
    "filepath = f\"{input_path}/{project}.rds\"\n",
    "\n",
    "from rpy2.robjects import pandas2ri\n",
    "pandas2ri.activate()\n",
    "\n",
    "readRDS = robjects.r['readRDS']\n",
    "df = readRDS(filepath)\n",
    "# df = pandas2ri.rpy2py_dataframe(df)\n",
    "data_keys = list(df.names)\n",
    "\n",
    "cell_ids = df[data_keys.index('cell_ids')]\n",
    "expression = df[data_keys.index('expression')]  \n",
    "genes = df[data_keys.index('feature_info')]['feature_id'].values\n",
    "\n",
    "N = len(cell_ids)  # number of cells \n",
    "G = len(genes)  # number of genes \n",
    "\n",
    "# load true traejctory and labels \n",
    "# true trajectory\n",
    "milestones_true = df[data_keys.index('milestone_ids')]\n",
    "network_true = df[data_keys.index('milestone_network')]\n",
    "M_true = len(milestones_true)\n",
    "\n",
    "# add node index; node index consistent with index in 'milestone_ids'\n",
    "# will use node index to present node from now on\n",
    "network_true['idx_from'] = [list(milestones_true).index(i) for i in network_true['from']]\n",
    "network_true['idx_to'] = [list(milestones_true).index(i) for i in network_true['to']]\n",
    "\n",
    "membership_true = df[data_keys.index('milestone_percentages')] \n",
    "# assign cells to the most probable node \n",
    "assignment_true = membership_true[membership_true.groupby(['cell_id'])['percentage'].transform(max) == membership_true['percentage']]\n",
    "assignment_true.set_index('cell_id', inplace=True)\n",
    "assignment_true = assignment_true.reindex(cell_ids)\n",
    "clusters_true = [list(milestones_true).index(c) for c in assignment_true['milestone_id'].values]\n",
    "\n",
    "# load trajectory inference results \n",
    "f = h5py.File(trajectory_filename, 'r')\n",
    "\n",
    "# # Check what keys are \n",
    "# for key in f.keys():\n",
    "#     print(key) \n",
    "\n",
    "key = 'data'\n",
    "\n",
    "# Get the HDF5 group\n",
    "group = f[key]\n",
    "\n",
    "# #Checkout what keys are inside that group.\n",
    "# for key in group.keys():\n",
    "#     print(key)\n",
    "\n",
    "_percentages = group['milestone_percentages']\n",
    "_network = group['milestone_network']\n",
    "_progressions = group['progressions']\n",
    "\n",
    "# # Check what keys are  \n",
    "#  data.keys()\n",
    "#  data['data'].keys()\n",
    "\n",
    "_cell_ids = list(_percentages['data']['cell_id'])\n",
    "_cell_ids = [i.decode('utf-8') for i in _cell_ids]\n",
    "estimated_percentages = pd.DataFrame(zip(_cell_ids, list(_percentages['data']['milestone_id']), list(_percentages['data']['percentage'])))\n",
    "estimated_percentages.columns = ['cell_id', 'milestone_id', 'percentage']\n",
    "\n",
    "_cell_ids = list(_progressions['data']['cell_id'])\n",
    "_cell_ids = [i.decode('utf-8') for i in _cell_ids]\n",
    "estimated_progressions = pd.DataFrame(zip(_cell_ids, list(_progressions['data']['from']), list(_progressions['data']['to']), list(_progressions['data']['percentage'])))\n",
    "estimated_progressions.columns = ['cell_id', 'from', 'to' ,'percentage']\n",
    "estimated_progressions = estimated_progressions.set_index(\"cell_id\")  \n",
    "estimated_progressions = estimated_progressions.reindex(assignment_true.index.values)  # assignment_true already reindexed by cell_ids\n",
    "\n",
    "estimated_network = pd.DataFrame(pd.DataFrame(zip(list(_network['data']['from']), list(_network['data']['to']), list(_network['data']['length']))))\n",
    "estimated_clusters = estimated_percentages.loc[estimated_percentages.groupby([\"cell_id\"])[\"percentage\"].idxmax()].set_index('cell_id').reindex(cell_ids)\n",
    "estimated_clusters['milestone_id'] = [_c.decode(\"utf-8\") for _c in estimated_clusters['milestone_id']]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 2. Prepare pseudotime results per cell for TraSig\n",
    "<a id=2></a>\n",
    "\n",
    "In this part, we prepare the following for TraSig:\n",
    "\n",
    "    1. assigned path (edge)\n",
    "    2. assigned time / progression on the edge \n",
    "    3. cell type labels (ground truth)\n",
    "\n",
    "* You may need to customize your code in this part to load the corresponding trajectory results from the specific pseudotime inference tool you used. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "estimated_progressions['from'] = [i.decode('utf-8') for i in estimated_progressions['from']] \n",
    "estimated_progressions['to'] = [i.decode('utf-8') for i in estimated_progressions['to']] \n",
    "estimated_progressions['edge'] = estimated_progressions['from'] + '_' + estimated_progressions['to'] \n",
    "\n",
    "# assign unique label (integer) to each edge \n",
    "\n",
    "edges = np.unique(estimated_progressions['edge'])\n",
    "\n",
    "edge2idx = {}\n",
    "for i, v in enumerate(edges):\n",
    "    edge2idx[v] = i\n",
    "    \n",
    "estimated_progressions['idx_edge'] = estimated_progressions['edge'].replace(edge2idx)\n",
    "\n",
    "hid_var = {'cell_path': estimated_progressions['idx_edge'].values,\n",
    "          'cell_time': estimated_progressions['percentage'].values,\n",
    "          'cell_labels':assignment_true['milestone_id'].values} "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Input file format\n",
    "This is dictionary with three keys, each correponding to a numpy array. \n",
    "\n",
    "    * \"cell_path\" (integer or strings): cluster / branch / edge label of each cell \n",
    "    * \"cell_time\" (float): progression on the branch assigned to each cell, a value in [0, 1]\n",
    "    * \"cell_labels\" (integer or strings): cell type label assigned to each cell \n",
    " \n",
    "See following for an example."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{'cell_path': array([0, 1, 1, ..., 1, 1, 1]),\n",
       " 'cell_time': array([0.        , 0.        , 0.        , ..., 0.87017241, 0.        ,\n",
       "        0.97737995]),\n",
       " 'cell_labels': array(['OPC', 'OPC', 'OPC', ..., 'NFOL', 'OPC', 'NFOL'], dtype=object)}"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "hid_var"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Filename format\n",
    "\n",
    "Remember to save this dictionary using the following file format."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "# save \n",
    "filename = f\"{project}{_preprocess}_{model_name}_it2_hid_var.pickle\"\n",
    "with open(os.path.join(output_path, filename), 'wb') as handle:\n",
    "    pickle.dump(hid_var, handle, protocol=pickle.HIGHEST_PROTOCOL)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 3. Subsetting expression data \n",
    "<a id=3></a>\n",
    "\n",
    "1. filter expression data to keep only ligands and receptors (and genes used for alignment, if applicable)\n",
    "2. the following take expression and ligand-receptor list (database) as input; alignment genes list is also required if using alignment and using genes other than ligands and receptors for alignment\n",
    "3. make sure your expression matrix is a numpy array (row: cells; column: genes)\n",
    "4. interaction list should be a pickle file wiht a python list of tuple: (ligand, receptor)\n",
    "    - see ../ligand_receptor_lists/ligand_receptor_FANTOM.pickle for an example list \n",
    "5. alignment genes list should be a pickle file wiht a python list of string: gene \n",
    "    - see ../alignment_genes_lists/cell_cycle_markers_Seurat.pickle for an example list "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "lr_list_path = \"../ligand_receptor_lists/ligand_receptor_FANTOM.pickle\"\n",
    "list_type = 'ligand_receptor'\n",
    "\n",
    "# use \"interaction\" if not using alignment or using 'interaction' for alignment\n",
    "gene_pair_type = \"cell_cycle\" \n",
    "# set as \"None\" if not using alignment or using 'interaction' for alignment\n",
    "align_gene_path = \"../alignment_genes_lists/cell_cycle_markers_Seurat.pickle\"  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "# get interaction file (list of (ligand, receptor))\n",
    "with open(os.path.join(lr_list_path), 'rb') as handle:\n",
    "    interaction_list = pickle.load(handle)\n",
    "    \n",
    "ligands_receptors = np.unique([i[0] for i in interaction_list] + [i[1] for i in interaction_list])\n",
    "\n",
    "# get align genes \n",
    "if gene_pair_type != \"None\" and align_gene_path != \"None\":\n",
    "    # get align genes (list of genes)\n",
    "    with open(os.path.join(align_gene_path), 'rb') as handle:\n",
    "        align_genes = pickle.load(handle)\n",
    "    \n",
    "    selected_genes = list(np.unique(list(ligands_receptors) + align_genes))\n",
    "else:\n",
    "    align_genes = None \n",
    "    selected_genes = ligands_receptors\n",
    "\n",
    "# get list of genes identified as ligand or receptor \n",
    "genes_upper = [g.upper() for g in genes]\n",
    "kepted_genes = list(set(genes_upper).intersection(set(selected_genes)))\n",
    "\n",
    "df = pd.DataFrame(expression)\n",
    "df.columns = genes_upper\n",
    "df.index = cell_ids\n",
    "\n",
    "df_sub = df[kepted_genes]\n",
    "\n",
    "# save filtered expression \n",
    "if align_genes == None:\n",
    "    filename = f\"{project}{_preprocess}_{list_type}.txt\"\n",
    "else:\n",
    "    filename = f\"{project}{_preprocess}_{list_type}_{gene_pair_type}.txt\"\n",
    "    \n",
    "data_file = os.path.join(output_path, filename)\n",
    "df_sub.to_csv(data_file)\n",
    "\n",
    "# save filtered interactions (list of (ligand, receptor) that are expressed)\n",
    "filtered_interactions = []\n",
    "for i, j in interaction_list:\n",
    "    if i in kepted_genes and j in kepted_genes:\n",
    "        filtered_interactions.append((i, j))\n",
    " \n",
    "filename = f\"{list_type}_{project}{_preprocess}.pickle\"\n",
    "with open(os.path.join(output_path, filename), 'wb') as handle:\n",
    "    pickle.dump(filtered_interactions, handle, protocol=pickle.HIGHEST_PROTOCOL)\n",
    "    \n",
    "# save filtered align genes \n",
    "if align_genes != None:\n",
    "    align_genes_filtered = list(set(kepted_genes).intersection(set(align_genes)))\n",
    "    if len(align_genes_filtered) == 0:\n",
    "        raise ValueError(f\"No align genes expressed in the dataset! Consider using another set of genes for alignment!\")\n",
    "    else:\n",
    "        filename = f\"align_{gene_pair_type}_filtered_{project}{_preprocess}.pickle\"\n",
    "        with open(os.path.join(output_path, filename), 'wb') as handle:\n",
    "            pickle.dump(align_genes_filtered, handle, protocol=pickle.HIGHEST_PROTOCOL)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 4. Save correspondence from sampling time to paths\n",
    "<a id=4></a>\n",
    "\n",
    "1. Note here cell_path refers to the edge where the cell is assigned to \n",
    "2. We will only find interactions between cells from the same sampling time and those from consecutive sampling times:\n",
    "    - i.e., between the ones from the same time, the ones from 1 sampling time before the ones from 1 sampling time after\n",
    "3. Given we don't know the sampling time for the example data, we set all sampling time as 0. For your own data, if you are not certain about sampling time, just assign the time for all cells as 0.\n",
    "4. If sampling time is known, rank the real time (e.g. day 0, day 17) first and assign the rank to the cell_ori_time variable below. \n",
    "    - e.g., for cells from two sampling time day 0 and day 17, assign those from day 0 as 0 and those from day 17 as 1. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy import stats"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### If known sampling time, then set the following variable = the sampling time of cells"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "cell_ori_time = np.repeat(0, N)  # put all cells at time 0 if sampling time unknow "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### The following is trying to assign each cluster / branch / edge a sampling time, determined by the majority of cells"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "unique_days = np.unique(cell_ori_time)\n",
    "sorted_days = list(np.sort(unique_days)) \n",
    "cell_paths = np.unique(hid_var[\"cell_path\"])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "sampleT2path = dict.fromkeys(range(len(sorted_days)))  # use index of sorted sampling time as key\n",
    "for k, v in sampleT2path.items():\n",
    "    sampleT2path[k] = []"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "current path (edge) 0\n",
      "Sampling time for the majority of cells: 0, making 1.0% percent\n",
      "current path (edge) 1\n",
      "Sampling time for the majority of cells: 0, making 1.0% percent\n",
      "current path (edge) 2\n",
      "Sampling time for the majority of cells: 0, making 1.0% percent\n",
      "current path (edge) 3\n",
      "Sampling time for the majority of cells: 0, making 1.0% percent\n"
     ]
    }
   ],
   "source": [
    "for i, cur_path in enumerate(cell_paths):\n",
    "    print(\"current path (edge)\", cur_path)\n",
    "    \n",
    "    # get data corresponding to a path\n",
    "    condition = hid_var[\"cell_path\"] == cur_path\n",
    "    cur_days = np.array(cell_ori_time)[condition]\n",
    "    \n",
    "    # get the sampling time for the majority cells \n",
    "    mode, count = stats.mode(cur_days)\n",
    "    print(f\"Sampling time for the majority of cells: {mode[0]}, making {round(float(count[0])/len(cur_days), 2)}% percent\")\n",
    "    cur_sampleT = mode[0]\n",
    "    \n",
    "    # will use index instead of input time \n",
    "    sampleT2path[sorted_days.index(cur_sampleT)].append(cur_path)\n",
    "    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Input file format\n",
    "This is dictionary with sampling time as keyes, each correponding to a list of clusters that are deterimined to be sampled from the sampling time. \n",
    " \n",
    "See following for an example."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{0: [0, 1, 2, 3]}"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sampleT2path"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "# save the dictionary\n",
    "filename = 'sampling_time_per_path_' + project + suffix + '.pickle'\n",
    "\n",
    "with open(os.path.join(output_path, filename), 'wb') as handle:\n",
    "    pickle.dump(sampleT2path, handle, protocol=pickle.HIGHEST_PROTOCOL)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python [conda env:spatem]",
   "language": "python",
   "name": "conda-env-spatem-py"
  },
  "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.6.11"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}