run_OnClass_example.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "from anndata import read_h5ad\n",
    "from scipy import stats, sparse\n",
    "import numpy as np\n",
    "import sys\n",
    "from collections import Counter\n",
    "from OnClass.OnClassModel import OnClassModel\n",
    "from utils import read_ontology_file, read_data, run_scanorama_multiply_datasets\n",
    "from config import ontology_data_dir, scrna_data_dir, model_dir, Run_scanorama_batch_correction, NHIDDEN, MAX_ITER"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## read data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "#Change train_file, test_file, train_label, test_label according to your datasets\n",
    "#train_label is the key of labels in .h5ad train_file.\n",
    "train_file = scrna_data_dir + '/Lemur/microcebusBernard.h5ad'\n",
    "test_file = scrna_data_dir + '/Lemur/microcebusAntoine.h5ad'\n",
    "\n",
    "train_label = 'cell_ontology_id'\n",
    "test_label = 'cell_ontology_id'\n",
    "model_path = model_dir + 'example_file_model'\n",
    "\n",
    "print ('read ontology data and initialize training model...')\n",
    "cell_type_nlp_emb_file, cell_type_network_file, cl_obo_file = read_ontology_file('cell ontology', ontology_data_dir)\n",
    "OnClass_train_obj = OnClassModel(cell_type_nlp_emb_file = cell_type_nlp_emb_file, cell_type_network_file = cell_type_network_file)\n",
    "\n",
    "print ('read training single cell data...')\n",
    "train_feature, train_genes, train_label, _, _ = read_data(train_file, cell_ontology_ids = OnClass_train_obj.cell_ontology_ids,\n",
    "\texclude_non_leaf_ontology = False, tissue_key = 'tissue', AnnData_label_key = train_label, filter_key = {},\n",
    "\tnlp_mapping = False, cl_obo_file = cl_obo_file, cell_ontology_file = cell_type_network_file, co2emb = OnClass_train_obj.co2vec_nlp)\n",
    "#you can also replace it with your own data and make sure that:\n",
    "#train_feature is a ncell by ngene matrix\n",
    "#train_genes is a ngene long vector of gene names\n",
    "#train_label is a ncell long vector\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Embed cell types based on the Cell Ontology Graph"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print ('embed cell types using the cell ontology...')\n",
    "OnClass_train_obj.EmbedCellTypes(train_label)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Training"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# print ('read test single cell data...')\n",
    "x = read_h5ad(test_file)\n",
    "test_label = x.obs[test_label].tolist()\n",
    "test_feature = x.X.toarray()\n",
    "test_genes = np.array([x.upper() for x in x.var.index])\n",
    "\n",
    "# optional batch correction\n",
    "if Run_scanorama_batch_correction:\n",
    "\ttrain_feature, test_feature = run_scanorama_multiply_datasets([train_feature, test_feature], [train_genes, test_genes], scan_dim = 10)[1]\n",
    "\tprint (np.shape(train_feature), np.shape(test_feature))\n",
    "\n",
    "print ('generate pretrain model. Save the model to $model_path...')\n",
    "cor_train_feature, cor_test_feature, cor_train_genes, cor_test_genes = OnClass_train_obj.ProcessTrainFeature(train_feature, train_label, train_genes, test_feature = test_feature, test_genes = test_genes)\n",
    "OnClass_train_obj.BuildModel(ngene = len(cor_train_genes), nhidden = NHIDDEN)\n",
    "OnClass_train_obj.Train(cor_train_feature, train_label, save_model = model_path, max_iter = MAX_ITER)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Classify test cells"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print ('initialize test model. Load the model from $model_path...')\n",
    "OnClass_test_obj = OnClassModel(cell_type_nlp_emb_file = cell_type_nlp_emb_file, cell_type_network_file = cell_type_network_file)\n",
    "cor_test_feature = OnClass_train_obj.ProcessTestFeature(cor_test_feature, cor_test_genes, use_pretrain = model_path, log_transform = False)\n",
    "OnClass_test_obj.BuildModel(ngene = None, use_pretrain = model_path)\n",
    "\n",
    "#use_normalize=False will return a tree-based prediction, where parent node often has higher score than child node. use_normalize=True will normalize among child nodes and parent nodes\n",
    "pred_Y_seen, pred_Y_all, pred_label = OnClass_test_obj.Predict(cor_test_feature, test_genes = cor_test_genes, use_normalize=True)\n",
    "pred_label_str = [OnClass_test_obj.i2co[l] for l in pred_label]\n",
    "#x.obs['OnClass_annotation_flat_based_ontology_ID'] = pred_label_str\n",
    "\n",
    "pred_Y_seen, pred_Y_all, pred_label = OnClass_test_obj.Predict(cor_test_feature, test_genes = cor_test_genes, use_normalize=False)\n",
    "pred_label_str = [OnClass_test_obj.i2co[l] for l in pred_label]\n",
    "#x.obs['OnClass_annotation_tree_based_ontology_ID'] = pred_label_str\n",
    "#x.write(scrna_data_dir + 'Pilot12.annotated.h5ad')"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python (tabula-muris-env)",
   "language": "python",
   "name": "tabula-muris-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.6.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}