{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "a5fd5ea5-5545-4717-a109-d8ba2d14f3d5",
   "metadata": {},
   "source": [
    "# Summarize CAV classifiers"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 53,
   "id": "6ba258e4-b037-4c8c-8278-b19986b49cd8",
   "metadata": {},
   "outputs": [],
   "source": [
    "import torch\n",
    "import os\n",
    "import logomaker\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import seaborn as sns\n",
    "\n",
    "from Bio import motifs\n",
    "from collections import defaultdict\n",
    "from glob import glob\n",
    "from matplotlib import gridspec\n",
    "from tqdm import tqdm\n",
    "from sklearn.metrics import precision_recall_fscore_support as score"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "id": "e406da46-e862-4105-a5ad-b98e8478c317",
   "metadata": {
    "editable": true,
    "papermill": {
     "duration": 0.024606,
     "end_time": "2025-04-01T16:59:29.576178",
     "exception": false,
     "start_time": "2025-04-01T16:59:29.551572",
     "status": "completed"
    },
    "slideshow": {
     "slide_type": ""
    },
    "tags": [
     "parameters"
    ]
   },
   "outputs": [],
   "source": [
    "# parameters\n",
    "cav_root_dir = \"output/512bp_regions_insert_3_motifs/cavs/\"\n",
    "\n",
    "motif_meme = \"0_data/motif-clustering-v2.1beta_consensus_pwms.meme\"\n",
    "\n",
    "topN = 30 # how many motifs to plot, ranked by fscore"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "id": "670e20d7-7a19-4e63-9fc4-2382375f5310",
   "metadata": {
    "papermill": {
     "duration": 0.027998,
     "end_time": "2025-04-01T16:59:30.388033",
     "exception": false,
     "start_time": "2025-04-01T16:59:30.360035",
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "tcav_expts = []\n",
    "for cav_dir in os.listdir(cav_root_dir):\n",
    "    concept, control = cav_dir.split(\"_control_\")\n",
    "    tcav_expts.append((concept, control, os.path.join(cav_root_dir, cav_dir )))\n",
    "tcav_expts = pd.DataFrame(tcav_expts, columns=['concept', 'control', 'cav_dir'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "id": "d228f8db-eb69-47eb-8d6f-857668042f3d",
   "metadata": {
    "papermill": {
     "duration": 0.04468,
     "end_time": "2025-04-01T16:59:30.446445",
     "exception": false,
     "start_time": "2025-04-01T16:59:30.401765",
     "status": "completed"
    },
    "scrolled": true,
    "tags": []
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>concept</th>\n",
       "      <th>control</th>\n",
       "      <th>cav_dir</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>AC0001:GATA-PROP:GATA</td>\n",
       "      <td>random_regions</td>\n",
       "      <td>output/512bp_regions_insert_3_motifs/cavs/AC00...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>AC0002:PROP-ALX:Homeodomain</td>\n",
       "      <td>random_regions</td>\n",
       "      <td>output/512bp_regions_insert_3_motifs/cavs/AC00...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>AC0003:HNF1A-HNF1B:Homeodomain</td>\n",
       "      <td>random_regions</td>\n",
       "      <td>output/512bp_regions_insert_3_motifs/cavs/AC00...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>AC0004:ZSCAN:C2H2_ZF</td>\n",
       "      <td>random_regions</td>\n",
       "      <td>output/512bp_regions_insert_3_motifs/cavs/AC00...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>AC0005:POU3F-POU1F:Homeodomain,POU</td>\n",
       "      <td>random_regions</td>\n",
       "      <td>output/512bp_regions_insert_3_motifs/cavs/AC00...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>...</th>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>632</th>\n",
       "      <td>AC0633:TBX-MGA::EVX:T-box</td>\n",
       "      <td>random_regions</td>\n",
       "      <td>output/512bp_regions_insert_3_motifs/cavs/AC06...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>633</th>\n",
       "      <td>AC0634:ZNF:C2H2_ZF</td>\n",
       "      <td>random_regions</td>\n",
       "      <td>output/512bp_regions_insert_3_motifs/cavs/AC06...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>634</th>\n",
       "      <td>AC0635:TBX-EOMES:T-box</td>\n",
       "      <td>random_regions</td>\n",
       "      <td>output/512bp_regions_insert_3_motifs/cavs/AC06...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>635</th>\n",
       "      <td>AC0636:SREBF:bHLH</td>\n",
       "      <td>random_regions</td>\n",
       "      <td>output/512bp_regions_insert_3_motifs/cavs/AC06...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>636</th>\n",
       "      <td>AC0637:AHR:bHLH</td>\n",
       "      <td>random_regions</td>\n",
       "      <td>output/512bp_regions_insert_3_motifs/cavs/AC06...</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>637 rows × 3 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "                                concept         control  \\\n",
       "0                 AC0001:GATA-PROP:GATA  random_regions   \n",
       "1           AC0002:PROP-ALX:Homeodomain  random_regions   \n",
       "2        AC0003:HNF1A-HNF1B:Homeodomain  random_regions   \n",
       "3                  AC0004:ZSCAN:C2H2_ZF  random_regions   \n",
       "4    AC0005:POU3F-POU1F:Homeodomain,POU  random_regions   \n",
       "..                                  ...             ...   \n",
       "632           AC0633:TBX-MGA::EVX:T-box  random_regions   \n",
       "633                  AC0634:ZNF:C2H2_ZF  random_regions   \n",
       "634              AC0635:TBX-EOMES:T-box  random_regions   \n",
       "635                   AC0636:SREBF:bHLH  random_regions   \n",
       "636                     AC0637:AHR:bHLH  random_regions   \n",
       "\n",
       "                                               cav_dir  \n",
       "0    output/512bp_regions_insert_3_motifs/cavs/AC00...  \n",
       "1    output/512bp_regions_insert_3_motifs/cavs/AC00...  \n",
       "2    output/512bp_regions_insert_3_motifs/cavs/AC00...  \n",
       "3    output/512bp_regions_insert_3_motifs/cavs/AC00...  \n",
       "4    output/512bp_regions_insert_3_motifs/cavs/AC00...  \n",
       "..                                                 ...  \n",
       "632  output/512bp_regions_insert_3_motifs/cavs/AC06...  \n",
       "633  output/512bp_regions_insert_3_motifs/cavs/AC06...  \n",
       "634  output/512bp_regions_insert_3_motifs/cavs/AC06...  \n",
       "635  output/512bp_regions_insert_3_motifs/cavs/AC06...  \n",
       "636  output/512bp_regions_insert_3_motifs/cavs/AC06...  \n",
       "\n",
       "[637 rows x 3 columns]"
      ]
     },
     "execution_count": 32,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "tcav_expts"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "id": "2c796cee-f4d2-4a46-8078-df5282c672b1",
   "metadata": {
    "editable": true,
    "papermill": {
     "duration": 1.318247,
     "end_time": "2025-04-01T16:59:31.780317",
     "exception": false,
     "start_time": "2025-04-01T16:59:30.462070",
     "status": "completed"
    },
    "scrolled": true,
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "637it [00:11, 53.18it/s]\n"
     ]
    }
   ],
   "source": [
    "dfs = []\n",
    "cavs = defaultdict(list)\n",
    "tcavs = defaultdict(list)\n",
    "\n",
    "fscores = []\n",
    "for item in tqdm(tcav_expts.itertuples()):\n",
    "\n",
    "    prefix_dir = item.cav_dir\n",
    "    \n",
    "    perform = pd.read_table(f\"{prefix_dir}/classifier_perform_on_test.txt\", \n",
    "                            comment=\"#\")\n",
    "\n",
    "    precision, recall, fscore, support = score(perform.Truth, perform.Pred)\n",
    "\n",
    "    fscores.append(np.mean(fscore))\n",
    "\n",
    "    # here predictions are model outputs, so it's usually expected that concepts 0 has a higher probability of being bound\n",
    "    #model_predictions_0 = np.loadtxt(f\"{prefix_dir}/predictions_concept_0.txt\")\n",
    "    #model_predictions_1 = np.loadtxt(f\"{prefix_dir}/predictions_concept_1.txt\")\n",
    "    \n",
    "    classifier_scores = pd.DataFrame({'precision': precision,\n",
    "                  'recall': recall,\n",
    "                  'fscore': fscore,\n",
    "                  'support': support})\n",
    "\n",
    "    dfs.append((np.mean(fscore), support, item.concept, item.control))\n",
    "\n",
    "dfs_merge = pd.DataFrame(dfs, columns=['fscore', 'support', 'concept', 'rep'])\n",
    "\n",
    "#dfs_merge['concept_name'] = dfs_merge.apply(lambda x: f'Random genome with\\n{x[\"concept\"]} inserted' if x['concept_id']==0 else 'Random genome', axis=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "id": "5a9e0a62-86dd-4bad-8a26-c8748742bd0a",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>fscore</th>\n",
       "      <th>support</th>\n",
       "      <th>concept</th>\n",
       "      <th>rep</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>0.570113</td>\n",
       "      <td>[469, 448]</td>\n",
       "      <td>AC0001:GATA-PROP:GATA</td>\n",
       "      <td>random_regions</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>0.627836</td>\n",
       "      <td>[469, 456]</td>\n",
       "      <td>AC0002:PROP-ALX:Homeodomain</td>\n",
       "      <td>random_regions</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>0.624060</td>\n",
       "      <td>[469, 462]</td>\n",
       "      <td>AC0003:HNF1A-HNF1B:Homeodomain</td>\n",
       "      <td>random_regions</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>0.882251</td>\n",
       "      <td>[469, 483]</td>\n",
       "      <td>AC0004:ZSCAN:C2H2_ZF</td>\n",
       "      <td>random_regions</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>0.686810</td>\n",
       "      <td>[469, 432]</td>\n",
       "      <td>AC0005:POU3F-POU1F:Homeodomain,POU</td>\n",
       "      <td>random_regions</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>...</th>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>632</th>\n",
       "      <td>0.683243</td>\n",
       "      <td>[469, 415]</td>\n",
       "      <td>AC0633:TBX-MGA::EVX:T-box</td>\n",
       "      <td>random_regions</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>633</th>\n",
       "      <td>0.637607</td>\n",
       "      <td>[469, 450]</td>\n",
       "      <td>AC0634:ZNF:C2H2_ZF</td>\n",
       "      <td>random_regions</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>634</th>\n",
       "      <td>0.502615</td>\n",
       "      <td>[469, 440]</td>\n",
       "      <td>AC0635:TBX-EOMES:T-box</td>\n",
       "      <td>random_regions</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>635</th>\n",
       "      <td>0.645897</td>\n",
       "      <td>[469, 463]</td>\n",
       "      <td>AC0636:SREBF:bHLH</td>\n",
       "      <td>random_regions</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>636</th>\n",
       "      <td>0.842541</td>\n",
       "      <td>[469, 415]</td>\n",
       "      <td>AC0637:AHR:bHLH</td>\n",
       "      <td>random_regions</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>637 rows × 4 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "       fscore     support                             concept             rep\n",
       "0    0.570113  [469, 448]               AC0001:GATA-PROP:GATA  random_regions\n",
       "1    0.627836  [469, 456]         AC0002:PROP-ALX:Homeodomain  random_regions\n",
       "2    0.624060  [469, 462]      AC0003:HNF1A-HNF1B:Homeodomain  random_regions\n",
       "3    0.882251  [469, 483]                AC0004:ZSCAN:C2H2_ZF  random_regions\n",
       "4    0.686810  [469, 432]  AC0005:POU3F-POU1F:Homeodomain,POU  random_regions\n",
       "..        ...         ...                                 ...             ...\n",
       "632  0.683243  [469, 415]           AC0633:TBX-MGA::EVX:T-box  random_regions\n",
       "633  0.637607  [469, 450]                  AC0634:ZNF:C2H2_ZF  random_regions\n",
       "634  0.502615  [469, 440]              AC0635:TBX-EOMES:T-box  random_regions\n",
       "635  0.645897  [469, 463]                   AC0636:SREBF:bHLH  random_regions\n",
       "636  0.842541  [469, 415]                     AC0637:AHR:bHLH  random_regions\n",
       "\n",
       "[637 rows x 4 columns]"
      ]
     },
     "execution_count": 34,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "dfs_merge"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "id": "26c409db-78ea-4336-829f-08fefa01ff03",
   "metadata": {},
   "outputs": [],
   "source": [
    "if topN is not None:\n",
    "    df_plot = dfs_merge.sort_values('fscore', ascending=False).iloc[:topN]\n",
    "else:\n",
    "    df_plot = dfs_merge.sort_values('fscore', ascending=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 73,
   "id": "3b862b20-5cf1-4c7c-aa1e-d192c3ecb742",
   "metadata": {},
   "outputs": [],
   "source": [
    "motif_pwms = motifs.parse(open(motif_meme), fmt='MINIMAL')\n",
    "motif_pwms = {m.name: m.pwm for m in motif_pwms} # index pwm by name\n",
    "def seq_logo(concept, ax=None, max_len=20):\n",
    "    pwm = motif_pwms.get(concept)\n",
    "\n",
    "    if pwm is not None:    \n",
    "        pwm_df = pd.DataFrame(pwm)\n",
    "        motif_len = len(pwm_df)\n",
    "        \n",
    "        # Logomaker expects columns as A,C,G,T, so ensure correct order\n",
    "        pwm_df = pwm_df[['A', 'C', 'G', 'T']]\n",
    "    \n",
    "        # Compute information content at each position (2 - entropy)\n",
    "        def compute_ic(pwm_row):\n",
    "            entropy = -sum([p * np.log2(p) if p > 0 else 0 for p in pwm_row])\n",
    "            return 2 - entropy\n",
    "        \n",
    "        # Compute IC matrix: IC_letter = p * IC_total\n",
    "        ic_df = pwm_df.copy()\n",
    "        for i in range(len(pwm_df)):\n",
    "            ic_total = compute_ic(pwm_df.iloc[i])\n",
    "            ic_df.iloc[i] = pwm_df.iloc[i] * ic_total\n",
    "    \n",
    "        \n",
    "        # Plot with logomaker\n",
    "        x0, y0, width, height = ax.get_position().bounds\n",
    "        ax.set_position([x0, y0, width * min(1., motif_len/max_len), height])\n",
    "        logo = logomaker.Logo(ic_df, color_scheme={'A': 'red', 'C': 'blue',\n",
    "                                                  'G': 'orange', 'T': 'green'},\n",
    "                              ax=ax)\n",
    "        logo.ax.axis('off')\n",
    "    else:\n",
    "        print(f\"Can't find PWM for concept {concept}!\")\n",
    "        ax.axis('off')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 74,
   "id": "fdc7a9f8-8b7c-4bb2-93c5-fd61cd77a43a",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Can't find PWM for concept AC0161:E2F-ZBTB:E2F!\n"
     ]
    },
    {
     "data": {
      "image/png": "[encoded data, 230400 characters]",
      "text/plain": [
       "<Figure size 640x480 with 31 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig = plt.figure()\n",
    "\n",
    "ax_bar = sns.barplot(data=df_plot, y='concept', x='fscore', orient=\"y\")\n",
    "ax_bar.set_position([0, 0, 0.7, 1])\n",
    "ax_bar_bbox = ax_bar.get_position()\n",
    "\n",
    "gs3 = gridspec.GridSpec(len(df_plot), 1)\n",
    "\n",
    "logo_height = ax_bar_bbox.height/len(df_plot)\n",
    "logo_width = 0.25\n",
    "for i, (concept, g) in enumerate(zip(reversed(df_plot.concept.tolist()), gs3)):\n",
    "    ax3 = plt.subplot(g)\n",
    "    ax3.set_position([0.71, ax_bar_bbox.y0+i*logo_height, logo_width+0.01, logo_height])\n",
    "    if concept is not None:\n",
    "        seq_logo(concept, ax=ax3)\n",
    "    else:\n",
    "        ax3.axis('off')"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "torch_latest",
   "language": "python",
   "name": "torch_latest"
  },
  "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.12.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
