MouseCellAtlas/combined_analyses.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "using Pandas\n",
    "using DataFrames\n",
    "using CSV\n",
    "using PyCall\n",
    "using PyPlot\n",
    "using LsqFit\n",
    "using Printf\n",
    "using KernelDensity\n",
    "using Distributions\n",
    "\n",
    "include(\"tacos.jl\") # Tools to Analyse COmponents Systems\n",
    "\n",
    "scipy_integrate = pyimport(\"scipy.integrate\")\n",
    "scipy_stats = pyimport(\"scipy.stats\")\n",
    "pushfirst!(PyVector(pyimport(\"sys\").\"path\"), \"\")\n",
    "frontiers = pyimport(\"frontiers_colors\");"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "wm = pyimport(\"watermark\").WaterMark()\n",
    "wm.watermark(\"-v -m  -u -n -p scipy -a Filippo_Valle -g -r -b -w\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Global Zipf "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "f_gtex = CSV.read(\"A_gtex.dat\", DataFrames.DataFrame).abundance\n",
    "f_gtex = sort(f_gtex./sum(f_gtex), rev=true);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "f_mca = Tacos.load_pickle(\"mca/data_all.pkl\")[\"freq\"]\n",
    "f_mca = sort(f_mca./sum(f_mca), rev=true);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "f_tm = Tacos.load_pickle(\"tm/data_all.pkl\")[\"freq\"]\n",
    "f_tm = sort(f_tm./sum(f_tm), rev=true);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "f_ss3 = CSV.read(\"A_ss3.dat\", DataFrames.DataFrame).abundance\n",
    "f_ss3 = sort(f_ss3./sum(f_ss3), rev=true);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "f_fibro = CSV.read(\"A_ss3_fibro.dat\", DataFrames.DataFrame).abundance\n",
    "f_fibro = sort(f_fibro./sum(f_fibro), rev=true);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = PyPlot.subplots(figsize=(18,15))\n",
    "\n",
    "X = collect(1:1e5)\n",
    "\n",
    "ax.plot(X,[1e-1*x^-1 for x in X], lw=20, ls=\"--\", color=\"gray\")\n",
    "\n",
    "ax.plot(f_gtex, label=\"GTEx\", lw=15, color=\"gray\", alpha=0.7)\n",
    "ax.plot(f_mca, label=\"Mouse Cell Atlas\", lw=15, color=\"blue\", alpha=0.6)\n",
    "ax.plot(f_tm, label=\"Tabula Muris\", lw=15, color=\"red\", alpha=0.5)\n",
    "\n",
    "\n",
    "ax.annotate(L\"$k*i^{-1}$\", (1e2,1e-2), fontsize=35)\n",
    "\n",
    "\n",
    "ax.legend(fontsize=35)\n",
    "ax.set_yscale(\"log\")\n",
    "ax.set_xscale(\"log\")\n",
    "\n",
    "ax.tick_params(labelsize=35, width=8, length=20)\n",
    "ax.set_xlabel(L\"Rank, $i$\", fontsize=35)\n",
    "ax.set_ylabel(L\"Frequency, $f_i$\", fontsize=35)\n",
    "\n",
    "ax.set_xlim(1,1e5)\n",
    "ax.set_ylim(1e-9,1e0)\n",
    "show()\n",
    "\n",
    "fig.savefig(\"zipf_alldatasets.pdf\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Fit ranges and regimes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# fit the power law\n",
    "fit_func = function(X, p)\n",
    "    return p[1]*[x^(-p[2]) for x in X]\n",
    "end\n",
    "\n",
    "fit_func_exp = function(X, p)\n",
    "    return p[1]*[exp(-p[2]*x) for x in X]\n",
    "end"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As we describe in the paper, here we search the bounds of the power-law like regime.\n",
    "\n",
    "Note that in Tabula Muris, we're looking for a local maximum, so *bounds[\"sx\"][\"max\"]* should be *82* instead of *122*"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "bounds = Dict()\n",
    "\n",
    "bounds[\"dx\"] = Dict()\n",
    "bounds[\"sx\"] = Dict()\n",
    "\n",
    "bounds[\"sx\"][\"min\"] = 2\n",
    "bounds[\"sx\"][\"max\"] = 82\n",
    "bounds[\"sx\"][\"step\"] = 10\n",
    "\n",
    "bounds[\"dx\"][\"min\"] = 5e2\n",
    "bounds[\"dx\"][\"max\"] = 2e4\n",
    "bounds[\"dx\"][\"step\"] = 100\n",
    "\n",
    "optimize_fit = function(f)\n",
    "    maxmimum_r2 = -Inf\n",
    "    optimal_fit = nothing\n",
    "    optimal_limits = nothing\n",
    "    fit_exp = nothing\n",
    "    \n",
    "    Rs = zeros(Float64, Int64((bounds[\"sx\"][\"max\"] - bounds[\"sx\"][\"min\"])/ bounds[\"sx\"][\"step\"]+1), Int64(( bounds[\"dx\"][\"max\"]- bounds[\"dx\"][\"min\"])/ bounds[\"dx\"][\"step\"]+1))\n",
    "    for left_side::Int32=range(bounds[\"sx\"][\"min\"],bounds[\"sx\"][\"max\"],step=bounds[\"sx\"][\"step\"]), right_side::Int32=range(bounds[\"dx\"][\"min\"],bounds[\"dx\"][\"max\"],step=bounds[\"dx\"][\"step\"])\n",
    "        #println(left_side,\" \", right_side)\n",
    "        \n",
    "        x_data = left_side:right_side\n",
    "        y_data = f[left_side:right_side]\n",
    "        fit = curve_fit(fit_func, x_data, y_data, [1e-1, 0.8])\n",
    "        popt, pcov = (coef(fit), estimate_covar(fit))\n",
    "        r_squared = Tacos.estimate_r2(fit_func, popt, x_data, y_data)\n",
    "        Rs[Int64((left_side-bounds[\"sx\"][\"min\"])/bounds[\"sx\"][\"step\"]+1),Int64((right_side-bounds[\"dx\"][\"min\"])/bounds[\"dx\"][\"step\"]+1)]=r_squared\n",
    "        if r_squared > maxmimum_r2\n",
    "            maxmimum_r2 = r_squared\n",
    "            optimal_fit = fit\n",
    "            optimal_limits = (left_side, right_side)\n",
    "        end    \n",
    "    end\n",
    "    return maxmimum_r2, optimal_fit, optimal_limits, Rs, bounds\n",
    "end"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Search the **bounds** of the regimes\n",
    "\n",
    "In the next cell you can chose from *f_gtex* (the one presented in the paper), *f_mca*, *f_tm* or *f_ss3*"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "r2, fit,  limits, Rs, windows = optimize_fit(f_mca);\n",
    "#CSV.write(\"rs_mca.csv\",DataFrames.DataFrame(Rs), writeheader=false);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = PyPlot.subplots(figsize=(18,15))\n",
    "\n",
    "popt = coef(fit)\n",
    "\n",
    "ax.plot(f_gtex, label=\"GTEx\", lw=15, color=\"gray\", alpha=0.7)\n",
    "ax.plot(f_mca, label=\"Mouse Cell Atlas\", lw=15, color=\"blue\", alpha=0.6)\n",
    "ax.plot(f_tm, label=\"Tabula Muris\", lw=15, color=\"red\", alpha=0.5)\n",
    "ax.plot(f_ss3, label=\"SmartSeq3\", lw=15, color=\"green\", alpha=0.5)\n",
    "ax.plot(f_fibro, label=\"Fibro\", lw=15, color=\"orange\", alpha=0.5)\n",
    "\n",
    "\n",
    "\n",
    "X = collect(1:1e5)\n",
    "\n",
    "#ax.plot(X,[1e-1*x^-0.8 for x in X], ls=\"--\", color=\"gray\")\n",
    "ax.plot(X,[popt[1]*x^-popt[2] for x in X], ls=\":\", lw=5, color=\"green\", label=\"$(@sprintf(\"%.2f\", popt[1]))*x^-$(@sprintf(\"%.2f\", popt[2]))\")\n",
    "\n",
    "#ax.annotate(L\"$k*i^{-0.8}$\", (1e2,1e-2), fontsize=22)\n",
    "ax.annotate(\"R^2=$(@sprintf(\"%.4f\", r2))\", (10^((log10(limits[1])+log10(limits[2]))/2),1e-6), fontsize=35)\n",
    "\n",
    "ax.vlines([limits[1] limits[2]], 1e-9, 1e-1, colors=\"gray\", ls=\"-.\", lw=15, alpha =0.5)\n",
    "ax.annotate(\"$(limits[1])\", (limits[1],2e-9), fontsize=35, alpha=0.8)\n",
    "ax.annotate(\"$(limits[2])\", (limits[2],2e-9), fontsize=35, alpha=0.8)\n",
    "\n",
    "\n",
    "ax.legend(fontsize=35)\n",
    "ax.set_yscale(\"log\")\n",
    "ax.set_xscale(\"log\")\n",
    "\n",
    "ax.set_xlabel(L\"Rank, $i$\", fontsize=35)\n",
    "ax.set_ylabel(L\"Frequency, $f_i$\", fontsize=35)\n",
    "\n",
    "ax.tick_params(labelsize=35, width=8, length=20)\n",
    "\n",
    "\n",
    "ax.set_xlim(1,1e5)\n",
    "ax.set_ylim(1e-9,1e0)\n",
    "show()\n",
    "\n",
    "fig.savefig(\"zipf_alldatasets_fit_ss3.pdf\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig,ax = PyPlot.subplots(5, 2, figsize=(18,40))\n",
    "\n",
    "for ((idataset,dataset), limits, sx_max, color) in zip(enumerate([\"tm\", \"mca\", \"gtex\", \"ss3\", \"ss3_fibro\"]),[(12,1000); (72,4100); (62,5100); (62,7900); (72,1300)],[82 121 121 72 82], [\"red\" \"blue\" \"gray\" \"green\" \"orange\"])\n",
    "    Rs = Matrix(CSV.read(\"rs_$(dataset).csv\",DataFrames.DataFrame))\n",
    "    \n",
    "    lax = idataset\n",
    "    rax = idataset+5\n",
    "    \n",
    "    ax[lax].plot(range(bounds[\"sx\"][\"min\"],sx_max-1,step=bounds[\"sx\"][\"step\"]), Rs[:,Int64((limits[2]-bounds[\"dx\"][\"min\"])/bounds[\"dx\"][\"step\"]+1)], lw=15, color=color)\n",
    "    ax[rax].plot(range(bounds[\"dx\"][\"min\"],bounds[\"dx\"][\"max\"],step=bounds[\"dx\"][\"step\"]), Rs[Int64((limits[1]-bounds[\"sx\"][\"min\"])/bounds[\"sx\"][\"step\"]),:], lw=15, color=color)\n",
    "    \n",
    "    ax[lax].vlines([limits[1]], 0.95, 1, colors=color, ls=\"-.\", lw=10, alpha =0.5)\n",
    "    ax[rax].vlines([limits[2]], 0.98, 1, colors=color, ls=\"-.\", lw=10, alpha =0.5)\n",
    "    \n",
    "    ax[lax].set_ylabel(L\"R^2\", fontsize=35)\n",
    "\n",
    "    ax[rax].set_xscale(\"log\")\n",
    "    \n",
    "    ax[lax].tick_params(labelsize=35, width=8, length=20)\n",
    "    ax[lax].tick_params(axis=\"both\", which=\"minor\", width=5, length=10)\n",
    "    \n",
    "    ax[rax].tick_params(labelsize=25, width=8, length=20)\n",
    "    ax[rax].tick_params(axis=\"both\", which=\"minor\", width=5, length=10)\n",
    "\n",
    "end\n",
    "\n",
    "ax[5].set_xlabel(\"left bound\", fontsize=35)\n",
    "ax[10].set_xlabel(\"right bound\", fontsize=35)\n",
    "\n",
    "ax[1].set_title(\"Tabula Muris\", fontsize=35)\n",
    "ax[2].set_title(\"Mouse Cell Atlas\", fontsize=35)\n",
    "ax[3].set_title(\"GTEx\", fontsize=35)\n",
    "ax[4].set_title(\"SS3\", fontsize=35)\n",
    "ax[5].set_title(\"SS3_fibro\", fontsize=35)\n",
    "\n",
    "\n",
    "ax[1].set_ylim(0.99,1)\n",
    "ax[4].set_ylim(0.98,1)\n",
    "ax[6].set_ylim(0.9950,1)\n",
    "ax[5].set_ylim(0.950,1)\n",
    "\n",
    "\n",
    "PyPlot.show()\n",
    "fig.savefig(\"bounds_panel.pdf\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for (f_data, name) in zip([f_mca, f_tm, f_gtex, f_ss3, f_fibro], [\"mca\" \"tm\" \"GTEx\" \"SSQ3\" \"fibro\"])\n",
    "    \n",
    "    print(\"\\\\multicolumn{2}{|c|}{$(name)} \\\\\\\\\",\"\\n\")\n",
    "    fit = curve_fit(fit_func, X[1:limits[1]], f_data[1:limits[1]], [1e-1, 0.5])\n",
    "    popt, pcov = (coef(fit), estimate_covar(fit))\n",
    "    print(round(popt[2], sigdigits=2),\" & \", round(sqrt(pcov[2*2]), sigdigits=1),\" \\\\\\\\ \\n\")\n",
    "    \n",
    "    fit = curve_fit(fit_func, X[limits[1]:limits[2]], f_data[limits[1]:limits[2]], [1e-1, 1])\n",
    "    popt, pcov = (coef(fit), estimate_covar(fit))\n",
    "    print(round(popt[2], sigdigits=2),\" & \", round(sqrt(pcov[2*2]), sigdigits=1),\"\\\\\\\\ \\n\")\n",
    "    \n",
    "    fit = curve_fit(fit_func_exp, X[limits[2]:length(f_data)], f_data[limits[2]:end], [1e-1, 1e-3])\n",
    "    popt_exp, pcov_exp = (coef(fit), estimate_covar(fit))\n",
    "    print(round(popt_exp[2], sigdigits=2),\" & \",round(sqrt(pcov_exp[2*2]), sigdigits=1),\"\\\\\\\\ \\\\hline \\n\")\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "limits = (Int64(round(mean([12 72 62]))), Int64(round(mean([1000 4100 5100])))) #tm mca gtex\n",
    "\n",
    "fig, ax = PyPlot.subplots(figsize=(18,15))\n",
    "\n",
    "ax.plot(f_gtex, label=\"GTEx\", lw=15, color=\"gray\", alpha=0.7)\n",
    "ax.plot(f_mca, label=\"Mouse Cell Atlas\", lw=15, color=\"blue\", alpha=0.6)\n",
    "ax.plot(f_tm, label=\"Tabula Muris\", lw=15, color=\"red\", alpha=0.5)\n",
    "\n",
    "X = collect(1:1e5)\n",
    "\n",
    "ax.vlines([limits[1] limits[2]], 1e-9, 1e-1, colors=\"gray\", ls=\"-.\", lw=15, alpha =0.5)\n",
    "ax.annotate(\"$(limits[1])\", (limits[1],2e-9), fontsize=35, alpha=0.8)\n",
    "ax.annotate(\"$(limits[2])\", (limits[2],2e-9), fontsize=35, alpha=0.8)\n",
    "\n",
    "\n",
    "ax.legend(fontsize=35)\n",
    "ax.set_yscale(\"log\")\n",
    "ax.set_xscale(\"log\")\n",
    "\n",
    "ax.set_xlabel(L\"Rank, $i$\", fontsize=35)\n",
    "ax.set_ylabel(L\"Frequency, $f_i$\", fontsize=35)\n",
    "\n",
    "ax.tick_params(labelsize=35, width=8, length=20)\n",
    "\n",
    "ax.set_xlim(1,1e5)\n",
    "ax.set_ylim(1e-9,1e0)\n",
    "show()\n",
    "\n",
    "fig.savefig(\"zipf_alldatasets_fit_avg.pdf\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Save first regime for GO"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The results of the GeneOntologies analyses are well described in the paper"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "first_regime = DataFrames.DataFrame()\n",
    "first_regime.mca = sort(CSV.read(\"A_mca.dat\",DataFrames.DataFrame), 2, rev=true).Column1[1:limits[1]]\n",
    "first_regime.tm = sort(CSV.read(\"A_tm.dat\",DataFrames.DataFrame), 2, rev=true).Column1[1:limits[1]]\n",
    "first_regime.gtex = sort(CSV.read(\"A_gtex.dat\",DataFrames.DataFrame), 2, rev=true).ensg[1:limits[1]]\n",
    "CSV.write(\"first_regime.csv\",first_regime)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fit = curve_fit(fit_func_exp, X[4100:length(f_mca)], f_mca[4100:end], [1e-1, 1e-3])\n",
    "popt_exp, pcov_exp = (coef(fit), estimate_covar(fit))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "df_nc = CSV.read(\"MGImarkerQuery_20200914_050053.txt\", DataFrames.DataFrame, delim=\"\\t\");\n",
    "names(df_nc)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "unique(df_nc[:,\"Feature Type\"])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pc = df_nc[df_nc[:,\"Feature Type\"] .== \"protein coding gene\",\"Symbol\"];\n",
    "nc = df_nc[collect(map(x->x in [\"lncRNA gene\" \"miRNA gene\" \"snoRNA gene\" \"rRNA gene\" \"tRNA gene\" \"non-coding RNA gene\"], df_nc[:,\"Feature Type\"])), \"Symbol\"];"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "f_mca = Tacos.load_pickle(\"mca/data_all.pkl\")[\"freq\"]\n",
    "mask = reshape(f_mca.index.isin(pc),(length(f_mca)))\n",
    "f_mca_pc = [v for v in f_mca][mask .== 1];\n",
    "f_mca_pc = sort(f_mca_pc./sum(f_mca_pc), rev=true);\n",
    "f_mca_nc = [v for v in f_mca][mask .!= 1];\n",
    "f_mca_nc = sort(f_mca_nc./sum(f_mca_nc), rev=true);\n",
    "f_mca = sort(f_mca./sum(f_mca), rev=true);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "limits = [122 4100]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = PyPlot.subplots(figsize=(18,15))\n",
    "\n",
    "X = collect(1:1e5)\n",
    "\n",
    "#ax.plot(X,[1e-1*x^-0.86 for x in X], lw=20, ls=\"--\", color=\"gray\")\n",
    "ax.plot(X[1:122],[1e-2*x^-0.5 for x in X[1:122]], lw=20, ls=\"--\", color=\"darkgray\")\n",
    "ax.plot(X[122:4100],[1e-1*x^-1 for x in X[122:4100]], lw=20, ls=\"--\", color=\"gray\")\n",
    "ax.plot(X[4100:end],fit_func_exp(X[4100:end], popt_exp), lw=20, ls=\"--\", color=\"gray\")\n",
    "\n",
    "\n",
    "#ax.plot(f_gtex, label=\"GTEx\", lw=15, color=\"gray\", alpha=0.7)\n",
    "ax.plot(f_mca_pc, label=\"Mouse Cell Atlas\", lw=15, color=\"red\", alpha=0.6)\n",
    "ax.plot(f_mca, label=\"Mouse Cell Atlas\", lw=15, color=\"blue\", alpha=0.6)\n",
    "#ax.plot(f_tm, label=\"Tabula Muris\", lw=15, color=\"red\", alpha=0.5)\n",
    "\n",
    "ax.vlines([limits[1]], 1e-9, 1e-1, colors=\"gray\", ls=\"-.\", lw=10, alpha =0.5)\n",
    "ax.vlines([limits[2]], 1e-9, 1e-1, colors=\"gray\", ls=\"-.\", lw=10, alpha =0.5)\n",
    "\n",
    "ax.annotate(L\"$k*i^{-0.5}$\", (3e0,0.5e-4), fontsize=35)\n",
    "ax.annotate(L\"$k*i^{-1}$\", (2e2,0.5e-4), fontsize=35)\n",
    "ax.annotate(L\"$k*e^{-0.002*i}$\", (0.8e4,0.5e-4), fontsize=35)\n",
    "\n",
    "ax.annotate(\"$(limits[1])\", (limits[1],2e-9), fontsize=35, alpha=0.8)\n",
    "ax.annotate(\"$(limits[2])\", (limits[2],2e-9), fontsize=35, alpha=0.8)\n",
    "\n",
    "ax.legend(fontsize=35)\n",
    "ax.set_yscale(\"log\")\n",
    "ax.set_xscale(\"log\")\n",
    "\n",
    "ax.tick_params(labelsize=35, width=8, length=20)\n",
    "ax.tick_params(axis=\"both\",which=\"minor\", width=5, length=10)\n",
    "\n",
    "ax.set_xlabel(L\"Rank, $i$\", fontsize=35)\n",
    "ax.set_ylabel(L\"Frequency, $f_i$\", fontsize=35)\n",
    "\n",
    "ax.set_xlim(1,1e5)\n",
    "ax.set_ylim(1e-9,1e0)\n",
    "show()\n",
    "\n",
    "fig.savefig(\"zipf_mca_nc.pdf\")\n",
    "fig.savefig(\"zipf_mca_nc.eps\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = PyPlot.subplots(figsize=(18,15))\n",
    "\n",
    "X = collect(1:1e5)\n",
    "\n",
    "#ax.plot(X,[1e-1*x^-0.86 for x in X], lw=20, ls=\"--\", color=\"gray\")\n",
    "ax.plot(X,[.5e-1*x^-0.8 for x in X], lw=20, ls=\"--\", color=\"darkgray\")\n",
    "\n",
    "ax.plot(f_mca_nc, label=\"Mouse Cell Atlas\", lw=15, color=\"blue\", alpha=0.6)\n",
    "\n",
    "ax.legend(fontsize=35)\n",
    "ax.set_yscale(\"log\")\n",
    "ax.set_xscale(\"log\")\n",
    "\n",
    "ax.tick_params(labelsize=35, width=8, length=20)\n",
    "ax.tick_params(axis=\"both\",which=\"minor\", width=5, length=10)\n",
    "\n",
    "ax.set_xlabel(L\"Rank, $i$\", fontsize=35)\n",
    "ax.set_ylabel(L\"Frequency, $f_i$\", fontsize=35)\n",
    "\n",
    "ax.set_xlim(1,1e5)\n",
    "ax.set_ylim(1e-7,1e0)\n",
    "show()\n",
    "fig.savefig(\"zipf_mca_only_nc.pdf\")\n",
    "fig.savefig(\"zipf_mca_only_nc.eps\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## CV^2 "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data_source = \"mca\"\n",
    "data = Tacos.load_pickle(\"$(data_source)/data_all.pkl\");"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = PyPlot.subplots(figsize=(18,15))\n",
    "\n",
    "means = [m for (m,v) in zip(data[\"means\"], data[\"var\"]) if v>1e-6]\n",
    "vars = [v for v in data[\"var\"] if v>1e-6]\n",
    "\n",
    "ax.scatter(means, [v/m/m for (m,v) in zip(means, vars)], c=\"gray\", alpha=0.3, s=225, label=\"data\")\n",
    "\n",
    "x = range(1e-5,10, length=10)\n",
    "ax.plot(x,[1/x_i for x_i in x], lw=5, color=\"red\", ls=\"--\", label=L\"$m_g^{-1}$ (Poisson)\")\n",
    "\n",
    "x = range(0.5,5e5, length=10)\n",
    "ax.plot(x,[1 for _ in x], lw=5, color=\"blue\", ls=\"--\", label=L\"$1$ (Taylor)\")\n",
    "\n",
    "ax.set_ylabel(L\"Coefficient of variation squared, $CV^2_g$\", fontsize=35)\n",
    "ax.set_xlabel(L\"Mean expression level, $m_g$\", fontsize=35)\n",
    "\n",
    "ax.set_xscale(\"log\")\n",
    "ax.set_yscale(\"log\")\n",
    "\n",
    "ax.tick_params(labelsize=35, width=8, length=20)\n",
    "\n",
    "ax.set_xlim(1e-5, maximum(means)*1.5)\n",
    "ax.set_ylim(5e-2, 1e5)\n",
    "\n",
    "ax.legend(fontsize=35)\n",
    "\n",
    "PyPlot.show()\n",
    "fig.savefig(\"cv2_mean_$(data_source).pdf\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Occurrency distribution of all the data. Aka the global U"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data = Tacos.load_pickle(\"mca/data_all.pkl\");"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = PyPlot.subplots(figsize=(18,15))\n",
    "\n",
    "bins = -0.05:0.1:1.05\n",
    "\n",
    "bin_counts, bin_edges = ax.hist(data[\"O\"].values, color=\"blue\", lw=10, histtype=\"step\", bins=bins)\n",
    "\n",
    "x_bins = (bin_edges[1:length(bin_edges)-1]+bin_edges[2:length(bin_edges)])/2\n",
    "ax.hist(x_bins, weights=bin_counts, bins=bins, histtype=\"stepfilled\", lw=0, alpha=0.8)\n",
    "\n",
    "ax.set_xlabel(L\"Occurrence, O_i\", fontsize=35)\n",
    "ax.set_ylabel(\"Number of genes\", fontsize=35)\n",
    "\n",
    "ax.tick_params(labelsize=35, width=8, length=20)\n",
    "ax.tick_params(which=\"minor\", width=6, length=10)\n",
    "\n",
    "PyPlot.show()\n",
    "fig.savefig(\"U_mca.pdf\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# UMI count histogram"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data = Tacos.load_pickle(\"mca/data_all.pkl\");"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = PyPlot.subplots(figsize=(18,15))\n",
    "\n",
    "x_bins = sort!([10^i for i in log10(500):0.01:5])\n",
    "bin_counts, bin_edges = ax.hist(data[\"M\"], density=true, color=\"blue\", histtype=\"step\", lw=10, bins=x_bins)\n",
    "\n",
    "m = mean(data[\"M\"])\n",
    "ax.vlines([m], 0,1e-3, lw=15, alpha=0.8,  ls=\"-.\")\n",
    "ax.annotate(\"Average: $(@sprintf(\"%.0f\", m))\", (m*1.1,5e-4), fontsize=35)\n",
    "\n",
    "m = median(data[\"M\"])\n",
    "ax.vlines([m], 0,2e-3, lw=15, alpha=0.8,  ls=\"-.\")\n",
    "ax.annotate(\"Median: $(@sprintf(\"%.0f\", m))\", (m*1.1,10e-4), fontsize=35)\n",
    "\n",
    "x_bins=(bin_edges[1:length(bin_edges)-1]+bin_edges[2:length(bin_edges)])/2\n",
    "\n",
    "scaled_pois = function(X, par)\n",
    "    mu = par[1]\n",
    "    C = par[2]\n",
    "    a = par[3]\n",
    "    d = Poisson(mu)\n",
    "    return C*[pdf(d,round(Int, x/a::Float64)) for x in X]\n",
    "end\n",
    "\n",
    "fit_pois = curve_fit(scaled_pois, x_bins, bin_counts, [0.5, 0.1, 1000])\n",
    "popt = coef(fit_pois)\n",
    "println(popt)\n",
    "\n",
    "#ax.plot(x_bins, scaled_pois(x_bins, popt), lw=10, color=\"orange\")\n",
    "\n",
    "ax.set_xlabel(L\"UMI per cell, $M$\", fontsize=35)\n",
    "ax.set_ylabel(\"Probability density function\", fontsize=35)\n",
    "\n",
    "ax.tick_params(labelsize=35, width=8, length=20)\n",
    "\n",
    "ax.set_xscale(\"log\")\n",
    "#ax.set_yscale(\"log\")\n",
    "\n",
    "#ax.set_xlim(1,3e8)\n",
    "\n",
    "#ax.legend(fontsize=22)\n",
    "\n",
    "PyPlot.show()\n",
    "fig.savefig(\"M_hist_mca.pdf\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Pdf of Frequencies"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data = Tacos.load_pickle(\"mca/data_all.pkl\");"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = PyPlot.subplots(figsize=(18,15))\n",
    "\n",
    "\n",
    "#ax.plot(sort!(data[\"freq\"].values,rev=true), label=\"nan=0\", lw=10)\n",
    "ax.plot(sort!(data[\"means\"].values/sum(filter(!isnan, data[\"means\"].values)),rev=true), label=\"w/o nan\", lw=10)\n",
    "#ax.plot(sort!(data[\"means_nonzero\"].values/sum(filter(!isnan, data[\"means_nonzero\"].values)),rev=true), ms=50, label=\"m0\", lw=10, color=\"gray\")\n",
    "\n",
    "ax.set_xscale(\"log\")\n",
    "ax.set_yscale(\"log\")\n",
    "\n",
    "ax.legend(fontsize=35)\n",
    "ax.set_yscale(\"log\")\n",
    "ax.set_xscale(\"log\")\n",
    "\n",
    "ax.set_ylim(1e-9,1e-1)\n",
    "\n",
    "ax.tick_params(labelsize=35, width=8, length=20)\n",
    "ax.set_xlabel(L\"Rank, $i$\", fontsize=35)\n",
    "ax.set_ylabel(L\"Frequency, $f_i$\", fontsize=35)\n",
    "fig.savefig(\"pdf_freq_mca.pdf\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = PyPlot.subplots(figsize=(18,15))\n",
    "\n",
    "f = filter(!isnan, data[\"means\"].values)\n",
    "f = f / sum(f)\n",
    "#f = data[\"freq\"].values\n",
    "f = f[f .> 0]\n",
    "#f = f[f .< 5e-3]\n",
    "f = sort(f, rev=true)\n",
    "\n",
    "\n",
    "#x_bins = sort!([10. ^ -x for x in 1:0.2:8])\n",
    "x_bins = sort!([10. ^ x for x in log10(minimum(f)):0.2:log10(maximum(f))])\n",
    "bin_counts, bin_edges = ax.hist(f, bins=x_bins, density=false, color=\"gray\", histtype=\"step\", lw=25)\n",
    "\n",
    "ax.hist(f[1:123], bins=x_bins, density=false, color=\"red\", histtype=\"step\", lw=25)\n",
    "ax.hist(f[123:15000], bins=x_bins, density=false, color=\"green\", histtype=\"step\", lw=25)\n",
    "ax.hist(f[15001:end], bins=x_bins, density=false, color=\"blue\", histtype=\"step\", lw=25)\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "ax.set_xlabel(L\"Frequency, f_i\", fontsize=45)\n",
    "ax.set_ylabel(\"Probability density function\", fontsize=45)\n",
    "\n",
    "ax.tick_params(labelsize=35, width=8, length=20)\n",
    "ax.tick_params(which=\"minor\", width=6, length=10)\n",
    "ax.tick_params(axis=\"x\", which=\"both\", rotation=45)\n",
    "\n",
    "#kde_values = pdf(kde(f), x_bins)\n",
    "#ax.plot(x_bins,kde_values/sum(kde_values)*30000)\n",
    "\n",
    "ax.set_xscale(\"log\")\n",
    "#ax.set_yscale(\"log\")\n",
    "\n",
    "ax.set_ylim(5e-1,4e3)\n",
    "\n",
    "PyPlot.tight_layout()\n",
    "PyPlot.show()\n",
    "fig.savefig(\"freq_pdf_mca.pdf\")\n",
    "fig.savefig(\"freq_pdf_mca.eps\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = PyPlot.subplots(figsize=(18,15))\n",
    "\n",
    "X = collect(1:1e5)\n",
    "\n",
    "ax.plot(X,[1e-1*x^-1 for x in X], lw=20, ls=\"--\", color=\"gray\")\n",
    "\n",
    "ax.plot(range(1,123, step=1), f[1:123], label=\"MCA\", lw=15, color=\"red\", alpha=0.7)\n",
    "ax.plot(range(123, 15001, step=1),f[123:15001], label=\"MCA\", lw=15, color=\"green\", alpha=0.7)\n",
    "ax.plot(range(15001,length(f), step=1),f[15001:end], label=\"MCA\", lw=15, color=\"blue\", alpha=0.7)\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "ax.annotate(L\"$k*i^{-1}$\", (1e2,1e-2), fontsize=35)\n",
    "\n",
    "\n",
    "ax.legend(fontsize=35)\n",
    "ax.set_yscale(\"log\")\n",
    "ax.set_xscale(\"log\")\n",
    "\n",
    "ax.tick_params(labelsize=35, width=8, length=20)\n",
    "ax.set_xlabel(L\"Rank, $i$\", fontsize=35)\n",
    "ax.set_ylabel(L\"Frequency, $f_i$\", fontsize=35)\n",
    "\n",
    "ax.set_xlim(1,1e5)\n",
    "ax.set_ylim(1e-9,1e0)\n",
    "show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Analyse single organ"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tissue = \"Bone-Marrow_c-kit\"\n",
    "\n",
    "data = Dict()\n",
    "data[\"data\"]=Tacos.load_pickle(\"mca/data_$(tissue)_data.pkl\")\n",
    "data[\"mazzolini\"]=Tacos.load_pickle(\"mca/data_$(tissue)_mazzolini.pkl\");\n",
    "#data[\"null_1\"]=load_pickle(\"data_$(tissue)_null_1.pkl\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This is not shown in the paper. It is a check that sampled and true frequency distributios are actually the same"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = PyPlot.subplots(figsize=(18,15))\n",
    "ax.plot(sort(data[\"data\"][\"freq\"].values, rev=true), lw=15, label=\"data\")\n",
    "ax.plot(sort(filter(!isnan,data[\"mazzolini\"][\"freq\"]), rev=true), lw=15, ls=\"--\", label=\"mazzolini\")\n",
    "\n",
    "ax.legend(fontsize=35)\n",
    "ax.set_yscale(\"log\")\n",
    "ax.set_xscale(\"log\")\n",
    "\n",
    "ax.tick_params(labelsize=35, width=8, length=20)\n",
    "ax.tick_params(which=\"minor\", labelsize=35, width=5, length=10)\n",
    "\n",
    "ax.set_xlabel(\"Rank, i\", fontsize=35)\n",
    "ax.set_ylabel(\"f_i\", fontsize=35)\n",
    "show()\n",
    "\n",
    "fig.savefig(\"null_zipf.pdf\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Heaps"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "loglog=true\n",
    "\n",
    "fig, ax = PyPlot.subplots(figsize=(18,15))\n",
    "\n",
    "ax.scatter(data[\"data\"][\"M\"], data[\"data\"][\"diffWords\"], c=\"gray\", alpha=0.3, s=225, label=\"data\")\n",
    "bin_means, bin_edges = scipy_stats.binned_statistic(data[\"data\"][\"M\"], data[\"data\"][\"diffWords\"], bins=15)\n",
    "\n",
    "x_bins = (bin_edges[1:10] + bin_edges[2:11])/2\n",
    "ax.plot(x_bins, bin_means[1:10], c=\"blue\", alpha=0.9, lw=5, ms=25, marker=\"o\", label=\"binned_data\")\n",
    "\n",
    "\n",
    "ax.set_xlabel(\"UMI per cell\", fontsize=35)\n",
    "ax.set_ylabel(\"Number of expressed genes\", fontsize=35)\n",
    "\n",
    "if loglog\n",
    "    ax.set_xscale(\"log\")\n",
    "    ax.set_yscale(\"log\")\n",
    "end\n",
    "\n",
    "ax.legend(fontsize=35)\n",
    "\n",
    "ax.tick_params(labelsize=35, width=8, length=35)\n",
    "\n",
    "PyPlot.show()\n",
    "if loglog\n",
    "    fig.savefig(\"heaps_boneMarrow_MCA_loglog.pdf\")\n",
    "else\n",
    "    fig.savefig(\"heaps_boneMarrow_MCA.pdf\")\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = PyPlot.subplots(figsize=(18,15))\n",
    "\n",
    "ax.scatter(data[\"data\"][\"M\"], data[\"data\"][\"diffWords\"], c=\"gray\", alpha=0.3, s=225, label=\"data\")\n",
    "ax.scatter(data[\"mazzolini\"][\"M\"], data[\"mazzolini\"][\"diffWords\"], c=\"red\", alpha=0.3, s=225, label=\"sampling model\")\n",
    "\n",
    "x_bins = (bin_edges[1:10] + bin_edges[2:11])/2\n",
    "ax.plot(x_bins, bin_means[1:10], c=\"blue\", alpha=0.9, ms=25, marker=\"o\", label=\"binned_data\")\n",
    "\n",
    "\n",
    "ax.set_xlabel(\"UMI per cell\", fontsize=35)\n",
    "ax.set_ylabel(\"Number of expressed genes\", fontsize=35)\n",
    "\n",
    "ax.tick_params(labelsize=35, width=8, length=35)\n",
    "ax.tick_params(which=\"minor\", labelsize=35, width=5, length=25)\n",
    "\n",
    "\n",
    "if loglog\n",
    "    ax.set_xscale(\"log\")\n",
    "    ax.set_yscale(\"log\")\n",
    "end\n",
    "\n",
    "ax.legend(fontsize=35)\n",
    "\n",
    "PyPlot.show()\n",
    "if loglog\n",
    "    fig.savefig(\"heaps_boneMarrow_sampling_MCA_loglog.pdf\")\n",
    "else\n",
    "    fig.savefig(\"heaps_boneMarrow_sampling_MCA.pdf\")\n",
    "end"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### U"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = PyPlot.subplots(figsize=(18,15))\n",
    "\n",
    "density=true\n",
    "\n",
    "bins = [10^x for x in range(log10(3e-5), log10(1.05), length=10)]\n",
    "#bins = 0.003:0.05:1.05\n",
    "bin_counts, bin_edges = ax.hist(data[\"data\"][\"O\"].values, color=\"gray\", lw=10, histtype=\"step\", bins=bins, density=density)\n",
    "\n",
    "x_bins = (bin_edges[1:length(bin_edges)-1]+bin_edges[2:length(bin_edges)])/2\n",
    "ax.hist(x_bins, weights=bin_counts, bins=bins, histtype=\"stepfilled\", color=\"gray\", lw=0, alpha=0.3)\n",
    "\n",
    "bin_counts, bin_edges = ax.hist(data[\"mazzolini\"][\"O\"], color=\"red\", lw=10, histtype=\"step\", bins=bins, density=density)\n",
    "\n",
    "\n",
    "ax.set_xlabel(L\"Occurrence, O_i\", fontsize=35)\n",
    "ax.set_ylabel(\"Number of genes\", fontsize=35)\n",
    "\n",
    "ax.tick_params(labelsize=35, width=8, length=20)\n",
    "ax.tick_params(which=\"minor\", width=6, length=10)\n",
    "\n",
    "ax.set_yscale(\"log\")\n",
    "ax.set_xscale(\"log\")\n",
    "\n",
    "PyPlot.show()\n",
    "fig.savefig(\"U_$(tissue)_mca.pdf\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Julia 1.5.1",
   "language": "julia",
   "name": "julia-1.5"
  },
  "language_info": {
   "file_extension": ".jl",
   "mimetype": "application/julia",
   "name": "julia",
   "version": "1.5.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}