Fig3_TF_Enrichment.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f1e0b609",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import math\n",
    "import matplotlib.pyplot as plt\n",
    "import csv\n",
    "import pandas as pd\n",
    "from scipy import stats\n",
    "import re\n",
    "from scipy.stats import pearsonr,spearmanr,fisher_exact,binom_test\n",
    "#import rpy2.robjects as robjects\n",
    "import random\n",
    "from statsmodels.stats.multitest import fdrcorrection\n",
    "import os\n",
    "import seaborn as sns\n",
    "#import gseapy as gs\n",
    "from scipy.stats import norm\n",
    "import gseapy as gs\n",
    "from collections import Counter\n",
    "import scanpy as sc"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fed1add3",
   "metadata": {},
   "outputs": [],
   "source": [
    "#Colors\n",
    "mouse = \"#F2C911\"\n",
    "rat = \"#65B0AC\"\n",
    "intrinsic = \"#F55F00\"\n",
    "extrinsic = \"#7D9AF4\"\n",
    "reinforcing = \"#9B00F5\"\n",
    "opposing = \"#F50901\"\n",
    "interaction = \"#1E771A\"\n",
    "\n",
    "sns.set(font_scale=1.5)\n",
    "sns.set_style(\"white\")\n",
    "d_ct_abrev = {\"Chondrocyte\":\"chondrocyte\", \"Forebrain glutamatergic progenitors\":\"brain.glut.prog\", \"Forebrain GABAergic progenitors\":\"brain.GABA.prog\", \"Intermediate progenitors\":\"inter.prog\", \"Forebrain glutamatergic neurons\":\"brain.glut.neu\", \"Forebrain GABAergic neurons\":\"brain.GABA.neu\", \"Spinal GABAergic neurons\":\"spine.GABA.neu\", \"Spinal glutamatergic neurons\":\"spine.glut.neu\", \"Chondrocytes\":\"chondrocyte\", \"Mesenchyme 0\":\"mesen.0\", \"Mesenchyme 2\":\"mesen.2\", \"Mesenchyme cycling\":\"mesen.cyc\", \"Forebrain GABAergic neurons 0\":\"Toss\"}\n",
    "\n",
    "#Function to convert file to cell type name\n",
    "def file_to_celltype(x):\n",
    "    x = x.split(\"/\")[-1]\n",
    "    x = x.replace(\"_GO_Molecular_Function_2023\", \"\").replace(\"_GO_Biological_Process_2023\", \"\").replace(\"Mesechyme\", \"Mesenchyme\").replace(\"Brain_Or_2010_Div_New4_NewNorm_\", \"\")\n",
    "    x = x.replace(\"Mesenchymal_Or_2010_Div_New4_NewNorm_\", \"\").replace(\"Chondrocyte_Or_2010_Div_New4_NewNorm_\", \"\").replace(\".csv\", \"\").replace(\"_\", \" \").replace(\" all\", \"\")\n",
    "    x = x.replace(\"Glutamatergic\", \"Forebrain glutamatergic\").replace(\"GABAergic\", \"Forebrain GABAergic\").replace(\"Spinal Forebrain\", \"Spinal\").replace(\"Chondrocytes\", \"Chondrocyte\")\n",
    "    return x\n",
    "\n",
    "plt.rcParams[\"font.family\"] = \"Arial\"\n",
    "\n",
    "#For more comments, see Enrichment_Figures3-4_SFigs.ipynb\n",
    "def summarize_enrichment(fold, out_file):\n",
    "    #Summarize the enrichments\n",
    "    out = []\n",
    "    for folder in os.listdir(\"Final/\" + fold):\n",
    "        v = pd.read_csv(\"Final/\" + fold + \"/\" + folder + \"/\" + \"gseapy.gene_set.prerank.report.csv\")\n",
    "        for index, row in v.iterrows():\n",
    "            if row[\"FDR q-val\"] < 0.25:\n",
    "                out.append([folder, row[\"Term\"], row[\"FDR q-val\"], row[\"NES\"], row[\"Lead_genes\"]])\n",
    "\n",
    "    df = pd.DataFrame(out)\n",
    "    df.columns = [\"Cluster_Category\", \"Term\", \"FDR\", \"NES\", \"Ledge genes\"]\n",
    "    df = df.sort_values(\"FDR\")\n",
    "    df.to_csv(\"Final/Summarized_Enrichments_\" + out_file + \".csv\", index = False)\n",
    "    \n",
    "vals = [\"Extrinsic\", \"Intrinsic\", \"Interaction\"]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e4113c94",
   "metadata": {},
   "outputs": [],
   "source": [
    "m2h = pd.read_csv(\"Mouse_To_Human_Gene_Conversions.csv\").dropna()\n",
    "d = {}\n",
    "for index, row in m2h.iterrows():\n",
    "    d[row[\"external_gene_name\"]] = row[\"hsapiens_homolog_associated_gene_name\"]\n",
    "gene_sets = ['TF_Perturbations_Followed_by_Expression']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "438218f5",
   "metadata": {},
   "outputs": [],
   "source": [
    "#Aggregating by tissue and doing enrichment analysis for signed things\n",
    "vals = [\"Extrinsic\", \"Intrinsic\", \"Interaction\"]\n",
    "for keep in [\"All\", \"Brain\", \"Conn\"]:\n",
    "    for divergence in vals:\n",
    "        df = pd.DataFrame()\n",
    "        ind = 1\n",
    "        #Iterate through the files\n",
    "        for file in os.listdir(\"Final/Div\"):\n",
    "            #If we want to keep all of them regardless of tissue\n",
    "            if keep == \"All\":\n",
    "                v = pd.read_csv(\"Final/Div/\" + file).set_index(\"Unnamed: 0\")\n",
    "                v[\"Signed \" + divergence.lower()] = np.sign(v[divergence])*v[\"Proportion \" + divergence.lower()]\n",
    "                v = v[[\"Signed \" + divergence.lower()]]\n",
    "                v.columns = [d_ct_abrev[file_to_celltype(file)].replace(\".\", \"_\")]\n",
    "                if ind:\n",
    "                    df = v\n",
    "                    ind = 0\n",
    "                else:\n",
    "                    df = df.join(v, how=\"outer\")\n",
    "            \n",
    "            else:\n",
    "                #If we only want neuronal cell types\n",
    "                if keep == \"Brain\":\n",
    "                    if \"Brain\" in file:\n",
    "                        v = pd.read_csv(\"Final/Div/\" + file).set_index(\"Unnamed: 0\")\n",
    "                        v[\"Signed \" + divergence.lower()] = np.sign(v[divergence])*v[\"Proportion \" + divergence.lower()]\n",
    "                        v = v[[\"Signed \" + divergence.lower()]]\n",
    "                        v.columns = [d_ct_abrev[file_to_celltype(file)].replace(\".\", \"_\")]\n",
    "                        if ind:\n",
    "                            df = v\n",
    "                            ind = 0\n",
    "                        else:\n",
    "                            df = df.join(v, how=\"outer\")\n",
    "                #If we only want connective tissue cell types\n",
    "                elif keep == \"Conn\":\n",
    "                    if \"Brain\" not in file:\n",
    "                        v = pd.read_csv(\"Final/Div/\" + file).set_index(\"Unnamed: 0\")\n",
    "                        v[\"Signed \" + divergence.lower()] = np.sign(v[divergence])*v[\"Proportion \" + divergence.lower()]\n",
    "                        v = v[[\"Signed \" + divergence.lower()]]\n",
    "                        v.columns = [d_ct_abrev[file_to_celltype(file)].replace(\".\", \"_\")]\n",
    "                        if ind:\n",
    "                            df = v\n",
    "                            ind = 0\n",
    "                        else:\n",
    "                            df = df.join(v, how=\"outer\")\n",
    "        out = []\n",
    "        indices = []\n",
    "        #Require a gene be measured in at least 5 cell types if all, at least 3 cell types if connective tissue or nervous only\n",
    "        for index, row in df.iterrows():\n",
    "            if keep == \"All\":\n",
    "                if len(row.dropna()) > 4:\n",
    "                    out.append(np.mean(row.dropna()))\n",
    "                    indices.append(index)\n",
    "            else:\n",
    "                if len(row.dropna()) > 2:\n",
    "                    out.append(np.mean(row.dropna()))\n",
    "                    indices.append(index)\n",
    "        df2 = pd.DataFrame(out)\n",
    "        df2.index = indices\n",
    "        df2.columns = [\"Mean signed \" + divergence.lower()]\n",
    "        df2.sort_values(\"Mean signed \" + divergence.lower())\n",
    "        out = []\n",
    "        \n",
    "        #Same enrichment strategy as above\n",
    "        for index, row in df2.iterrows():\n",
    "            try:\n",
    "                out.append([d[index], row[\"Mean signed \" + divergence.lower()]])\n",
    "            except:\n",
    "                pass\n",
    "        ranking = pd.DataFrame(out)\n",
    "        ranking.columns = [0, 1]\n",
    "        for gene_set in gene_sets:\n",
    "            try:\n",
    "                gs.prerank(rnk=ranking, gene_sets=gene_set, threads=4, permutation_num=1000, outdir= 'C:/Users/astar/Chimerism/Final/GSEAPY_New4_Signed_Proportion_Mean_Signed_' + divergence.lower() + '/' + keep + \"_\" + gene_set, format='png', seed=6, min_size = 10, max_size = 300)\n",
    "            except:\n",
    "                pass"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6498529c",
   "metadata": {},
   "outputs": [],
   "source": [
    "#Deduplicate the output to retain only one TF + manipulation + direction\n",
    "divergence = \"extrinsic\"\n",
    "\n",
    "v = pd.read_csv(\"Final/GSEAPY_Signed_Proportion_Mean_Signed_TF_\" + divergence + \"/All_TF_TF_Perturbations_Followed_by_Expression/gseapy.gene_set.prerank.report.csv\")\n",
    "v = v.sort_values(\"FDR q-val\")\n",
    "dedup = []\n",
    "c = 0\n",
    "\n",
    "#For each entry\n",
    "for index, row in v.iterrows():\n",
    "    to_add = \"\"\n",
    "    \n",
    "    #Split the term to only retain the name of the TF and whether it is down or up\n",
    "    l = row[\"Term\"].split(\" \")\n",
    "    for i in range(len(l)):\n",
    "        if i == 0:\n",
    "            to_add = l[i]\n",
    "        elif i == 1:\n",
    "            to_add = to_add + \" \" + l[i]\n",
    "        elif \"GSE\" in l[i]:\n",
    "            to_add = to_add + \" \" + l[i]\n",
    "        elif l[i] == \"DOWN\" or l[i] == \"UP\":\n",
    "            to_add = to_add + \" \" + l[i]\n",
    "            c += 1\n",
    "    dedup.append(to_add)\n",
    "v[\"Dedup\"] = dedup\n",
    "\n",
    "#Remove duplicates.  E.g. if there are multiple Jun UP, there will only be one Jun UP after this.\n",
    "v = v.drop_duplicates(\"Dedup\")\n",
    "\n",
    "#Replace 0 q-values with actual minimum based on number of permutations\n",
    "v[\"-Log$_{10}$ FDR\"] = -np.log10(np.maximum(v[\"FDR q-val\"], 0.0001))\n",
    "\n",
    "#Rename columns and make plot/write out\n",
    "v[\"Enrichment score for TF targets\"] = v[\"ES\"]\n",
    "v[\"Significance\"] = list(np.repeat(\"FDR < 0.1\", v[v[\"FDR q-val\"] < 0.1].shape[0])) + list(np.repeat(\"Not significant\", v[v[\"FDR q-val\"] >= 0.1].shape[0]))\n",
    "sns.scatterplot(data = v, x = \"Enrichment score for TF targets\", y = \"-Log$_{10}$ FDR\", hue = \"Significance\", palette = {\"FDR < 0.1\":\"red\", \"Not significant\":\"grey\"})\n",
    "plt.title(\"Signed proportion \" + divergence + \" divergence\")\n",
    "plt.legend([],[], frameon=False)\n",
    "v.to_csv(\"Final/GSEAPY_Signed_Proportion_Mean_Signed_TF_\" + divergence + \"/All_TF_TF_Perturbations_Followed_by_Expression/Dedupped_table.csv\", index = False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9fab2347",
   "metadata": {},
   "outputs": [],
   "source": [
    "#Recomputed dataframe that was used for enrichment\n",
    "#See cell 4 for detailed comments on the procedure\n",
    "for keep in [\"All\"]:\n",
    "    divergence = \"Intrinsic\"\n",
    "    df = pd.DataFrame()\n",
    "    ind = 1\n",
    "    for file in os.listdir(\"Final/Div\"):\n",
    "\n",
    "        if keep == \"All\":\n",
    "            v = pd.read_csv(\"Final/Div/\" + file).set_index(\"Unnamed: 0\")\n",
    "            v[\"Signed \" + divergence.lower()] = np.sign(v[divergence])*v[\"Proportion \" + divergence.lower()]\n",
    "            v = v[[\"Signed \" + divergence.lower()]]\n",
    "            v.columns = [d_ct_abrev[file_to_celltype(file)].replace(\".\", \"_\")]\n",
    "            if ind:\n",
    "                df = v\n",
    "                ind = 0\n",
    "            else:\n",
    "                df = df.join(v, how=\"outer\")\n",
    "        else:\n",
    "            if keep == \"Brain\":\n",
    "                if \"Brain\" in file:\n",
    "                    v = pd.read_csv(\"Final/Div/\" + file).set_index(\"Unnamed: 0\")\n",
    "                    v[\"Signed \" + divergence.lower()] = np.sign(v[divergence])*v[\"Proportion \" + divergence.lower()]\n",
    "                    v = v[[\"Signed \" + divergence.lower()]]\n",
    "                    v.columns = [d_ct_abrev[file_to_celltype(file)].replace(\".\", \"_\")]\n",
    "                    if ind:\n",
    "                        df = v\n",
    "                        ind = 0\n",
    "                    else:\n",
    "                        df = df.join(v, how=\"outer\")\n",
    "            elif keep == \"Conn\":\n",
    "                if \"Brain\" not in file:\n",
    "                    v = pd.read_csv(\"Final/Div/\" + file).set_index(\"Unnamed: 0\")\n",
    "                    v[\"Signed \" + divergence.lower()] = np.sign(v[divergence])*v[\"Proportion \" + divergence.lower()]\n",
    "                    v = v[[\"Signed \" + divergence.lower()]]\n",
    "                    v.columns = [d_ct_abrev[file_to_celltype(file)].replace(\".\", \"_\")]\n",
    "                    if ind:\n",
    "                        df = v\n",
    "                        ind = 0\n",
    "                    else:\n",
    "                        df = df.join(v, how=\"outer\")\n",
    "    out = []\n",
    "    indices = []\n",
    "    for index, row in df.iterrows():\n",
    "        if keep == \"All\":\n",
    "            if len(row.dropna()) > 4:\n",
    "                out.append(np.mean(row.dropna()))\n",
    "                indices.append(index)\n",
    "        else:\n",
    "            if len(row.dropna()) > 2:\n",
    "                out.append(np.mean(row.dropna()))\n",
    "                indices.append(index)\n",
    "    df2 = pd.DataFrame(out)\n",
    "    df2.index = indices\n",
    "    df2.columns = [\"Mean signed \" + divergence.lower()]\n",
    "    df2.sort_values(\"Mean signed \" + divergence.lower())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "93e2fb7c",
   "metadata": {},
   "outputs": [],
   "source": [
    "#New style volcano plots\n",
    "#Reran for intrinsic and extrinsic\n",
    "\n",
    "fdr = []\n",
    "enrichment = []\n",
    "ct_list = []\n",
    "term_list = []\n",
    "vv = pd.read_csv(\"Final/GSEAPY_Signed_Proportion_Mean_Signed_TF_intrinsic/\" + \"All_TF_TF_Perturbations_Followed_by_Expression\" + \"/Dedupped_table.csv\", sep = \",\")\n",
    "out = []\n",
    "fdr_cat = []\n",
    "for index, row in df2.iterrows():\n",
    "    try:\n",
    "        out.append([d[index], row[\"Mean signed intrinsic\"]])\n",
    "    except:\n",
    "        pass\n",
    "v = pd.DataFrame(out)\n",
    "v.columns = [\"Gene\", \"Mean signed intrinsic\"]\n",
    "for index, row in vv.iterrows():\n",
    "    term = row[\"Term\"]\n",
    "    term_list.append(row[\"Term\"])\n",
    "    gene_set = d_BP[term]\n",
    "    genes = row[\"Lead_genes\"].split(\";\")\n",
    "    fdr.append(row[\"-Log$_{10}$ FDR\"])\n",
    "    fdr_cat.append(row[\"Significance\"])\n",
    "    if row[\"ES\"] < 0:\n",
    "        asc = True\n",
    "    else:\n",
    "        asc = False\n",
    "    v = v.sort_values(\"Mean signed intrinsic\", ascending = asc)\n",
    "\n",
    "    cut = 0\n",
    "    for index, row in v.iterrows():\n",
    "        cut += 1\n",
    "        if genes[-1] == row[\"Gene\"]:\n",
    "            break\n",
    "    genes_agree = len(np.intersect1d(list(v[\"Gene\"])[:cut], gene_set))\n",
    "    genes_disagree = len(np.intersect1d(list(v[\"Gene\"])[::-1][:cut], gene_set))\n",
    "    ct_list.append(d_ct_abrev[file_to_celltype(file)])\n",
    "    if asc:\n",
    "        enrichment.append(-np.log2((genes_agree + 1)/(genes_disagree + 1)))\n",
    "    else:\n",
    "        enrichment.append(np.log2((genes_agree + 1)/(genes_disagree + 1)))\n",
    "\n",
    "df_plot_new = pd.DataFrame([term_list, ct_list, enrichment, fdr, fdr_cat]).T\n",
    "df_plot_new.columns = [\"Gene set\", \"Cell type\", \"Log$_{2}$ fold-enrichment\", \"-Log$_{10}$ FDR\", \"Significance\"]\n",
    "df_plot_new[\"Log$_{2}$ fold-enrichment\"] = df_plot_new[\"Log$_{2}$ fold-enrichment\"].astype(float)\n",
    "\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "# Create subplots\n",
    "fig, ax = plt.subplots(sharex=True, sharey=True)\n",
    "\n",
    "# Separate data by color\n",
    "grey_data = df_plot_new[df_plot_new['Significance'] != 'FDR < 0.1']\n",
    "blue_data = df_plot_new[df_plot_new['Significance'] == 'FDR < 0.1']\n",
    "\n",
    "# Scatter plots with zorder\n",
    "sns.scatterplot(x=grey_data['Log$_{2}$ fold-enrichment'], y=grey_data['-Log$_{10}$ FDR'], c='#BBBBBC', label='Not significant', zorder=1)\n",
    "sns.scatterplot(x=blue_data['Log$_{2}$ fold-enrichment'], y=blue_data['-Log$_{10}$ FDR'], c='red', label='FDR < 0.1', zorder=2)\n",
    "\n",
    "# Customize plot settings\n",
    "plt.title(\"Signed proportion extrinsic divergence\")\n",
    "ax.set_xlabel('Log$_{2}$ fold-enrichment for TF targets')\n",
    "ax.set_ylabel('-Log$_{10}$ FDR')\n",
    "\n",
    "# Repeat similar steps for other subplots (7 Day, 14 Day, 28 Day)\n",
    "\n",
    "# Show the legend\n",
    "ax.legend(bbox_to_anchor= (1.1, 1))\n",
    "\n",
    "# Display the plot\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d80bec57",
   "metadata": {},
   "outputs": [],
   "source": [
    "out = []\n",
    "df = pd.DataFrame()\n",
    "divergence = \"Extrinsic\"\n",
    "for file in os.listdir(\"Final/Div\"):\n",
    "    v = pd.read_csv(\"Final/Div/\" + file).set_index(\"Unnamed: 0\")\n",
    "    #vv = pd.read_csv(\"Final/NoDiv/\" + file.replace(\"Div\", \"NoDiv\")).set_index(\"Unnamed: 0\")\n",
    "    #v = pd.concat([v, vv])\n",
    "    v[\"Signed \" + divergence.lower()] = np.sign(v[divergence])*v[\"Proportion \" + divergence.lower()]\n",
    "    v = v[[\"Signed \" + divergence.lower()]]\n",
    "    v.columns = [d_ct_abrev[file_to_celltype(file)].replace(\".\", \"_\")]\n",
    "    if ind:\n",
    "        df = v\n",
    "        ind = 0\n",
    "    else:\n",
    "        df = df.join(v, how=\"outer\")\n",
    "\n",
    "df = df.loc[[\"Xbp1\", \"Jun\", \"Atf3\", \"Atf6\"]]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9fed5d4b",
   "metadata": {},
   "outputs": [],
   "source": [
    "#Checking whether Xbp1's divergence correlates with that of these other genes\n",
    "out = []\n",
    "\n",
    "#Compute all-by-all correlation\n",
    "for index, row in df.iterrows():\n",
    "    for index2, row2 in df.iterrows():\n",
    "        if index != index2:\n",
    "            df2 = pd.DataFrame([row, row2]).T.dropna().T\n",
    "            out.append([spearmanr(df2.loc[index], df2.loc[index2])[0], spearmanr(df2.loc[index], df2.loc[index2])[1], index, index2, df2.shape[1]])\n",
    "        else:\n",
    "            df2 = pd.DataFrame([row]).T.dropna().T\n",
    "            out.append([spearmanr(df2.loc[index], df2.loc[index2])[0], spearmanr(df2.loc[index], df2.loc[index2])[1], index, index2, df2.shape[1]])\n",
    "\n",
    "dff = pd.DataFrame(out).dropna()\n",
    "dff.sort_values(1)\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "75fd9569",
   "metadata": {},
   "outputs": [],
   "source": [
    "#Reformat into something that can be plotted easily with seaborn\n",
    "\n",
    "cols  = []\n",
    "out = []\n",
    "rows = []\n",
    "for index, row in dff.iterrows():\n",
    "    if row[\"Gene 1\"] not in cols:\n",
    "        cols.append(row[\"Gene 1\"])\n",
    "    if row[\"Gene 2\"] not in rows:\n",
    "        rows.append(row[\"Gene 2\"])\n",
    "out = []\n",
    "for i in range(len(list(set(dff[\"Gene 1\"])))):\n",
    "    out.append(list(np.repeat(0, len(list(set(dff[\"Gene 1\"]))))))\n",
    "for index, row in dff.iterrows():\n",
    "    xind = cols.index(row[\"Gene 1\"])\n",
    "    yind = rows.index(row[\"Gene 2\"])\n",
    "    if xind != yind:\n",
    "        out[xind][yind] = row[\"Spearman's rho\"]\n",
    "    #Added the diagonal correlations that were missing\n",
    "    else:\n",
    "        out[xind][yind] = 1\n",
    "new_df = pd.DataFrame(out)\n",
    "new_df.columns = cols\n",
    "new_df.index = rows"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "701fd573",
   "metadata": {},
   "outputs": [],
   "source": [
    "#Create clustered heat map\n",
    "sns.set(font_scale = 2)\n",
    "sns.set_style(\"white\")\n",
    "g=sns.clustermap(new_df, metric = \"Euclidean\", tree_kws=dict(linewidths=2))\n",
    "colorbar = g.ax_heatmap.collections[0].colorbar\n",
    "colorbar.set_ticks([-0.4, 0, 0.4, 0.8])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "951a1a06",
   "metadata": {},
   "outputs": [],
   "source": [
    "#Explained in more detail in Fig_2A-F.ipynb\n",
    "def bp(gene, file, extra_title = \"\"):\n",
    "    vv = pd.read_csv(file).set_index(\"Unnamed: 0\")\n",
    "    vp = vv.loc[gene]\n",
    "    r1 = [\"Rat-like env.\", \"Mouse\", vp[\"MR1_mi Norm CPM\"]]\n",
    "    r2 = [\"Rat-like env.\", \"Rat\", vp[\"MR1_ri Norm CPM\"]]\n",
    "    r3 = [\"Mouse-like env.\", \"Mouse\", vp[\"RM1_mi Norm CPM\"]]\n",
    "    r4 = [\"Mouse-like env.\", \"Rat\", vp[\"RM1_ri Norm CPM\"]]\n",
    "    r5 = [\"Mouse-like env.\", \"Mouse\", vp[\"RM2_mi Norm CPM\"]]\n",
    "    r6 = [\"Mouse-like env.\", \"Rat\", vp[\"RM2_ri Norm CPM\"]]\n",
    "    r7 = [\"Wildtype\", \"Mouse\", vp[\"WT_mi Norm CPM\"]]\n",
    "    r8 = [\"Wildype\", \"Rat\", vp[\"WT_ri Norm CPM\"]]\n",
    "    if \"Brain\" in file:\n",
    "        r21 = [\"Host\", \"Rat\", vp[\"MR1_ri Norm CPM\"]]\n",
    "        r31 = [\"Host\", \"Mouse\", vp[\"RM1_mi Norm CPM\"]]\n",
    "        tp = pd.DataFrame([r1, r2, r3, r4, r31, r21])\n",
    "    else:\n",
    "        r21 = [\"Host\", \"Rat\", vp[\"MR1_ri Norm CPM\"]]\n",
    "        r31 = [\"Host\", \"Mouse\", vp[\"RM2_mi Norm CPM\"]]\n",
    "        tp = pd.DataFrame([r31, r21, r1, r2, r6, r5])\n",
    "    tp.columns = [\"Chimera (Donor-Host)\", \"Cell species\", \"Norm CPM\"]\n",
    "    #fig, ax = plt.subplots(figsize = (6, 4.5))\n",
    "    sns.set(font_scale=1.5)\n",
    "    sns.set_style(\"white\")\n",
    "    sns.barplot(data = tp, y = \"Norm CPM\", x = \"Chimera (Donor-Host)\", hue = \"Cell species\", palette = {\"Mouse\":mouse, \"Rat\":rat}, alpha = 1, order = [\"Host\", \"Rat-like env.\", \"Mouse-like env.\"])\n",
    "    plt.ylabel(\"Pseudobulked normalized counts\")\n",
    "    #plt.xlabel(\"Chimera (donor\" + r'$\\rightarrow$' + \"host)\")\n",
    "    plt.xlabel(\"Extrinsic environment\")\n",
    "    if not extra_title:\n",
    "        plt.title(\"Expression of $\\it{\" + gene + \"}$ in \" + d_ct_abrev[file_to_celltype(file)])\n",
    "    else:\n",
    "        plt.title(extra_title + \"\\n$\\it{\" + gene + \"}$ in \" + d_ct_abrev[file_to_celltype(file)])\n",
    "    xmin, xmax, ymin, ymax = plt.axis()\n",
    "    #plt.ylim(0, 1)\n",
    "    plt.axvline(0.5, color = \"black\", linewidth = 2.5, alpha = 1, linestyle=\"dashed\")\n",
    "    #plt.legend(bbox_to_anchor=(1.375, 1.05))\n",
    "    plt.legend([],[], frameon=False)\n",
    "    plt.show()\n",
    "    return vp\n",
    "\n",
    "file = \"Final/Div/Brain_Or_2010_Div_New4_NewNorm_Glutamatergic_progenitors.csv\"\n",
    "gene = \"Xbp1\"\n",
    "bp(gene, file)\n",
    "file = \"Final/Div/Mesenchymal_Or_2010_Div_New4_NewNorm_Mesenchyme_0.csv\"\n",
    "gene = \"Xbp1\"\n",
    "bp(gene, file)\n",
    "file = \"Final/Div/Mesenchymal_Or_2010_Div_New4_NewNorm_Mesenchyme_0.csv\"\n",
    "gene = \"Pomp\"\n",
    "bp(gene, file)\n",
    "file = \"Final/Div/Brain_Or_2010_Div_New4_NewNorm_Spinal_glutamatergic_neurons.csv\"\n",
    "gene = \"Xbp1\"\n",
    "bp(gene, file)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cedbfb66",
   "metadata": {},
   "outputs": [],
   "source": [
    "#Plot correlations between targets and TFs\n",
    "#Can change divergence and \"find\" variables as needed\n",
    "#Can replace scatterplot with regplot for XBP1\n",
    "\n",
    "find_xbp1_up = \"XBP1 OE MOUSE GSE46178 CREEDSID GENE 1437 UP\"\n",
    "find_xbp1_down = \"XBP1 KO MOUSE GSE40273 CREEDSID GENE 1060 DOWN\"\n",
    "find_stat3 = \"STAT3 DEFICIENCY MOUSE GSE6846 CREEDSID GENE 85 UP\"\n",
    "find_nrf1 = \"NRF1 KO MOUSE GSE35124 CREEDSID GENE 2150 DOWN\"\n",
    "find_nfe2l2 = \"NFE2L2 KO MOUSE GSE18344 CREEDSID GENE 966 UP\"\n",
    "find_atf6 = \"ATF6 INDUCED MOUSE GSE8322 CREEDSID GENE 1473 UP\"\n",
    "\n",
    "#Set find\n",
    "find = find_nrf1\n",
    "\n",
    "#Set gene of interest\n",
    "gene = \"Nfe2l1\"\n",
    "gene_caps = gene.upper()\n",
    "\n",
    "#Set type of divergence to plot\n",
    "divergence = \"Intrinsic\"\n",
    "df = pd.DataFrame()\n",
    "sns.set(font_scale = 1.3)\n",
    "sns.set_style(\"white\")\n",
    "ind = 1\n",
    "\n",
    "#Creates table of all estimates across all cell types\n",
    "for file in os.listdir(\"Final/Div\"):\n",
    "    #Read in normalized counts and divergence estimates\n",
    "    v = pd.read_csv(\"Final/Div/\" + file).set_index(\"Unnamed: 0\")\n",
    "    #vv = pd.read_csv(\"Final/NoDiv/\" + file.replace(\"Div\", \"NoDiv\")).set_index(\"Unnamed: 0\")\n",
    "    #v = pd.concat([v, vv])\n",
    "    \n",
    "    #Compute the siggned divergence\n",
    "    v[\"Signed \" + divergence.lower()] = np.sign(v[divergence])*v[\"Proportion \" + divergence.lower()]\n",
    "    v = v[[\"Signed \" + divergence.lower()]]\n",
    "    v.columns = [d_ct_abrev[file_to_celltype(file)].replace(\".\", \"_\")]\n",
    "    if ind:\n",
    "        df = v\n",
    "        ind = 0\n",
    "    else:\n",
    "        df = df.join(v, how=\"outer\")\n",
    "\n",
    "out = []\n",
    "\n",
    "#Get the enrichment score for the TF category of interest\n",
    "for folder in os.listdir(\"Final/GSEAPY_New4_Signed_Proportion_TF_\" + divergence.lower()):\n",
    "    v = pd.read_csv(\"Final/GSEAPY_New4_Signed_Proportion_TF_\" + divergence.lower() + \"/\" + folder + \"/\" + \"gseapy.gene_set.prerank.report.csv\").set_index(\"Term\")\n",
    "    #val = -np.log10(v.loc[find][\"FDR q-val\"])\n",
    "    val = v.loc[find][\"ES\"]\n",
    "    out.append([folder, val])\n",
    "    \n",
    "#Create new data frame with divergence estimate for TF and TF category enrichment\n",
    "df_new = pd.DataFrame(out)\n",
    "df_new = df_new.set_index(0)\n",
    "df_new.index = [x.replace(\"Brain_Or_2010_Div_New4_NewNorm_\", \"\").replace(\"Mesenchymal_Or_2010_Div_New4_NewNorm_\", \"\").replace(\"Chondrocyte_Or_2010_Div_New4_NewNorm_\", \"\").replace(\".csv\", \"\").replace(\"_TF_Perturbations_Followed_by_Expression\", \"\") for x in list(df_new.index)]\n",
    "df_new.columns = [gene_caps + \" enrichment\"]\n",
    "\n",
    "#Create scatterplot\n",
    "from scipy.stats import spearmanr\n",
    "sns.scatterplot(x=pd.DataFrame(df.loc[gene]).join(df_new).dropna()[gene], y=pd.DataFrame(df.loc[gene]).join(df_new).dropna()[gene_caps + \" enrichment\"], color = intrinsic)\n",
    "plt.xlabel(\"$\\it{\" + gene + \"}$\" + \" signed \" + divergence.lower() + \" proportion\")\n",
    "plt.ylabel(\"$\\it{\" + gene + \"}$\" + \" target \" + divergence.lower() + \" enrichment score\")\n",
    "plt.ylim(-0.6, 0.6)\n",
    "print(spearmanr(pd.DataFrame(df.loc[gene]).join(df_new).dropna()[gene], pd.DataFrame(df.loc[gene]).join(df_new).dropna()[gene_caps + \" enrichment\"]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4f6d44f5",
   "metadata": {},
   "outputs": [],
   "source": [
    "#Correlate with ER stress response gene expression\n",
    "#Nearly identical to the above\n",
    "find_xbp1_up = \"XBP1 OE MOUSE GSE46178 CREEDSID GENE 1437 UP\"\n",
    "find_xbp1_down = \"XBP1 KO MOUSE GSE40273 CREEDSID GENE 1060 DOWN\"\n",
    "term = \"Response To Endoplasmic Reticulum Stress (GO:0034976)\"\n",
    "find = find_xbp1_up\n",
    "gene = \"Xbp1\"\n",
    "gene_caps = gene.upper()\n",
    "divergence = \"Extrinsic\"\n",
    "df = pd.DataFrame()\n",
    "sns.set(font_scale = 1.3)\n",
    "sns.set_style(\"white\")\n",
    "ind = 1\n",
    "\n",
    "for file in os.listdir(\"Final/Div\"):\n",
    "    v = pd.read_csv(\"Final/Div/\" + file).set_index(\"Unnamed: 0\")\n",
    "    #vv = pd.read_csv(\"Final/NoDiv/\" + file.replace(\"Div\", \"NoDiv\")).set_index(\"Unnamed: 0\")\n",
    "    #v = pd.concat([v, vv])\n",
    "    v[\"Signed \" + divergence.lower()] = np.sign(v[divergence])*v[\"Proportion \" + divergence.lower()]\n",
    "    v = v[[\"Signed \" + divergence.lower()]]\n",
    "    v.columns = [d_ct_abrev[file_to_celltype(file)].replace(\".\", \"_\")]\n",
    "    if ind:\n",
    "        df = v\n",
    "        ind = 0\n",
    "    else:\n",
    "        df = df.join(v, how=\"outer\")\n",
    "\n",
    "out = []\n",
    "for folder in os.listdir(\"Final/GSEAPY_New4_Signed_Proportion_NoXbp1\" + divergence.lower()):\n",
    "    if \"Biological_Process\" in folder:\n",
    "        v = pd.read_csv(\"Final/GSEAPY_New4_Signed_Proportion_NoXbp1\" + divergence.lower() + \"/\" + folder + \"/\" + \"gseapy.gene_set.prerank.report.csv\").set_index(\"Term\")\n",
    "        #val = -np.log10(v.loc[find][\"FDR q-val\"])\n",
    "        val = v.loc[term][\"ES\"]\n",
    "        out.append([folder, val])\n",
    "df_new = pd.DataFrame(out)\n",
    "df_new = df_new.set_index(0)\n",
    "df_new.index = [x.replace(\"Brain_Or_2010_Div_New4_NewNorm_\", \"\").replace(\"Mesenchymal_Or_2010_Div_New4_NewNorm_\", \"\").replace(\"Chondrocyte_Or_2010_Div_New4_NewNorm_\", \"\").replace(\".csv\", \"\").replace(\"_TF_Perturbations_Followed_by_Expression\", \"\").replace(\"_GO_Biological_Process_2023\", \"\") for x in list(df_new.index)]\n",
    "df_new.columns = [gene_caps + \" enrichment\"]\n",
    "\n",
    "from scipy.stats import spearmanr\n",
    "sns.regplot(x=pd.DataFrame(df.loc[gene]).join(df_new).dropna()[gene], y=pd.DataFrame(df.loc[gene]).join(df_new).dropna()[gene_caps + \" enrichment\"], color = extrinsic)\n",
    "plt.xlabel(\"$\\it{\" + gene + \"}$\" + \" signed \" + divergence.lower() + \" proportion\")\n",
    "plt.ylabel(\"Response to ER stress enrichment score\")\n",
    "plt.title(\"Signed propotion \" + divergence.lower() + \" divergence\")\n",
    "print(spearmanr(pd.DataFrame(df.loc[gene]).join(df_new).dropna()[gene], pd.DataFrame(df.loc[gene]).join(df_new).dropna()[gene_caps + \" enrichment\"]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "60eb5a24",
   "metadata": {},
   "outputs": [],
   "source": [
    "#Count number of significant TFs\n",
    "v_dedup = pd.read_csv(\"Final/GSEAPY_Signed_Proportion_Mean_Signed_TF_\" + divergence + \"/All_TF_TF_Perturbations_Followed_by_Expression/Dedupped_table.csv\")\n",
    "\n",
    "#Filter based on q-value\n",
    "v_sig = v_dedup[v_dedup[\"FDR q-val\"] < 0.1]\n",
    "count_this = []\n",
    "for index, row in v_sig.iterrows():\n",
    "    count_this.append(row[\"Term\"].split(\" \")[0])\n",
    "\n",
    "#Remove duplicate TFs\n",
    "count_this = list(set(count_this))\n",
    "print(len(count_this))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "15257a85",
   "metadata": {},
   "outputs": [],
   "source": [
    "#Doing enrichement analysis for signed extrinsic\n",
    "#Negative is higher in mouse, positive is higher in rat\n",
    "#Across cell types, not used\n",
    "for divergence in vals:\n",
    "    try:\n",
    "        os.mkdir(\"Final/GSEAPY_\" + \"Signed_Proportion_\" + divergence.lower())\n",
    "    except:\n",
    "        pass\n",
    "    for gene_set in gene_sets:\n",
    "        for file in os.listdir(\"Final/Div\"):\n",
    "            v = pd.read_csv(\"Final/Div/\" + file, sep = \",\")\n",
    "            v[\"Signed \" + divergence.lower()] = np.sign(v[divergence])*v[\"Proportion \" + divergence.lower()]\n",
    "            v = v.sort_values(\"Signed \" + divergence.lower(), ascending = False)\n",
    "            ranking = v[[\"Unnamed: 0\", \"Signed \" + divergence.lower()]]\n",
    "            out = []\n",
    "            for index, row in ranking.iterrows():\n",
    "                try:\n",
    "                    out.append([d[row[\"Unnamed: 0\"]], row[\"Signed \" + divergence.lower()]])\n",
    "                except:\n",
    "                    pass\n",
    "            ranking = pd.DataFrame(out)\n",
    "            try:\n",
    "                gs.prerank(rnk=ranking, gene_sets=gene_set, threads=4, permutation_num=1000, outdir= 'C:/Users/astar/Chimerism/Final/GSEAPY_New4_Signed_Proportion_TF_' + divergence.lower() + '/' + d_ct_abrev[file_to_celltype(file)].replace(\".\", \"_\") + \"_\" + gene_set, format='png', seed=6, min_size = 10, max_size = 300)\n",
    "            except:\n",
    "                pass"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.11.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}