0_Data_download_and_processing.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "tracked-cartridge",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "from bokeh.io import show, output_notebook, reset_output\n",
    "from IPython.display import display, HTML\n",
    "import matplotlib.pyplot as plt\n",
    "import matplotlib as mpl\n",
    "import numpy as np\n",
    "import os\n",
    "import pandas as pd\n",
    "import scanpy as sc\n",
    "import scanpy.external as sce\n",
    "import seaborn as sns\n",
    "from scipy.sparse import csr_matrix\n",
    "import subprocess\n",
    "from tqdm import tqdm\n",
    "import triku as tk\n",
    "import bbknn\n",
    "\n",
    "reset_output()\n",
    "output_notebook()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5fe3d0db-510b-4a4b-b5a3-cda0371474a6",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "from cellassign import assign_cats"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "yellow-fishing",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "data_dir = 'data/'\n",
    "fig_dir = 'figures/'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "statutory-founder",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# Palettes for UMAP gene expression\n",
    "\n",
    "magma = [plt.get_cmap('magma')(i) for i in np.linspace(0,1, 80)]\n",
    "magma[0] = (0.88, 0.88, 0.88, 1)\n",
    "magma = mpl.colors.LinearSegmentedColormap.from_list(\"\", magma[:65])\n",
    "\n",
    "seed = 0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8c2d0ef6-93a2-4dc0-88b3-57da52c538a0",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "dict_cats_general = {'FAP Lum+': ['Apod', 'Lum', 'Ly6a', 'Pdgfra', 'Mfap5', 'Dcn'], \n",
    "                     'FAP Prg4+': ['Prg4', 'Fbn1', 'Ly6a', 'Pdgfra', 'Mfap5', 'Dcn'], \n",
    "                     'Endothelial': ['Pecam1', 'Kdr', 'Fabp4', 'Cav1', 'Cdh5', 'Tek'], \n",
    "                     'Pericyte': ['Rgs5', 'Notch3', 'Myl9', 'Ndufa4l2', 'Itga7', 'Myh11', 'Pln', 'Abcc9'], \n",
    "                     'Satellite cell': ['Pax7', 'Myod1', 'Chodl', 'Vcam1', 'Sdc4', 'Myf5',], \n",
    "                     'Myonuclei': ['Tnnc2', 'Myh4', 'Acta1', 'Ckm', 'Tpm2', 'Eno3', 'Slc25a4'], \n",
    "                     'Tenocyte': ['Scx', 'Tnmd', 'Mkx', 'Col12a1', 'Col1a1', 'Tnc', 'Fmod', 'Comp'], \n",
    "                     'Neural cell': ['Mpz', 'Ptn', 'S100b'], \n",
    "                     'Glial cell': ['Plp1', 'Kcna1', 'S100b', 'Mbp', 'Mpz',],\n",
    "                     'Guide cell': ['Ncam2'],\n",
    "                     'Immune': ['H2-Aa', 'Cd74'], \n",
    "                     'APC': ['H2-Eb1', 'H2-Ab1'],\n",
    "                     'B cell': ['Cd19', 'Cd22', 'Ms4a1', 'Ptprc'], \n",
    "                     'T cell': ['Cd3d', 'Cd3e', 'Cd3g', 'Cd8a', 'Cd4', 'Ptprc', 'Cd28'], \n",
    "                     'Monocyte': ['Csf1r', 'Adgre1'], \n",
    "                     'Macrophage': ['Itgam', 'Csf1r', 'Adgre1', 'Itgb1', 'Cd68'],\n",
    "                     'Neutrophil': ['S100a8', 'S100a9', 'Itgam', 'Cd14', ], \n",
    "                     }"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "alien-recognition",
   "metadata": {},
   "source": [
    "# Oprescu 2020"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "dependent-classics",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "oprescu_dir = data_dir + '/oprescu'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "handmade-headset",
   "metadata": {},
   "outputs": [],
   "source": [
    "link = 'https://ftp.ncbi.nlm.nih.gov/geo/series/GSE138nnn/GSE138826/suppl/GSE138826%5Fexpression%5Fmatrix%2Etxt%2Egz'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "korean-freeze",
   "metadata": {},
   "outputs": [],
   "source": [
    "!aria2c -x 16 {link} -d {oprescu_dir}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "baking-software",
   "metadata": {},
   "outputs": [],
   "source": [
    "adata_oprescu = sc.read_text(oprescu_dir+'/GSE138826_expression_matrix.txt.gz').transpose()\n",
    "adata_oprescu.write_loom(oprescu_dir+'/adata_oprescu.loom')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "measured-entry",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "adata_oprescu = sc.read_loom(oprescu_dir+'/adata_oprescu.loom')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "708b4ac1-88bf-463d-9c4b-c322c6c1314b",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "adata_oprescu.obs.index = adata_oprescu.obs['obs_names']\n",
    "adata_oprescu.var.index = adata_oprescu.var['var_names']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "increased-pride",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "adata_oprescu.obs['batch'] = [i.split('_')[0] for i in adata_oprescu.obs_names]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "dietary-parade",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# Basic QC filtering\n",
    "adata_oprescu.var['mt'] = adata_oprescu.var_names.str.startswith('mt-')  # annotate the group of mitochondrial genes as 'mt'\n",
    "sc.pp.calculate_qc_metrics(adata_oprescu, qc_vars=['mt'], percent_top=None, inplace=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "vocational-trouble",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "sc.pl.violin(adata_oprescu, ['n_genes_by_counts', 'total_counts', 'pct_counts_mt'],\n",
    "             jitter=0.4, multi_panel=True)\n",
    "\n",
    "sc.pl.scatter(adata_oprescu, x='total_counts', y='pct_counts_mt')\n",
    "sc.pl.scatter(adata_oprescu, x='total_counts', y='n_genes_by_counts', color='batch')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "broad-league",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "sc.pp.filter_cells(adata_oprescu, min_genes=150)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d4a7ae4c-67e3-4c61-92cf-fdad7a65d5b7",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "adata_oprescu.obs['batch']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fifteen-tennessee",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "adata_oprescu_d0 = adata_oprescu[adata_oprescu.obs['batch'] == 'Noninjured'].copy()\n",
    "adata_oprescu_d05 = adata_oprescu[adata_oprescu.obs['batch'] == 'X0.5.DPI'].copy()\n",
    "adata_oprescu_d2 = adata_oprescu[adata_oprescu.obs['batch'] == 'X2.DPI'].copy()\n",
    "adata_oprescu_d35 = adata_oprescu[adata_oprescu.obs['batch'] == 'X3.5.DPI'].copy()\n",
    "adata_oprescu_d5 = adata_oprescu[adata_oprescu.obs['batch'] == 'X5.DPI'].copy()\n",
    "adata_oprescu_d10 = adata_oprescu[adata_oprescu.obs['batch'] == 'X10.DPI'].copy()\n",
    "adata_oprescu_d21 = adata_oprescu[adata_oprescu.obs['batch'] == 'X21.DPI'].copy()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "above-hampton",
   "metadata": {},
   "source": [
    "## Individual adata processing"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5cfd6414-a672-43dd-9275-a241bc5bc008",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "for adata in [adata_oprescu_d0, adata_oprescu_d05, adata_oprescu_d2, adata_oprescu_d35, adata_oprescu_d5, adata_oprescu_d10, adata_oprescu_d21]:\n",
    "    adata.raw= adata\n",
    "    sc.pp.filter_genes(adata, min_counts=1)\n",
    "    sc.pp.normalize_per_cell(adata)\n",
    "    sc.pp.log1p(adata)\n",
    "    \n",
    "    sc.pp.pca(adata, random_state=seed, n_comps=30)\n",
    "    sc.pp.neighbors(adata, random_state=seed, n_neighbors=int(len(adata) ** 0.5 // 5), metric='cosine')\n",
    "    tk.tl.triku(adata, use_raw=False)\n",
    "    \n",
    "    sc.tl.umap(adata, min_dist=0.25, random_state=seed)\n",
    "    sc.tl.leiden(adata, resolution=3, random_state=seed)\n",
    "    sc.pl.umap(adata, color=['leiden', 'n_counts'], legend_loc='on data')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "humanitarian-alexander",
   "metadata": {},
   "source": [
    "### Subpopulation characterization (based on paper markers)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b09ed53f-3a92-4e3d-8c3e-6b29f250dcc9",
   "metadata": {
    "scrolled": true,
    "tags": []
   },
   "outputs": [],
   "source": [
    "for adata in [adata_oprescu_d0, adata_oprescu_d05, adata_oprescu_d2, adata_oprescu_d35, adata_oprescu_d5, adata_oprescu_d10, adata_oprescu_d21]:\n",
    "    assign_cats(adata, dict_cats_general, column_groupby='leiden', quantile_gene_sel=0.7, do_return=False, intermediate_states=False, diff=0.05, \n",
    "                key_added='cell_type', min_score=0.45, others_name='unassigned')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e05f08f6-bc81-41c4-bde3-f44d278e3902",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "for adata in [adata_oprescu_d0, adata_oprescu_d05, adata_oprescu_d2, adata_oprescu_d35, adata_oprescu_d5, adata_oprescu_d10, adata_oprescu_d21]:\n",
    "    sc.pl.umap(adata, color=['cell_type'])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5c872d78-e1aa-41ca-8939-b2243cfa9ba5",
   "metadata": {},
   "source": [
    "### Adata saving"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "362618e7-d2ad-4653-aae6-bbd43d5bb0c2",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "for adata, day in zip([adata_oprescu_d0, adata_oprescu_d05, adata_oprescu_d2, adata_oprescu_d35, adata_oprescu_d5, adata_oprescu_d10, adata_oprescu_d21], ['D0', 'D0.5', 'D2', 'D3.5', 'D5', 'D10', 'D21']):\n",
    "    adata.obs['dataset'] = f'Oprescu {day}'\n",
    "    day_save = day.replace('.', '').lower()\n",
    "    adata.write_h5ad(f'data/processed/oprescu_{day_save}.h5')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "material-village",
   "metadata": {},
   "source": [
    "# Scott 2019"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "front-onion",
   "metadata": {},
   "outputs": [],
   "source": [
    "!wget https://ftp.ncbi.nlm.nih.gov/geo/samples/GSM2976nnn/GSM2976778/suppl/GSM2976778%5Fqsnt%5Fbarcodes%2Etsv%2Egz -P {data_dir}/scott\n",
    "!wget https://ftp.ncbi.nlm.nih.gov/geo/samples/GSM2976nnn/GSM2976778/suppl/GSM2976778%5Fqsnt%5Fgenes%2Etsv%2Egz -P {data_dir}/scott\n",
    "!wget https://ftp.ncbi.nlm.nih.gov/geo/samples/GSM2976nnn/GSM2976778/suppl/GSM2976778%5Fqsnt%5Fmatrix%2Emtx%2Egz -P {data_dir}/scott"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "composed-solomon",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "adata_scott_d0 = sc.read_mtx(data_dir+'scott/GSM2976778_qsnt_matrix.mtx.gz').transpose()\n",
    "barcodes = pd.read_csv(data_dir+'scott/GSM2976778_qsnt_barcodes.tsv.gz', sep='\\t', header=None)[0].values\n",
    "features = pd.read_csv(data_dir+'scott/GSM2976778_qsnt_genes.tsv.gz', sep='\\t', header=None)[1].values\n",
    "\n",
    "adata_scott_d0.var_names = features\n",
    "adata_scott_d0.obs_names = barcodes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "racial-webster",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "adata_scott_d0.var_names_make_unique()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "signed-material",
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "!wget https://ftp.ncbi.nlm.nih.gov/geo/samples/GSM2976nnn/GSM2976779/suppl/GSM2976779%5Fd1%5Fbarcodes%2Etsv%2Egz -P {data_dir}/scott\n",
    "!wget https://ftp.ncbi.nlm.nih.gov/geo/samples/GSM2976nnn/GSM2976779/suppl/GSM2976779%5Fd1%5Fgenes%2Etsv%2Egz -P {data_dir}/scott\n",
    "!wget https://ftp.ncbi.nlm.nih.gov/geo/samples/GSM2976nnn/GSM2976779/suppl/GSM2976779%5Fd1%5Fmatrix%2Emtx%2Egz -P {data_dir}/scott"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "conventional-humidity",
   "metadata": {},
   "outputs": [],
   "source": [
    "adata_scott_d1 = sc.read_mtx(data_dir+'scott/GSM2976779_d1_matrix.mtx.gz').transpose()\n",
    "barcodes = pd.read_csv(data_dir+'scott/GSM2976779_d1_barcodes.tsv.gz', sep='\\t', header=None)[0].values\n",
    "features = pd.read_csv(data_dir+'scott/GSM2976779_d1_genes.tsv.gz', sep='\\t', header=None)[1].values\n",
    "\n",
    "adata_scott_d1.var_names = features\n",
    "adata_scott_d1.obs_names = barcodes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "downtown-northwest",
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "!wget https://ftp.ncbi.nlm.nih.gov/geo/samples/GSM2976nnn/GSM2976780/suppl/GSM2976780%5Fd2%5Fbarcodes%2Etsv%2Egz -P {data_dir}/scott\n",
    "!wget https://ftp.ncbi.nlm.nih.gov/geo/samples/GSM2976nnn/GSM2976780/suppl/GSM2976780%5Fd2%5Fgenes%2Etsv%2Egz -P {data_dir}/scott\n",
    "!wget https://ftp.ncbi.nlm.nih.gov/geo/samples/GSM2976nnn/GSM2976780/suppl/GSM2976780%5Fd2%5Fmatrix%2Emtx%2Egz -P {data_dir}/scott"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "packed-cleanup",
   "metadata": {},
   "outputs": [],
   "source": [
    "adata_scott_d2 = sc.read_mtx(data_dir+'scott/GSM2976780_d2_matrix.mtx.gz').transpose()\n",
    "barcodes = pd.read_csv(data_dir+'scott/GSM2976780_d2_barcodes.tsv.gz', sep='\\t', header=None)[0].values\n",
    "features = pd.read_csv(data_dir+'scott/GSM2976780_d2_genes.tsv.gz', sep='\\t', header=None)[1].values\n",
    "\n",
    "adata_scott_d2.var_names = features\n",
    "adata_scott_d2.obs_names = barcodes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "killing-creator",
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "!wget https://ftp.ncbi.nlm.nih.gov/geo/samples/GSM2976nnn/GSM2976781/suppl/GSM2976781%5Fd4%5Fbarcodes%2Etsv%2Egz -P {data_dir}/scott\n",
    "!wget https://ftp.ncbi.nlm.nih.gov/geo/samples/GSM2976nnn/GSM2976781/suppl/GSM2976781%5Fd4%5Fgenes%2Etsv%2Egz -P {data_dir}/scott\n",
    "!wget https://ftp.ncbi.nlm.nih.gov/geo/samples/GSM2976nnn/GSM2976781/suppl/GSM2976781%5Fd4%5Fmatrix%2Emtx%2Egz -P {data_dir}/scott"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "republican-novel",
   "metadata": {},
   "outputs": [],
   "source": [
    "adata_scott_d4 = sc.read_mtx(data_dir+'scott/GSM2976781_d4_matrix.mtx.gz').transpose()\n",
    "barcodes = pd.read_csv(data_dir+'scott/GSM2976781_d4_barcodes.tsv.gz', sep='\\t', header=None)[0].values\n",
    "features = pd.read_csv(data_dir+'scott/GSM2976781_d4_genes.tsv.gz', sep='\\t', header=None)[1].values\n",
    "\n",
    "adata_scott_d4.var_names = features\n",
    "adata_scott_d4.obs_names = barcodes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "sunset-notion",
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "!wget https://ftp.ncbi.nlm.nih.gov/geo/samples/GSM2976nnn/GSM2976782/suppl/GSM2976782%5Fd14%5Fbarcodes%2Etsv%2Egz -P {data_dir}/scott\n",
    "!wget https://ftp.ncbi.nlm.nih.gov/geo/samples/GSM2976nnn/GSM2976782/suppl/GSM2976782%5Fd14%5Fgenes%2Etsv%2Egz -P {data_dir}/scott\n",
    "!wget https://ftp.ncbi.nlm.nih.gov/geo/samples/GSM2976nnn/GSM2976782/suppl/GSM2976782%5Fd14%5Fmatrix%2Emtx%2Egz -P {data_dir}/scott"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "settled-discussion",
   "metadata": {},
   "outputs": [],
   "source": [
    "adata_scott_d14 = sc.read_mtx(data_dir+'scott/GSM2976782_d14_matrix.mtx.gz').transpose()\n",
    "barcodes = pd.read_csv(data_dir+'scott/GSM2976782_d14_barcodes.tsv.gz', sep='\\t', header=None)[0].values\n",
    "features = pd.read_csv(data_dir+'scott/GSM2976782_d14_genes.tsv.gz', sep='\\t', header=None)[1].values\n",
    "\n",
    "adata_scott_d14.var_names = features\n",
    "adata_scott_d14.obs_names = barcodes"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "standard-bahrain",
   "metadata": {},
   "source": [
    "## Adata Scott D0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "musical-strike",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "adata_scott_d0.X = np.asarray(adata_scott_d0.X.todense())\n",
    "adata_scott_d0.raw = adata_scott_d0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "sized-steps",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# Basic QC filtering\n",
    "adata_scott_d0.var['mt'] = adata_scott_d0.var_names.str.startswith('mt-')  # annotate the group of mitochondrial genes as 'mt'\n",
    "sc.pp.calculate_qc_metrics(adata_scott_d0, qc_vars=['mt'], percent_top=None, inplace=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ethical-apple",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "sc.pl.violin(adata_scott_d0, ['n_genes_by_counts', 'total_counts', 'pct_counts_mt'],\n",
    "             jitter=0.4, multi_panel=True)\n",
    "\n",
    "sc.pl.scatter(adata_scott_d0, x='total_counts', y='pct_counts_mt')\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "filled-confusion",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "adata_scott_d0 = adata_scott_d0[adata_scott_d0.obs.n_genes_by_counts < 2500, :]\n",
    "adata_scott_d0 = adata_scott_d0[adata_scott_d0.obs.pct_counts_mt < 10, :]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "polar-works",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "sc.pp.filter_cells(adata_scott_d0, min_genes=250)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "double-infection",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "sc.pp.filter_genes(adata_scott_d0, min_counts=1)\n",
    "sc.pp.normalize_per_cell(adata_scott_d0)\n",
    "sc.pp.log1p(adata_scott_d0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "welcome-schema",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "sc.pp.pca(adata_scott_d0, random_state=seed, n_comps=30)\n",
    "sc.pp.neighbors(adata_scott_d0, random_state=seed, n_neighbors=int(len(adata_scott_d0) ** 0.5 // 4), metric='cosine')\n",
    "tk.tl.triku(adata_scott_d0, use_raw=False)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "negative-piece",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "sc.tl.umap(adata_scott_d0, min_dist=0.7, random_state=seed)\n",
    "sc.tl.leiden(adata_scott_d0, resolution=3, random_state=seed)\n",
    "sc.pl.umap(adata_scott_d0, color=['leiden', 'n_counts'], legend_loc='on data')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cross-commons",
   "metadata": {},
   "source": [
    "### Subpopulation characterization"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "20bc5445-79d4-4c00-a67e-b8383d98ac51",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "assign_cats(adata_scott_d0, dict_cats_general, column_groupby='leiden', quantile_gene_sel=0.7, do_return=False, intermediate_states=False, diff=0.05, \n",
    "                key_added='cell_type', min_score=0.4, others_name='unassigned')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8532dc43-7392-4553-b56d-a89404e04411",
   "metadata": {
    "scrolled": true,
    "tags": []
   },
   "outputs": [],
   "source": [
    "for cat in dict_cats_general.keys():\n",
    "    print(cat)\n",
    "    sc.pl.umap(adata_scott_d0, color=['cell_type'] + [i for i in dict_cats_general[cat] if i in adata_scott_d0.var_names], ncols=3, cmap=magma)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "aa978487-de8c-4c2c-bbce-99be9a562a46",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "sc.pl.umap(adata_scott_d0, color=['leiden'] + ['Bricd5', 'Col9a2', 'Dlk1', 'Grin2b', 'Mpzl2', 'Saa1', 'Shisa3', 'Tenm2'], legend_loc='on data', cmap=magma, ncols=3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "front-fancy",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "sc.pl.umap(adata_scott_d0, color=['Nipal1', 'Trpm6', 'S100b', 'Gpld1', 'Plxdc1', 'Gfra2',\n",
    "                                  'Cd38', 'Cd300lg'], legend_loc='on data', cmap=magma, ncols=3)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0c47aa80-6562-4019-8bae-7f558881ba33",
   "metadata": {},
   "source": [
    "### Adata saving"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a1167f71-fc94-45ec-9907-1d459e1103c1",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "adata_scott_d0.obs['dataset'] = 'Scott D0'\n",
    "adata_scott_d0.write_h5ad('data/processed/scott_d0.h5')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ancient-oregon",
   "metadata": {},
   "source": [
    "# De Micheli et al 2020 (mouse)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "senior-ontario",
   "metadata": {},
   "source": [
    "## Adata download and preprocessing"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "speaking-defensive",
   "metadata": {},
   "outputs": [],
   "source": [
    "!wget https://ftp.ncbi.nlm.nih.gov/geo/series/GSE143nnn/GSE143437/suppl/GSE143437%5FDeMicheli%5FMuSCatlas%5Frawdata%2Etxt%2Egz -P {data_dir}/demicheli_mouse\n",
    "!wget https://ftp.ncbi.nlm.nih.gov/geo/series/GSE143nnn/GSE143437/suppl/GSE143437%5FDeMicheli%5FMuSCatlas%5Fmetadata%2Etxt%2Egz -P {data_dir}/demicheli_mouse"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "freelance-proceeding",
   "metadata": {},
   "outputs": [],
   "source": [
    "adata_de_micheli_mouse = sc.read_text(data_dir+'demicheli_mouse/GSE143437_DeMicheli_MuSCatlas_rawdata.txt.gz').transpose()\n",
    "adata_de_micheli_mouse.obs['batch'] = ['_'.join(i.split('_')[:2]) for i in adata_de_micheli_mouse.obs_names]\n",
    "adata_de_micheli_mouse.obs_names = [i.split('_')[-1] for i in adata_de_micheli_mouse.obs_names]\n",
    "adata_de_micheli_mouse.obs_names_make_unique()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "registered-premises",
   "metadata": {},
   "outputs": [],
   "source": [
    "adata_de_micheli_mouse.write_loom(data_dir+'demicheli_mouse/adata_demicheli_mouse.loom')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "important-darkness",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "adata_de_micheli_mouse = sc.read_loom(data_dir+'demicheli_mouse/adata_demicheli_mouse.loom')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0faf116a-2c24-4c83-a41b-4b03a83c6d87",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "adata_de_micheli_mouse.var.index = adata_de_micheli_mouse.var['var_names'].values"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "tropical-surface",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# We'll remove some RPSs and RPLs\n",
    "not_RPS = [i for i in adata_de_micheli_mouse.var_names if not (('Rps' in i) | ('Rpl' in i))]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "da5e1df7-93e7-4784-a020-e9c058b8242c",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "adata_de_micheli_mouse.var"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "least-indianapolis",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "adata_de_micheli_mouse.raw = adata_de_micheli_mouse\n",
    "adata_de_micheli_mouse = adata_de_micheli_mouse[:, not_RPS]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "standard-sugar",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# Basic QC filtering\n",
    "adata_de_micheli_mouse.var['mt'] = adata_de_micheli_mouse.var_names.str.startswith('mt-')  # annotate the group of mitochondrial genes as 'mt'\n",
    "sc.pp.calculate_qc_metrics(adata_de_micheli_mouse, qc_vars=['mt'], percent_top=None, inplace=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "standard-craps",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "sc.pl.violin(adata_de_micheli_mouse, ['n_genes_by_counts', 'total_counts', 'pct_counts_mt'],\n",
    "             jitter=0.4, multi_panel=True)\n",
    "\n",
    "sc.pl.scatter(adata_de_micheli_mouse, x='n_genes_by_counts', y='pct_counts_mt')\n",
    "sc.pl.scatter(adata_de_micheli_mouse, x='total_counts', y='n_genes_by_counts')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "greater-lobby",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "adata_de_micheli_mouse = adata_de_micheli_mouse[adata_de_micheli_mouse.obs.n_genes_by_counts < 6000, :]\n",
    "adata_de_micheli_mouse = adata_de_micheli_mouse[adata_de_micheli_mouse.obs.total_counts < 30000, :]\n",
    "adata_de_micheli_mouse = adata_de_micheli_mouse[adata_de_micheli_mouse.obs.pct_counts_mt < 10, :]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "attempted-factor",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "sc.pl.violin(adata_de_micheli_mouse, ['n_genes_by_counts', 'total_counts', 'pct_counts_mt'],\n",
    "             jitter=0.4, multi_panel=True)\n",
    "\n",
    "sc.pl.scatter(adata_de_micheli_mouse, x='n_genes_by_counts', y='pct_counts_mt')\n",
    "sc.pl.scatter(adata_de_micheli_mouse, x='total_counts', y='n_genes_by_counts')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "clinical-yorkshire",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "sc.pp.filter_genes(adata_de_micheli_mouse, min_counts=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "eight-disney",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "adata_de_micheli_mouse.X = np.array(adata_de_micheli_mouse.X.todense())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "mobile-furniture",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "set(adata_de_micheli_mouse.obs.batch)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "active-traveler",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "adata_de_micheli_mouse_d0 = adata_de_micheli_mouse[adata_de_micheli_mouse.obs['batch'].isin(['D0_A', 'D0_B', 'D0_Cv3'])].copy()\n",
    "adata_de_micheli_mouse_d2 = adata_de_micheli_mouse[adata_de_micheli_mouse.obs['batch'].isin(['D2_C', 'D2_D'])].copy()\n",
    "adata_de_micheli_mouse_d5 = adata_de_micheli_mouse[adata_de_micheli_mouse.obs['batch'].isin(['D5_A', 'D5_B', 'D5_C'])].copy()\n",
    "adata_de_micheli_mouse_d7 = adata_de_micheli_mouse[adata_de_micheli_mouse.obs['batch'].isin(['D7_C', 'D7_D'])].copy()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "assisted-yesterday",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "for adata in [adata_de_micheli_mouse_d0, adata_de_micheli_mouse_d2, \n",
    "              adata_de_micheli_mouse_d5, adata_de_micheli_mouse_d7]:\n",
    "\n",
    "    sc.pp.filter_genes(adata, min_counts=1)\n",
    "    sc.pp.normalize_total(adata)\n",
    "    sc.pp.log1p(adata)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "subject-terror",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "for adata in [adata_de_micheli_mouse_d0, adata_de_micheli_mouse_d2, \n",
    "              adata_de_micheli_mouse_d5, adata_de_micheli_mouse_d7]:\n",
    "    \n",
    "    sc.pp.pca(adata, random_state=seed)\n",
    "    bbknn.bbknn(adata, metric='angular')\n",
    "    tk.tl.triku(adata, use_raw=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "floating-spirituality",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "for adata in [adata_de_micheli_mouse_d0, adata_de_micheli_mouse_d2, \n",
    "              adata_de_micheli_mouse_d5, adata_de_micheli_mouse_d7]:\n",
    "\n",
    "    sc.tl.umap(adata, min_dist=0.4,  random_state=seed)\n",
    "    sc.tl.leiden(adata, resolution=2.5, random_state=seed)\n",
    "    sc.pl.umap(adata, color=['leiden', 'batch', 'total_counts'], legend_loc='on data', ncols=2)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "needed-pride",
   "metadata": {},
   "source": [
    "## Population characterization"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "strategic-addiction",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "for adata in [adata_de_micheli_mouse_d0, adata_de_micheli_mouse_d2, \n",
    "              adata_de_micheli_mouse_d5, adata_de_micheli_mouse_d7]:\n",
    "    sc.tl.rank_genes_groups(adata, groupby='leiden', method='wilcoxon')\n",
    "    assign_cats(adata, dict_cats_general, column_groupby='leiden', quantile_gene_sel=0.7, do_return=False, intermediate_states=False, diff=0.05, \n",
    "                key_added='cell_type', min_score=0.4, others_name='unassigned')\n",
    "    sc.pl.umap(adata, color=['leiden', 'cell_type',], legend_loc='on data', ncols=2)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ace0278e-b2c7-440e-b31a-eb86477bb753",
   "metadata": {},
   "source": [
    "### Adata saving"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8a2c4749-d40e-4ce1-8b30-9f59f77ba5a5",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "adata_de_micheli_mouse_d0.obs['dataset'] = 'De Micheli mouse D0'\n",
    "adata_de_micheli_mouse_d0.write_h5ad('data/processed/de_micheli_mouse_d0.h5')\n",
    "\n",
    "adata_de_micheli_mouse_d2.obs['dataset'] = 'De Micheli mouse D2'\n",
    "adata_de_micheli_mouse_d2.write_h5ad('data/processed/de_micheli_mouse_d2.h5')\n",
    "\n",
    "adata_de_micheli_mouse_d5.obs['dataset'] = 'De Micheli mouse D5'\n",
    "adata_de_micheli_mouse_d5.write_h5ad('data/processed/de_micheli_mouse_d5.h5')\n",
    "\n",
    "adata_de_micheli_mouse_d7.obs['dataset'] = 'De Micheli mouse D7'\n",
    "adata_de_micheli_mouse_d7.write_h5ad('data/processed/de_micheli_mouse_d7.h5')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "offensive-anaheim",
   "metadata": {},
   "source": [
    "# De Micheli et al 2020 (human)\n",
    "\n",
    "This is a human dataset, consisting of 10 patients. To remove batch effects, we will directly run kallisto in all samples at once, and then run the analysis on them."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "referenced-union",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "de_micheli_dir = os.getcwd() + '/data/demicheli'\n",
    "os.makedirs(de_micheli_dir, exist_ok=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "romance-explanation",
   "metadata": {},
   "outputs": [],
   "source": [
    "SRA_list = \"\"\"\n",
    "SRR10897760\n",
    "SRR10897761\n",
    "SRR10897762\n",
    "SRR10897763\n",
    "SRR10897764\n",
    "SRR10897765\n",
    "SRR10897766\n",
    "SRR10897767\n",
    "SRR10897768\n",
    "SRR10897769\n",
    "\"\"\"\n",
    "\n",
    "with open(de_micheli_dir + '/accession.txt', 'w') as f:\n",
    "    f.write(SRA_list)\n",
    "    \n",
    "df = pd.DataFrame({'name': ['de_micheli'], 'technology': ['10xv3'], 'targetnumcells': [5000]})\n",
    "\n",
    "df.to_csv(de_micheli_dir + '/metadata.tab', sep='\\t', index=None)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "controversial-grill",
   "metadata": {},
   "outputs": [],
   "source": [
    "!fastq-dump SRR10897768 --split-files --gzip  # This is just one. Run from 60 to 69"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "selected-secondary",
   "metadata": {},
   "outputs": [],
   "source": [
    "from time import sleep"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "enabling-explanation",
   "metadata": {},
   "outputs": [],
   "source": [
    "sleep(3600 * 4)\n",
    "for SRA_idx, SRA in enumerate(SRA_list.split('\\n')[1:-1]):\n",
    "    os.rename(de_micheli_dir + f'/{SRA}_2.fastq.gz', \n",
    "              de_micheli_dir + f'/de_micheli_L00{SRA_idx+1}_R1_001.fastq.gz')\n",
    "    os.rename(de_micheli_dir + f'/{SRA}_3.fastq.gz', \n",
    "              de_micheli_dir + f'/de_micheli_L00{SRA_idx+1}_R2_001.fastq.gz')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "acknowledged-portrait",
   "metadata": {},
   "outputs": [],
   "source": [
    "!loompy fromfq {de_micheli_dir}/de_micheli.loom de_micheli /media/seth/SETH_DATA/SETH_Alex/Programs/human_GRCh38_gencode.v31.600 {de_micheli_dir}/metadata.tab \\\n",
    "{de_micheli_dir}/de_micheli_L001_R1_001.fastq.gz {de_micheli_dir}/de_micheli_L001_R2_001.fastq.gz \\\n",
    "{de_micheli_dir}/de_micheli_L002_R1_001.fastq.gz {de_micheli_dir}/de_micheli_L002_R2_001.fastq.gz \\\n",
    "{de_micheli_dir}/de_micheli_L003_R1_001.fastq.gz {de_micheli_dir}/de_micheli_L003_R2_001.fastq.gz \\\n",
    "{de_micheli_dir}/de_micheli_L004_R1_001.fastq.gz {de_micheli_dir}/de_micheli_L004_R2_001.fastq.gz \\\n",
    "{de_micheli_dir}/de_micheli_L005_R1_001.fastq.gz {de_micheli_dir}/de_micheli_L005_R2_001.fastq.gz \\\n",
    "{de_micheli_dir}/de_micheli_L006_R1_001.fastq.gz {de_micheli_dir}/de_micheli_L006_R2_001.fastq.gz \\\n",
    "{de_micheli_dir}/de_micheli_L007_R1_001.fastq.gz {de_micheli_dir}/de_micheli_L007_R2_001.fastq.gz \\\n",
    "{de_micheli_dir}/de_micheli_L008_R1_001.fastq.gz {de_micheli_dir}/de_micheli_L008_R2_001.fastq.gz \\\n",
    "{de_micheli_dir}/de_micheli_L009_R1_001.fastq.gz {de_micheli_dir}/de_micheli_L009_R2_001.fastq.gz \\\n",
    "{de_micheli_dir}/de_micheli_L0010_R1_001.fastq.gz {de_micheli_dir}/de_micheli_L0010_R2_001.fastq.gz "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "appreciated-gentleman",
   "metadata": {},
   "source": [
    "We will also load the matrix file because it contains the batches, and it will be useful for us to \n",
    "label some of the batches and map them later on. This will help us remove possible small spurious clusters."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "serial-chart",
   "metadata": {},
   "outputs": [],
   "source": [
    "link = 'https://ftp.ncbi.nlm.nih.gov/geo/series/GSE143nnn/GSE143704/suppl/GSE143704%5FDeMicheli%5FHumanMuscleAtlas%5Frawdata%2Etxt%2Egz'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "obvious-seafood",
   "metadata": {},
   "outputs": [],
   "source": [
    "!wget {link} -P {data_dir}/demicheli"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "enormous-moment",
   "metadata": {},
   "source": [
    "## Adata load and preprocessing"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "worst-bandwidth",
   "metadata": {},
   "outputs": [],
   "source": [
    "adata_de_micheli_batch = sc.read_text(data_dir+'/demicheli/GSE143704_DeMicheli_HumanMuscleAtlas_rawdata.txt.gz').transpose()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "beginning-bachelor",
   "metadata": {},
   "outputs": [],
   "source": [
    "adata_de_micheli_batch.obs['batch'] = [i.split('_')[1] for i in adata_de_micheli_batch.obs_names]\n",
    "adata_de_micheli_batch.var_names_make_unique()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c8dae16d-5321-4dcc-948d-8624e7d06c79",
   "metadata": {},
   "outputs": [],
   "source": [
    "adata_de_micheli_batch.write_loom(de_micheli_dir+'/de_micheli_human.loom')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "alike-input",
   "metadata": {},
   "outputs": [],
   "source": [
    "adata_de_micheli_human = sc.read_loom(de_micheli_dir + '/de_micheli_human.loom')\n",
    "adata_de_micheli_human.ra"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "least-return",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Basic QC filtering\n",
    "adata_de_micheli_human.var['mt'] = adata_de_micheli_human.var_names.str.startswith('MT-')  # annotate the group of mitochondrial genes as 'mt'\n",
    "sc.pp.calculate_qc_metrics(adata_de_micheli_human, qc_vars=['mt'], percent_top=None, inplace=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "loaded-silver",
   "metadata": {},
   "outputs": [],
   "source": [
    "sc.pl.violin(adata_de_micheli_human, ['n_genes_by_counts', 'total_counts', 'pct_counts_mt'],\n",
    "             jitter=0.4, multi_panel=True)\n",
    "\n",
    "sc.pl.scatter(adata_de_micheli_human, x='n_genes_by_counts', y='pct_counts_mt')\n",
    "sc.pl.scatter(adata_de_micheli_human, x='total_counts', y='n_genes_by_counts')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "dynamic-horizon",
   "metadata": {},
   "outputs": [],
   "source": [
    "adata_de_micheli_human = adata_de_micheli_human[adata_de_micheli_human.obs.n_genes_by_counts < 3500, :]\n",
    "adata_de_micheli_human = adata_de_micheli_human[adata_de_micheli_human.obs.total_counts < 15000, :]\n",
    "adata_de_micheli_human = adata_de_micheli_human[adata_de_micheli_human.obs.pct_counts_mt < 80, :]\n",
    "adata_de_micheli_human = adata_de_micheli_human[~ ((adata_de_micheli_human.obs.pct_counts_mt < 10) & \n",
    "                                       (adata_de_micheli_human.obs.pct_counts_mt < 400)), :]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "distinguished-possible",
   "metadata": {},
   "outputs": [],
   "source": [
    "sc.pl.violin(adata_de_micheli_human, ['n_genes_by_counts', 'total_counts', 'pct_counts_mt'],\n",
    "             jitter=0.4, multi_panel=True)\n",
    "\n",
    "sc.pl.scatter(adata_de_micheli_human, x='n_genes_by_counts', y='pct_counts_mt')\n",
    "sc.pl.scatter(adata_de_micheli_human, x='total_counts', y='n_genes_by_counts')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "exciting-shell",
   "metadata": {},
   "outputs": [],
   "source": [
    "sc.pp.filter_genes(adata_de_micheli_human, min_counts=1)\n",
    "sc.pp.normalize_total(adata_de_micheli_human)\n",
    "sc.pp.log1p(adata_de_micheli_human)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "declared-depth",
   "metadata": {},
   "outputs": [],
   "source": [
    "sc.pp.pca(adata_de_micheli_human, random_state=seed)\n",
    "bbknn.bbknn(adata_de_micheli_human, neighbors_within_batch=3)\n",
    "tk.tl.triku(adata_de_micheli_human, use_raw=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ecological-nursery",
   "metadata": {},
   "outputs": [],
   "source": [
    "sc.tl.umap(adata_de_micheli_human, min_dist=0.3, random_state=seed)\n",
    "sc.tl.leiden(adata_de_micheli_human, resolution=2, random_state=seed)\n",
    "sc.pl.umap(adata_de_micheli_human, color=['leiden', 'batch', 'total_counts'], legend_loc='on data', ncols=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "administrative-coalition",
   "metadata": {},
   "outputs": [],
   "source": [
    "sc.pl.umap(adata_de_micheli_human, color=['SHISA3', 'NIPAL1', 'SAA1', 'COL9A2', 'S100B', 'GFRA2', 'CD300LG'], \n",
    "           legend_loc='on data', ncols=2, cmap=magma)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "862db39e-86c2-4dc6-921a-a2999946ac9c",
   "metadata": {},
   "outputs": [],
   "source": [
    "# ESTO HAY QUE PONERLE UPPER\n",
    "assign_cats(adata_de_micheli_human, dict_cats_general, column_groupby='leiden', quantile_gene_sel=0.7, do_return=False, intermediate_states=False, diff=0.05, \n",
    "                key_added='cell_type', min_score=0.4, others_name='unassigned')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6d68b4d2-3913-445b-a6e2-cfb3ae78df5b",
   "metadata": {},
   "outputs": [],
   "source": [
    "for cat in dict_cats_general.keys():\n",
    "    print(cat)\n",
    "    sc.pl.umap(adata_de_micheli_human, color=['cell_type'] + dict_cats_general[cat], ncols=3, cmap=magma)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "dc4b474c-e94c-4fb2-9414-067b9c755952",
   "metadata": {},
   "source": [
    "### Adata saving"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c6802217-99f0-41a9-be69-1c18e26b5fda",
   "metadata": {},
   "outputs": [],
   "source": [
    "adata_de_micheli_human.obs['dataset'] = 'De Micheli human'\n",
    "adata_de_micheli_human.write_h5ad('data/processed/de_micheli_human.h5')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "structural-audit",
   "metadata": {},
   "source": [
    "# Giordani 2020"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "treated-reservoir",
   "metadata": {},
   "source": [
    "## Adata download and preprocessing"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cognitive-supervisor",
   "metadata": {},
   "outputs": [],
   "source": [
    "giordani_dir = os.getcwd() + '/data/giordani'\n",
    "os.makedirs(giordani_dir, exist_ok=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "meaningful-ethnic",
   "metadata": {},
   "outputs": [],
   "source": [
    "SRA_list = \"\"\"\n",
    "SRR8352705\n",
    "SRR8352706\n",
    "\"\"\"\n",
    "\n",
    "with open(giordani_dir + '/accession.txt', 'w') as f:\n",
    "    f.write(SRA_list)\n",
    "\n",
    "df = pd.DataFrame({'name': ['bamtofastq'], 'technology': ['10xv2'], 'targetnumcells': [5000]})\n",
    "\n",
    "df.to_csv(giordani_dir + '/metadata.tab', sep='\\t', index=None)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "pleased-blond",
   "metadata": {},
   "outputs": [],
   "source": [
    "!cd {giordani_dir} && aria2c -x 16 https://sra-pub-src-1.s3.amazonaws.com/SRR8352706/20180917_uninjured_wt_possorted_genome_bam.bam.1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "sapphire-myanmar",
   "metadata": {},
   "outputs": [],
   "source": [
    "!cd {giordani_dir} && aria2c -x 16 https://sra-pub-src-1.s3.amazonaws.com/SRR8352705/20171018_uninjured_wt_possorted_genome_bam.bam.1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "optional-compiler",
   "metadata": {},
   "outputs": [],
   "source": [
    "program_dir = '/media/seth/SETH_DATA/SETH_Alex/Programs/'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "breeding-return",
   "metadata": {},
   "outputs": [],
   "source": [
    "!cd {program_dir} && bamtofastq-1.3.2 {giordani_dir}/20171018_uninjured_wt_possorted_genome_bam.bam.1 {giordani_dir}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "unique-renaissance",
   "metadata": {},
   "outputs": [],
   "source": [
    "!cd {program_dir} && bamtofastq-1.3.2 {giordani_dir}/20180917_uninjured_wt_possorted_genome_bam.bam.1.1 {giordani_dir}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "alert-insert",
   "metadata": {},
   "outputs": [],
   "source": [
    "bamdata_str_20171018 = ''\n",
    "for ROOT, DIR, FILES in os.walk(giordani_dir + f'/bam_20171018'):\n",
    "    for R1_file in FILES:\n",
    "        if \"_R1_\" in R1_file:\n",
    "            R2_file = R1_file.replace('_R1_', '_R2_')\n",
    "            bamdata_str_20171018 += f\"{ROOT}/{R1_file} {ROOT}/{R2_file} \"\n",
    "\n",
    "bamdata_str_20180917 = ''\n",
    "for ROOT, DIR, FILES in os.walk(giordani_dir + f'/bam_20180917'):\n",
    "    for R1_file in FILES:\n",
    "        if \"_R1_\" in R1_file:\n",
    "            R2_file = R1_file.replace('_R1_', '_R2_')\n",
    "            bamdata_str_20180917 += f\"{ROOT}/{R1_file} {ROOT}/{R2_file} \"\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "casual-tucson",
   "metadata": {},
   "outputs": [],
   "source": [
    "!loompy fromfq {giordani_dir}/giordani_20171018.loom bamtofastq /media/seth/SETH_DATA/SETH_Alex/Programs/mouse_GRCm38_gencode.v31 {giordani_dir}/metadata.tab {bamdata_str_20171018}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "manual-johnson",
   "metadata": {},
   "outputs": [],
   "source": [
    "!loompy fromfq {giordani_dir}/giordani_20180917.loom bamtofastq /media/seth/SETH_DATA/SETH_Alex/Programs/mouse_GRCm38_gencode.v31 {giordani_dir}/metadata.tab {bamdata_str_20180917}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "working-calgary",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "adata_giordani_20180917 = sc.read_loom(data_dir+'/giordani/giordani_20180917.loom')\n",
    "adata_giordani_20180917.var_names_make_unique()\n",
    "adata_giordani_20171018 = sc.read_loom(data_dir+'/giordani/giordani_20171018.loom')\n",
    "adata_giordani_20171018.var_names_make_unique()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "afraid-timber",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "adata_giordani = sc.AnnData.concatenate(adata_giordani_20180917, adata_giordani_20171018)\n",
    "adata_giordani.raw = adata_giordani"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "respiratory-freedom",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# Basic QC filtering\n",
    "adata_giordani.var['mt'] = adata_giordani.var_names.str.startswith('mt-')  # annotate the group of mitochondrial genes as 'mt'\n",
    "sc.pp.calculate_qc_metrics(adata_giordani, qc_vars=['mt'], percent_top=None, inplace=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "governing-rocket",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "len(np.sum(adata_giordani.X.todense(), 1).ravel().tolist()[0])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "powerful-buffalo",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "h = plt.hist(np.sum(adata_giordani.X.todense(), 1).ravel().tolist()[0], bins=100)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "individual-census",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "sc.pl.violin(adata_giordani, ['n_genes_by_counts', 'total_counts', 'pct_counts_mt'],\n",
    "             jitter=0.4, multi_panel=True)\n",
    "\n",
    "sc.pl.scatter(adata_giordani, x='total_counts', y='pct_counts_mt')\n",
    "sc.pl.scatter(adata_giordani, x='total_counts', y='n_genes_by_counts', color='batch')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "natural-richardson",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "sc.pp.filter_cells(adata_giordani, min_genes=250)\n",
    "sc.pp.filter_genes(adata_giordani, min_cells=1)\n",
    "adata_giordani = adata_giordani[adata_giordani.obs['pct_counts_mt'] < 20]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "direct-appreciation",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "sc.pp.filter_genes(adata_giordani, min_cells=1)\n",
    "sc.pp.normalize_per_cell(adata_giordani)\n",
    "sc.pp.log1p(adata_giordani)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "meaning-drill",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "sc.pp.pca(adata_giordani, random_state=seed)\n",
    "bbknn.bbknn(adata_giordani)\n",
    "tk.tl.triku(adata_giordani, use_raw=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "graduate-content",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "sc.tl.umap(adata_giordani, min_dist=0.3, random_state=seed)\n",
    "sc.tl.leiden(adata_giordani, resolution=1.3, random_state=seed)\n",
    "sc.pl.umap(adata_giordani, color=['leiden', 'batch', 'n_genes_by_counts'], legend_loc='on data')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "amateur-fifth",
   "metadata": {},
   "source": [
    "## Population characterization"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b1d5efc8-0995-4503-b74d-41e711e8f3f4",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "assign_cats(adata_giordani, dict_cats_general, column_groupby='leiden', quantile_gene_sel=0.7, do_return=False, intermediate_states=False, diff=0.05, \n",
    "                key_added='cell_type', min_score=0.3, others_name='unassigned')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6eb88156-c8c9-4f3e-8f14-90462980e31f",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "for cat in dict_cats_general.keys():\n",
    "    print(cat)\n",
    "    sc.pl.umap(adata_giordani, color=['cell_type'] + dict_cats_general[cat], ncols=3, cmap=magma)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8f3e4a9f-6ae3-4fdc-a2ed-a8ea9ec467fa",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "sc.pl.umap(adata_giordani, color=['cell_type'] + ['Saa1', 'Shisa3', 'Col9a2'], ncols=3, cmap=magma)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2f75c5e6-c93f-47bf-9d25-8bc155aec511",
   "metadata": {},
   "source": [
    "### Adata saving"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e28ebc31-2e0b-4521-a461-8f1cfe3f4128",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "adata_giordani.obs['dataset'] = 'Giordani'\n",
    "adata_giordani.write_h5ad('data/processed/giordani.h5')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cathedral-medication",
   "metadata": {},
   "source": [
    "# Proietti 2021"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "handy-shell",
   "metadata": {},
   "source": [
    "## Adata download and preprocessing"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "pointed-dressing",
   "metadata": {},
   "outputs": [],
   "source": [
    "proietti_dir = os.getcwd() + '/data/proietti_2021'\n",
    "os.makedirs(proietti_dir, exist_ok=True)\n",
    "\n",
    "mouse_gencode_dir = '/media/seth/SETH_DATA/SETH_Alex/Programs/mouse_GRCm38_gencode.v31'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cultural-ceremony",
   "metadata": {},
   "outputs": [],
   "source": [
    "SRA_list_ITGA_muscle_PRJNA626530_DEN = ['SRR11574458', 'SRR11574462', 'SRR11574463', 'SRR11574464']\n",
    "SRA_list_ITGA_muscle_PRJNA626530_CTRL = ['SRR11574459', 'SRR11574460', 'SRR11574461', 'SRR11574465']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "mature-eating",
   "metadata": {},
   "outputs": [],
   "source": [
    "df = pd.DataFrame({'name': ['PRJNA626530_DEN', 'PRJNA626530_CTRL'], \n",
    "                   'technology': ['10xv3'] * 2, \n",
    "                   'targetnumcells': [6000] * 2})\n",
    "\n",
    "df.to_csv(proietti_dir + '/metadata.tab', sep='\\t', index=None)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "unlikely-berry",
   "metadata": {},
   "outputs": [],
   "source": [
    "for idx, name in enumerate(SRA_list_ITGA_muscle_PRJNA626530_DEN):\n",
    "    os.system(f'cd {proietti_dir} && parallel-fastq-dump -s {name} -t {16} --split-files -O . --tmpdir tmp --gzip')\n",
    "    os.system(f'cd {proietti_dir} && mv {name}_1.fastq.gz PRJNA626530_DEN_L00{idx+1}_R1_001.fastq.gz')\n",
    "    os.system(f'cd {proietti_dir} && mv {name}_2.fastq.gz PRJNA626530_DEN_L00{idx+1}_R2_001.fastq.gz')\n",
    "\n",
    "for idx, name in enumerate(SRA_list_ITGA_muscle_PRJNA626530_CTRL):\n",
    "    os.system(f'cd {proietti_dir} && parallel-fastq-dump -s {name} -t {16} --split-files -O . --tmpdir tmp --gzip')\n",
    "    os.system(f'cd {proietti_dir} && mv {name}_1.fastq.gz PRJNA626530_CTRL_L00{idx+1}_R1_001.fastq.gz')\n",
    "    os.system(f'cd {proietti_dir} && mv {name}_2.fastq.gz PRJNA626530_CTRL_L00{idx+1}_R2_001.fastq.gz')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "moral-vehicle",
   "metadata": {},
   "outputs": [],
   "source": [
    "print(f'''cd {proietti_dir} && loompy fromfq PRJNA626530_DEN.loom PRJNA626530_DEN {mouse_gencode_dir} metadata.tab \n",
    "      PRJNA626530_DEN_L001_R1_001.fastq.gz PRJNA626530_DEN_L001_R2_001.fastq.gz \n",
    "      PRJNA626530_DEN_L002_R1_001.fastq.gz PRJNA626530_DEN_L002_R2_001.fastq.gz \n",
    "      PRJNA626530_DEN_L003_R1_001.fastq.gz PRJNA626530_DEN_L003_R2_001.fastq.gz \n",
    "      PRJNA626530_DEN_L004_R1_001.fastq.gz PRJNA626530_DEN_L004_R2_001.fastq.gz''')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "lined-cowboy",
   "metadata": {},
   "outputs": [],
   "source": [
    "os.system(f'''cd {proietti_dir} && loompy fromfq PRJNA626530_DEN.loom PRJNA626530_DEN {mouse_gencode_dir} metadata.tab \n",
    "      PRJNA626530_DEN_L001_R1_001.fastq.gz PRJNA626530_DEN_L001_R2_001.fastq.gz \n",
    "      PRJNA626530_DEN_L002_R1_001.fastq.gz PRJNA626530_DEN_L002_R2_001.fastq.gz \n",
    "      PRJNA626530_DEN_L003_R1_001.fastq.gz PRJNA626530_DEN_L003_R2_001.fastq.gz \n",
    "      PRJNA626530_DEN_L004_R1_001.fastq.gz PRJNA626530_DEN_L004_R2_001.fastq.gz''')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "behind-analysis",
   "metadata": {},
   "outputs": [],
   "source": [
    "print(f'''cd {proietti_dir} && loompy fromfq PRJNA626530_CTRL.loom PRJNA626530_CTRL {mouse_gencode_dir} metadata.tab \n",
    "      PRJNA626530_CTRL_L001_R1_001.fastq.gz PRJNA626530_CTRL_L001_R2_001.fastq.gz \n",
    "      PRJNA626530_CTRL_L002_R1_001.fastq.gz PRJNA626530_CTRL_L002_R2_001.fastq.gz \n",
    "      PRJNA626530_CTRL_L003_R1_001.fastq.gz PRJNA626530_CTRL_L003_R2_001.fastq.gz \n",
    "      PRJNA626530_CTRL_L004_R1_001.fastq.gz PRJNA626530_CTRL_L004_R2_001.fastq.gz''')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "increased-bristol",
   "metadata": {},
   "outputs": [],
   "source": [
    "os.system(f'''cd {proietti_dir} && loompy fromfq PRJNA626530_CTRL.loom PRJNA626530_CTRL {mouse_gencode_dir} metadata.tab \n",
    "      PRJNA626530_CTRL_L001_R1_001.fastq.gz PRJNA626530_CTRL_L001_R2_001.fastq.gz \n",
    "      PRJNA626530_CTRL_L002_R1_001.fastq.gz PRJNA626530_CTRL_L002_R2_001.fastq.gz \n",
    "      PRJNA626530_CTRL_L003_R1_001.fastq.gz PRJNA626530_CTRL_L003_R2_001.fastq.gz \n",
    "      PRJNA626530_CTRL_L004_R1_001.fastq.gz PRJNA626530_CTRL_L004_R2_001.fastq.gz''')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "informational-rolling",
   "metadata": {},
   "source": [
    "## Preprocess dataset"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "endangered-photography",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "adata_proietti_DEN = sc.read(data_dir+'/proietti_2021/PRJNA626530_DEN.loom')\n",
    "adata_proietti_DEN.var_names_make_unique()\n",
    "\n",
    "adata_proietti_CTRL = sc.read(data_dir+'/proietti_2021/PRJNA626530_CTRL.loom')\n",
    "adata_proietti_CTRL.var_names_make_unique()\n",
    "\n",
    "adata_proietti = sc.AnnData.concatenate(adata_proietti_DEN, adata_proietti_CTRL, batch_key='condition', batch_categories=['DEN', 'CTRL'])\n",
    "adata_proietti.raw = adata_proietti"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "improving-dictionary",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# Basic QC filtering\n",
    "adata_proietti.var['mt'] = adata_proietti.var_names.str.startswith('mt-')  # annotate the group of mitochondrial genes as 'mt'\n",
    "sc.pp.calculate_qc_metrics(adata_proietti, qc_vars=['mt'], percent_top=None, inplace=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "conscious-plumbing",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "sc.pl.violin(adata_proietti, ['n_genes_by_counts', 'total_counts', 'pct_counts_mt'],\n",
    "             jitter=0.4, multi_panel=True)\n",
    "\n",
    "sc.pl.scatter(adata_proietti, x='total_counts', y='pct_counts_mt')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "superior-jones",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "adata_proietti = adata_proietti[adata_proietti.obs['pct_counts_mt'] < 25]\n",
    "adata_proietti = adata_proietti[adata_proietti.obs['n_genes_by_counts'] < 3000]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "thrown-metadata",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "sc.pp.filter_cells(adata_proietti, min_genes=250)\n",
    "sc.pp.filter_genes(adata_proietti, min_cells=5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "indie-progress",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "sc.pp.normalize_per_cell(adata_proietti)\n",
    "sc.pp.log1p(adata_proietti)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "first-mistress",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "sc.pp.pca(adata_proietti, random_state=seed)\n",
    "bbknn.bbknn(adata_proietti, metric='angular', neighbors_within_batch=3, batch_key='condition')\n",
    "tk.tl.triku(adata_proietti, use_raw=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "isolated-islam",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "sc.tl.umap(adata_proietti, min_dist=0.3, random_state=seed)\n",
    "sc.tl.leiden(adata_proietti, resolution=1.5, random_state=seed)\n",
    "sc.pl.umap(adata_proietti, color=['leiden', 'n_counts', 'condition', 'pct_counts_mt'], legend_loc='on data', cmap=magma, ncols=2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3b678b5b-6367-43bd-91b6-1e009953f2da",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "assign_cats(adata_proietti, dict_cats_general, column_groupby='leiden', quantile_gene_sel=0.7, do_return=False, \n",
    "            intermediate_states=False, diff=0.05, \n",
    "                key_added='cell_type', min_score=0.35,  others_name='unassigned')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9500b1c9-2402-4366-82aa-b06bb88af4a7",
   "metadata": {
    "scrolled": true,
    "tags": []
   },
   "outputs": [],
   "source": [
    "for cat in dict_cats_general.keys():\n",
    "    print(cat)\n",
    "    sc.pl.umap(adata_proietti, color=['cell_type'] + [i for i in dict_cats_general[cat] if i in adata_proietti.var_names], ncols=3, cmap=magma)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "324ff9b1-43c8-49c0-bf2a-275e05153c0c",
   "metadata": {},
   "source": [
    "### Adata saving"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fcd86ac6-9ae4-45e7-8ba0-125818e2c82e",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "adata_proietti.obs['dataset'] = 'Proietti'\n",
    "adata_proietti.write_h5ad('data/processed/proietti.h5')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9945f7cf-975a-49f8-ad41-558061c67b6a",
   "metadata": {},
   "source": [
    "# Ronzoni 2021"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ba34a68d-6c49-4ad6-a10f-2b1b1b0e1f20",
   "metadata": {},
   "source": [
    "## Adata download and preprocessing"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5b1632da-57d1-4f81-96c1-e78a1854e111",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "ronzoni_dir = os.getcwd() + '/data/ronzoni'\n",
    "os.makedirs(ronzoni_dir, exist_ok=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "afc29d74-6c4f-4817-9abe-2722f703dd30",
   "metadata": {},
   "outputs": [],
   "source": [
    "!aria2c -x 16 -d {ronzoni_dir} https://ftp.ncbi.nlm.nih.gov/geo/samples/GSM4447nnn/GSM4447776/suppl/GSM4447776%5Fraw%5Fcounts%2Ecsv%2Egz\n",
    "!aria2c -x 16 -d {ronzoni_dir} https://ftp.ncbi.nlm.nih.gov/geo/samples/GSM4447nnn/GSM4447776/suppl/GSM4447776%5Fmetadata%2Ecsv%2Egz"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e7a6b184-211b-4e34-8fb7-04b2d665ba29",
   "metadata": {},
   "source": [
    "## Preprocess dataset"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b33780d8-fd3e-4163-81c5-f19cda8a6c5c",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "gene2symbol = pd.read_csv('data/mouse_gene_to_symbol.txt', sep='\\t')\n",
    "dictgene2symbol = dict(zip(gene2symbol['Gene stable ID'].values, gene2symbol['Gene name'].values))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "283a47a5-2c43-4b46-8b4c-d73d28fddaaa",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "adata_ronzoni = sc.read(ronzoni_dir+'/GSM4447776_raw_counts.csv.gz').transpose()\n",
    "\n",
    "list_names = []\n",
    "for gene in adata_ronzoni.var_names:\n",
    "    try:\n",
    "        list_names.append(dictgene2symbol[gene])\n",
    "    except:\n",
    "        list_names.append(gene)\n",
    "\n",
    "adata_ronzoni.var_names = list_names\n",
    "adata_ronzoni.var_names_make_unique()\n",
    "\n",
    "adata_ronzoni.raw = adata_ronzoni"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "81e919d7-7422-4b05-abd7-58450407f8e5",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# Basic QC filtering\n",
    "adata_ronzoni.var['mt'] = adata_ronzoni.var_names.str.startswith('mt-')  # annotate the group of mitochondrial genes as 'mt'\n",
    "sc.pp.calculate_qc_metrics(adata_ronzoni, qc_vars=['mt'], percent_top=None, inplace=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "872c3883-c954-4aa9-8529-fe74f83a6dbc",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "sc.pl.violin(adata_ronzoni, ['n_genes_by_counts', 'total_counts', 'pct_counts_mt'],\n",
    "             jitter=0.4, multi_panel=True)\n",
    "\n",
    "sc.pl.scatter(adata_ronzoni, x='total_counts', y='pct_counts_mt')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a16e4f06-07b1-4e81-ab91-c64a2be5f084",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "adata_ronzoni = adata_ronzoni[adata_ronzoni.obs['pct_counts_mt'] < 12]\n",
    "adata_ronzoni = adata_ronzoni[adata_ronzoni.obs['n_genes_by_counts'] < 6000]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "98b3d509-5bcf-41b2-8ff4-9517390dfa45",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "sc.pp.filter_cells(adata_ronzoni, min_genes=250)\n",
    "sc.pp.filter_genes(adata_ronzoni, min_cells=5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "dfe744db-c8ee-4a19-9438-deccb20890fe",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "sc.pp.normalize_per_cell(adata_ronzoni)\n",
    "sc.pp.log1p(adata_ronzoni, base=10)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f010023c-b6bb-4ff1-95ac-503a8b27be43",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "sc.pp.pca(adata_ronzoni, random_state=seed)\n",
    "sc.pp.neighbors(adata_ronzoni, random_state=seed)\n",
    "tk.tl.triku(adata_ronzoni, use_raw=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d3571753-442f-404b-89f7-d1381ae10017",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "sc.tl.umap(adata_ronzoni, min_dist=0.7, random_state=seed)\n",
    "sc.tl.leiden(adata_ronzoni, resolution=0.5, random_state=seed)\n",
    "sc.pl.umap(adata_ronzoni, color=['leiden', 'n_counts'], legend_loc='on data', cmap=magma, ncols=2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ce86a2ad-3874-4111-9ae4-15b61aeb6ea4",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "assign_cats(adata_ronzoni, dict_cats_general, column_groupby='leiden', quantile_gene_sel=0.7, do_return=False, intermediate_states=False, diff=0.05, \n",
    "                key_added='cell_type', min_score=0.5,  others_name='unassigned')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e254f4f0-39b7-420a-b1af-ecabafcde8ff",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "sc.pl.umap(adata_ronzoni, color=['cell_type'] +  ['Rgs5', 'Ndufa4l2', 'Pax7', 'Ncam2'], ncols=3, cmap=magma)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8373fe39-c3c1-488b-be23-0125436662a0",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "for cat in dict_cats_general.keys():\n",
    "    print(cat)\n",
    "    sc.pl.umap(adata_ronzoni, color=['cell_type'] + [i for i in dict_cats_general[cat] if i in adata_ronzoni.var_names], ncols=3, cmap=magma)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "af54e502-bc2a-44b3-8790-2da66fe84a66",
   "metadata": {},
   "source": [
    "### Adata saving"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "eda93128-198d-4add-be27-6cf51fb21015",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "adata_ronzoni.obs['dataset'] = 'Ronzoni'\n",
    "adata_ronzoni.write_h5ad('data/processed/ronzoni.h5')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f0b5accc-57c2-4d54-9703-b96c25f91ba2",
   "metadata": {},
   "source": [
    "# Dell'Orso 2019"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "397e5605-c0a9-4bb3-b39c-264e85bc6798",
   "metadata": {},
   "source": [
    "## Adata download and preprocessing"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d9bd144a-3a5f-4447-9bfe-8baff64b106a",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "dellorso_dir = os.getcwd() + '/data/dellorso'\n",
    "os.makedirs(dellorso_dir, exist_ok=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a4cff03d-47a4-4bff-bfcd-aceccc8d7289",
   "metadata": {},
   "outputs": [],
   "source": [
    "!aria2c -x 16 -d {dellorso_dir} https://ftp.ncbi.nlm.nih.gov/geo/samples/GSM3614nnn/GSM3614992/suppl/GSM3614992%5Ftotal%5Fmuscle%5Fwt%5Frep1%5Fbarcodes%2Etsv%2Egz\n",
    "!aria2c -x 16 -d {dellorso_dir} https://ftp.ncbi.nlm.nih.gov/geo/samples/GSM3614nnn/GSM3614992/suppl/GSM3614992%5Ftotal%5Fmuscle%5Fwt%5Frep1%5Fgenes%2Etsv%2Egz\n",
    "!aria2c -x 16 -d {dellorso_dir} https://ftp.ncbi.nlm.nih.gov/geo/samples/GSM3614nnn/GSM3614992/suppl/GSM3614992%5Ftotal%5Fmuscle%5Fwt%5Frep1%5Fmatrix%2Emtx%2Egz\n",
    "\n",
    "!aria2c -x 16 -d {dellorso_dir} https://ftp.ncbi.nlm.nih.gov/geo/samples/GSM3614nnn/GSM3614993/suppl/GSM3614993%5Ftotal%5Fmuscle%5Fwt%5Frep2%5Fbarcodes%2Etsv%2Egz\n",
    "!aria2c -x 16 -d {dellorso_dir} https://ftp.ncbi.nlm.nih.gov/geo/samples/GSM3614nnn/GSM3614993/suppl/GSM3614993%5Ftotal%5Fmuscle%5Fwt%5Frep2%5Fgenes%2Etsv%2Egz\n",
    "!aria2c -x 16 -d {dellorso_dir} https://ftp.ncbi.nlm.nih.gov/geo/samples/GSM3614nnn/GSM3614993/suppl/GSM3614993%5Ftotal%5Fmuscle%5Fwt%5Frep2%5Fmatrix%2Emtx%2Egz"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4310ce99-6b79-44dc-9b99-84af2fac9863",
   "metadata": {},
   "source": [
    "## Preprocess dataset"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1183570b-4e08-4773-9477-f21456f9c88c",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "gene2symbol = pd.read_csv('data/mouse_gene_to_symbol.txt', sep='\\t')\n",
    "dictgene2symbol = dict(zip(gene2symbol['Gene stable ID'].values, gene2symbol['Gene name'].values))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b9ad2cf2-4372-4633-afa4-8284a410f13c",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "adata_dellorso_1 = sc.read(dellorso_dir+'/GSM3614992_total_muscle_wt_rep1_matrix.mtx.gz').transpose()\n",
    "adata_dellorso_1.obs_names = pd.read_csv(dellorso_dir+'/GSM3614992_total_muscle_wt_rep1_barcodes.tsv.gz', header=None)[0].values\n",
    "adata_dellorso_1.var_names = pd.read_csv(dellorso_dir+'/GSM3614992_total_muscle_wt_rep1_genes.tsv.gz', header=None, sep='\\t')[1].values\n",
    "adata_dellorso_1.var_names_make_unique()\n",
    "\n",
    "adata_dellorso_2 = sc.read(dellorso_dir+'/GSM3614993_total_muscle_wt_rep2_matrix.mtx.gz').transpose()\n",
    "adata_dellorso_2.obs_names = pd.read_csv(dellorso_dir+'/GSM3614993_total_muscle_wt_rep2_barcodes.tsv.gz', header=None)[0].values\n",
    "adata_dellorso_2.var_names = pd.read_csv(dellorso_dir+'/GSM3614993_total_muscle_wt_rep2_genes.tsv.gz', header=None, sep='\\t')[1].values\n",
    "adata_dellorso_2.var_names_make_unique()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "70330b39-2dcc-4ec1-8825-3745f586af15",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "adata_dellorso = sc.AnnData.concatenate(adata_dellorso_1, adata_dellorso_2)\n",
    "\n",
    "adata_dellorso.raw = adata_dellorso"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2f885749-297b-40d0-967b-598d7d525015",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# Basic QC filtering\n",
    "adata_dellorso.var['mt'] = adata_dellorso.var_names.str.startswith('mt-')  # annotate the group of mitochondrial genes as 'mt'\n",
    "sc.pp.calculate_qc_metrics(adata_dellorso, qc_vars=['mt'], percent_top=None, inplace=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6bef931b-da74-48f0-88e2-a22a80c77994",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "sc.pl.violin(adata_dellorso, ['n_genes_by_counts', 'total_counts', 'pct_counts_mt'],\n",
    "             jitter=0.4, multi_panel=True)\n",
    "\n",
    "sc.pl.scatter(adata_dellorso, x='total_counts', y='pct_counts_mt')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0a6ee949-d2ff-4da9-8b57-84e34befe7e0",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "adata_dellorso = adata_dellorso[adata_dellorso.obs['pct_counts_mt'] < 2.5]\n",
    "adata_dellorso = adata_dellorso[adata_dellorso.obs['n_genes_by_counts'] < 3000]\n",
    "adata_dellorso = adata_dellorso[adata_dellorso.obs['n_genes_by_counts'] > 500]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "22325b67-a8e4-49ab-a2f8-e394a4ec42cb",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "sc.pp.filter_cells(adata_dellorso, min_genes=250)\n",
    "sc.pp.filter_genes(adata_dellorso, min_cells=5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d5e22b38-242a-4eb0-af81-adfd975f87fc",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "sc.pp.normalize_per_cell(adata_dellorso)\n",
    "sc.pp.log1p(adata_dellorso)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7e15cce0-d510-492f-9c8a-ad7b59d94107",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "sc.pp.pca(adata_dellorso, random_state=seed)\n",
    "bbknn.bbknn(adata_dellorso, neighbors_within_batch=3)\n",
    "tk.tl.triku(adata_dellorso, use_raw=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1427dc28-8968-4688-b99b-110980791ce6",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "sc.tl.umap(adata_dellorso, min_dist=0.7, random_state=seed)\n",
    "sc.tl.leiden(adata_dellorso, resolution=0.8, random_state=seed)\n",
    "sc.pl.umap(adata_dellorso, color=['leiden', 'batch', 'total_counts'], legend_loc='on data', cmap=magma, ncols=2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "aa12c088-5b9a-40a7-af2b-802d3f9f63c2",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "assign_cats(adata_dellorso, dict_cats_general, column_groupby='leiden', quantile_gene_sel=0.7, do_return=False, intermediate_states=False, diff=0.05, \n",
    "                key_added='cell_type', min_score=0.5,  others_name='unassigned')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6ec26429-a875-4fb4-89c3-98df7fc58db3",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "sc.pl.umap(adata_dellorso, color=['cell_type'] + ['Mpz', 'S100b'], ncols=3, cmap=magma)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "495b26ac-f217-4ca1-9569-15602987289e",
   "metadata": {
    "scrolled": true,
    "tags": []
   },
   "outputs": [],
   "source": [
    "for cat in dict_cats_general.keys():\n",
    "    print(cat)\n",
    "    sc.pl.umap(adata_dellorso, color=['cell_type'] + [i for i in dict_cats_general[cat] if i in adata_dellorso.var_names], ncols=3, cmap=magma)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6cb5983a-afe3-446e-8e4d-d630d0bbb1a7",
   "metadata": {},
   "source": [
    "### Adata saving"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "de2fd336-14e3-4c70-8ae3-2db4df06910f",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "adata_dellorso.obs['dataset'] = \"Dell'Orso\"\n",
    "adata_dellorso.write_h5ad('data/processed/dellorso.h5')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1bf3d2be-751f-4ac2-b018-b452083f0369",
   "metadata": {},
   "source": [
    "# Fan et al 2021 (mouse) [it is only endothelial!!!]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f5327a60-ada0-4a68-a474-ac9a2d15768a",
   "metadata": {},
   "outputs": [],
   "source": [
    "fan_dir = os.getcwd() + '/data/fan'\n",
    "os.makedirs(fan_dir, exist_ok=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ff473714-3b80-4018-a23d-19cf612384a8",
   "metadata": {},
   "outputs": [],
   "source": [
    "SRR_dict = {'SRR12769482': 'Sedentary_1', 'SRR12769483': 'Sedentary_2', 'SRR12769484': '2weeks_run_1', 'SRR12769485': '2weeks_run_2'}\n",
    "\n",
    "df = pd.DataFrame({'name': list(SRR_dict.values()), 'technology': ['10xv2'] * len(SRR_dict), 'targetnumcells': [5000] * len(SRR_dict)})\n",
    "df.to_csv(fan_dir + '/metadata.tab', sep='\\t', index=None)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "be39c1fc-6b16-4f79-890b-3567e289e585",
   "metadata": {},
   "outputs": [],
   "source": [
    "for SRR, name in SRR_dict.items():\n",
    "    print(SRR, name)\n",
    "    os.system(f\"cd {fan_dir} && parallel-fastq-dump  -t 32 -s {SRR} --gzip --split-files\")\n",
    "    os.system(f\"loompy fromfq {fan_dir}/{name}.loom {name} /media/seth/SETH_DATA/SETH_Alex/Programs/mouse_GRCm38_gencode.v31 {fan_dir}/metadata.tab {fan_dir}/{SRR}_2.fastq.gz {fan_dir}/{SRR}_3.fastq.gz\") "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "be3224e0-3dba-4d7d-b015-d91e6d644fc5",
   "metadata": {},
   "source": [
    "## Adata load and preprocessing"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3cb835dc-b0be-44ad-9458-96c66b77fa8e",
   "metadata": {},
   "outputs": [],
   "source": [
    "adata_fan_sedentary1 = sc.read(fan_dir + '/Sedentary_1.loom')\n",
    "adata_fan_sedentary1.var_names_make_unique()\n",
    "adata_fan_sedentary2 = sc.read(fan_dir + '/Sedentary_2.loom')\n",
    "adata_fan_sedentary2.var_names_make_unique()\n",
    "\n",
    "\n",
    "adata_fan_sedentary = sc.AnnData.concatenate(adata_fan_sedentary1, adata_fan_sedentary2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d593f770-bd83-454e-a35d-d8d10bc28b58",
   "metadata": {},
   "outputs": [],
   "source": [
    "adata_fan_2weeks_run_1 = sc.read(fan_dir + '/2weeks_run_1.loom')\n",
    "adata_fan_2weeks_run_1.var_names_make_unique()\n",
    "adata_fan_2weeks_run_2 = sc.read(fan_dir + '/2weeks_run_2.loom')\n",
    "adata_fan_2weeks_run_2.var_names_make_unique()\n",
    "\n",
    "\n",
    "adata_fan_2weeks_run = sc.AnnData.concatenate(adata_fan_2weeks_run_1, adata_fan_2weeks_run_2)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "059c34a9-c8a0-48c0-865d-8cf7b5b3d027",
   "metadata": {},
   "source": [
    "### sedentary processing"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a59b8bfd-69d8-4475-9d32-6615a9fdb41f",
   "metadata": {},
   "outputs": [],
   "source": [
    "adata_fan_sedentary.var['mt'] = adata_fan_sedentary.var_names.str.startswith('mt-')  # annotate the group of mitochondrial genes as 'mt'\n",
    "sc.pp.calculate_qc_metrics(adata_fan_sedentary, qc_vars=['mt'], percent_top=None, inplace=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9a99e87c-30cc-443d-ad68-dba3802e5f17",
   "metadata": {},
   "outputs": [],
   "source": [
    "sc.pl.violin(adata_fan_sedentary, ['n_genes_by_counts', 'total_counts', 'pct_counts_mt'],\n",
    "             jitter=0.4, multi_panel=True)\n",
    "\n",
    "sc.pl.scatter(adata_fan_sedentary, x='n_genes_by_counts', y='pct_counts_mt')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4841b889-15fb-4946-83cf-8697d9ebab8b",
   "metadata": {},
   "outputs": [],
   "source": [
    "batches = sorted(list(set(adata_fan_sedentary.obs['batch'].values)))\n",
    "for batch in batches:\n",
    "    counts = adata_fan_sedentary.obs['n_genes_by_counts'].loc[adata_fan_sedentary.obs['batch'] == batch].values\n",
    "    sns.distplot(counts)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f3ff9e56-6912-4940-aa1e-ce85f8314dbc",
   "metadata": {},
   "outputs": [],
   "source": [
    "adata_fan_sedentary = adata_fan_sedentary[((adata_fan_sedentary.obs.n_genes_by_counts < 1500) & \n",
    "                                    (adata_fan_sedentary.obs.n_genes_by_counts > 500)).values, :]\n",
    "adata_fan_sedentary = adata_fan_sedentary[adata_fan_sedentary.obs.pct_counts_mt < 2.5, :]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f49fdd59-382f-4714-80d0-4cc157c11b4c",
   "metadata": {},
   "outputs": [],
   "source": [
    "sc.pp.filter_genes(adata_fan_sedentary, min_counts=1)\n",
    "sc.pp.normalize_total(adata_fan_sedentary)\n",
    "sc.pp.log1p(adata_fan_sedentary)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f92e737f-1ea4-4957-96d3-94ff19ec5ace",
   "metadata": {},
   "outputs": [],
   "source": [
    "sc.pp.pca(adata_fan_sedentary, random_state=seed, n_comps=30)\n",
    "sce.pp.harmony_integrate(adata_fan_sedentary, key='batch', max_iter_harmony=50)\n",
    "sc.pp.neighbors(adata_fan_sedentary, use_rep='X_pca_harmony',  n_neighbors=int(0.5 * len(adata_fan_sedentary) ** 0.5 // 2), metric='cosine')\n",
    "tk.tl.triku(adata_fan_sedentary)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4911278a-9807-4d03-9e9c-50c8b3c5bc66",
   "metadata": {},
   "outputs": [],
   "source": [
    "sc.tl.umap(adata_fan_sedentary, min_dist=0.3, random_state=seed)\n",
    "sc.tl.leiden(adata_fan_sedentary, resolution=2, random_state=seed)\n",
    "sc.pl.umap(adata_fan_sedentary, color=['leiden', 'batch', 'total_counts'], legend_loc='on data', ncols=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b09a509b-d74c-4266-9036-33fa36a77a0b",
   "metadata": {},
   "outputs": [],
   "source": [
    "assign_cats(adata_fan_sedentary, dict_cats_general, column_groupby='leiden', quantile_gene_sel=0.7, do_return=False, intermediate_states=False, diff=0.05, \n",
    "                key_added='cell_type', min_score=0.4, others_name='unassigned')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7668b4e6-d9ef-4473-83b7-6849855c45d9",
   "metadata": {},
   "outputs": [],
   "source": [
    "sc.pl.umap(adata_fan_sedentary, color=['Vcam1', 'Clu', 'Car4', 'Itm2a'], ncols=3, cmap=magma)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fab2a626-85a7-4259-b818-b7d6e3bdad26",
   "metadata": {},
   "outputs": [],
   "source": [
    "for cat in dict_cats_general.keys():\n",
    "    print(cat)\n",
    "    sc.pl.umap(adata_fan_sedentary, color=['cell_type'] + dict_cats_general[cat], ncols=3, cmap=magma)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "70b33159-7c99-45ae-9826-5d346e999012",
   "metadata": {},
   "source": [
    "### 2 weeks processing"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8a05e07d-376d-4ae3-bebd-11bcd0dd7e2d",
   "metadata": {},
   "outputs": [],
   "source": [
    "adata_fan_2weeks_run.var['mt'] = adata_fan_2weeks_run.var_names.str.startswith('mt-')  # annotate the group of mitochondrial genes as 'mt'\n",
    "sc.pp.calculate_qc_metrics(adata_fan_2weeks_run, qc_vars=['mt'], percent_top=None, inplace=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "92d1693f-0bf5-4869-bee3-480bb225c540",
   "metadata": {},
   "outputs": [],
   "source": [
    "sc.pl.violin(adata_fan_2weeks_run, ['n_genes_by_counts', 'total_counts', 'pct_counts_mt'],\n",
    "             jitter=0.4, multi_panel=True)\n",
    "\n",
    "sc.pl.scatter(adata_fan_2weeks_run, x='n_genes_by_counts', y='pct_counts_mt')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fdc68c41-4287-4ae4-bef6-a3e19344e291",
   "metadata": {},
   "outputs": [],
   "source": [
    "batches = sorted(list(set(adata_fan_sedentary.obs['batch'].values)))\n",
    "for batch in batches:\n",
    "    counts = adata_fan_2weeks_run.obs['n_genes_by_counts'].loc[adata_fan_2weeks_run.obs['batch'] == batch].values\n",
    "    sns.distplot(counts)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5fe3d87d-79a4-4867-85cd-42df27a81f88",
   "metadata": {},
   "outputs": [],
   "source": [
    "adata_fan_2weeks_run = adata_fan_2weeks_run[((adata_fan_2weeks_run.obs.n_genes_by_counts < 2000) & \n",
    "                                    (adata_fan_2weeks_run.obs.n_genes_by_counts > 300)).values, :]\n",
    "adata_fan_2weeks_run = adata_fan_2weeks_run[adata_fan_2weeks_run.obs.pct_counts_mt < 4, :]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "46610ca7-c3d7-46f6-9d04-593f91907007",
   "metadata": {},
   "outputs": [],
   "source": [
    "sc.pp.filter_genes(adata_fan_2weeks_run, min_counts=1)\n",
    "sc.pp.normalize_total(adata_fan_2weeks_run)\n",
    "sc.pp.log1p(adata_fan_2weeks_run)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0324315b-967f-4e79-81c9-2f6d9d574b94",
   "metadata": {},
   "outputs": [],
   "source": [
    "sc.pp.pca(adata_fan_2weeks_run, random_state=seed, n_comps=30)\n",
    "sce.pp.harmony_integrate(adata_fan_2weeks_run, key='batch', max_iter_harmony=50)\n",
    "sc.pp.neighbors(adata_fan_2weeks_run, use_rep='X_pca_harmony',  n_neighbors=int(0.5 * len(adata_fan_2weeks_run) ** 0.5 // 2), metric='cosine')\n",
    "tk.tl.triku(adata_fan_2weeks_run)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c4ff3db6-0b4e-4dc7-b517-1cd028d618b5",
   "metadata": {},
   "outputs": [],
   "source": [
    "sc.tl.umap(adata_fan_2weeks_run, min_dist=0.3, random_state=seed)\n",
    "sc.tl.leiden(adata_fan_2weeks_run, resolution=2, random_state=seed)\n",
    "sc.pl.umap(adata_fan_2weeks_run, color=['leiden', 'batch', 'total_counts'], legend_loc='on data', ncols=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fec16055-9ac3-46c9-93d5-b6fc400c410a",
   "metadata": {},
   "outputs": [],
   "source": [
    "assign_cats(adata_fan_2weeks_run, dict_cats_general, column_groupby='leiden', quantile_gene_sel=0.7, do_return=False, intermediate_states=False, diff=0.05, \n",
    "                key_added='cell_type', min_score=0.4, others_name='unassigned')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c7d88a04-72b5-498c-9838-c42136fc73dd",
   "metadata": {},
   "outputs": [],
   "source": [
    "sc.pl.umap(adata_fan_2weeks_run, color=['Vcam1', 'Clu', 'Car4', 'Itm2a'], ncols=3, cmap=magma)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ecc04f18-e542-4b77-ab00-d387b62885a7",
   "metadata": {},
   "outputs": [],
   "source": [
    "for cat in dict_cats_general.keys():\n",
    "    print(cat)\n",
    "    sc.pl.umap(adata_fan_2weeks_run, color=['cell_type'] + dict_cats_general[cat], ncols=3, cmap=magma)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6124b075-6e11-4075-95ab-72c8b8329173",
   "metadata": {},
   "source": [
    "### Adata saving"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "33638363-62cf-476d-8003-48bf167426bf",
   "metadata": {},
   "outputs": [],
   "source": [
    "adata_fan_sedentary.obs['dataset'] = 'Fan - Sedentary'\n",
    "adata_fan_sedentary.write_h5ad('data/processed/fan_sedentary.h5')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c5a1eb91-f0db-4a72-a976-0cfb66c46a71",
   "metadata": {},
   "outputs": [],
   "source": [
    "adata_fan_2weeks_run_2.obs['dataset'] = 'Fan - 2 weeks run'\n",
    "adata_fan_2weeks_run_2.write_h5ad('data/processed/fan_2weeks_run.h5')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e45f265e-0b1a-48ee-9517-9bfd22543273",
   "metadata": {},
   "source": [
    "# De Leinroth et al 2022 (mouse)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "66f897eb-2a9c-422d-9550-7ab26bfe6160",
   "metadata": {},
   "source": [
    "## Adata download and preprocessing"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9ea86ff0-19e0-4475-b9dc-63fedd1be305",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# The data was downloaded as an RDS file. I have the script in R to change it to a csv matrix"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "dd7afc61-8235-405a-a702-f0f774ee8e2c",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "df = pd.read_csv(data_dir+'leinroth/leinroth.csv', sep=',', index_col=0).transpose()\n",
    "df_metadata = pd.read_csv(data_dir+'leinroth/leinroth_metadata.csv', sep=',', index_col=0)\n",
    "df_metadata['seurat_clusters'] = df_metadata['seurat_clusters'].astype(str)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c91731c8-8e32-4c20-9feb-8fe0b03ab955",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "adata_leinroth = sc.AnnData(df)\n",
    "adata_leinroth.obs = df_metadata\n",
    "\n",
    "adata_leinroth.var_names_make_unique()\n",
    "\n",
    "adata_leinroth.raw = adata_leinroth"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "75bbb077-6e89-4c56-85c8-a3d7c076270c",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# Basic QC filtering\n",
    "adata_leinroth.var['mt'] = adata_leinroth.var_names.str.startswith('mt-')  # annotate the group of mitochondrial genes as 'mt'\n",
    "sc.pp.calculate_qc_metrics(adata_leinroth, qc_vars=['mt'], percent_top=None, inplace=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7ef14c72-4105-4107-9cef-9f996507cd82",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "sc.pl.violin(adata_leinroth, ['n_genes_by_counts', 'total_counts', 'pct_counts_mt'],\n",
    "             jitter=0.4, multi_panel=True)\n",
    "\n",
    "sc.pl.scatter(adata_leinroth, x='n_genes_by_counts', y='pct_counts_mt')\n",
    "sc.pl.scatter(adata_leinroth, x='total_counts', y='n_genes_by_counts')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "14060fd5-3c0d-4bfc-9adb-4e80cfb1b4c9",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "adata_leinroth = adata_leinroth[adata_leinroth.obs.n_genes_by_counts < 5500, :]\n",
    "adata_leinroth = adata_leinroth[adata_leinroth.obs.n_genes_by_counts > 2000, :]\n",
    "\n",
    "adata_leinroth = adata_leinroth[adata_leinroth.obs.total_counts < 25000, :]\n",
    "adata_leinroth = adata_leinroth[adata_leinroth.obs.pct_counts_mt < 10, :]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "560fac4d-66ab-4c35-b697-c1c2a1def346",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "sc.pl.violin(adata_leinroth, ['n_genes_by_counts', 'total_counts', 'pct_counts_mt'],\n",
    "             jitter=0.4, multi_panel=True)\n",
    "\n",
    "sc.pl.scatter(adata_leinroth, x='n_genes_by_counts', y='pct_counts_mt')\n",
    "sc.pl.scatter(adata_leinroth, x='total_counts', y='n_genes_by_counts')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b6c7ca32-16ac-42e4-afcc-98e458663a27",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "sc.pp.filter_genes(adata_leinroth, min_counts=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7f2bd48e-3d74-4800-aeb9-c24de2c2d11e",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "sc.pp.filter_genes(adata_leinroth, min_counts=1)\n",
    "sc.pp.normalize_total(adata_leinroth)\n",
    "sc.pp.log1p(adata_leinroth)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6addf95e-25fe-44cd-9d0b-88d46060cf53",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "sc.pp.pca(adata_leinroth, random_state=seed)\n",
    "sc.pp.neighbors(adata_leinroth, metric='cosine')\n",
    "tk.tl.triku(adata_leinroth, use_raw=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "851edc67-d2d2-4ea5-a421-6d74f4d2a5dd",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "sc.tl.umap(adata_leinroth, min_dist=0.4,  random_state=seed)\n",
    "sc.tl.leiden(adata_leinroth, resolution=3.5, random_state=seed)\n",
    "sc.pl.umap(adata_leinroth, color=['leiden', 'seurat_clusters', 'total_counts'], legend_loc='on data', ncols=2)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a70c2d2b-b44e-44c6-a437-3a99becac9ae",
   "metadata": {},
   "source": [
    "## Population characterization"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ca0e385b-99f0-4f1a-93ac-0e9d2e9984ce",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "sc.tl.rank_genes_groups(adata_leinroth, groupby='leiden', method='wilcoxon')\n",
    "assign_cats(adata_leinroth, dict_cats_general, column_groupby='leiden', quantile_gene_sel=0.7, do_return=False, intermediate_states=False, diff=0.05, \n",
    "            key_added='cell_type', min_score=0.4, others_name='unassigned')\n",
    "sc.pl.umap(adata_leinroth, color=['leiden', 'seurat_clusters', 'cell_type',], legend_loc='on data', ncols=2)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e04d4d72-f3df-4d76-93fa-6ed07c95af0e",
   "metadata": {},
   "source": [
    "### Adata saving"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1178f334-d8b1-4c81-a16d-f4484e9ff6ca",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "adata_leinroth.obs['dataset'] = 'Leinroth'\n",
    "adata_leinroth.write_h5ad('data/processed/leinroth.h5')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "66cb6eef-fd6b-4351-8524-8af8026ede42",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python [conda env:alex-base] *",
   "language": "python",
   "name": "conda-env-alex-base-py"
  },
  "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.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}