notebooks/00_data_proc.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "758da1a7-f1fa-4001-8f35-9a3ad32560af",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "import anndata as ad\n",
    "import scanpy as sc\n",
    "from scipy.sparse import csr_matrix"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "406f7632",
   "metadata": {},
   "source": [
    "**Count matrix and metadata for VISp dataset**\n",
    " - Download the count data from Allen Institute portal\n",
    " - Convert to AnnData format - see [this getting started with AnnData tutorial](https://anndata-tutorials.readthedocs.io/en/latest/getting-started.html)\n",
    " - Save the resulting object for later use"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "a3a398cd",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "  % Total    % Received % Xferd  Average Speed   Time    Time     Time  Current\n",
      "                                 Dload  Upload   Total   Spent    Left  Speed\n",
      "100  291M  100  291M    0     0  54.3M      0  0:00:05  0:00:05 --:--:-- 55.5M     0  53.6M      0  0:00:05  0:00:03  0:00:02 53.6M\n",
      "Archive:  ../data/VISp.zip\n",
      "  inflating: ../data/VISp/mouse_VISp_2018-06-14_exon-matrix.csv  \n",
      "  inflating: ../data/VISp/mouse_VISp_2018-06-14_genes-rows.csv  \n",
      "  inflating: ../data/VISp/mouse_VISp_2018-06-14_intron-matrix.csv  \n",
      "  inflating: ../data/VISp/mouse_VISp_2018-06-14_readme.txt  \n",
      "  inflating: ../data/VISp/mouse_VISp_2018-06-14_samples-columns.csv  \n"
     ]
    }
   ],
   "source": [
    "# Download count matrices from https://portal.brain-map.org/atlases-and-data/rnaseq/mouse-v1-and-alm-smart-seq or use the shell commands below.\n",
    "!curl -o ../data/VISp.zip https://celltypes.brain-map.org/api/v2/well_known_file_download/694413985\n",
    "!unzip -d ../data/VISp ../data/VISp.zip\n",
    "!rm ../data/VISp.zip"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "c577e965",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>seq_name</th>\n",
       "      <th>class</th>\n",
       "      <th>subclass</th>\n",
       "      <th>cluster</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>sample_id</th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>F1S4_160108_001_A01</th>\n",
       "      <td>LS-15006_S09_E1-50</td>\n",
       "      <td>GABAergic</td>\n",
       "      <td>Vip</td>\n",
       "      <td>Vip Arhgap36 Hmcn1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>F1S4_160108_001_B01</th>\n",
       "      <td>LS-15006_S10_E1-50</td>\n",
       "      <td>GABAergic</td>\n",
       "      <td>Lamp5</td>\n",
       "      <td>Lamp5 Lsp1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>F1S4_160108_001_C01</th>\n",
       "      <td>LS-15006_S11_E1-50</td>\n",
       "      <td>GABAergic</td>\n",
       "      <td>Lamp5</td>\n",
       "      <td>Lamp5 Lsp1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>F1S4_160108_001_D01</th>\n",
       "      <td>LS-15006_S12_E1-50</td>\n",
       "      <td>GABAergic</td>\n",
       "      <td>Vip</td>\n",
       "      <td>Vip Crispld2 Htr2c</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>F1S4_160108_001_E01</th>\n",
       "      <td>LS-15006_S13_E1-50</td>\n",
       "      <td>GABAergic</td>\n",
       "      <td>Lamp5</td>\n",
       "      <td>Lamp5 Plch2 Dock5</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                               seq_name      class subclass  \\\n",
       "sample_id                                                     \n",
       "F1S4_160108_001_A01  LS-15006_S09_E1-50  GABAergic      Vip   \n",
       "F1S4_160108_001_B01  LS-15006_S10_E1-50  GABAergic    Lamp5   \n",
       "F1S4_160108_001_C01  LS-15006_S11_E1-50  GABAergic    Lamp5   \n",
       "F1S4_160108_001_D01  LS-15006_S12_E1-50  GABAergic      Vip   \n",
       "F1S4_160108_001_E01  LS-15006_S13_E1-50  GABAergic    Lamp5   \n",
       "\n",
       "                                cluster  \n",
       "sample_id                                \n",
       "F1S4_160108_001_A01  Vip Arhgap36 Hmcn1  \n",
       "F1S4_160108_001_B01          Lamp5 Lsp1  \n",
       "F1S4_160108_001_C01          Lamp5 Lsp1  \n",
       "F1S4_160108_001_D01  Vip Crispld2 Htr2c  \n",
       "F1S4_160108_001_E01   Lamp5 Plch2 Dock5  "
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Load VISp dataset\n",
    "filename = '../data/VISp/mouse_VISp_2018-06-14_exon-matrix.csv'\n",
    "expr_df = pd.read_csv(filename, header=0, index_col=0, delimiter=',').transpose()\n",
    "expr = expr_df.values\n",
    "\n",
    "# Find gene names\n",
    "filename = '../data/VISp/mouse_VISp_2018-06-14_genes-rows.csv'\n",
    "genes_df = pd.read_csv(filename, header=0, index_col=0, delimiter=',')\n",
    "gene_symbol = genes_df.index.values\n",
    "gene_ids = genes_df['gene_entrez_id'].values\n",
    "gene_names = np.array([gene_symbol[np.where(gene_ids == name)[0][0]] for name in expr_df.columns])\n",
    "\n",
    "# Get metadata and save restrict to relevant fields\n",
    "filename = '../data/VISp/mouse_VISp_2018-06-14_samples-columns.csv'\n",
    "obs = pd.read_csv(filename, header=0, index_col=0, delimiter=',', encoding='iso-8859-1')\n",
    "\n",
    "obs = obs.reset_index()\n",
    "obs = obs[['sample_name','seq_name','class','subclass','cluster']]\n",
    "obs = obs.rename(columns={'sample_name':'sample_id'})\n",
    "obs = obs.set_index('sample_id')\n",
    "obs.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "2443ffbd",
   "metadata": {},
   "outputs": [],
   "source": [
    "# compose and store anndata object for efficient read/write\n",
    "adata = ad.AnnData(X=csr_matrix(expr), dtype=expr.dtype)\n",
    "adata.var_names = gene_names\n",
    "adata.var.index.set_names('genes', inplace=True)\n",
    "adata.obs = obs\n",
    "adata.write('../data/VISp.h5ad')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c5a11dd0",
   "metadata": {},
   "source": [
    "**Filtering samples**\n",
    "\n",
    "The next code block is optional, and requires `VISp_PERSIST_metadata.csv` which contains:\n",
    "- cell type labels at different resolutions of the taxonomy (see manuscript for details)\n",
    "- sample ids to filter out non-neuronal cells\n",
    "\n",
    "In the following, we will\n",
    "1. restrict cells only to those samples specified in `VISp_PERSIST_metadata.csv`\n",
    "2. append metadata from `VISp_PERSIST_metadata.csv` to the AnnData object\n",
    "3. normalize counts, determine highly variable genes using scanpy functions\n",
    "3. save a filtered AnnData object into a .h5ad file for subsequent use"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "10b09db4",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "             seq_name       cell_types_98 cell_types_50 cell_types_25\n",
      "0  LS-15006_S09_E1-50  Vip Arhgap36 Hmcn1           n70           n66\n",
      "1  LS-15006_S10_E1-50          Lamp5 Lsp1    Lamp5 Lsp1           n78\n",
      "\n",
      "old shape: (15413, 45768)\n",
      "new shape: (13349, 45768)\n"
     ]
    }
   ],
   "source": [
    "adata = ad.read_h5ad('../data/VISp.h5ad')\n",
    "persist_df = pd.read_csv('../data/VISp_PERSIST_metadata.csv')\n",
    "print(persist_df.head(2))\n",
    "print(f'\\nold shape: {adata.shape}')\n",
    "\n",
    "adata = adata[adata.obs['seq_name'].isin(persist_df['seq_name']), :]\n",
    "obs = adata.obs.copy().reset_index()\n",
    "obs = obs.merge(right=persist_df, how='left', left_on='seq_name', right_on='seq_name')\n",
    "obs = obs.set_index('sample_id')\n",
    "\n",
    "adata = ad.AnnData(X=adata.X, dtype=adata.X.dtype, obs=obs, var=adata.var)\n",
    "print(f'new shape: {adata.shape}')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3b1398f1",
   "metadata": {},
   "source": [
    "**Normalization and preliminary gene selection**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "11401b0c",
   "metadata": {},
   "outputs": [],
   "source": [
    "# transforms data in adata.X\n",
    "adata.layers['log1pcpm'] = sc.pp.normalize_total(adata, target_sum=1e6, inplace=False)['X']\n",
    "\n",
    "# transforms data in layers['lognorm'] inplace\n",
    "sc.pp.log1p(adata, layer='log1pcpm')\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "2b8cd61f",
   "metadata": {},
   "outputs": [],
   "source": [
    "import toml\n",
    "dat = toml.load('../data/VISp_markers.toml')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "c42e06e6",
   "metadata": {},
   "outputs": [],
   "source": [
    "# introduces \"highly_variable\" column to adata.var\n",
    "sc.pp.highly_variable_genes(adata, layer='log1pcpm', flavor='cell_ranger',\n",
    "                            n_top_genes=10000, inplace=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "a5cdf670",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "AnnData object with n_obs × n_vars = 13349 × 45768\n",
       "    obs: 'seq_name', 'class', 'subclass', 'cluster', 'cell_types_98', 'cell_types_50', 'cell_types_25'\n",
       "    var: 'highly_variable', 'means', 'dispersions', 'dispersions_norm', 'markers'\n",
       "    uns: 'log1p', 'hvg'\n",
       "    layers: 'log1pcpm'"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Create new field with marker genes\n",
    "adata.var['markers'] = np.isin(adata.var.index.values,dat['markers'])\n",
    "adata"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "bb0bb3c2",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<13349x45768 sparse matrix of type '<class 'numpy.float32'>'\n",
       "\twith 120801360 stored elements in Compressed Sparse Row format>"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# This is a sparse matrix\n",
    "adata.layers['log1pcpm']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "dc70f25f",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[0.        , 0.        , 3.8425822 , 4.405392  , 0.        ],\n",
       "       [0.        , 0.        , 4.164545  , 4.5287347 , 0.42111826],\n",
       "       [0.        , 0.        , 3.8250968 , 3.562683  , 0.        ],\n",
       "       [0.        , 0.        , 3.935433  , 0.        , 0.        ],\n",
       "       [0.        , 0.        , 5.4067717 , 4.6221585 , 0.        ]],\n",
       "      dtype=float32)"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# For sparse matrix `M`, one can use `M.A` to convert it to dense\n",
    "adata.layers['log1pcpm'].A[:5,:5]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f8b21c66",
   "metadata": {},
   "source": [
    "**Parting notes:**\n",
    "\n",
    "The anndata object created in this way has a few different fields that we will end up using with PERSIST\n",
    "1. The raw counts are in `adata.X`\n",
    "2. The normalized counts (log1p of CPM values) are in `adata.layers['log1pcpm']`\n",
    "3. All metadata (cell type labels etc. for supervised mode in PERSIST) is in `adata.obs`\n",
    "4. A coarse selection of genes is in `adata.var['highly_variable']`\n",
    "5. Marker genes defined by Tasic et al. are indicated by `adata.var['markers']`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "53129c25",
   "metadata": {},
   "outputs": [],
   "source": [
    "# adata_hvg is a view. We'll convert it to a new AnnData object and write it out. \n",
    "adata_hvg = ad.AnnData(X=adata.X, dtype=adata.X.dtype,\n",
    "                       obs=adata.obs, var=adata.var[['highly_variable']],\n",
    "                       layers=adata.layers, uns=adata.uns)\n",
    "adata_hvg.write('../data/VISp_filtered_cells.h5ad')"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3.8.15 ('persist')",
   "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.8.15"
  },
  "vscode": {
   "interpreter": {
    "hash": "660691d0bb1e24e4a68343475da76b12317d18d7509dab1fd0158534dd4eebe4"
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}