09_simulate_tasic_like_data.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "69589703-921d-45f7-87af-cceaaf41ad23",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import anndata\n",
    "from scipy import sparse"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "3a1fac73-e7e5-4623-a8ec-bd64af4af845",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/usr/local/lib/python3.8/dist-packages/anndata/_core/anndata.py:121: ImplicitModificationWarning: Transforming to str index.\n",
      "  warnings.warn(\"Transforming to str index.\", ImplicitModificationWarning)\n",
      "/usr/local/lib/python3.8/dist-packages/anndata/_core/anndata.py:1830: UserWarning: Variable names are not unique. To make them unique, call `.var_names_make_unique`.\n",
      "  utils.warn_names_duplicates(\"var\")\n"
     ]
    }
   ],
   "source": [
    "basepath = 'data/tasic/'\n",
    "adata_tasic = anndata.read_h5ad(f'{basepath}adata.h5ad')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "3aa328a0-754d-4a5a-850b-c22f58ed6be4",
   "metadata": {},
   "outputs": [],
   "source": [
    "def simulate(adata,only_marker_genes_DE,read_umi_factor = 100, theta_umi = 100, theta_z = 1, seed = 42):\n",
    "\n",
    "    cluster_ids = np.unique(adata.obs['clusters'])\n",
    "\n",
    "    readcounts_sim = []\n",
    "    umi_counts_sim = []\n",
    "    clusterlabels_sim = []\n",
    "    clustercolors_sim = []\n",
    "\n",
    "    housekeeping_idx = ~adata.var['marker_idx']\n",
    "    across_cluster_ps = np.array(adata.X.sum(axis=0) / adata.X.sum()).flatten()\n",
    "\n",
    "    for cluster_i in cluster_ids:\n",
    "\n",
    "        print(f'simulating cluster {cluster_i}..')\n",
    "\n",
    "        cluster_idx = adata.obs['clusters'] == cluster_i\n",
    "        ad_cl = adata[cluster_idx]\n",
    "\n",
    "        ###get UMIs\n",
    "        #get seqdepths per cell\n",
    "        readcount_ns = np.array(ad_cl.X.sum(axis=1)).flatten()\n",
    "        umi_ns = readcount_ns / read_umi_factor\n",
    "        #get proportions per gene\n",
    "        ps = np.array(ad_cl.X.sum(axis=0) / np.sum(readcount_ns)).flatten()\n",
    "        \n",
    "        if only_marker_genes_DE:\n",
    "            #set ps to across-cluster ps for all non-marker genes\n",
    "            ps[housekeeping_idx] = across_cluster_ps[housekeeping_idx]\n",
    "\n",
    "        #sample UMIs from NB(n*p)\n",
    "        mus = np.outer(umi_ns,ps)\n",
    "        nb_p = theta_umi / (theta_umi + mus)\n",
    "        np.random.seed(seed)\n",
    "        umi_counts = np.random.negative_binomial(np.ones_like(nb_p)*theta_umi, nb_p)\n",
    "\n",
    "        ###amplify UMIs to get readcounts\n",
    "        #intialize readcounts\n",
    "        readcounts=np.zeros(umi_counts.shape)\n",
    "        #only nonzero counts need amplification observations \n",
    "        umi_nonzero_idx = umi_counts>0\n",
    "        umi_nonzero_counts = umi_counts[umi_nonzero_idx]\n",
    "        #index for which amplification sample(s) get to amplify which UMIs\n",
    "        split_idx = np.cumsum(umi_nonzero_counts).astype(int)\n",
    "        #one large sample instead of separate ones is more efficient (one sample per molecule observed)\n",
    "        \n",
    "        mean_z = read_umi_factor\n",
    "        np.random.seed(seed)        \n",
    "        zs = np.random.geometric(1/mean_z,size=sum(umi_nonzero_counts))\n",
    "\n",
    "        #splitting up into separate groups of samples for each (gene-x-cell)-observation\n",
    "        zs_per_cell_x_gene=np.split(zs,split_idx[:-1])\n",
    "        #summing reads for each (gene-x-cell)-observation\n",
    "        zs_sums =[sum(z) for z in zs_per_cell_x_gene]\n",
    "        #mapping back to count matrix\n",
    "        readcounts[umi_nonzero_idx]=zs_sums\n",
    "\n",
    "        n_cells = sum(cluster_idx)\n",
    "        clusterlabels_sim += [cluster_i]*n_cells\n",
    "        clustercolors_sim += [adata.uns['clustercolors'][cluster_i]]*n_cells\n",
    "        readcounts_sim += [readcounts]\n",
    "        umi_counts_sim += [umi_counts]\n",
    "    \n",
    "    adata_sim = anndata.AnnData(X=sparse.csc_matrix(np.concatenate(readcounts_sim)))\n",
    "    adata_sim.obs['clustercolor'] = clustercolors_sim\n",
    "    adata_sim.obs['clusters'] = clusterlabels_sim\n",
    "    adata_sim.layers['umis_sim'] = sparse.csc_matrix(np.concatenate(umi_counts_sim))\n",
    "    adata_sim.uns['theta_umi']=theta_umi\n",
    "    adata_sim.uns['theta_z']=theta_z\n",
    "    adata_sim.uns['seed']=seed\n",
    "    adata_sim.uns['read_umi_factor']=read_umi_factor\n",
    "    adata_sim.var_names = adata.var_names.copy()\n",
    "    adata_sim.var['marker_idx'] = adata.var['marker_idx'].copy()\n",
    "    \n",
    "    return adata_sim\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "905c29d2-8d85-4634-9909-5fa8556dc675",
   "metadata": {},
   "outputs": [],
   "source": [
    "theta_umi = 100\n",
    "theta_z = 1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "27026c62-3f79-4499-9ce9-4777bbf8c6ad",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "simulating cluster 0..\n",
      "simulating cluster 1..\n",
      "simulating cluster 2..\n",
      "simulating cluster 3..\n",
      "simulating cluster 4..\n",
      "simulating cluster 5..\n",
      "simulating cluster 6..\n",
      "simulating cluster 7..\n",
      "simulating cluster 8..\n",
      "simulating cluster 9..\n",
      "simulating cluster 10..\n",
      "simulating cluster 11..\n",
      "simulating cluster 12..\n",
      "simulating cluster 13..\n",
      "simulating cluster 14..\n",
      "simulating cluster 15..\n",
      "simulating cluster 16..\n",
      "simulating cluster 17..\n",
      "simulating cluster 18..\n",
      "simulating cluster 19..\n",
      "simulating cluster 20..\n",
      "simulating cluster 21..\n",
      "simulating cluster 22..\n",
      "simulating cluster 23..\n",
      "simulating cluster 24..\n",
      "simulating cluster 25..\n",
      "simulating cluster 26..\n",
      "simulating cluster 27..\n",
      "simulating cluster 28..\n",
      "simulating cluster 29..\n",
      "simulating cluster 30..\n",
      "simulating cluster 31..\n",
      "simulating cluster 32..\n",
      "simulating cluster 33..\n",
      "simulating cluster 34..\n",
      "simulating cluster 35..\n",
      "simulating cluster 36..\n",
      "simulating cluster 37..\n",
      "simulating cluster 38..\n",
      "simulating cluster 39..\n",
      "simulating cluster 40..\n",
      "simulating cluster 41..\n",
      "simulating cluster 42..\n",
      "simulating cluster 43..\n",
      "simulating cluster 44..\n",
      "simulating cluster 45..\n",
      "simulating cluster 46..\n",
      "simulating cluster 47..\n",
      "simulating cluster 48..\n",
      "simulating cluster 49..\n",
      "simulating cluster 50..\n",
      "simulating cluster 51..\n",
      "simulating cluster 52..\n",
      "simulating cluster 53..\n",
      "simulating cluster 54..\n",
      "simulating cluster 55..\n",
      "simulating cluster 56..\n",
      "simulating cluster 57..\n",
      "simulating cluster 58..\n",
      "simulating cluster 59..\n",
      "simulating cluster 60..\n",
      "simulating cluster 61..\n",
      "simulating cluster 62..\n",
      "simulating cluster 63..\n",
      "simulating cluster 64..\n",
      "simulating cluster 65..\n",
      "simulating cluster 66..\n",
      "simulating cluster 67..\n",
      "simulating cluster 68..\n",
      "simulating cluster 69..\n",
      "simulating cluster 70..\n",
      "simulating cluster 71..\n",
      "simulating cluster 72..\n",
      "simulating cluster 73..\n",
      "simulating cluster 74..\n",
      "simulating cluster 75..\n",
      "simulating cluster 76..\n",
      "simulating cluster 77..\n",
      "simulating cluster 78..\n",
      "simulating cluster 79..\n",
      "simulating cluster 80..\n",
      "simulating cluster 81..\n",
      "simulating cluster 82..\n",
      "simulating cluster 83..\n",
      "simulating cluster 84..\n",
      "simulating cluster 85..\n",
      "simulating cluster 86..\n",
      "simulating cluster 87..\n",
      "simulating cluster 88..\n",
      "simulating cluster 89..\n",
      "simulating cluster 90..\n",
      "simulating cluster 91..\n",
      "simulating cluster 92..\n",
      "simulating cluster 93..\n",
      "simulating cluster 94..\n",
      "simulating cluster 95..\n",
      "simulating cluster 96..\n",
      "simulating cluster 97..\n",
      "simulating cluster 98..\n",
      "simulating cluster 99..\n",
      "simulating cluster 100..\n",
      "simulating cluster 101..\n",
      "simulating cluster 102..\n",
      "simulating cluster 103..\n",
      "simulating cluster 104..\n",
      "simulating cluster 105..\n",
      "simulating cluster 106..\n",
      "simulating cluster 107..\n",
      "simulating cluster 108..\n",
      "simulating cluster 109..\n",
      "simulating cluster 110..\n",
      "simulating cluster 111..\n",
      "simulating cluster 112..\n",
      "simulating cluster 113..\n",
      "simulating cluster 114..\n",
      "simulating cluster 115..\n",
      "simulating cluster 116..\n",
      "simulating cluster 117..\n",
      "simulating cluster 118..\n",
      "simulating cluster 119..\n",
      "simulating cluster 120..\n",
      "simulating cluster 121..\n",
      "simulating cluster 122..\n",
      "simulating cluster 123..\n",
      "simulating cluster 124..\n",
      "simulating cluster 125..\n",
      "simulating cluster 126..\n",
      "simulating cluster 127..\n",
      "simulating cluster 128..\n",
      "simulating cluster 129..\n",
      "simulating cluster 130..\n",
      "simulating cluster 131..\n",
      "simulating cluster 132..\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "tcmalloc: large alloc 8152080384 bytes == 0x4505e4000 @ \n",
      "tcmalloc: large alloc 1985470464 bytes == 0x63b02e000 @ \n",
      "/tmp/ipykernel_28877/474217884.py:64: FutureWarning: X.dtype being converted to np.float32 from float64. In the next version of anndata (0.9) conversion will not be automatic. Pass dtype explicitly to avoid this warning. Pass `AnnData(X, dtype=X.dtype, ...)` to get the future behavour.\n",
      "  adata_sim = anndata.AnnData(X=sparse.csc_matrix(np.concatenate(readcounts_sim)))\n",
      "tcmalloc: large alloc 8152080384 bytes == 0x4505e4000 @ \n",
      "tcmalloc: large alloc 1985470464 bytes == 0x72832c000 @ \n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "simulating cluster 0..\n",
      "simulating cluster 1..\n",
      "simulating cluster 2..\n",
      "simulating cluster 3..\n",
      "simulating cluster 4..\n",
      "simulating cluster 5..\n",
      "simulating cluster 6..\n",
      "simulating cluster 7..\n",
      "simulating cluster 8..\n",
      "simulating cluster 9..\n",
      "simulating cluster 10..\n",
      "simulating cluster 11..\n",
      "simulating cluster 12..\n",
      "simulating cluster 13..\n",
      "simulating cluster 14..\n",
      "simulating cluster 15..\n",
      "simulating cluster 16..\n",
      "simulating cluster 17..\n",
      "simulating cluster 18..\n",
      "simulating cluster 19..\n",
      "simulating cluster 20..\n",
      "simulating cluster 21..\n",
      "simulating cluster 22..\n",
      "simulating cluster 23..\n",
      "simulating cluster 24..\n",
      "simulating cluster 25..\n",
      "simulating cluster 26..\n",
      "simulating cluster 27..\n",
      "simulating cluster 28..\n",
      "simulating cluster 29..\n",
      "simulating cluster 30..\n",
      "simulating cluster 31..\n",
      "simulating cluster 32..\n",
      "simulating cluster 33..\n",
      "simulating cluster 34..\n",
      "simulating cluster 35..\n",
      "simulating cluster 36..\n",
      "simulating cluster 37..\n",
      "simulating cluster 38..\n",
      "simulating cluster 39..\n",
      "simulating cluster 40..\n",
      "simulating cluster 41..\n",
      "simulating cluster 42..\n",
      "simulating cluster 43..\n",
      "simulating cluster 44..\n",
      "simulating cluster 45..\n",
      "simulating cluster 46..\n",
      "simulating cluster 47..\n",
      "simulating cluster 48..\n",
      "simulating cluster 49..\n",
      "simulating cluster 50..\n",
      "simulating cluster 51..\n",
      "simulating cluster 52..\n",
      "simulating cluster 53..\n",
      "simulating cluster 54..\n",
      "simulating cluster 55..\n",
      "simulating cluster 56..\n",
      "simulating cluster 57..\n",
      "simulating cluster 58..\n",
      "simulating cluster 59..\n",
      "simulating cluster 60..\n",
      "simulating cluster 61..\n",
      "simulating cluster 62..\n",
      "simulating cluster 63..\n",
      "simulating cluster 64..\n",
      "simulating cluster 65..\n",
      "simulating cluster 66..\n",
      "simulating cluster 67..\n",
      "simulating cluster 68..\n",
      "simulating cluster 69..\n",
      "simulating cluster 70..\n",
      "simulating cluster 71..\n",
      "simulating cluster 72..\n",
      "simulating cluster 73..\n",
      "simulating cluster 74..\n",
      "simulating cluster 75..\n",
      "simulating cluster 76..\n",
      "simulating cluster 77..\n",
      "simulating cluster 78..\n",
      "simulating cluster 79..\n",
      "simulating cluster 80..\n",
      "simulating cluster 81..\n",
      "simulating cluster 82..\n",
      "simulating cluster 83..\n",
      "simulating cluster 84..\n",
      "simulating cluster 85..\n",
      "simulating cluster 86..\n",
      "simulating cluster 87..\n",
      "simulating cluster 88..\n",
      "simulating cluster 89..\n",
      "simulating cluster 90..\n",
      "simulating cluster 91..\n",
      "simulating cluster 92..\n",
      "simulating cluster 93..\n",
      "simulating cluster 94..\n",
      "simulating cluster 95..\n",
      "simulating cluster 96..\n",
      "simulating cluster 97..\n",
      "simulating cluster 98..\n",
      "simulating cluster 99..\n",
      "simulating cluster 100..\n",
      "simulating cluster 101..\n",
      "simulating cluster 102..\n",
      "simulating cluster 103..\n",
      "simulating cluster 104..\n",
      "simulating cluster 105..\n",
      "simulating cluster 106..\n",
      "simulating cluster 107..\n",
      "simulating cluster 108..\n",
      "simulating cluster 109..\n",
      "simulating cluster 110..\n",
      "simulating cluster 111..\n",
      "simulating cluster 112..\n",
      "simulating cluster 113..\n",
      "simulating cluster 114..\n",
      "simulating cluster 115..\n",
      "simulating cluster 116..\n",
      "simulating cluster 117..\n",
      "simulating cluster 118..\n",
      "simulating cluster 119..\n",
      "simulating cluster 120..\n",
      "simulating cluster 121..\n",
      "simulating cluster 122..\n",
      "simulating cluster 123..\n",
      "simulating cluster 124..\n",
      "simulating cluster 125..\n",
      "simulating cluster 126..\n",
      "simulating cluster 127..\n",
      "simulating cluster 128..\n",
      "simulating cluster 129..\n",
      "simulating cluster 130..\n",
      "simulating cluster 131..\n",
      "simulating cluster 132..\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "tcmalloc: large alloc 8152080384 bytes == 0x3e4c00000 @ \n",
      "tcmalloc: large alloc 2127740928 bytes == 0x79ddc8000 @ \n",
      "/tmp/ipykernel_28877/474217884.py:64: FutureWarning: X.dtype being converted to np.float32 from float64. In the next version of anndata (0.9) conversion will not be automatic. Pass dtype explicitly to avoid this warning. Pass `AnnData(X, dtype=X.dtype, ...)` to get the future behavour.\n",
      "  adata_sim = anndata.AnnData(X=sparse.csc_matrix(np.concatenate(readcounts_sim)))\n",
      "tcmalloc: large alloc 8152080384 bytes == 0x3e4c00000 @ \n"
     ]
    }
   ],
   "source": [
    "adata_sim = simulate(adata_tasic,theta_umi=theta_umi,theta_z=theta_z,only_marker_genes_DE=False)\n",
    "adata_sim_only_marker_genes_DE = simulate(adata_tasic,theta_umi=theta_umi,theta_z=theta_z,only_marker_genes_DE=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "4071c2d7-f4af-4d76-aaa3-14a00718e0d2",
   "metadata": {},
   "outputs": [],
   "source": [
    "adata_sim.write_h5ad(f'data/tasic/simulations/adata_sim_thetaUMI_{theta_umi}_GeomZ.h5ad')\n",
    "adata_sim_only_marker_genes_DE.write_h5ad(f'data/tasic/simulations/adata_sim_only_marker_genes_DE_thetaUMI_{theta_umi}_GeomZ.h5ad')"
   ]
  }
 ],
 "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"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}