01_preprocess_PMBC_and_negative_controls.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pickle\n",
    "import pandas as pd\n",
    "import numpy as np\n",
    "import anndata\n",
    "from scipy.io import mmread\n",
    "from scipy import sparse\n",
    "\n",
    "from tools import remove_rare_genes, rnaseqTools"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 33k PMBC data \n",
    "### (10X v1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "def load_PMBC_dataset(folder,genename_file='genes.tsv'):\n",
    "    counts = mmread(folder + 'matrix.mtx')\n",
    "    counts = sparse.csr_matrix(counts).T\n",
    "    genes = pd.read_csv(folder + genename_file, header=None, sep='\\t')\n",
    "    genes = genes.values[:,1].squeeze().astype(str)\n",
    "    return counts,genes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Of 20678 total genes, returning 16809 genes that are detected in 5 or more cells.\n",
      "Output shape: (33148, 16809)\n",
      "{'counts': <33148x16809 sparse matrix of type '<class 'numpy.float64'>'\n",
      "\twith 25374746 stored elements in Compressed Sparse Row format>, 'genes': array(['AL627309.1', 'AP006222.2', 'RP11-206L10.2', ..., 'KIR2DL2',\n",
      "       'PNRC2', 'SRSF10'], dtype='<U19'), 'label': '33k PBMC (v1)', 'folder': 'datasets/33k_pbmc/'}\n",
      "CPU times: user 45.2 s, sys: 847 ms, total: 46.1 s\n",
      "Wall time: 46.3 s\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "folder = 'datasets/33k_pbmc/'\n",
    "label = '33k PBMC (v1)'\n",
    "counts,genes=load_PMBC_dataset(folder)\n",
    "counts,genes=remove_rare_genes(counts,genes,minimum_detected_cells_per_gene=5)\n",
    "\n",
    "counts = sparse.csr_matrix(counts)\n",
    "\n",
    "dataset = dict(counts=counts,genes=genes,label=label,folder=folder)\n",
    "print(dataset)\n",
    "pickle.dump(dataset, open(folder+'preprocessed.pickle', \"wb\"), protocol=4)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## Write dense data to csv for R to read\n",
    "counts = dataset['counts'].toarray()\n",
    "df = pd.DataFrame(counts,columns=dataset['genes'])\n",
    "ad = anndata.AnnData(df)\n",
    "ad.var_names_make_unique()\n",
    "df_unique_genes_T = ad.to_df().T #transpose as R/sctransform expects rows as genes\n",
    "df_unique_genes_T.to_csv('datasets/33k_pbmc/preproc_for_R.csv')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Negative control data\n",
    "### 10X v2 - Svensson 2017 (technical control)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "code_folding": []
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Of 20647 total genes, returning 13025 genes that are detected in 5 or more cells.\n",
      "Output shape: (2000, 13025)\n",
      "{'counts': <2000x13025 sparse matrix of type '<class 'numpy.float32'>'\n",
      "\twith 1239772 stored elements in Compressed Sparse Row format>, 'genes': array(['ENSG00000000003', 'ENSG00000000419', 'ENSG00000000457', ...,\n",
      "       'LRG_759', 'LRG_788', 'LRG_92'], dtype=object), 'label': '2k Chromium control Svensson 2017 (10x v2) - sample 1', 'folder': 'datasets/10x/', 'protocol': '10X'}\n",
      "Of 21411 total genes, returning 13239 genes that are detected in 5 or more cells.\n",
      "Output shape: (2000, 13239)\n",
      "{'counts': <2000x13239 sparse matrix of type '<class 'numpy.float32'>'\n",
      "\twith 1135300 stored elements in Compressed Sparse Row format>, 'genes': array(['ENSG00000000003', 'ENSG00000000419', 'ENSG00000000971', ...,\n",
      "       'LRG_81', 'LRG_89', 'LRG_92'], dtype=object), 'label': '2k Chromium control Svensson 2017 (10x v2) - sample 2', 'folder': 'datasets/10x/', 'protocol': '10X'}\n",
      "CPU times: user 1.63 s, sys: 35.7 ms, total: 1.67 s\n",
      "Wall time: 1.67 s\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "folder = 'datasets/10x/'\n",
    "file = folder + 'svensson_chromium_control.h5ad'\n",
    "\n",
    "raw_input = anndata.read(file)\n",
    "\n",
    "for i,sample_id in enumerate([20311, 20312]):\n",
    "    sample_input = raw_input[raw_input.obs.query('sample == \"%u\"' % (sample_id)).index]\n",
    "    label = '2k Chromium control Svensson 2017 (10x v2) - sample %u' % (i+1)\n",
    "    \n",
    "    counts_pd = sample_input.to_df()\n",
    "    genes = np.array(counts_pd.columns)\n",
    "    counts = counts_pd.to_numpy()\n",
    "    counts,genes=remove_rare_genes(counts,genes,minimum_detected_cells_per_gene=5)\n",
    "\n",
    "    counts = sparse.csr_matrix(counts)\n",
    "\n",
    "    dataset = dict(counts=counts,genes=genes,label=label,folder=folder,protocol='10X')\n",
    "    print(dataset)\n",
    "    pickle.dump(dataset, open(folder+'sample%u_preprocessed.pickle' % (i+1), \"wb\"), protocol=4)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# inDrop"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Klein 2015 (one technial and various biological controls)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "paths =['datasets/indrop/GSM1599501_K562_pure_RNA.csv']\n",
    "ns = ['1k']\n",
    "labelsshort = ['K562_pure_RNA']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "def pickle_klein_data(path,labelshort,n):\n",
    "    label = '%s %s Klein 2015 (inDrop)' % (n,labelshort)\n",
    "    df = pd.read_csv(path,index_col=0)\n",
    "    genes = df.index\n",
    "    counts = df.values.T\n",
    "    counts,genes = remove_rare_genes(counts,genes,minimum_detected_cells_per_gene=5)\n",
    "    counts = sparse.csr_matrix(counts)\n",
    "    dataset = dict(counts=counts,genes=genes,folder=path,label=label,protocol='inDrop')\n",
    "    picklefile = '%s_preprocessed.pickle' % (path)\n",
    "    print(dataset)\n",
    "    pickle.dump(dataset, open(picklefile, \"wb\"), protocol=4)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Of 25266 total genes, returning 25025 genes that are detected in 5 or more cells.\n",
      "Output shape: (953, 25025)\n",
      "{'counts': <953x25025 sparse matrix of type '<class 'numpy.longlong'>'\n",
      "\twith 9261429 stored elements in Compressed Sparse Row format>, 'genes': Index(['A1BG', 'A1BG-AS1', 'A1CF', 'A2M', 'A2M-AS1', 'A2ML1', 'A2MP1',\n",
      "       'A3GALT2', 'A4GALT', 'A4GNT',\n",
      "       ...\n",
      "       'ZWILCH', 'ZWINT', 'ZXDA', 'ZXDB', 'ZXDC', 'ZYG11A', 'ZYG11B', 'ZYX',\n",
      "       'ZZEF1', 'ZZZ3'],\n",
      "      dtype='object', length=25025), 'folder': 'datasets/indrop/GSM1599501_K562_pure_RNA.csv', 'label': '1k K562_pure_RNA Klein 2015 (inDrop)', 'protocol': 'inDrop'}\n",
      "CPU times: user 2.98 s, sys: 200 ms, total: 3.18 s\n",
      "Wall time: 3.38 s\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "## pickle-prepares all klein datasets\n",
    "for p,n,l in zip(paths,ns,labelsshort):\n",
    "    pickle_klein_data(p,l,n)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# MicrowellSeq"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Han 2018 (biological controls)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "..................... done\n",
      "GSM2906413_EmbryonicStemCell_dge.txt.gz (9994, 20682)\n",
      "Of 9994 cells, returning 9994 cells that have a depth larger than 500\n",
      "New shape: (9994, 20682)\n",
      "Of 20677 total genes, returning 15069 genes that are detected in 5 or more cells.\n",
      "Output shape: (9994, 15069)\n",
      "{'counts': <9994x15069 sparse matrix of type '<class 'numpy.float64'>'\n",
      "\twith 6762689 stored elements in Compressed Sparse Row format>, 'genes': array(['0610007P14Rik', '0610009B22Rik', '0610009L18Rik', ..., 'n-R5s151',\n",
      "       'n-R5s2', 'n-R5s88'], dtype='<U14'), 'cells': array(['EmbryonicStemCells_1.CTCGCAATCAACAACCTA',\n",
      "       'EmbryonicStemCells_1.CCGCTAATTCCAGTTGCC',\n",
      "       'EmbryonicStemCells_1.CTCGCACACAAGACACCC', ...,\n",
      "       'EmbryonicStemCells_1.CGAGTACCAGACAGATGG',\n",
      "       'EmbryonicStemCells_1.CGGCAGCTGAAAGGCTGC',\n",
      "       'EmbryonicStemCells_1.CTGTGTGAGATCGTCCCG'], dtype=object), 'folder': 'datasets/microwellseq/GSM2906413_EmbryonicStemCell_dge.txt.gz', 'label': 'Han 2018 EmbryonicStemCell control (MicrowellSeq)', 'protocol': 'microwellSeq'}\n",
      "CPU times: user 1min 7s, sys: 2.36 s, total: 1min 9s\n",
      "Wall time: 1min 9s\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "## Preproc files for ES cells\n",
    "files_raw = ['datasets/microwellseq/GSM2906413_EmbryonicStemCell_dge.txt.gz']\n",
    "separators = [' ']\n",
    "\n",
    "### iterate over files and make pickles for each\n",
    "for i,(file_raw,sep) in enumerate(zip(files_raw,separators)):\n",
    "\n",
    "    file = file_raw.split('/')[-1]\n",
    "    cellline = file.split('_')[1]\n",
    "    \n",
    "    counts, genes, cells = rnaseqTools.sparseload(file_raw, sep=sep)\n",
    "\n",
    "    counts = counts.toarray()\n",
    "    print(file,counts.shape)\n",
    "    \n",
    "    #remove low depth cells\n",
    "    depths = np.array(np.sum(counts,axis=1)).flatten()\n",
    "    minimum_depth = 500\n",
    "    cells = cells[depths>minimum_depth]\n",
    "    counts = counts[depths>minimum_depth,:]\n",
    "    print('Of',len(depths),'cells, returning',sum(depths>minimum_depth),'cells that have a depth larger than', minimum_depth)\n",
    "    print('New shape:', counts.shape)    \n",
    "    \n",
    "    counts,genes = remove_rare_genes(counts,genes,minimum_detected_cells_per_gene=5)\n",
    "    counts = sparse.csr_matrix(counts)\n",
    "\n",
    "    dataset = dict(counts=counts,genes=genes,cells=cells,folder=file_raw,\n",
    "                   label='Han 2018 %s control (MicrowellSeq)' % (cellline),\n",
    "                   protocol='microwellSeq')\n",
    "    \n",
    "    print(dataset)\n",
    "    picklefile = '%s_preprocessed.pickle' % (file_raw)\n",
    "    pickle.dump(dataset, open(picklefile, \"wb\"), protocol=4)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The watermark extension is already loaded. To reload it, use:\n",
      "  %reload_ext watermark\n"
     ]
    }
   ],
   "source": [
    "%load_ext watermark"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Last updated: 2021-05-21T18:31:57.975558+02:00\n",
      "\n",
      "Python implementation: CPython\n",
      "Python version       : 3.8.0\n",
      "IPython version      : 7.21.0\n",
      "\n",
      "Compiler    : GCC 8.3.0\n",
      "OS          : Linux\n",
      "Release     : 3.10.0-957.el7.x86_64\n",
      "Machine     : x86_64\n",
      "Processor   : x86_64\n",
      "CPU cores   : 40\n",
      "Architecture: 64bit\n",
      "\n"
     ]
    }
   ],
   "source": [
    "watermark"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "pandas : 1.2.0\n",
      "numpy  : 1.20.1\n",
      "anndata: 0.7.5\n",
      "scipy  : 1.6.0\n",
      "\n"
     ]
    }
   ],
   "source": [
    "watermark --iversions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "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.8.10"
  },
  "varInspector": {
   "cols": {
    "lenName": 16,
    "lenType": 16,
    "lenVar": 40
   },
   "kernels_config": {
    "python": {
     "delete_cmd_postfix": "",
     "delete_cmd_prefix": "del ",
     "library": "var_list.py",
     "varRefreshCmd": "print(var_dic_list())"
    },
    "r": {
     "delete_cmd_postfix": ") ",
     "delete_cmd_prefix": "rm(",
     "library": "var_list.r",
     "varRefreshCmd": "cat(var_dic_list()) "
    }
   },
   "types_to_exclude": [
    "module",
    "function",
    "builtin_function_or_method",
    "instance",
    "_Feature"
   ],
   "window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}