final_analysis.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "import pylcs as LCS"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Import cell data\n",
    "cell_data = pd.read_csv('data/cells_no_repeats.csv', index_col=0)\n",
    "print(cell_data.shape)\n",
    "display(cell_data.head())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Import sequence data\n",
    "sequences_df = pd.read_csv('data/dendritic_sequences/dendritic_sequences_ids.csv', index_col=0)\n",
    "print(sequences_df.shape)\n",
    "display(sequences_df.iloc[:10])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import math\n",
    "import string\n",
    "\n",
    "# Convert a pt_root_id index to a string of printable characters\n",
    "def map_to_printable_string(number, length):\n",
    "    if 0 <= number <= 56208:\n",
    "        printable_chars = string.printable[:-6]  # Exclude non-printable characters\n",
    "        base = len(printable_chars)\n",
    "        \n",
    "        encoded = \"\"\n",
    "        while number > 0:\n",
    "            number, index = divmod(number, base)\n",
    "            encoded = printable_chars[index] + encoded\n",
    "        \n",
    "        # Pad with leading zeros if needed\n",
    "        padding = length - len(encoded)\n",
    "        encoded = printable_chars[0] * padding + encoded\n",
    "        \n",
    "        return encoded\n",
    "    else:\n",
    "        raise ValueError(\"Number must be in the range 0 to 56208\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_printable(encoded_str, encoding_map):\n",
    "    return \"\".join([encoding_map[c] for c in encoded_str])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pt_root_ids = cell_data['pt_root_id'].unique()\n",
    "\n",
    "max_value = len(pt_root_ids) - 1\n",
    "num_unique_chars = len(string.printable[:-6])  # Excluding non-printable characters\n",
    "min_string_length = math.ceil(math.log(max_value + 1, num_unique_chars))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create dictionary mapping pt_root_id to a single unique unicode character\n",
    "pt_root_id_to_char = {pt_root_id: chr(i) for i, pt_root_id in enumerate(pt_root_ids)}\n",
    "pt_root_id_to_char[-1] = '' # Add empty string to replace padding\n",
    "\n",
    "char_to_pt_root_id = {v: k for k, v in pt_root_id_to_char.items()}\n",
    "print(f'Number of unique pt_root_ids: {len(pt_root_ids)}')\n",
    "\n",
    "# Create dictionary mapping pt_root_id char encoding to printable string\n",
    "char_encoding_to_printable_string = {pt_root_id_to_char[pt_root_id]: '<' + map_to_printable_string(i, min_string_length) + '>' for i, pt_root_id in enumerate(pt_root_ids)}\n",
    "printable_string_to_pt_root_id = {v: char_to_pt_root_id[k] for k, v in char_encoding_to_printable_string.items()}\n",
    "print(f'Minimum string length: {min_string_length}')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Test encoding\n",
    "test_id = 864691134884741370\n",
    "print(f'pt_root_id: {test_id}')\n",
    "print(f'Size of unique Encoding: {len(pt_root_id_to_char[test_id])}')\n",
    "print(f'Encoded: {char_encoding_to_printable_string[pt_root_id_to_char[test_id]]}')\n",
    "print(f'Decoded: {printable_string_to_pt_root_id[char_encoding_to_printable_string[pt_root_id_to_char[test_id]]]}')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Pre-process sequences: convert from pandas dataframe to list of strings\n",
    "sequences = []\n",
    "sequences_printable = []\n",
    "for _, sequence_row in sequences_df.iterrows():\n",
    "    cur_sequence = ''\n",
    "    for rank in sequences_df.columns[1:]:\n",
    "        cur_sequence += pt_root_id_to_char[sequence_row[rank]]\n",
    "    sequences.append(cur_sequence)\n",
    "    sequences_printable.append(get_printable(cur_sequence, char_encoding_to_printable_string))\n",
    "    print(_, len(cur_sequence), get_printable(cur_sequence, char_encoding_to_printable_string))\n",
    "\n",
    "sequences = np.array(sequences)\n",
    "sequences_printable = np.array(sequences_printable)\n",
    "print(len(sequences))\n",
    "        "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Sort sequences by sequence length\n",
    "sequences_sorted = sorted(sequences, key=len)\n",
    "sequences_printable_sorted = sorted(sequences_printable, key=len)\n",
    "print(sequences_printable_sorted[:10])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Find the index of each sequence length increase\n",
    "# Each element of sequence_length_indexes is the index in sequences where the first sequence of a new length starts\n",
    "sequence_length_indexes = np.array([0])\n",
    "remaining_indexes = np.where(np.array([len(s) - len(sequences_sorted[i-1]) for i, s in enumerate(sequences_sorted)])[1:] > 0)[0] + 1\n",
    "sequence_length_indexes = np.append(sequence_length_indexes, remaining_indexes)\n",
    "\n",
    "print(sequence_length_indexes)\n",
    "print(len(sequence_length_indexes))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "test_a = 'abcd'\n",
    "test_b = ['abc', 'abcd', 'abcde', 'abcdef']\n",
    "print(max(LCS.lcs_string_of_list(test_a, test_b)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "all_distributions = []\n",
    "\n",
    "for k, length_index in enumerate(sequence_length_indexes):\n",
    "    kth_distribution = []\n",
    "\n",
    "    # Disregard sequences shorter than the current length\n",
    "    cur_sequences = sequences_sorted[length_index:]\n",
    "\n",
    "    # Single out the sequences of exactly the current length\n",
    "    if k != len(sequence_length_indexes) - 1:\n",
    "        short_sequences = sequences_sorted[length_index:sequence_length_indexes[k+1]]\n",
    "    else:\n",
    "        short_sequences = sequences_sorted[length_index:]\n",
    "    \n",
    "    # Get the LCS for each sequence of the current length\n",
    "    for short_sequence in short_sequences:\n",
    "        kth_distribution.append(max(LCS.lcs_string_of_list(short_sequence, cur_sequences)))\n",
    "        print(kth_distribution[-1])\n",
    "        break\n",
    "    break\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "microns",
   "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.8.10"
  },
  "orig_nbformat": 4
 },
 "nbformat": 4,
 "nbformat_minor": 2
}