12_prepare_ziegenhain.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "import gzip\n",
    "import pickle\n",
    "import anndata\n",
    "from scipy import sparse\n",
    "import mygene"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "basepath = 'data/reads_per_umi_tables/'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "filename_multiexp = 'Ziegenhain2017.hd1.txt.gz'\n",
    "filename_smartseq3 = 'Hagemann-Jensen2020_Smartseq3_SE.hd1.txt.gz' \n",
    "filename_smartseq3_PE ='Johnsson2022_Smartseq3_PE.hd1.txt.gz'\n",
    "filename_smartseq3xpress = 'Hagemann-Jensen2022_Smartseq3xpress.hd1.txt.gz'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "def preload_long_data(path,filename):\n",
    "    \n",
    "    with gzip.open(f'{path}{filename}') as f:\n",
    "        \n",
    "        lines=[]\n",
    "        for i,line in enumerate(f):\n",
    "\n",
    "            line_decoded = line.decode(\"utf-8\").replace('\\n','').split('\\t')\n",
    "\n",
    "            if i == 0:\n",
    "                header = line_decoded\n",
    "            else:\n",
    "                lines.append(line_decoded)\n",
    "    df_long = pd.DataFrame(np.array(lines),columns=header)\n",
    "    df_long['N'] = df_long['N'].astype(int)\n",
    "    \n",
    "    return df_long"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "tcmalloc: large alloc 10691682304 bytes == 0x344bc000 @ \n",
      "tcmalloc: large alloc 8619941888 bytes == 0x7b1ae000 @ \n",
      "tcmalloc: large alloc 22113992704 bytes == 0x2ece26000 @ \n",
      "tcmalloc: large alloc 1922957312 bytes == 0xb1540000 @ \n",
      "tcmalloc: large alloc 2221088768 bytes == 0x93344000 @ \n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 9min 13s, sys: 1min 37s, total: 10min 51s\n",
      "Wall time: 10min 54s\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "# 17min\n",
    "df_multiexp = preload_long_data(basepath,filename_multiexp)\n",
    "df_smartseq3 = preload_long_data(basepath,filename_smartseq3)  \n",
    "df_smartseq3_PE = preload_long_data(basepath,filename_smartseq3_PE)\n",
    "df_smartseq3xpress = preload_long_data(basepath,filename_smartseq3xpress)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "tcmalloc: large alloc 2550333440 bytes == 0x17fa3a000 @ \n"
     ]
    }
   ],
   "source": [
    "df_multiexp.rename(columns={\"experiment\": \"Experiment\"},inplace=True)\n",
    "\n",
    "df_smartseq3xpress['Experiment'] = 'Smartseq3xpressA'\n",
    "\n",
    "df_multiexp.set_index('Experiment',inplace=True)\n",
    "df_smartseq3.set_index('Experiment',inplace=True)\n",
    "df_smartseq3_PE.set_index('Experiment',inplace=True)\n",
    "df_smartseq3xpress.set_index('Experiment',inplace=True)\n",
    "\n",
    "df_all = pd.concat([df_multiexp,df_smartseq3,df_smartseq3_PE,df_smartseq3xpress],join='inner')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "tcmalloc: large alloc 2550333440 bytes == 0x31f8e2000 @ \n"
     ]
    }
   ],
   "source": [
    "dfs = [exp_df_tuple[1] for exp_df_tuple in df_all.groupby('Experiment')]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "tcmalloc: large alloc 4294967296 bytes == 0x49ce8a000 @ \n",
      "tcmalloc: large alloc 8589934592 bytes == 0x59ce8a000 @ \n",
      "tcmalloc: large alloc 17179869184 bytes == 0x827e90000 @ \n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 3min 24s, sys: 54.8 s, total: 4min 19s\n",
      "Wall time: 4min 20s\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "#4min\n",
    "with open(f'{basepath}reads_per_umi_dfs_hd1_preproc.pickle','wb') as f:\n",
    "    pickle.dump(dfs,f) "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Prepare anndata objects for each experiment"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "experiments = [df.index[0] for df in dfs]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "['CELseq2A',\n",
       " 'CELseq2B',\n",
       " 'DropSeqA',\n",
       " 'DropSeqB',\n",
       " 'MARSseqA',\n",
       " 'MARSseqB',\n",
       " 'SCRBseqA',\n",
       " 'SCRBseqB',\n",
       " 'Smartseq3_Fibroblast',\n",
       " 'Smartseq3_Fibroblast_PE',\n",
       " 'Smartseq3xpressA']"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "experiments"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CELseq2A\n",
      "unpacking..\n",
      "extracting cells.................................."
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "<timed exec>:42: FutureWarning: X.dtype being converted to np.float32 from int64. 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"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "read/umi counts: AnnData object with n_obs × n_vars = 34 × 23555\n",
      "    obs: 'cells'\n",
      "    uns: 'experiment', 'type'\n",
      "    layers: 'reads'\n",
      "CELseq2B\n",
      "unpacking..\n",
      "extracting cells....................................."
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "<timed exec>:42: FutureWarning: X.dtype being converted to np.float32 from int64. 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"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "read/umi counts: AnnData object with n_obs × n_vars = 37 × 25479\n",
      "    obs: 'cells'\n",
      "    uns: 'experiment', 'type'\n",
      "    layers: 'reads'\n",
      "DropSeqA\n",
      "unpacking..\n",
      "extracting cells.........................................."
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "<timed exec>:42: FutureWarning: X.dtype being converted to np.float32 from int64. 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"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "read/umi counts: AnnData object with n_obs × n_vars = 42 × 23530\n",
      "    obs: 'cells'\n",
      "    uns: 'experiment', 'type'\n",
      "    layers: 'reads'\n",
      "DropSeqB\n",
      "unpacking..\n",
      "extracting cells.................................."
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "<timed exec>:42: FutureWarning: X.dtype being converted to np.float32 from int64. 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"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "read/umi counts: AnnData object with n_obs × n_vars = 34 × 22104\n",
      "    obs: 'cells'\n",
      "    uns: 'experiment', 'type'\n",
      "    layers: 'reads'\n",
      "MARSseqA\n",
      "unpacking..\n",
      "extracting cells............................."
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "<timed exec>:42: FutureWarning: X.dtype being converted to np.float32 from int64. 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"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "read/umi counts: AnnData object with n_obs × n_vars = 29 × 20307\n",
      "    obs: 'cells'\n",
      "    uns: 'experiment', 'type'\n",
      "    layers: 'reads'\n",
      "MARSseqB\n",
      "unpacking..\n",
      "extracting cells...................................."
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "<timed exec>:42: FutureWarning: X.dtype being converted to np.float32 from int64. 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"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "read/umi counts: AnnData object with n_obs × n_vars = 36 × 21150\n",
      "    obs: 'cells'\n",
      "    uns: 'experiment', 'type'\n",
      "    layers: 'reads'\n",
      "SCRBseqA\n",
      "unpacking..\n",
      "extracting cells......................................."
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "<timed exec>:42: FutureWarning: X.dtype being converted to np.float32 from int64. 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"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "read/umi counts: AnnData object with n_obs × n_vars = 39 × 24054\n",
      "    obs: 'cells'\n",
      "    uns: 'experiment', 'type'\n",
      "    layers: 'reads'\n",
      "SCRBseqB\n",
      "unpacking..\n",
      "extracting cells............................................."
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "<timed exec>:42: FutureWarning: X.dtype being converted to np.float32 from int64. 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"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "read/umi counts: AnnData object with n_obs × n_vars = 45 × 24281\n",
      "    obs: 'cells'\n",
      "    uns: 'experiment', 'type'\n",
      "    layers: 'reads'\n",
      "Smartseq3_Fibroblast\n",
      "unpacking..\n",
      "extracting cells................................................................................................................................................."
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "<timed exec>:42: FutureWarning: X.dtype being converted to np.float32 from int64. 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"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "read/umi counts: AnnData object with n_obs × n_vars = 145 × 23161\n",
      "    obs: 'cells'\n",
      "    uns: 'experiment', 'type'\n",
      "    layers: 'reads'\n",
      "Smartseq3_Fibroblast_PE\n",
      "unpacking..\n",
      "extracting cells..............................................................................................................................................................................................................................................................................................................................."
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "<timed exec>:42: FutureWarning: X.dtype being converted to np.float32 from int64. 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"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "read/umi counts: AnnData object with n_obs × n_vars = 319 × 25427\n",
      "    obs: 'cells'\n",
      "    uns: 'experiment', 'type'\n",
      "    layers: 'reads'\n",
      "Smartseq3xpressA\n",
      "unpacking..\n",
      "extracting cells................................................................................................read/umi counts: AnnData object with n_obs × n_vars = 96 × 25595\n",
      "    obs: 'cells'\n",
      "    uns: 'experiment', 'type'\n",
      "    layers: 'reads'\n",
      "CPU times: user 10h 2min 44s, sys: 19min 33s, total: 10h 22min 17s\n",
      "Wall time: 10h 21min 58s\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "<timed exec>:42: FutureWarning: X.dtype being converted to np.float32 from int64. 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"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "ads = []\n",
    "for exp,df_exp in zip(experiments,dfs):\n",
    "    df_exp.reset_index(inplace=True)\n",
    "    print(exp)\n",
    "   \n",
    "    #set index as multilevel index\n",
    "    df_idxd = df_exp.set_index(keys = ['Experiment','RG','GE','UB'])\n",
    "    #change remaining count column to sparse datatype\n",
    "    series_sparse = df_idxd['N'].astype('Sparse[int32]')\n",
    "    #unpack long raw data into sparse coo format (which is implicitly wide) \n",
    "    print('unpacking..')\n",
    "    mtx_sparse, rows_idxs, column_idxs = series_sparse.sparse.to_coo(row_levels=[\"Experiment\", \"RG\", 'UB'], column_levels=[\"GE\"],sort_labels=True)\n",
    "            \n",
    "    #extract row and column annotations\n",
    "    genes = column_idxs\n",
    "    reads_per_umi_cells = [r[1] for r in rows_idxs]\n",
    "    reads_per_umi_umis =  [r[2] for r in rows_idxs]\n",
    "    \n",
    "    #make sparse dataframe to extract read/UMI counts per cell\n",
    "    df_sparse = pd.DataFrame.sparse.from_spmatrix(mtx_sparse,index=[reads_per_umi_cells,reads_per_umi_umis],columns=column_idxs)\n",
    "    df_sparse.index.set_names(['cells','umis'],inplace=True)\n",
    "    \n",
    "    #group by cell\n",
    "    df_grouped = df_sparse.groupby(by='cells')\n",
    "    #extract read and UMI counts per cell\n",
    "    UMIs = []\n",
    "    reads = []\n",
    "    cells = []\n",
    "    print('extracting cells',end='')\n",
    "    for idx,group in df_grouped:\n",
    "        print('.',end='')\n",
    "        UMI = np.sum(group>0,axis=0) #count number of UMIs per gene\n",
    "        read = np.sum(group,axis=0)  #sum number of reads per gene\n",
    "        UMIs.append(UMI.values)\n",
    "        reads.append(read.values)\n",
    "        cells.append(idx)\n",
    "    UMIs = np.vstack(UMIs)\n",
    "    reads= np.vstack(reads)\n",
    "\n",
    "\n",
    "    #make read/umi count adata\n",
    "    ad_readcounts_umicounts = anndata.AnnData(sparse.csr_matrix(UMIs),\n",
    "                             var=dict(genes=genes),\n",
    "                             obs=dict(cells=cells),\n",
    "                             uns=dict(experiment=exp,\n",
    "                                      type='UMI counts in X, read counts in `layers`'),\n",
    "                             layers=dict(reads=sparse.csr_matrix(reads)))\n",
    "    \n",
    "    ad_readcounts_umicounts.var.set_index('genes',inplace=True)\n",
    "\n",
    "    print('read/umi counts:',ad_readcounts_umicounts)\n",
    "    ad_readcounts_umicounts.write_h5ad(f'{basepath}ad_readcounts_umicounts_hd1_{exp}')\n",
    "    \n",
    "    ads.append(ad_readcounts_umicounts)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "def add_annotations(ad,loadfile='',savefile='',species='mouse'):\n",
    "    \n",
    "    mg = mygene.MyGeneInfo()\n",
    "    \n",
    "    if not loadfile:\n",
    "        print('Not loading, re-fetching annotations..')\n",
    "        gene_names = ad.var.index\n",
    "        gene_annotations = mg.querymany(list(gene_names), scopes='ensemblgene',species=species,fields=['name','symbol','type_of_gene','ensembl.type_of_gene'],as_dataframe=True)\n",
    "        \n",
    "        #if any genes gave duplicate hits, we dont assign the annotations and add the `not_uniquely_assignable` category\n",
    "        gene_annotation_names,n_hits = np.unique(gene_annotations.index,return_counts=True)\n",
    "        multi_hits = gene_annotation_names[n_hits>1]\n",
    "        gene_annotations_dedup = gene_annotations.drop(multi_hits)\n",
    "        #make and append one row for each of the multi hits\n",
    "        duplicates_df = pd.DataFrame({clm:{dup:'not_uniquely_assignable' for dup in multi_hits} for clm in gene_annotations.columns})\n",
    "        gene_annotations = pd.concat([gene_annotations_dedup,duplicates_df])\n",
    "        #after that the number of genes should match again\n",
    "        gene_annotations = gene_annotations.fillna('missing')\n",
    "\n",
    "    else:\n",
    "        print('Loading annotations from',loadfile)\n",
    "        with open(loadfile,'rb') as f:\n",
    "            gene_annotations = pickle.load(f)\n",
    "        gene_names = ad.var.index\n",
    "\n",
    "    #update ad.var\n",
    "    columns_to_keep = ['type_of_gene','symbol','name','ensembl.type_of_gene']\n",
    "    columns_to_drop = [c for c in gene_annotations.columns if c not in columns_to_keep]\n",
    "    gene_annotations_slim = gene_annotations.drop(columns=columns_to_drop)\n",
    "    \n",
    "    ad.var = gene_annotations_slim.loc[ad.var.index,:]\n",
    "    \n",
    "    #add ERCC info\n",
    "    ercc_idx = np.array([s[:5]=='gERCC' for s in ad.var.index])\n",
    "    ad.var['type_of_gene'][ercc_idx] = 'artificial/spike-in'\n",
    "    ad.var['ensembl.type_of_gene'][ercc_idx] = 'artificial/spike-in'\n",
    "    ad.var['name'][ercc_idx] = ad.var.index[ercc_idx]\n",
    "    ad.var['symbol'][ercc_idx] = [s.split('-')[1] for s in ad.var.index[ercc_idx]]\n",
    "\n",
    "    if savefile:\n",
    "        print('Saving annotations to',savefile)\n",
    "        with open(savefile,'wb') as f:\n",
    "            pickle.dump(gene_annotations,f)\n",
    "        \n",
    "def add_coarse_gene_annotations_single_exp(ad):\n",
    "    \n",
    "    types_master = np.unique(ad.var['ensembl.type_of_gene'])\n",
    "    coarse_types_pseudo_idx = ['pseudogene' in t for t in types_master]\n",
    "    \n",
    "    fine2coarse_map = {fine:'pseudogene' for fine in types_master[coarse_types_pseudo_idx]}\n",
    "    fine2coarse_map['artificial/spike-in']='artificial/spike-in'\n",
    "    fine2coarse_map['missing']='missing'\n",
    "    fine2coarse_map['not_uniquely_assignable']='missing'\n",
    "    fine2coarse_map['protein_coding']='protein_coding'\n",
    "\n",
    "    other_types = types_master[~np.isin(types_master,list(fine2coarse_map.keys()))]\n",
    "    for other_t in other_types:\n",
    "        fine2coarse_map[other_t] = 'other'\n",
    "\n",
    "    ad.var['coarse_types'] = ad.var['ensembl.type_of_gene'].copy()\n",
    "    ad.var.replace({'coarse_types':fine2coarse_map},inplace=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/usr/local/lib/python3.8/dist-packages/anndata/_core/anndata.py:1828: UserWarning: Observation names are not unique. To make them unique, call `.obs_names_make_unique`.\n",
      "  utils.warn_names_duplicates(\"obs\")\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Not loading, re-fetching annotations..\n",
      "querying 1-1000...done.\n",
      "querying 1001-2000...done.\n",
      "querying 2001-3000...done.\n",
      "querying 3001-4000...done.\n",
      "querying 4001-5000...done.\n",
      "querying 5001-6000...done.\n",
      "querying 6001-7000...done.\n",
      "querying 7001-8000...done.\n",
      "querying 8001-9000...done.\n",
      "querying 9001-10000...done.\n",
      "querying 10001-11000...done.\n",
      "querying 11001-12000...done.\n",
      "querying 12001-13000...done.\n",
      "querying 13001-14000...done.\n",
      "querying 14001-15000...done.\n",
      "querying 15001-16000...done.\n",
      "querying 16001-17000...done.\n",
      "querying 17001-18000...done.\n",
      "querying 18001-19000...done.\n",
      "querying 19001-20000...done.\n",
      "querying 20001-21000...done.\n",
      "querying 21001-22000...done.\n",
      "querying 22001-23000...done.\n",
      "querying 23001-24000...done.\n",
      "querying 24001-25000...done.\n",
      "querying 25001-26000...done.\n",
      "querying 26001-27000...done.\n",
      "querying 27001-28000...done.\n",
      "querying 28001-29000...done.\n",
      "querying 29001-30000...done.\n",
      "querying 30001-31000...done.\n",
      "querying 31001-32000...done.\n",
      "querying 32001-33000...done.\n",
      "querying 33001-34000...done.\n",
      "querying 34001-35000...done.\n",
      "querying 35001-36000...done.\n",
      "querying 36001-37000...done.\n",
      "querying 37001-38000...done.\n",
      "querying 38001-39000...done.\n",
      "querying 39001-40000...done.\n",
      "querying 40001-41000...done.\n",
      "querying 41001-42000...done.\n",
      "querying 42001-43000...done.\n",
      "querying 43001-44000...done.\n",
      "querying 44001-45000...done.\n",
      "querying 45001-46000...done.\n",
      "querying 46001-47000...done.\n",
      "querying 47001-48000...done.\n",
      "querying 48001-49000...done.\n",
      "querying 49001-50000...done.\n",
      "querying 50001-51000...done.\n",
      "querying 51001-52000...done.\n",
      "querying 52001-53000...done.\n",
      "querying 53001-54000...done.\n",
      "querying 54001-55000...done.\n",
      "querying 55001-56000...done.\n",
      "querying 56001-57000...done.\n",
      "querying 57001-58000...done.\n",
      "querying 58001-59000...done.\n",
      "querying 59001-60000...done.\n",
      "querying 60001-61000...done.\n",
      "querying 61001-62000...done.\n",
      "querying 62001-63000...done.\n",
      "querying 63001-64000...done.\n",
      "querying 64001-64298...done.\n",
      "Finished.\n",
      "6 input query terms found dup hits:\n",
      "\t[('ENSG00000249738', 2), ('ENSG00000278903', 3), ('ENSMUSG00000094002', 2), ('ENSMUSG00000095706', 2\n",
      "1626 input query terms found no hit:\n",
      "\t['ENSG00000112096', 'ENSG00000116883', 'ENSG00000130489', 'ENSG00000130723', 'ENSG00000154537', 'ENS\n",
      "Pass \"returnall=True\" to return complete lists of duplicate or missing query terms.\n",
      "Saving annotations to data/reads_per_umi_tables/mygene_annotations.pickle\n",
      "Loading annotations from data/reads_per_umi_tables/mygene_annotations.pickle\n",
      "read/umi counts: AnnData object with n_obs × n_vars = 34 × 23555\n",
      "    obs: 'cells'\n",
      "    var: 'name', 'symbol', 'type_of_gene', 'ensembl.type_of_gene', 'coarse_types'\n",
      "    uns: 'experiment', 'type'\n",
      "    layers: 'reads'\n",
      "Loading annotations from data/reads_per_umi_tables/mygene_annotations.pickle\n",
      "read/umi counts: AnnData object with n_obs × n_vars = 37 × 25479\n",
      "    obs: 'cells'\n",
      "    var: 'name', 'symbol', 'type_of_gene', 'ensembl.type_of_gene', 'coarse_types'\n",
      "    uns: 'experiment', 'type'\n",
      "    layers: 'reads'\n",
      "Loading annotations from data/reads_per_umi_tables/mygene_annotations.pickle\n",
      "read/umi counts: AnnData object with n_obs × n_vars = 42 × 23530\n",
      "    obs: 'cells'\n",
      "    var: 'name', 'symbol', 'type_of_gene', 'ensembl.type_of_gene', 'coarse_types'\n",
      "    uns: 'experiment', 'type'\n",
      "    layers: 'reads'\n",
      "Loading annotations from data/reads_per_umi_tables/mygene_annotations.pickle\n",
      "read/umi counts: AnnData object with n_obs × n_vars = 34 × 22104\n",
      "    obs: 'cells'\n",
      "    var: 'name', 'symbol', 'type_of_gene', 'ensembl.type_of_gene', 'coarse_types'\n",
      "    uns: 'experiment', 'type'\n",
      "    layers: 'reads'\n",
      "Loading annotations from data/reads_per_umi_tables/mygene_annotations.pickle\n",
      "read/umi counts: AnnData object with n_obs × n_vars = 29 × 20307\n",
      "    obs: 'cells'\n",
      "    var: 'name', 'symbol', 'type_of_gene', 'ensembl.type_of_gene', 'coarse_types'\n",
      "    uns: 'experiment', 'type'\n",
      "    layers: 'reads'\n",
      "Loading annotations from data/reads_per_umi_tables/mygene_annotations.pickle\n",
      "read/umi counts: AnnData object with n_obs × n_vars = 36 × 21150\n",
      "    obs: 'cells'\n",
      "    var: 'name', 'symbol', 'type_of_gene', 'ensembl.type_of_gene', 'coarse_types'\n",
      "    uns: 'experiment', 'type'\n",
      "    layers: 'reads'\n",
      "Loading annotations from data/reads_per_umi_tables/mygene_annotations.pickle\n",
      "read/umi counts: AnnData object with n_obs × n_vars = 39 × 24054\n",
      "    obs: 'cells'\n",
      "    var: 'name', 'symbol', 'type_of_gene', 'ensembl.type_of_gene', 'coarse_types'\n",
      "    uns: 'experiment', 'type'\n",
      "    layers: 'reads'\n",
      "Loading annotations from data/reads_per_umi_tables/mygene_annotations.pickle\n",
      "read/umi counts: AnnData object with n_obs × n_vars = 45 × 24281\n",
      "    obs: 'cells'\n",
      "    var: 'name', 'symbol', 'type_of_gene', 'ensembl.type_of_gene', 'coarse_types'\n",
      "    uns: 'experiment', 'type'\n",
      "    layers: 'reads'\n",
      "Loading annotations from data/reads_per_umi_tables/mygene_annotations.pickle\n",
      "read/umi counts: AnnData object with n_obs × n_vars = 145 × 23161\n",
      "    obs: 'cells'\n",
      "    var: 'name', 'symbol', 'type_of_gene', 'ensembl.type_of_gene', 'coarse_types'\n",
      "    uns: 'experiment', 'type'\n",
      "    layers: 'reads'\n",
      "Loading annotations from data/reads_per_umi_tables/mygene_annotations.pickle\n",
      "read/umi counts: AnnData object with n_obs × n_vars = 319 × 25427\n",
      "    obs: 'cells'\n",
      "    var: 'name', 'symbol', 'type_of_gene', 'ensembl.type_of_gene', 'coarse_types'\n",
      "    uns: 'experiment', 'type'\n",
      "    layers: 'reads'\n",
      "Loading annotations from data/reads_per_umi_tables/mygene_annotations.pickle\n",
      "read/umi counts: AnnData object with n_obs × n_vars = 96 × 25595\n",
      "    obs: 'cells'\n",
      "    var: 'name', 'symbol', 'type_of_gene', 'ensembl.type_of_gene', 'coarse_types'\n",
      "    uns: 'experiment', 'type'\n",
      "    layers: 'reads'\n"
     ]
    }
   ],
   "source": [
    "ad_all = anndata.concat(ads,join='outer')\n",
    "#annotate all ads at once\n",
    "annotation_filename = f'{basepath}mygene_annotations.pickle'\n",
    "add_annotations(ad_all,savefile=annotation_filename,species='human,mouse')\n",
    "\n",
    "for exp,ad in zip(experiments,ads):\n",
    "    add_annotations(ad,loadfile=annotation_filename)\n",
    "    add_coarse_gene_annotations_single_exp(ad)\n",
    "    print('read/umi counts:',ad)\n",
    "    ad.write_h5ad(f'{basepath}ad_readcounts_umicounts_hd1_{exp}')"
   ]
  },
  {
   "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"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}