4_Analysis_of_krano_dataset.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import anndata as ad\n",
    "\n",
    "from cellassign import assign_cats\n",
    "\n",
    "import numpy as np\n",
    "\n",
    "import matplotlib as mpl\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "import pandas as pd\n",
    "\n",
    "import seaborn as sns\n",
    "\n",
    "import scanpy as sc\n",
    "import scanpy.external as sce\n",
    "\n",
    "import scipy.sparse as sparse\n",
    "\n",
    "import triku as tk"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "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])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sc.settings.set_figure_params(dpi=100) \n",
    "seed = 0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Alevin outputs Ensembl IDs, and we will transform those to Gene symbols\n",
    "\n",
    "from pybiomart import Server\n",
    "\n",
    "server = Server(host='http://www.ensembl.org')\n",
    "\n",
    "df = server.marts['ENSEMBL_MART_ENSEMBL'].datasets['mmusculus_gene_ensembl'].query(attributes=['ensembl_gene_id', 'external_gene_name'])\n",
    "dict_ensemble_gene = dict(zip(df['Gene stable ID'], df['Gene name']))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Data loading and QC"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "adata = sc.read('data/ARAUZO_03/20230623/output_nfcore/alevin/mtx_conversions/combined_matrix.h5ad')\n",
    "\n",
    "adata.var['ensemble_ID'] = adata.var_names.copy()\n",
    "\n",
    "adata.var_names = [dict_ensemble_gene[i] if i in dict_ensemble_gene.keys() else 'NA' for i in adata.var_names]\n",
    "adata = adata[:, adata.var_names != 'NA']\n",
    "\n",
    "adata.var['gene_symbol'] = adata.var_names.copy()\n",
    "\n",
    "adata.var_names_make_unique()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "adata.raw = adata\n",
    "adata.X = adata.X.astype(np.float32) # If int64 it cannot be later normalised"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sc.pl.highest_expr_genes(adata, n_top=20, )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sc.pp.filter_cells(adata, min_counts=100) # based on a previous analysis\n",
    "sc.pp.filter_genes(adata, min_counts=20)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "adata.var['mt'] = adata.var_names.str.startswith('mt-')  # annotate the group of mitochondrial genes as 'mt'\n",
    "sc.pp.calculate_qc_metrics(adata, qc_vars=['mt'], percent_top=None, log1p=True, inplace=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sc.pl.violin(adata, ['log1p_n_genes_by_counts', 'log1p_total_counts', 'pct_counts_mt'], jitter=0.4, multi_panel=True)\n",
    "sc.pl.scatter(adata, x='log1p_total_counts', y='pct_counts_mt')\n",
    "sc.pl.scatter(adata, x='log1p_total_counts', y='log1p_n_genes_by_counts')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(1, 1, figsize=(10, 4))\n",
    "df = pd.DataFrame({'x': adata.obs['sample'], 'y': adata.obs['log1p_n_genes_by_counts']})\n",
    "sns.violinplot(x='x', y='y', data=df, ax=ax)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "adata = adata[(((adata.obs['sample'] == 'AZ7845') & (adata.obs.log1p_n_genes_by_counts < 8.5) & \n",
    "                                                (adata.obs.log1p_n_genes_by_counts > 7.2)) | \n",
    "                                              ((adata.obs['sample'] == 'AZ7846') & (adata.obs.log1p_n_genes_by_counts < 8.5) & \n",
    "                                                (adata.obs.log1p_n_genes_by_counts > 7))  \n",
    "                                             ).values, :]\n",
    "adata = adata[adata.obs.pct_counts_mt < 12, :]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Doublet flagging"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sce.pp.scrublet(adata, batch_key=None, expected_doublet_rate=0.05, random_state=seed, knn_dist_metric='cosine', log_transform=False)\n",
    "\n",
    "print(f\"Number of doublets: {adata.obs['predicted_doublet'].sum()}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(f'N cells in AZ7845: {len(adata[adata.obs[\"sample\"] == \"AZ7845\"])}')\n",
    "print(f'N cells in AZ7846: {len(adata[adata.obs[\"sample\"] == \"AZ7846\"])}')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Data processing"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "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,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Feature selection\n",
    "sc.pp.pca(adata, random_state=seed)\n",
    "sc.pp.neighbors(adata, n_neighbors=int(0.5 * len(adata) ** 0.5), random_state=seed, metric='correlation')\n",
    "tk.tl.triku(adata, use_raw=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sc.pp.pca(adata, random_state=seed, use_highly_variable=True)\n",
    "sc.pp.neighbors(adata, n_neighbors=int(0.5 * len(adata) ** 0.5), random_state=seed, metric='correlation')\n",
    "\n",
    "sc.tl.umap(adata)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sc.tl.leiden(adata, resolution=0.07, key_added='leiden')\n",
    "sc.tl.leiden(adata, resolution=1, key_added='leiden_sub')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sc.pl.umap(adata, color=['leiden', 'leiden_sub', 'sample'], ncols=3, alpha=0.4, legend_loc='on data')\n",
    "sc.pl.umap(adata, color=['log1p_n_genes_by_counts', 'pct_counts_mt'])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Kranocyte characterisation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "A_markers = ['Smim41', 'Col9a2', 'Dlk1', 'Shisa3',  'Saa1',  'Nipal1']\n",
    "B_markers = ['Lypd2', 'Wnt6', 'Cldn1', 'Moxd1', 'Mansc4', 'Dleu7', 'Efnb3', 'Stra6', 'Sbspon',\n",
    "              'Hcn4', 'Cldn22']  \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sc.pl.umap(adata, color=A_markers, cmap=magma, use_raw=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sc.pl.umap(adata, color=B_markers, cmap=magma, use_raw=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# The FACs was done as PDPN(+) CD31(-)\n",
    "\n",
    "sc.pl.umap(adata, color=['Pdpn', 'Pecam1', 'Pdgfra', 'Tnmd', 'Lum', 'Prg4'], cmap=magma, use_raw=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Analysis of major populations"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "krano_12 = ['Nr2f2', 'Gfra1', 'Vwa1', 'Cdh19', 'Rasgrp2', 'Sipa1l1', 'Sorcs2', 'Smim41', 'Ptgfr', 'Plxdc1', 'Nrp2', 'Septin9', \n",
    "            'Gng2', 'Ngfr', 'Malt1', 'Col9a2', 'Plxnc1', 'Rcsd1', 'Nipal1', 'Greb1', 'Pla2g7', 'Sox9', 'Matn2', 'Ogfrl1', \n",
    "            'Gfra2', 'Cp', 'Vwa1', 'P2ry1', 'Hlf', 'Sdc3', 'Mylk', 'Trf', 'Gli1', 'Cst6', 'Pcp4l1', 'Malt1', 'Kcnk2', \n",
    "            'Greb1', 'Nipal1']\n",
    "\n",
    "FAP = ['Cd248', 'Dpt', 'Dpep1', 'Adamts5', 'Emilin2', 'Pcsk6', 'Ifi205', 'Col6a6', 'Procr', 'Osr1', 'Ifi207', \n",
    "                        'Hspb8', 'Wnt11', 'Tmeff2']\n",
    "\n",
    "tnmd = ['Col12a1', 'Cilp2', 'Fmod', 'Col11a1', 'Chad', 'Col8a2', 'Comp', 'Tnmd', 'Pdgfa', 'Actn1', 'Tns3', 'Scx', 'Gas2', \n",
    "        'Ccn4', 'Col11a2', 'Kctd1', 'Myo1d', 'Ptpn14', 'Runx1', 'Tnc', 'Ptx4', 'Col13a1', 'Zfp185', 'Gnai1', 'Cdh2', 'P3h2', \n",
    "        'Emb', 'Atp6v0a4', 'Rflnb', 'Kcnma1', 'Edil3', 'Ppfibp2']\n",
    "\n",
    "satellite = ['Chodl', 'Pax7', 'Edn3', 'Notch3', 'Cavin2', 'Peg3', 'Fgfr4', 'Myf5', 'Cdh15', 'Chrdl2', \n",
    "                          'Heyl', 'Prox1', 'Gal', 'Tanc2', 'Des', 'Rapsn', 'Ank3', 'Jsrp1', 'Notch1', 'Cd82', 'Hs6st2', \n",
    "                          'Traf3ip3']\n",
    "\n",
    "dict_markers = {'FAP_1': FAP + ['Cd248', 'Efemp1', 'Efhd1', 'Dpp4', 'Sema3c', 'Pi16', 'Emilin2', 'Anxa3', 'Pla1a', 'Pcolce2', \n",
    "                                'Gfpt2', 'Procr', 'Cadm3', 'Stmn4', 'Limch1', 'Il18', 'Cd55', 'Cmah', 'Lurap1l', 'Mustn1', \n",
    "                                'Arl4d', 'Car8', 'Aif1l', 'Pde8a', 'Uchl1', 'Smpd3', 'Dact2', 'Rorb', 'Sytl2', 'Chst1', 'Tmem158', \n",
    "                                'Fam167a', 'Tubb4a', 'Sbsn', 'Islr2', 'Ackr2', 'Adamts16', 'Duoxa1', 'Aldh1a3', 'Gap43', 'Kcnk5'],\n",
    "                'FAP_2': FAP + ['Smoc2', 'Col4a1', 'Col15a1', 'Lamb1', 'Hsd11b1', 'Mme', 'Col4a2', 'Vwa1', 'Cxcl14', 'Atp1a2',\n",
    "                                'Lifr', 'Rgma', 'Ret', 'Fbln7', 'Gpm6b', 'Cyria', 'P2ry1', 'Hlf', 'Cldn15', 'Colgalt2', \n",
    "                                'Sorl1', 'Rem1', 'Vtn', 'Prex2', 'Alpl', 'Ces1d', 'Tmem64', 'Cbfa2t3', 'Clec14a', \n",
    "                                'P2ry14', 'Tspan9', 'Fscn1', 'Spon1', 'Csgalnact1', 'Plxna2', 'Kcna2', 'Adcyap1r1', 'Hunk'],\n",
    "                'FAP_3': FAP + ['Srpx', 'Fbln1', 'Gdf10', 'Igfbp7', 'Gpnmb', 'C2', 'Sfrp1', 'C7', 'Dkk2', 'Serpina3n', \n",
    "                                'Armh4', 'Clec11a', 'Inmt', 'C4b', 'Lox', 'Chrdl1', 'Ism1', 'Emb', 'Abcc9', 'Grem2'],\n",
    "                'FAP_4': FAP + ['Mgp', 'Meox2', 'Meox1', 'Hmcn1', 'Clec1a', 'Serpine2', 'Prg4', 'Kctd12', 'Clu', 'Cdh11', \n",
    "                                'Ap1s2', 'Daam2', 'Tmem204', 'Ezr', 'Plekha6', 'Etl4', 'Myo10', 'Ptgir', 'Tenm3', 'Cgnl1', \n",
    "                                'Etv1', 'Myo1b', 'Arhgdib', 'Trib2', 'Entrep1', 'Piezo2', 'Tbx2', 'Ptn', 'Enpp1', 'Col8a2', \n",
    "                                'Rspo3', 'Klf5', 'Mettl24', 'Tmem176a', 'Ddit4l', 'Stmnd1', 'Dchs2', 'Fam180a'], \n",
    "                'Krano_1': krano_12 + ['Cpe', 'Slc1a3', 'Tec', 'Tenm2', 'Piezo2', 'Kif21a', 'Foxd1', 'Mamdc2', \n",
    "                                       'Foxd2os', 'Unc13c', 'Cldn1', \n",
    "                                       'Col28a1', # Found in terminally differentiated Schwann and Merkel cells\n",
    "                                       'Kcnk2', 'Klf15', 'Prxl2a', 'Spp1', 'Trabd2b', 'Pilra', 'Lrp4', 'Slc27a1', \n",
    "                                       'Foxd2', 'Foxs1', 'Grin2b', \n",
    "                                       'Shisa2', 'Sox8', 'Homer2',  # Positive in Tnmc+ population\n",
    "                                       'Acot1', 'Clic6', \n",
    "                                       'Col26a1', 'Cdkn2b', # Partially in Krano_2\n",
    "                                       'Ndnf', 'Asns', 'Gjb5', 'Celf4', 'Gdf11', 'Sfxn5', 'Cerox1', \n",
    "                                       'Rab20', 'Sbspon', 'Lypd6', 'Ccl9', 'Rida', \n",
    "                                       'Ifitm1'], \n",
    "                'Krano_2': krano_12 + ['Sorl1', 'Spon1', 'Alpl', 'Gpld1', 'Rgs17', 'Psat1', 'Trpm6', 'Nkd1', 'Sphkap', \n",
    "                                       'Dlk1', 'Thrsp', 'Susd5', 'Fetub', 'Nrk', 'Saa1', 'Itga8', 'Usp35', 'Lbh', \n",
    "                                       'Grm8', 'Emilin3', 'Necab1', 'Gria1', 'Il12a', 'Col6a6'],\n",
    "                'TNMD_1': tnmd + ['Sparcl1', 'Col22a1', 'Chodl', 'Rbp1', 'Naalad2', 'Cdh2', 'Col18a1', 'Tnfrsf21', 'Scube2', \n",
    "                                  'Postn', 'Rab31', 'Metrnl', 'Plcb4', 'Crispld1', 'Lrrn2', 'Cdh5', 'Tubb2b', \n",
    "                                  'Reln', 'Fbn2', 'Cotl1', 'Olfml2a', 'Pla2g7', 'Cdk5rap2', 'Shisal2b', 'Man2a2', 'Adam23', \n",
    "                                  'Prph', 'Gab2', 'Prex2', 'Frem1', 'Tyms', 'Rapgef4', 'Epas1', 'Atp1a2', 'Egflam'],\n",
    "                'TNMD_2': tnmd + ['Fmod', 'Serping1', 'Ccdc3', 'Fxyd6', 'Cav1', 'Sema3b', 'Angptl7', 'Mylk', 'C3', 'Egfl6', \n",
    "                                  'Wif1', 'Pcolce2', 'Grem2', 'Sema3a', 'Ntn1', 'Cgref1', 'Itga2', 'Ggta1', 'Bmp3', 'Hpgd', \n",
    "                                  'Uts2r', 'Fgf9', 'Cav2', 'Bmpr1b', 'Chrdl1', 'Fndc5', 'Sned1', 'Thy1', 'Il33', 'Fez1', \n",
    "                                  'Loxl4'],\n",
    "                'TNMD_3': tnmd + ['Itm2a', 'Gpx3', 'Bicc1', 'Dlx5', 'Fgl2', 'Gpc1', 'Cemip2', 'Ptn', 'Cdkn1c', 'Wnt16', \n",
    "                                  'Tmem100', 'Flrt2', 'Crabp2', 'Fgfr2', 'Slc20a2', 'Epha3', 'Cdc42ep3', 'Gfra1', 'S100b', \n",
    "                                  'Ostn', 'Shisa2', 'Megf6', 'H19', 'Prss23', 'Slc1a3', 'Spon1', 'Nr4a2', 'Spon1', 'Gdf10', \n",
    "                                  'Pthlh', 'Igsf3', 'Kcns1', 'Igf2', 'Ror2', 'Aqp1', 'Serpine2', 'Nppc'],\n",
    "                'Sat Cav1+': satellite + ['Cav1', 'Kitl', 'Cd36', 'Emcn', 'Itga1', 'Sparcl1', 'Cxcl12', 'Myct1', 'Tspan13', 'Fabp4', \n",
    "                              'Flt1', 'Kdr', 'Cdh5', 'Esam', 'Tie1', 'S1pr1', 'Cped1', 'Adgrf5', 'Apbb2', 'Nrp1', 'Rasgrp3', \n",
    "                              'Adgrl4', 'Elk3', 'Slfn5', 'Epas1', 'Tns1', 'Clic5', 'Pecam1', 'Cd93', 'Egfl7', 'Car8', 'F11r', \n",
    "                              'Ptprm', 'Gimap6', 'Tmem204', 'Ptprb', 'Sult1a1', 'Rgcc', 'Arhgef15', 'Nos3', 'Meox2', 'Fgd5', \n",
    "                              'Aqp7', 'Ushbp1', 'Sox18', 'Arhgap31'],\n",
    "                'Sat Pax7+': satellite + ['Chrdl2', 'Myf5', 'Pax7', 'Edn3', 'Olfml2a', 'Heyl', 'Erfe', 'Fgfr4', 'Msc', \n",
    "                              'Hs6st2', 'Rapsn', 'Tac4', 'Tenm4', 'Serinc2', 'Gm13703', 'Golm2', 'Megf10', 'Nppc', 'Pde1c', \n",
    "                              'Flnc', 'Iqgap2', 'Slc7a2'],\n",
    "                'Sat fibro': satellite + ['Rnase4', 'S100a6', 'Dcn', 'Col1a2', 'Clec3b', 'Gsn', 'Serpinf1', 'C1s1', 'Col6a2', \n",
    "                             'Mfap5', 'Col5a2', 'Col14a1', 'Ifi207', 'Pcolce', 'Lum', 'Igfbp6', 'Axl', 'Myoc', 'Adamts2', 'Ecm1', \n",
    "                             'Sulf2', 'Lgi2', 'Pdlim2', 'Cd248', 'C3', 'Slc1a5', 'Dpep1', 'Emilin2', 'Dpysl3'],\n",
    "                    \n",
    "                }\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for key, val in dict_markers.items():\n",
    "    print(key, [i for i in val if i not in adata.var_names])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sc.tl.rank_genes_groups(adata, groupby='leiden')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sc.pl.umap(adata, color=adata.uns['rank_genes_groups']['names']['3'][0:150], cmap=magma, use_raw=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Analysing Tnmd+ populations"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "adata_tnmd = adata[adata.obs['leiden'] == '2']\n",
    "\n",
    "sc.pp.pca(adata_tnmd, random_state=seed, use_highly_variable=True)\n",
    "sc.pp.neighbors(adata_tnmd, n_neighbors=int(0.5 * len(adata_tnmd) ** 0.5), random_state=seed, metric='correlation')\n",
    "\n",
    "sc.tl.umap(adata_tnmd)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sc.tl.leiden(adata_tnmd, resolution=0.2, key_added='leiden_tnmd')\n",
    "sc.pl.umap(adata_tnmd, color=['leiden_tnmd'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sc.tl.rank_genes_groups(adata_tnmd, groupby='leiden_tnmd')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sc.pl.umap(adata_tnmd, color=adata_tnmd.uns['rank_genes_groups']['names']['2'][0:170], cmap=magma, use_raw=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Analysing FAP populations"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "adata_FAP = adata[adata.obs['leiden'] == '0']\n",
    "\n",
    "sc.pp.pca(adata_FAP, random_state=seed, use_highly_variable=True)\n",
    "sc.pp.neighbors(adata_FAP, n_neighbors=int(0.5 * len(adata_FAP) ** 0.5), random_state=seed, metric='correlation')\n",
    "\n",
    "sc.tl.umap(adata_FAP)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sc.tl.leiden(adata_FAP, resolution=0.3, key_added='leiden_fap')\n",
    "sc.pl.umap(adata_FAP, color=['leiden_fap'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sc.tl.rank_genes_groups(adata_FAP, groupby='leiden_fap')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sc.pl.umap(adata_FAP, color=adata_FAP.uns['rank_genes_groups']['names']['3'][0:170], cmap=magma, use_raw=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Analysing Satellite populations\n",
    "\n",
    "Part of the characterisation is from https://elifesciences.org/articles/51576"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "adata_SAT = adata[adata.obs['leiden'] == '1']\n",
    "\n",
    "sc.pp.pca(adata_SAT, random_state=seed, use_highly_variable=True)\n",
    "sc.pp.neighbors(adata_SAT, n_neighbors=int(0.5 * len(adata_SAT) ** 0.5), random_state=seed, metric='correlation')\n",
    "\n",
    "sc.tl.umap(adata_SAT)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sc.tl.leiden(adata_SAT, resolution=0.2, key_added='leiden_sat')\n",
    "sc.pl.umap(adata_SAT, color=['leiden_sat'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sc.tl.rank_genes_groups(adata_SAT, groupby='leiden_sat')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sc.pl.umap(adata_SAT, color=adata_SAT.uns['rank_genes_groups']['names']['0'][0:170], cmap=magma, use_raw=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Applying cellasign"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sc.tl.leiden(adata, resolution=8, key_added='leiden_assigncats')\n",
    "sc.pl.umap(adata, color=['leiden_assigncats'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "assign_cats(adata, dict_cats=dict_markers, column_groupby='leiden_assigncats', quantile_gene_sel=0.9, diff=0.05, )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sc.pl.umap(adata, color=['assigned_cats', 'Pdgfra'], cmap=magma, use_raw=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Plotting markers of populations"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### FAPs (general)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sc.pl.umap(adata, color=FAP, cmap=magma, use_raw=False, ncols=5)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### FAPs (1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sc.pl.umap(adata, color=['assigned_cats'] + dict_markers['FAP_1'][len(FAP):], cmap=magma, use_raw=False, ncols=5)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### FAPs (2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sc.pl.umap(adata, color=['assigned_cats'] + dict_markers['FAP_2'][len(FAP):], cmap=magma, use_raw=False, ncols=5)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### FAPs (3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sc.pl.umap(adata, color=['assigned_cats'] + dict_markers['FAP_3'][len(FAP):], cmap=magma, use_raw=False, ncols=5)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### FAPs (4)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sc.pl.umap(adata, color=['assigned_cats'] + dict_markers['FAP_4'][len(FAP):], cmap=magma, use_raw=False, ncols=5)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Krano (general)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sc.pl.umap(adata, color=['assigned_cats'] + krano_12, cmap=magma, use_raw=False, ncols=5)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Krano (1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sc.pl.umap(adata, color=['assigned_cats'] + dict_markers['Krano_1'][len(krano_12):], cmap=magma, use_raw=False, ncols=5)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Krano (2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sc.pl.umap(adata, color=['assigned_cats'] + dict_markers['Krano_2'][len(krano_12):], cmap=magma, use_raw=False, ncols=5)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Satellite (general)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sc.pl.umap(adata, color=['assigned_cats'] + satellite, cmap=magma, use_raw=False, ncols=5)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Satellite (1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sc.pl.umap(adata, color=['assigned_cats'] + dict_markers['Sat Cav1+'][len(satellite):], cmap=magma, use_raw=False, ncols=5)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Satellite (2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sc.pl.umap(adata, color=['assigned_cats'] + dict_markers['Sat Pax7+'][len(satellite):], cmap=magma, use_raw=False, ncols=5)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### TNMD (general)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sc.pl.umap(adata, color=['assigned_cats'] + tnmd, cmap=magma, use_raw=False, ncols=5)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### TNMD (1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sc.pl.umap(adata, color=['assigned_cats'] + dict_markers['TNMD_1'][len(tnmd):], cmap=magma, use_raw=False, ncols=5)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### TNMD (2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sc.pl.umap(adata, color=['assigned_cats'] + dict_markers['TNMD_2'][len(tnmd):], cmap=magma, use_raw=False, ncols=5)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### TNMD (3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sc.pl.umap(adata, color=['assigned_cats'] + dict_markers['TNMD_3'][len(tnmd):], cmap=magma, use_raw=False, ncols=5)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Preparing adata for cellxgene"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "adata_cellxgene = adata.copy()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "adata_cellxgene.uns['schema_version'] = '3.0.0'\n",
    "adata_cellxgene.uns['title'] = 'Analysis of muscle fibroblast heterogeneity'\n",
    "adata_cellxgene.uns['batch_condition'] = ['sample']\n",
    "adata_cellxgene.uns['default_embedding'] = 'X_umap'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "columns_in = ['scrublet', 'log1p', 'pca', 'neighbors', 'triku_params', 'umap', 'cell_assign', 'schema_version', \n",
    "'title', 'batch_condition', 'default_embedding']\n",
    "columns_out = [i for i in adata_cellxgene.uns.keys() if i not in columns_in ]\n",
    "\n",
    "print('Columns out: ', columns_out)\n",
    "\n",
    "for i in columns_out:\n",
    "    del adata_cellxgene.uns[i]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "adata_cellxgene.obs['organism_ontology_term_id'] = 'NCBITaxon:10090'\n",
    "adata_cellxgene.obs['tissue_ontology_term_id'] = 'UBERON:0001134'\n",
    "adata_cellxgene.obs['assay_ontology_term_id'] = 'EFO:0009922' # 10x 3' v3\n",
    "adata_cellxgene.obs['disease_ontology_term_id'] = 'PATO:0000461'\n",
    "adata_cellxgene.obs['cell_type_ontology_term_id'] = 'CL:1001609'\n",
    "adata_cellxgene.obs['self_reported_ethnicity_ontology_term_id'] = 'na'\n",
    "adata_cellxgene.obs['development_stage_ontology_term_id'] = 'MmusDv:0000063'\n",
    "adata_cellxgene.obs['sex_ontology_term_id'] = 'PATO:0000384'\n",
    "adata_cellxgene.obs['donor_id'] = adata_cellxgene.obs['sample']\n",
    "adata_cellxgene.obs['suspension_type'] = 'cell'\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "adata_cellxgene.obs = adata_cellxgene.obs.rename(columns={'assigned_cats': 'author_cell_type', 'leiden_assignedcats': 'leiden'})\n",
    "\n",
    "accepted_cats = ['sample', 'log1p_n_genes_by_counts', 'log1p_total_counts', 'pct_counts_mt',\n",
    "                'author_cell_type', 'organism_ontology_term_id', \n",
    "                'tissue_ontology_term_id', 'assay_ontology_term_id', 'disease_ontology_term_id', 'cell_type_ontology_term_id', \n",
    "                'self_reported_ethnicity_ontology_term_id', 'development_stage_ontology_term_id', 'sex_ontology_term_id', 'donor_id', \n",
    "                'suspension_type']\n",
    "  \n",
    "adata_cellxgene.obs = adata_cellxgene.obs[accepted_cats]\n",
    "\n",
    "for cat in accepted_cats:\n",
    "    if cat not in ['log1p_n_genes_by_counts', 'log1p_total_counts', 'pct_counts_mt']:\n",
    "        adata_cellxgene.obs[cat] = adata_cellxgene.obs[cat].astype('category')\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for var in [\"highly_variable\", \"triku_highly_variable\", \"mt\"]:\n",
    "    del adata_cellxgene.var[var]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "del adata_cellxgene.obsm['X_triku']\n",
    "del adata_cellxgene.varm\n",
    "del adata_cellxgene.obsp"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "adata_cellxgene.var_names = adata.var['ensemble_ID']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "adata_cellxgene.X = adata_cellxgene.X.astype(np.float32)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "adata_cellxgene.raw.var.set_index('ensemble_ID', inplace=True)\n",
    "adata_cellxgene.raw.var.drop(columns=['gene_symbol'], inplace=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "genes_add = [e for e in adata_cellxgene.raw.var.index if e not in adata_cellxgene.var.index]\n",
    "\n",
    "new_matrix = sparse.csr_matrix((adata_cellxgene.X.data, adata_cellxgene.X.indices, adata_cellxgene.X.indptr), shape = adata_cellxgene.raw.shape)\n",
    "\n",
    "all_genes = adata_cellxgene.var.index.to_list()\n",
    "all_genes.extend(genes_add)\n",
    "\n",
    "new_var = pd.DataFrame(index=all_genes)\n",
    "new_var = pd.merge(new_var, adata_cellxgene.var, left_index=True, right_index=True, how='left')\n",
    "new_var.loc[genes_add, 'feature_is_filtered'] = True\n",
    "\n",
    "new_adata = ad.AnnData(X=new_matrix, dtype=new_matrix.dtype, \n",
    "                        obs=adata_cellxgene.obs, var=new_var, uns=adata_cellxgene.uns, obsm=adata_cellxgene.obsm, raw = adata_cellxgene.raw)\n",
    "new_adata = new_adata[:,adata_cellxgene.raw.var.index.to_list()]\n",
    "new_adata.var.loc[adata_cellxgene.var.index, 'feature_is_filtered'] = False\n",
    "new_adata.var['feature_is_filtered'] = new_adata.var['feature_is_filtered'].astype(bool)\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "adata_cellxgene = new_adata"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "adata_cellxgene.write_h5ad('data/ARAUZO_03/20230623/output_nfcore/alevin/mtx_conversions/combined_matrix_cellxgene.h5ad', \n",
    "                 compression='gzip')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Exporting HTML"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!jupyter nbconvert --to html /data/Proyectos/kranocito/4_Analysis_of_krano_dataset.ipynb"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "single-cell",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.11.4"
  },
  "orig_nbformat": 4
 },
 "nbformat": 4,
 "nbformat_minor": 2
}