00.convert_tree_data_to_event_matrices.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import ete3\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "from evcouplings.align import Alignment\n",
    "from collections import defaultdict\n",
    "from copy import deepcopy\n",
    "import time\n",
    "import random\n",
    "from scipy.stats import mode\n",
    "#import sparse\n",
    "from scipy import sparse\n",
    "import seaborn as sns\n",
    "from Bio import SeqIO, SeqRecord, Seq\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Converts phylogeny and ancestral sequence reconstruction data into matrices of mutation events "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "\"\"\"\n",
    "Classes and functions for determining pairs of homoplasies on a phylogenetic tree\n",
    "\n",
    "Authors:\n",
    "  Anna G. Green\n",
    "\"\"\"\n",
    "\n",
    "def name_branches(t1, lineage):\n",
    "    \"\"\"\n",
    "    Assigns uniquen ames to each branch in forumlation: {lineage}-{anc_node}-{des_node}\n",
    "\n",
    "    Parameters\n",
    "    ----------\n",
    "    t1: ete3.Tree\n",
    "        the phylogeny of strains\n",
    "    lineage: int or str\n",
    "        name to prepend to branch names\n",
    "\n",
    "    Returns\n",
    "    -------\n",
    "    list of str\n",
    "\n",
    "    \"\"\"\n",
    "    branches = []\n",
    "\n",
    "    for node in t1.traverse(\"levelorder\"):\n",
    "        if len(node.children) > 1:\n",
    "\n",
    "            for idx, child_node in enumerate(node.children):\n",
    "                name_anc = node.name.split(\"_\")[0]\n",
    "                name_des = child_node.name.split(\"_\")[0]\n",
    "\n",
    "                branch_name = f\"{lineage}-{name_anc}-{name_des}\"\n",
    "                branches.append(branch_name)\n",
    "\n",
    "    return branches\n",
    "\n",
    "def get_ancestral_sequence(ancestral_seq_file, all_positions):\n",
    "    \"\"\"\n",
    "    Fetches the ancestral sequence of H37Rv for all positions in input array\n",
    "    all_positions MUST be in H37Rv numbering (same as ancestral seq file)\n",
    "\n",
    "    Parameters\n",
    "    ----------\n",
    "    ancestral_seq_file: str\n",
    "        path to file to input\n",
    "    all_positions: list of int\n",
    "        list of positions to grab from input\n",
    "\n",
    "    Returns\n",
    "    np.array of str\n",
    "    \"\"\"\n",
    "    # subtract 1 or get wrecked. H37Rv positions are 1-indexed\n",
    "    all_positions = [int(x) - 1 for x in all_positions]\n",
    "    \n",
    "    ancestral_seq = Alignment.from_file(open(ancestral_seq_file, \"r\"))\n",
    "    \n",
    "    ancestral_seq_subset = np.array([ancestral_seq.matrix[0,x] for x in all_positions])\n",
    "    \n",
    "    return ancestral_seq_subset\n",
    "\n",
    "class HomoplasyAnalysis:\n",
    "    \"\"\"\n",
    "    Class for determining homoplasic events, and calculating matrices of \n",
    "    mutation and reversion events as well as genetic backgrounds\n",
    "    \"\"\"\n",
    "\n",
    "    def __init__(self,mutation_event_table, phylo, lineage, ancestral_seq_file, phylo_aln):\n",
    "        \"\"\"\n",
    "        Parameters\n",
    "        ----------\n",
    "        mutation_event_table: pd.DataFrame\n",
    "            has column Position which contains the 1-indexed genomic positions analyzed with SNPPar\n",
    "        phylo: ete3.Tree\n",
    "            tree output by SNPPar with internal nodes labelled and matching the names in \n",
    "            phylo_aln\n",
    "        lineage: str\n",
    "            str with lineage name, prepended to branch names in output files\n",
    "        ancestral_seq_file: str\n",
    "            path to ancestral MTB sequence\n",
    "        phylo_aln: Evcouplings.align.Alignment\n",
    "            Alignment object with names matching phylo tree\n",
    "        \"\"\"\n",
    "\n",
    "        self.mutation_table = mutation_event_table\n",
    "        self.positions = sorted(list(mutation_table.Position.unique()))\n",
    "        self.pos_to_index = {pos:idx for idx,pos in enumerate(self.positions)}\n",
    "\n",
    "        self.tree = phylo\n",
    "        self.lineage_label = lineage\n",
    "\n",
    "        self.branches = name_branches(self.tree, self.lineage_label)\n",
    "        self.idx_to_branch = {idx:x for idx, x in enumerate(self.branches)}\n",
    "        \n",
    "        self.aln = phylo_aln\n",
    "\n",
    "        self.ancestral_seq = get_ancestral_sequence(ancestral_seq_file, self.positions)\n",
    "        \n",
    "\n",
    "        self.mutation_event_matrix = None\n",
    "        self.reversion_event_matrix = None\n",
    "        self.background_matrix = None\n",
    "        \n",
    "        # Make sure that ancestral call listed in SNPPar table matches the inferred\n",
    "        # ancestral sequence produced by SNPPar\n",
    "        indices_to_keep = []\n",
    "        for idx,row in self.mutation_table.iterrows():\n",
    "            if row.Ancestor_Call == self.aln.matrix[self.aln.id_to_index[row.Ancestor_Node], \n",
    "                                                    self.pos_to_index[row.Position]]:\n",
    "                indices_to_keep.append(idx)\n",
    "                \n",
    "        if (len(self.mutation_table) - len(indices_to_keep)) > 0:\n",
    "            print(f\"had to eliminate {len(self.mutation_table) - len(indices_to_keep)} rows due to outside error\")\n",
    "        \n",
    "        self.mutation_table = self.mutation_table.loc[indices_to_keep,:]\n",
    "\n",
    "    def filter_positions(self, new_positions):\n",
    "        \"\"\"\n",
    "        Parameters\n",
    "        ----------\n",
    "        new_positions: list of int\n",
    "        \"\"\"\n",
    "\n",
    "        if not (set(new_positions).issubset(set(self.positions))):\n",
    "            raise ValueError(\"positions provided are not a subset of positions in the model\")\n",
    "\n",
    "        new_positions = sorted(new_positions)    \n",
    "        \n",
    "        new_pos_to_idx = {pos:idx for idx,pos in enumerate(new_positions)}\n",
    "        indices_to_keep = [self.pos_to_index[x] for x in new_positions]\n",
    "\n",
    "        self.positions = new_positions\n",
    "        self.pos_to_index = new_pos_to_idx\n",
    "\n",
    "        self.ancestral_seq = self.ancestral_seq[indices_to_keep]\n",
    "        \n",
    "        self.mutation_table = self.mutation_table.query(\"Position in @new_positions\")\n",
    "        \n",
    "        print(self.aln.matrix.shape)\n",
    "        self.aln = self.aln.select(columns=indices_to_keep)\n",
    "        print(self.aln.matrix.shape)\n",
    "        \n",
    "        if self.mutation_event_matrix is not None:\n",
    "            self.mutation_event_matrix = self.mutation_event_matrix[:,indices_to_keep]\n",
    "\n",
    "        if self.reversion_event_matrix is not None:\n",
    "            self.reversion_event_matrix = self.reversion_event_matrix[:,indices_to_keep]\n",
    "\n",
    "        if self.background_matrix is not None:\n",
    "            self.background_matrix = self.background_matrix[:,indices_to_keep]\n",
    "\n",
    "    def create_mutation_matrix(self, mode=\"mutation\"):\n",
    "        \"\"\"\n",
    "        Creates a numpy matrix where each branch gets a row, and each position gets a column.\n",
    "        Contains a 1 if a mutation (or reversion, in reversion mode) happened on that branch\n",
    "\n",
    "        \"\"\"\n",
    "\n",
    "        #initialize mutation event matrix\n",
    "        mutation_event_matrix = np.zeros((len(self.branches), len(self.positions)),dtype=int)\n",
    "\n",
    "        for _,row in self.mutation_table.iterrows():\n",
    "\n",
    "            # look up row in mutation event matrix based on branch name\n",
    "            branch_name = f\"{self.lineage_label}-{row.Ancestor_Node}-{row.Derived_Node}\"\n",
    "            branch_index = self.branches.index(branch_name)\n",
    "\n",
    "            # look up column in mutation eevn matrix based on position\n",
    "            position_index = self.pos_to_index[row.Position]\n",
    "\n",
    "            # function can either find all mutations AWAY FROM reference\n",
    "            # or all mutations TO reference (reversions)\n",
    "            use = True\n",
    "            if mode==\"mutation\":\n",
    "                if row.Derived_Call == self.ancestral_seq[position_index]:\n",
    "                    use=False\n",
    "            elif mode==\"reversion\":\n",
    "                if row.Derived_Call != self.ancestral_seq[position_index]:\n",
    "                    use=False\n",
    "            else:\n",
    "                print(f\"select either mutation or reversion for mode, you selected {mode}\")\n",
    "                return None\n",
    "            if use:\n",
    "                mutation_event_matrix[branch_index, position_index] = 1\n",
    "\n",
    "        return mutation_event_matrix\n",
    "\n",
    "    def create_background_matrix(self):\n",
    "\n",
    "            \"\"\"\n",
    "            Creates a numpy matrix where each branch gets a row, and each position gets a column.\n",
    "            Contains a 0 if the ancestor node of that branch had the ancestral state, and a 1 otherwise\n",
    "\n",
    "            \"\"\"\n",
    "\n",
    "            # create a dictionary of ancestral node to branches descending from that node\n",
    "            node_to_branch_names = defaultdict(list)\n",
    "            for x in self.branches:\n",
    "                node_to_branch_names[x.split(\"-\")[1]].append(x)\n",
    "\n",
    "            #initialize mutation event matrix\n",
    "            mutation_background_matrix = np.zeros((len(self.branches), len(self.positions)),dtype=int)\n",
    "\n",
    "            # For each node in the tree\n",
    "            for node in self.tree.traverse(\"levelorder\"):\n",
    "                # ignore if the node is a leaf\n",
    "                if len(node.children) < 1:\n",
    "                    continue\n",
    "\n",
    "                # get the name of the ancestral node and its sequence\n",
    "                name_anc = node.name.split(\"_\")[0]\n",
    "                seq_anc = self.aln.matrix[self.aln.id_to_index[name_anc], :]\n",
    "\n",
    "                # find positions where different from LUCA sequence\n",
    "                differences = (self.ancestral_seq != seq_anc).astype(int)\n",
    "\n",
    "                # for every branch descending from this node\n",
    "                for b in node_to_branch_names[name_anc]:\n",
    "\n",
    "                    # save positions where there's a difference as 1, else 0\n",
    "                    mutation_background_matrix[self.branches.index(b)] = differences\n",
    "\n",
    "            return mutation_background_matrix\n",
    "\n",
    "    def create_background_matrix_new(self):\n",
    "\n",
    "            \"\"\"\n",
    "            Creates a numpy matrix where each branch gets a row, and each position gets a column.\n",
    "            Contains a 0 if the DESCENDANT node of that branch had the ancestral state and no mutation\n",
    "            ocurred, and a 1 otherwise\n",
    "\n",
    "            NB this was corrected in February 2022 from previous versions - this is an attempt to patch SNPPar's \n",
    "            mislabelling of some nodes in the tree\n",
    "\n",
    "            \"\"\"\n",
    "\n",
    "            # create a dictionary of descendant nodes to ancestral branches leading to that node\n",
    "            node_to_branch_names = defaultdict(list)\n",
    "            for x in self.branches:\n",
    "                print(x)\n",
    "                node_to_branch_names[x.split(\"-\")[-1]].append(x)\n",
    "\n",
    "            #initialize mutation event matrix\n",
    "            mutation_background_matrix = np.zeros((len(self.branches), len(self.positions)),dtype=int)\n",
    "\n",
    "            # For each node in the tree\n",
    "            for node in self.tree.traverse(\"levelorder\"):\n",
    "\n",
    "                node_name = node.name.split(\"_\")[0]\n",
    "                # ignore if the node has no ancestors\n",
    "                if len(node_to_branch_names[node_name]) < 1:\n",
    "                    print(\"skipping root node\", node_name)\n",
    "                    continue\n",
    "\n",
    "                # get the name of the node and its sequence\n",
    "                print(node.name, node_to_branch_names[node.name])\n",
    "\n",
    "                seq_node = self.aln.matrix[self.aln.id_to_index[node_name], :]\n",
    "\n",
    "                # find positions where descendant is different from LUCA sequence\n",
    "                differences = (self.ancestral_seq != seq_node).astype(int)\n",
    "\n",
    "                # for every branch ancestral to this node\n",
    "                for b in node_to_branch_names[node_name]:\n",
    "                    print(b)\n",
    "                    # save positions where there's a difference as 1, else 0\n",
    "                    mutation_background_matrix[self.branches.index(b)] = differences\n",
    "                break\n",
    "\n",
    "            return mutation_background_matrix\n",
    "\n",
    "    def dump(self, output_path):\n",
    "        \"\"\"\n",
    "        Saves calculated matrices to a given output path\n",
    "        \"\"\"\n",
    "        np.save(output_path+\"_positions.npy\", self.positions)\n",
    "        np.save(output_path+\"_branches.npy\", self.branches)\n",
    "        np.save(output_path+\"_mutation_matrix.npy\", self.mutation_event_matrix)\n",
    "        np.save(output_path+\"_background_matrix.npy\", self.background_matrix)\n",
    "        np.save(output_path+\"_reversion_matrix.npy\", self.reversion_matrix)\n",
    "\n",
    "    def dump_sparse(self, output_path):\n",
    "        \"\"\"\n",
    "        Saves calculated matrices to a given output path using scipy.sparse module\n",
    "        \"\"\"\n",
    "        np.save(output_path+\"_positions.npy\", self.positions)        \n",
    "        np.save(output_path+\"_branches.npy\", self.branches)\n",
    "        \n",
    "        sparse.save_npz(output_path+\"_mutation_matrix.npz\", sparse.coo_matrix(self.mutation_event_matrix))\n",
    "        sparse.save_npz(output_path+\"_background_matrix.npz\", sparse.coo_matrix(self.background_matrix))\n",
    "        sparse.save_npz(output_path+\"_reversion_matrix.npz\", sparse.coo_matrix(self.reversion_matrix))\n",
    "        \n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Clean newick trees\n",
    "# Trees output by SNPpar have info which states the # of events per branch, this confuses ete3 and needs to be removed\n",
    "import re\n",
    "def clean_tree(file):\n",
    "    with open(file) as infile: \n",
    "        string = infile.read()\n",
    "        #print(string)\n",
    "        \n",
    "        substituted = re.sub(\"\\[.*?\\]\",\"\",string)\n",
    "        \n",
    "    outfile = file.split(\".tre\")[0] + \"_cleaned.tre\"\n",
    "    \n",
    "    with open(outfile, 'w') as of:\n",
    "        of.write(substituted)\n",
    "\n",
    "for lineage in [\"L1\", \"L2\", \"L3\", \"L4A\", \"L4B\", \"L4C\", \"L5\", \"L6\"]:\n",
    "    \n",
    "    clean_tree(f\"input/{lineage}node_labelled_newick.tre\")\n",
    "#clean_tree(f\"input/L5node_labelled_newick.tre\")\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "## Cleans mutation tables\n",
    "## re-written to be faster\n",
    "def mutation_events_check_against_seqs(global_lineage):\n",
    "\n",
    "    #Load mutation events DataFrame\n",
    "    mutation_events_DF = pd.read_csv(\n",
    "       f'/n/data1/hms/dbmi/farhat/anna/ABR_epistasis/snppar_repo/reconstruction_treetime07/lineage{global_lineage}/L{global_lineage}all_mutation_events.tsv', sep = '\\t') \n",
    "\n",
    "    #Load inferred sequences at the nodes (inferred)\n",
    "    print(\"reading node alignment\")\n",
    "    ali = Alignment.from_file(open(\n",
    "        f\"/n/data1/hms/dbmi/farhat/anna/ABR_epistasis/snppar_repo/reconstruction_treetime07/lineage{global_lineage}/L{global_lineage}node_sequences.fasta\", \n",
    "    ))\n",
    "    seqs_dict = {}\n",
    "    for id in ali.ids:\n",
    "        seqs_dict[id] = \"\".join(ali.matrix[ali.id_to_index[id],:])\n",
    "\n",
    "    #Load the Reference Position SNP annotation from the SNP table that was used as input\n",
    "    SNP_table = pd.read_csv('/n/data1/hms/dbmi/farhat/Roger/homoplasy_project/rolling_DB_scrape/Genotypes_Filtered_2/SNP_ancestral_reconstruction/lineage{0}/SNP_table.csv'.format(global_lineage) , sep = ',')\n",
    "\n",
    "    print(\"constructing isolate alignment\")\n",
    "    #convert columns of SNP table into a SeqRecord for each isolate (terminal)\n",
    "    for isolate_i in SNP_table.columns[1:]:\n",
    "\n",
    "        isolate_i_seq = SNP_table.loc[: , isolate_i]\n",
    "        isolate_i_seq[isolate_i_seq == '-'] = 'N'\n",
    "        isolate_i_seq = ''.join(isolate_i_seq)\n",
    "        isolate_i_seq_record = SeqRecord.SeqRecord(Seq.Seq(isolate_i_seq ) , id = isolate_i)\n",
    "\n",
    "        #add seq record to dict\n",
    "        seqs_dict[isolate_i] = isolate_i_seq_record\n",
    "\n",
    "    print(\"Concatenating alignment\")\n",
    "    ali = Alignment.from_dict(seqs_dict)     \n",
    "\n",
    "    #create a series from the SNP table that maps each SNP concatenate index to H37Rv Reference Position\n",
    "    SNP_table.loc[: , 'index'] = SNP_table.index\n",
    "\n",
    "    SNP_pos_to_index = {x:idx for idx,x in enumerate(SNP_table.pos)}\n",
    "\n",
    "    #iterate through each mutation in events table and check against fasta seqs\n",
    "    mutations_checked_against_seqs_filter = []\n",
    "\n",
    "    import time\n",
    "\n",
    "    time0=time.time()\n",
    "\n",
    "    #print(len(mutation_events_DF))\n",
    "    count=0\n",
    "    for mutation_event_i, row in mutation_events_DF.iterrows():\n",
    "\n",
    "        #time1=time.time()\n",
    "        #get the Ancestor Call & Derived Call from the inferred/terminal sequences\n",
    "        position = SNP_pos_to_index[row.Position]\n",
    "        #time15 = time.time()\n",
    "        Anc_call_from_seqs = ali.matrix[ali.id_to_index[row.Ancestor_Node], position]\n",
    "        Der_call_from_seqs = ali.matrix[ali.id_to_index[row.Derived_Node], position]\n",
    "        #time2=time.time()\n",
    "        #check calls from sequences against what was reported in the mutation table\n",
    "        if (row.Ancestor_Call == Anc_call_from_seqs) and (row.Derived_Call == Der_call_from_seqs):\n",
    "            mutations_checked_against_seqs_filter.append(True)\n",
    "        else:\n",
    "            mutations_checked_against_seqs_filter.append(False)\n",
    "        #time3=time.time()\n",
    "        if count % 10000 == 0:\n",
    "            print(count)\n",
    "        count +=1\n",
    "\n",
    "    #time4 = time.time()\n",
    "    mutations_checked_against_seqs_filter = np.array(mutations_checked_against_seqs_filter)\n",
    "    print('There were ' + str(len(mutations_checked_against_seqs_filter)) + ' reported events in the mutation table between the inferred node & isolate sequences.')\n",
    "    print('There were ' + str(sum(mutations_checked_against_seqs_filter == False)) + \\\n",
    "          ' reported events in the mutation table that were not located between the inferred node & isolate sequences.')\n",
    "\n",
    "    #filter out the mutation events that were not found among the inferred / isolate sequences and save\n",
    "    mutation_events_DF = mutation_events_DF[mutations_checked_against_seqs_filter]\n",
    "    mutation_events_DF.set_index('Position' , drop = True , inplace = True)\n",
    "    return mutation_events_DF\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 0. Modify MTB ancestor position 2030521 (EsxM) because no extant strain has the ancestral base"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "# aln = Alignment.from_file(open(\"input/MTB_ancestor_reference.fasta\", \"r\"))\n",
    "# aln.matrix[0, 2030521 - 1] # Current value is A, replace with T as this is hte most common genotype\n",
    "\n",
    "# aln.matrix[0, 2030521 - 1] = \"T\"\n",
    "\n",
    "# aln.write(open(\"input/MTB_ancestor_reference_modified.fasta\", \"w\"), format=\"fasta\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "mkdir: cannot create directory ‘figures’: File exists\n"
     ]
    }
   ],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 0. Compare branch lengths from all phylogenies"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "L1\n",
      "L2\n",
      "L3\n",
      "L4A\n",
      "L4B\n",
      "L4C\n",
      "L5\n",
      "L6\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZMAAAEKCAYAAADXdbjqAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAAsTAAALEwEAmpwYAAAjbElEQVR4nO3de5hU1Znv8e8vqOANBUQfpEX0PMSISFBbg1HnYNCB5GSEZCSiJpJoJDGYG5MYnZmMJqMjjg4z0SQYvEzjiYoMySje8EJEJwaVjjJyi0AihD4wcokaEEXB9/xRq7Horr7uququ5vd5nnpq16q99n6rdne/vdbae21FBGZmZll8qKMDMDOzyudkYmZmmTmZmJlZZk4mZmaWmZOJmZll5mRiZmaZ7dXRAZTK6NGjY+7cuR0dhplZpVF7KnXZlsmmTZs6OgQzsz1Gl00mZmZWPk4mZmaWmZOJmZll1mUH4M3MSu29996jrq6Od955p6NDabMePXpQVVXF3nvvXZTtOZmYmbVTXV0dBx54IAMHDkRq10lQHSIi2Lx5M3V1dRx11FFF2aa7uczM2umdd96hT58+FZVIACTRp0+foraonEzMzDKotERSr9hxO5mYmRXZAQccAMC6des499xzOzia8vCYSZGNn76g3XVnTjy1iJGYWUc7/PDDmT17dkeHURZumZiZlcjq1asZMmQIADU1NXz2s59l9OjRDBo0iCuuuGLXeo8//jinnnoqJ554IuPGjWPr1q0A/PCHP+Tkk09myJAhTJw4kfo74y5cuJChQ4dy6qmn8t3vfnfXPnbu3Ml3v/tdTj75ZIYOHcrPfvazsn1WJxMzszJZtGgR9913H4sXL+a+++5j7dq1bNq0iWuvvZYnn3ySF198kerqaqZOnQrA5ZdfzsKFC1myZAlvv/02Dz30EABf+tKXuPXWW1mwYAHdunXbtf077riDgw46iIULF7Jw4UJuu+02Xn311bJ8NndzmZmVyciRIznooIMAGDx4MGvWrOGNN95g2bJlnHbaaQC8++67nHpqrsv7qaee4p//+Z/Ztm0bf/rTnzjuuOM444wz2LJlCx//+McBuOCCC3Ylmccff5yXX355V9fam2++ycqVK4t2+m9znEy6EI/XmHVu3bt337XcrVs3duzYQURw9tlnc++99+627jvvvMPXvvY1amtrOeKII7jmmmt45513dnV1FRIR3HLLLYwaNapkn6Ep7uYyM+tAw4cP59lnn2XVqlUAbNu2jRUrVuy6BuSQQw5h69atu1obvXr14sADD+S5554DYObMmbu2NWrUKKZNm8Z7770HwIoVK3jrrbfK8jncMjEz60B9+/alpqaG888/n+3btwNw7bXX8uEPf5hLL72U448/noEDB3LyySfvqnPHHXdw6aWXsv/++zNixIhdXWdf/vKXWb16NSeeeCIRQd++fbn//vvL8jnUXJOpklVXV0dtbW3Z99uRXU3u5jIrr+XLl3PssceWfb9bt27ddS3LlClTWL9+PT/60Y/avJ0m4m/X1YxumZiZVZiHH36Y66+/nh07dnDkkUdSU1PT0SE5mZiZVZrzzjuP8847r6PD2I2TSZHN3Oe6DLUfKlocZmbl5LO5zMwss5IlE0nHSFqU9/izpG9J6i3pCUkr03OvvDpXSVol6RVJo/LKT5K0OL13syp1mk4zsy6qZMkkIl6JiGERMQw4CdgG/CdwJTAvIgYB89JrJA0GxgPHAaOBn0qqnydgGjARGJQeo0sVt5mZtV25urlGAr+PiDXAGGBGKp8BjE3LY4CZEbE9Il4FVgGnSOoH9IyIBZE7j/muvDpmZnu0+lOE8z3zzDOceOKJ7LXXXmWbtbhcA/Djgfq5Ag6LiPUAEbFe0qGpvD/wXF6dulT2XlpuWN6IpInkWjAMGDCgaMGbmbVGlmu9Cmnv9V8DBgygpqaGm266qajxNKfkyUTSPsA5wFUtrVqgLJopb1wYMR2YDrmLFtsQZpfgM8nMDGDgwIEAfOhD5TvHqhx7+iTwYkS8ll6/lrquSM8bUnkdcERevSpgXSqvKlBuZmadRDmSyfl80MUFMAeYkJYnAA/klY+X1F3SUeQG2l9IXWJbJA1PZ3FdlFfHzMw6gZJ2c0naDzgb+Epe8RRglqRLgD8C4wAiYqmkWcAyYAcwKSJ2pjqXATXAvsCj6WFmZp1ESZNJRGwD+jQo20zu7K5C618HNOr4j4haYEgpYjQzs+w8nUonkvVMkJn7FCkQM6sY27Zto6rqg2HlyZMnc8YZZ/CZz3yG119/nQcffJCrr76apUuXljQOJxMzsyLpiFs5vP/++wXL6+rqCpaXiufmMjOzzJxMzMwsMycTMzPLzMnEzMwyczIxM7PMnEzMzCwzJxMzswpWaAr6qVOnMnjwYIYOHcrIkSNZs2ZNyePwdSZmZsVS8+nibu+L7ZvN+4QTTqC2tpb99tuPadOmccUVV3DfffcVN7YG3DIxM+tizjzzTPbbbz8Ahg8fXpYLGN0y6USy3Y8kmyxTuXTEVb9m1jp33HEHn/zkJ0u+HycTM7Mu6uc//zm1tbU8/fTTJd+Xk4mZWRf05JNPct111/H000/TvXv3ku/PycTMrIt56aWX+MpXvsLcuXM59NBDy7JPJxMzswpWaAr6Rx55hK1btzJu3DgABgwYwJw5c0oah5OJmVmxtPNU3iwKTUE/efLkssfhU4PNzCyzkiYTSQdLmi3pd5KWSzpVUm9JT0hamZ575a1/laRVkl6RNCqv/CRJi9N7N0tSKeM2M7O2KXXL5EfA3Ij4CPBRYDlwJTAvIgYB89JrJA0GxgPHAaOBn0rqlrYzDZgIDEqP0SWO28zM2qBkyURST+AvgDsAIuLdiHgDGAPMSKvNAMam5THAzIjYHhGvAquAUyT1A3pGxIKICOCuvDpmZtYJlLJlcjSwEfh3SS9Jul3S/sBhEbEeID3Xn7fWH1ibV78ulfVPyw3LzcyskyhlMtkLOBGYFhEnAG+RurSaUGgcJJopb7wBaaKkWkm1GzdubGu8ZmbWTqVMJnVAXUQ8n17PJpdcXktdV6TnDXnrH5FXvwpYl8qrCpQ3EhHTI6I6Iqr79u1btA9iZtZZFZqC/tZbb+X4449n2LBhnH766SxbtqzkcZTsOpOI+B9JayUdExGvACOBZekxAZiSnh9IVeYA90iaChxObqD9hYjYKWmLpOHA88BFwC2lintPlW2SyfKfW2/WGV382MVF3d6do+5sV70LLriAr371qwDMmTOHyZMnM3fu3GKG1kipL1r8OnC3pH2APwBfItcamiXpEuCPwDiAiFgqaRa5ZLMDmBQRO9N2LgNqgH2BR9PDzMwK6Nmz567lt956i3JcTVHSZBIRi4DqAm+NbGL964BG/yJHRC0wpKjBmZl1YT/5yU+YOnUq7777Lr/61a9Kvj9fAW9m1gVNmjSJ3//+99xwww1ce+21Jd+fk4mZWRc2fvx47r///pLvx8nEzKyLWbly5a7lhx9+mEGDBpV8n541uMgu7tf+ewfcuX5DyyuZmeUpNAX9mjVrePLJJ9l7773p1asXM2bMaGYLxeFkYmZWJO09lTeLQlPQdwQnE8ts/PQF7a47c+KpRYzEzDqKx0zMzCwzJxMzM8vM3VwFZOm22e/IIgZiZlYh3DIxM7PMnEzMzCwzJxMzswpWaAr6erNnz0YStbW1JY/DYyZmZkWy5qIJRd3ekXe1/2LDLVu2cPPNN/Oxj32siBE1zS0TM7Mu6Pvf/z5XXHEFPXr0KMv+nEzMzLqYl156ibVr1/LpT3+6bPt0N5eZWRfy/vvv8+1vf5uampqy7tctEzOzLmTLli0sWbKEESNGMHDgQJ577jnOOeeckg/Cu2ViZtaFHHTQQWzatGnX6xEjRnDTTTdRXV3oprfF42RSwMx9Gt05uNUupv1T0JuZtVWhKegnT55c9jhKmkwkrQa2ADuBHRFRLak3cB8wEFgNfC4iXk/rXwVcktb/RkQ8lspPAmqAfYFHgG9GRJQy9va6+u6d7a/8ieLFUU5Zki88VLQ4zDpallN526ulKejnz59fljjKMWZyZkQMi4j6NtaVwLyIGATMS6+RNBgYDxwHjAZ+KqlbqjMNmAgMSo/RZYjbzMxaqSMG4McA9el7BjA2r3xmRGyPiFeBVcApkvoBPSNiQWqN3JVXx8zMOoFSJ5MAHpf0W0kTU9lhEbEeID3XDzL0B9bm1a1LZf3TcsPyRiRNlFQrqXbjxo1F/BhmZtacUg/AnxYR6yQdCjwh6XfNrKsCZdFMeePCiOnAdIDq6upOOaZiu/NdGs26hpK2TCJiXXreAPwncArwWuq6Ij1vSKvXAUfkVa8C1qXyqgLlZmbWSZQsmUjaX9KB9cvAXwJLgDlA/WxoE4AH0vIcYLyk7pKOIjfQ/kLqCtsiabgkARfl1TEzs06glC2Tw4BfS/pv4AXg4YiYC0wBzpa0Ejg7vSYilgKzgGXAXGBSRNSfZ3sZcDu5QfnfA4+WMG4zs4pRaAr6mpoa+vbty7Bhwxg2bBi33357yeMo2ZhJRPwB+GiB8s3AyCbqXAc0umghImqBIcWO0cysmO6f+mJRtzd28ontrnveeefx4x//uIjRNM9zc5mZWWZOJmZmXdAvfvELhg4dyrnnnsvatWtbrpCRk4mZWRfzV3/1V6xevZqXX36Zs846iwkTinsHyEKcTMzMupg+ffrQvXt3AC699FJ++9vflnyfnjW4E7m4X7YZh+9cv6Hllcysy1u/fj39+vUDYM6cORx77LEl36eTSReSJRk5EZlVpkJT0G/cuJE5c+aw11570bt377LcddHJxMysSLKcytteTU1Bf/3115c1Do+ZmJlZZk4mZmaWmbu5rEP5Lo1mXYOTiVWsLNPXg6ewNysmd3OZmVlmTiZmZpaZu7kKyHK9xtXsbHklM7MiOeCAA9i6dWuj8lmzZnHNNdcgiY9+9KPcc889JY3DycTMrEhm/eCqom7vc1e371qRlStXcv311/Pss8/Sq1cvNmwo/UXJ7uYyM+tibrvtNiZNmkSvXr0AOPTQbFM1tYaTiZlZF7NixQpWrFjBaaedxvDhw5k7d27J9+luLjOzLmbHjh2sXLmS+fPnU1dXxxlnnMGSJUs4+OCDS7bPkrdMJHWT9JKkh9Lr3pKekLQyPffKW/cqSaskvSJpVF75SZIWp/dulqRSx21mVqmqqqoYM2YMe++9N0cddRTHHHMMK1euLOk+y9HN9U1ged7rK4F5ETEImJdeI2kwMB44DhgN/FRSt1RnGjARGJQeo8sQt5lZRRo7dixPPfUUAJs2bWLFihUcffTRJd1nSZOJpCrg/wC35xWPAWak5RnA2LzymRGxPSJeBVYBp0jqB/SMiAUREcBdeXXMzPZo9VPQ1z+mTp3KqFGj6NOnD4MHD+bMM8/kxhtvpE+fPiWNo9RjJv8GXAEcmFd2WESsB4iI9ZLqTzPoDzyXt15dKnsvLTcsb0TSRHItGAYMGFCE8Mvr6ruzXaPygwu7tbySmZVMe0/lzaKpKeinTp3K1KlTyxZHs8lE0j8083ZExD82U/fTwIaI+K2kEa2IpdA4SDRTXiig6cB0gOrq6oLrmJlZ8bXUMnmrQNl+wJeBPkCTyQQ4DThH0qeAHkBPST8HXpPUL7VK+gH1V9PUAUfk1a8C1qXyqgLl1kCmls0niheHme15mk0mEfEv9cuSDiQ3mH4xMBP4l6bqpbpXAVeluiOA70TE5yXdCEwApqTnB1KVOcA9kqYCh5MbaH8hInZK2iJpOPA8cBFwS9s+Zttk7W6y8sg2fT2Mn/537d+3Zxw2202LYyaSegOTgQvJDZifGBGvZ9jnFGCWpEuAPwLjACJiqaRZwDJgBzApIur/ql8G1AD7Ao+mh5lZh4sIKvFqhdz5TMXT0pjJjcBnyY1DHB8RjWcTa4WImA/MT8ubgZFNrHcd0OjfzYioBYa0Z99mZqXSo0cPNm/eTJ8+fSoqoUQEmzdvpkePHkXbZkstk78BtgN/D/xd3pelXDzRs2iRmJlVmKqqKurq6ti4cWNHh9JmPXr0oKqqquUVW6mlMRPP3bWHyDLt/p3rSz8jqVlnVH+FuXmiRzMzKwInEzMzy8zJxMzMMnMyMTOzzJxMzMwsMycTMzPLzMnEzMwyczIxM7PMfA9422NlmyjyoaLFYdYVuGViZmaZOZmYmVlmTiZmZpaZk4mZmWXmAXizdhg/fUG76/oujdYVOZlYZp6+3szczWVmZpmVrGUiqQfwDNA97Wd2RFyd7il/HzAQWA18rv6e8pKuAi4BdgLfiIjHUvlJfHAP+EeAb0axb2C8h7v67p3trvuDC7sVMRIzq0SlbJlsBz4RER8FhgGjJQ0HrgTmRcQgYF56jaTBwHjgOGA08FNJ9X+lpgETgUHpMbqEcZuZWRuVLJlEztb0cu/0CGAMMCOVzwDGpuUxwMyI2B4RrwKrgFMk9QN6RsSC1Bq5K6+OmZl1AiUdM5HUTdIiYAPwREQ8DxwWEesB0nP96G1/YG1e9bpU1j8tNyw3M7NOoqTJJCJ2RsQwoIpcK2NIM6ur0CaaKW+8AWmipFpJtRs3bmxzvGZm1j5lOZsrIt4A5pMb63gtdV2RnuvPDa0DjsirVgWsS+VVBcoL7Wd6RFRHRHXfvn2L+RHMzKwZpTybqy/wXkS8IWlf4CzgBmAOMAGYkp4fSFXmAPdImgocTm6g/YWI2ClpSxq8fx64CLilVHGbtYZnHDbbXSkvWuwHzEhnZH0ImBURD0laAMySdAnwR2AcQEQslTQLWAbsACZFRP35qpfxwanBj6aHmZl1EiVLJhHxMnBCgfLNwMgm6lwHNPqXLyJqgebGW6xC+ep5s67BV8CbmVlmTiZmZpaZJ3o0KzPPOGxdkVsmZmaWmVsmlpkniTQzJxOrWFnOBAOfDWZWTO7mMjOzzNwysT2Wr3ExKx63TMzMLDMnEzMzy8zJxMzMMnMyMTOzzDwAb1Zmnr7euiK3TMzMLDMnEzMzy8zdXNahPBWLWdfglomZmWXmlolZO/jqebPduWViZmaZlSyZSDpC0lOSlktaKumbqby3pCckrUzPvfLqXCVplaRXJI3KKz9J0uL03s2SVKq4zcys7UrZMtkB/E1EHAsMByZJGgxcCcyLiEHAvPSa9N544DhgNPBTSfUjrNOAicCg9BhdwrjNzKyNSjZmEhHrgfVpeYuk5UB/YAwwIq02A5gPfC+Vz4yI7cCrklYBp0haDfSMiAUAku4CxgKPlip2qwxZzgQDnw1mVkxlGYCXNBA4AXgeOCwlGiJivaT6kcz+wHN51epS2XtpuWF5of1MJNeCYcCAAUX8BNYV+bRks+IpeTKRdADwC+BbEfHnZoY7Cr0RzZQ3LoyYDkwHqK6uLriOWSUbP31BpvozJ55apEjMdlfSs7kk7U0ukdwdEb9Mxa9J6pfe7wfUnydZBxyRV70KWJfKqwqUm5lZJ1HKs7kE3AEsj4ipeW/NASak5QnAA3nl4yV1l3QUuYH2F1KX2BZJw9M2L8qrY2ZmnUApu7lOA74ALJa0KJX9LTAFmCXpEuCPwDiAiFgqaRawjNyZYJMior5T+zKgBtiX3MC7B9+tYmW54JE1xYvDrJhKeTbXryk83gEwsok61wGN5ueOiFpgSPGiMzOzYvIV8GZmlpmTiZmZZeaJHs3awdeomO3OLRMzM8vMycTMzDJzN5dZBZm5T6OTHdvooaLEYdaQWyZmZpaZk4mZmWXmZGJmZpl5zMSszDLdh+UTxYvDrJjcMjEzs8ycTMzMLDN3c5ntQbLcXMs31rLmuGViZmaZuWViVkEy3QsFfD8UKxm3TMzMLDMnEzMzy8zJxMzMMitZMpF0p6QNkpbklfWW9ISklem5V957V0laJekVSaPyyk+StDi9d7Okpm4FbGZmHaSUA/A1wI+Bu/LKrgTmRcQUSVem19+TNBgYDxwHHA48KenDEbETmAZMBJ4DHgFGA4+WMG6zLivbrMOecdiaVrKWSUQ8A/ypQfEYYEZangGMzSufGRHbI+JVYBVwiqR+QM+IWBARQS4xjcXMzDqVco+ZHBYR6wHSc/15jv2BtXnr1aWy/mm5YbmZmXUineU6k0LjINFMeeGNSBPJdYkxYMCA4kRm1olkmiQSPFGklUy5Wyavpa4r0vOGVF4HHJG3XhWwLpVXFSgvKCKmR0R1RFT37du3qIGbmVnTyt0ymQNMAKak5wfyyu+RNJXcAPwg4IWI2Clpi6ThwPPARcAtZY7ZrMvIcgX9nUWMw7qekiUTSfcCI4BDJNUBV5NLIrMkXQL8ERgHEBFLJc0ClgE7gEnpTC6Ay8idGbYvubO4fCaXmVknU7JkEhHnN/HWyCbWvw5odN5iRNQCQ4oYmpmZFZmvgDczs8ycTMzMLDMnEzMzy6yzXGdiZp2c79JozXEyMduDZLno8XunFzEQ63LczWVmZpk5mZiZWWZOJmZmlpmTiZmZZeZkYmZmmflsLjNrFd+l0ZrjZGJmreIZh6057uYyM7PMnEzMzCwzJxMzM8vMYyZm1ipZpmIZv8bzenV1bpmYmVlmTiZmZpaZk4mZmWVWMWMmkkYDPwK6AbdHxJQODsnMWmm/I2/LUNtjJpWgIpKJpG7AT4CzgTpgoaQ5EbGsYyMzs9bw4H3XVxHJBDgFWBURfwCQNBMYA5Qkmbw07Jul2KwV2btb/qPD9r3PgeM6ZL8nLPpRh+y3I7lVUxkUER0dQ4sknQuMjogvp9dfAD4WEZc3WG8iMDG9PAZ4pcGmDgLebEXZIcCmIoTeXoViKue2WlunpfXa+76PU3HrtGa95tbxccq2rUo7TpsiYnQLcTQWEZ3+AYwjN05S//oLwC3t2M70VpbVdvDnbRRTObfV2jotrdfe932cynucWlrHx8nHqTWPSjmbqw44Iu91FbCuHdt5sJVlHa2YMbVnW62t09J67X3fx6m4dVqzXnPr+Dhl29YecZwqpZtrL2AFMBL4f8BC4IKIWFqi/dVGRHUptm3F4+NUGXycKkPW41QRA/ARsUPS5cBj5E4NvrNUiSSZXsJtW/H4OFUGH6fKkOk4VUTLxMzMOrdKGTMxM7NOzMnEzMwyczIxM7PMnExaQdL+kmZIuk3ShR0djxUm6WhJd0ia3dGxWNMkjU2/Sw9I+suOjscKk3SspFslzZZ0WUvr77HJRNKdkjZIWtKgfLSkVyStknRlKv4sMDsiLgXOKXuwe7C2HKeI+ENEXNIxke7Z2nic7k+/S18EzuuAcPdYbTxOyyPiq8DngBZPGd5jkwlQA+w2ZUDehJKfBAYD50saTO4iybVptfbPWGftUUPrj5N1nBrafpz+Pr1v5VNDG46TpHOAXwPzWtrwHptMIuIZ4E8NindNKBkR7wL1E0rWkUsosAd/Zx2hjcfJOkhbjpNybgAejYgXyx3rnqytv08RMSciPg602L3vP4y7688HLRDIJZH+wC+Bv5Y0jc45XcSepuBxktRH0q3ACZKu6pjQLE9Tv09fB84CzpX01Y4IzHbT1O/TCEk3S/oZ8EhLG6mIK+DLSAXKIiLeAr5U7mCsSU0dp82A/zh1Hk0dp5uBm8sdjDWpqeM0H5jf2o24ZbK7Yk0oaaXl41QZfJwqQ1GOk5PJ7hYCgyQdJWkfYDwwp4NjssZ8nCqDj1NlKMpx2mOTiaR7gQXAMZLqJF0SETuA+gkllwOzSjyhpLXAx6ky+DhVhlIeJ0/0aGZmme2xLRMzMyseJxMzM8vMycTMzDJzMjEzs8ycTMzMLDMnEzMzy8zJpEJJGthwGuky7HO+pGanoi5VXGmeoI/nva6RdG4r6u0r6WlJ3VJsIenree//WNIXixRji99PkfbzDUnLJd2dcTutjjd9dxfkva6WdHNa3u3YtGcbnYGkuZLekPRQg/KZkgZ1VFyVwsmki0vTS3cFI4AW/2AVcDHwy4iov3XABuCb6UrfTkNSW+bJ+xrwqYgo543aBgK7EkFE1EbEN9LLEbTu2DS3jZKStLoVq90IfKFA+TTgiqIG1AU5mVS2vZS7A+TL6W5o+0HuF0fSP0j6NTBO0qWSFkr6b0m/yFuvJs0K+htJf8j/T1/SFZIWpzpT8vY5TtILklZIOqO54FJr4Ma075clfSWVj0j/Fc+W9DtJd0tSeu9TqezXKbaHJA0kN4HjtyUtytvvXxSKvYELgQfyXm8kd2+GCQXi3fWfuqRD6v8ASfqipPslPSjpVUmXS5os6SVJz0nqnbeZz6eYlkg6JdXfX7mbEi1Mdcbkbfc/JD0IPF4gnslpO0skfSuV3QocDcyR9O0G6x+Xjs2i9H0PathSlPQdSde0EO//TttYlOI9EJgCnJHKvp2OYcFjowatRklb02LBbaR1eqfv+OX0nQ5N5dek725+Os7fyPtOH04/n0skZb7JVkTMA7YUeOu/gLPUtoS/54kIPyrwQe6/vABOS6/vBL6TllcDV+St2ydv+Vrg62m5BvgPcv9UDCZ3TwPI3STnN8B+6XXv9Dwf+Je0/CngySbiWpKWJwJ/n5a7A7XAUeT+k32T3IRyHyI3vcPpQA9yU2EflercCzyUlq+p/3zNxd4gln2A/2kYW4rhd0A34MfAF/M+X3VaPgRYnZa/CKwCDgT6pti/mt77V+BbefVvS8t/kfc9/BPw+bR8MLAC2D9tt67++20Q+0nA4rTeAcBS4IS843tIgTq3ABfmffZ9849HKv8OcE0L8T7IBz9XB5CbXXxE/bFI5bteN3Fszs17vbVhnQLbuAW4Oi1/AliUt+3fkPv5OQTYDOwN/HV97Gm9g1r4fVndyt+r3WLMK38COKmjf+8788Mtk8q2NiKeTcs/J/cHud59ectDJP2XpMXk/lM/Lu+9+yPi/YhYBhyWys4C/j0itgFERP7NdH6Znn9L7g9Vc/4SuEjSIuB5oA9Q3/f8QkTURcT7wKK0rY8Af4iIV9M697aw/UKx5zsEeKNhYdr+C+R1ubTCUxGxJSI2kksm9fe1Wczu38O9aR/PAD0lHUzue7gyfQ/zySXNAWn9Jxp8v/VOB/4zIt6KiK3kvvdmW4LkkvLfSvoecGREvN2Kz1Uo3meBqakVcHDk5m4qtdOB/5ti+RXQR9JB6b2HI2J7RGwi1015GLnv/SxJN0g6IyLebLhBSX9X38ICDs9rbbXn7o4bgMPbUW+P4WRS2RpOrJb/+q285Rrg8og4HvgBuT9m9bbnLSvvualJ2+rX30nL98MRuVbQsPQ4KiLqu3Py91u/rUL3VWhOodjzvc3unzXfPwHfY/ffgR15rxvWy9/X+3mv32f376HQMRHw13nfw4CIWJ7ef4vC2vpdEBH3AOeQ+9yPSfoEu38maPy5GsUbEVOAL5Nr2Twn6SNtDGXXPlP3ZWvGpwreUyM9N/pZiYgVfNB6u17SPzSqHHFd/XcOrMv7/ie1/qPs0oPc92pNcDKpbAMknZqWzyd3r+ZCDgTWS9qbVtx+k1z//cX6YGyldwvrN+Ux4LK0XyR9WNL+zaz/O+Do1A8PkN8PvoXc52i1iHgd6CapUUKJiN8By4BP5xWvJvcHCqDFM8WacB6ApNOBN9N/zI8BX88bFzqhFdt5Bhgrab/0nX2GXN99kyQdTa5ldzO5KcSHAq8Bhyp3F8ru7P55C8Yr6X9FxOKIuIFc1+RHaP77b/jeaj74HseQ65YqtF7Dz3thimUEsCki/tzMZz0c2BYRPwduAk5sat0i+TC5rkZrggeUKttyYIJyt9VcSe6sk0K+T66baQ25/+Sa/aMcEXMlDQNqJb1L7padf9uO+G4n1wX0YvpDuhEY28x+35b0NWCupE3kuqLqPQjMToPXXy+4gcIeJ9eF8mSB964DXsp7fRMwS9IXgF+1YR/5Xpf0G6AnuTPJAP4R+Dfg5fQ9rKbxH/XdRMSLkmr44Du4PSJeaqYK5BLD5yW9B/wP8MOIeE/SD8kd/1fJJeyW4v2WpDPJtQKWAY+Sa4HtkPTf5Fq6+bE0PDa3AQ9IeoHcyQ71ra+Xm9nGNcC/S3oZ2EaBEyQaOB64UdL7wHvAZS2s3yJJ/0UucR4gqQ64JCIek3QY8HZErM+6j67MU9BbpyLpgIjYmv7o/gRYGRH/mmF7JwCTI6LQKZ9mLUpnzf05Iu7o6Fg6M3dzWWdzaRowXQocBPwsy8bSf/NPqetcb2Pl9wYwo6OD6OzcMjEzs8zcMjEzs8ycTMzMLDMnEzMzy8zJxMzMMnMyMTOzzJxMzMwss/8Pw9i3MVld+gMAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "def extract_branch_lengths(t1):\n",
    "    \"\"\"\n",
    "    Extracts branch lengths for each branch in the tree\n",
    "    \n",
    "    Parameters\n",
    "    ----------\n",
    "    t1: ete3.Tree\n",
    "        the phylogeny of strains\n",
    "\n",
    "    Returns\n",
    "    -------\n",
    "    list of tuples of (branch_name, float)\n",
    "\n",
    "    \"\"\"\n",
    "    branch_names = []\n",
    "    nodes = []\n",
    "\n",
    "    for node in t1.traverse(\"levelorder\"):\n",
    "        if len(node.children) > 1:\n",
    "\n",
    "            for idx, child_node in enumerate(node.children):\n",
    "                name_anc = node.name.split(\"_\")[0]\n",
    "                name_des = child_node.name.split(\"_\")[0]\n",
    "\n",
    "                branch_name = f\"{name_anc}-{name_des}\"\n",
    "                \n",
    "                if name_des=='canettii':\n",
    "                    #print(\"excluding canetti\")\n",
    "                    continue\n",
    "                    \n",
    "                branch_names.append(branch_name)\n",
    "                nodes.append(node.dist)\n",
    "\n",
    "    return branch_names, nodes\n",
    "\n",
    "# Note: each tree is constructed from a different # of positions(the mutation concatenate for that lineage)\n",
    "# So, to find the number of mutations per branch we must multliply by the # of mutated positions used to construct the tree\n",
    "\n",
    "lineages = [\"L1\", \"L2\", \"L3\", \"L4A\", \"L4B\", \"L4C\", \"L5\", \"L6\"]\n",
    "all_branch_lens = []\n",
    "l_label = []\n",
    "# extracts the branch length for each lineage and multiply by total # of muts. considered\n",
    "for lineage in lineages: \n",
    "    \n",
    "    print(lineage)\n",
    "    # t1 = ete3.Tree(f\"../newinput/{lineage}node_labelled_newick_cleaned.tre\", format=1)\n",
    "    t1 = ete3.Tree(f\"input/{lineage}node_labelled_newick_cleaned.tre\", format=1)\n",
    "    b, n = extract_branch_lengths(t1)\n",
    "    \n",
    "    mutation_table = pd.read_csv(f\"input/{lineage}all_mutation_events.tsv\", sep=\"\\t\", low_memory=False)\n",
    "    \n",
    "    num_mutations = len(mutation_table.Position.unique())\n",
    "    #print(lineage, num_mutations)\n",
    "    n = np.array(n)\n",
    "    n = n * num_mutations\n",
    "    \n",
    "    # This normalizes to mutation rate per genome\n",
    "    #n = n / 4411532\n",
    "\n",
    "    all_branch_lens += list(n)\n",
    "    l_label += [lineage]*len(n)\n",
    "\n",
    "data = pd.DataFrame({\"length\":all_branch_lens, \"lineage\":l_label})\n",
    "# Add 1 to all branch lengths because some are length 0\n",
    "data.length = data.length+1\n",
    "\n",
    "## Plotting data with colors by lineage\n",
    "\n",
    "data[\"lineage\"] = [\"L4\" if \"4\" in x else x for x in data.lineage]\n",
    "sns.histplot(data, x=\"length\", hue=\"lineage\", stat=\"count\", \n",
    "             multiple=\"stack\",log_scale=True, linewidth=0, bins=20)\n",
    "plt.xlabel(\"branch length (Number of substitutions + 1)\")\n",
    "plt.ylabel(\"N\")\n",
    "sns.despine()\n",
    "plt.savefig(\"figures/substitutions_per_genome_by_lineage.pdf\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1: Create Mutation Matrices for each lineage independently based on the respective tree\n",
    "\n",
    "- Uses data previously computed by Roger Vargas using SNPPar\n",
    "- Each lineage phylogeny was constructed independently. Lineage 4 was partitioned into three groups due to its large size"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "mkdir: cannot create directory ‘output’: File exists\n",
      "mkdir: cannot create directory ‘output/00.mutation_data_by_lineage’: File exists\n",
      "mkdir: cannot create directory ‘output/00.mutation_data_by_lineage/Lineage1’: File exists\n",
      "mkdir: cannot create directory ‘output/00.mutation_data_by_lineage/Lineage2’: File exists\n",
      "mkdir: cannot create directory ‘output/00.mutation_data_by_lineage/Lineage3’: File exists\n",
      "mkdir: cannot create directory ‘output/00.mutation_data_by_lineage/Lineage4A’: File exists\n",
      "mkdir: cannot create directory ‘output/00.mutation_data_by_lineage/Lineage4B’: File exists\n",
      "mkdir: cannot create directory ‘output/00.mutation_data_by_lineage/Lineage4C’: File exists\n",
      "mkdir: cannot create directory ‘output/00.mutation_data_by_lineage/Lineage5’: File exists\n",
      "mkdir: cannot create directory ‘output/00.mutation_data_by_lineage/Lineage6’: File exists\n"
     ]
    }
   ],
   "source": [
    "!mkdir output\n",
    "!mkdir output/00.mutation_data_by_lineage\n",
    "!mkdir output/00.mutation_data_by_lineage/Lineage1\n",
    "!mkdir output/00.mutation_data_by_lineage/Lineage2\n",
    "!mkdir output/00.mutation_data_by_lineage/Lineage3\n",
    "!mkdir output/00.mutation_data_by_lineage/Lineage4A\n",
    "!mkdir output/00.mutation_data_by_lineage/Lineage4B\n",
    "!mkdir output/00.mutation_data_by_lineage/Lineage4C\n",
    "!mkdir output/00.mutation_data_by_lineage/Lineage5\n",
    "!mkdir output/00.mutation_data_by_lineage/Lineage6"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Number of branches: 5630\n",
      "Number of tips: 2816\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/tmp/ipykernel_30924/2934379494.py:15: DtypeWarning: Columns (7,8,9,16,17,19,20) have mixed types. Specify dtype option on import or set low_memory=False.\n",
      "  mutation_table = pd.read_csv(f\"input/L1all_mutation_events.tsv\", sep=\"\\t\")\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Number of unique mutations: 189040\n",
      "Number of mutations: 204695\n"
     ]
    }
   ],
   "source": [
    "prefix = \"output/00.mutation_data_by_lineage/Lineage1/L1\"\n",
    "\n",
    "ancestral_file=\"input/MTB_ancestor_reference_modified.fasta\"\n",
    "\n",
    "# Read in the phylogenetic tree\n",
    "t1 = ete3.Tree(f\"input/L1node_labelled_newick_cleaned.tre\", format=1)\n",
    "\n",
    "print(\"Number of branches:\", len(t1.get_edges()))\n",
    "print(\"Number of tips:\", len(t1.get_leaves()))\n",
    "\n",
    "# Read in concatenated SNP sequences used to do ASR on phylogeny\n",
    "aln = Alignment.from_file(open(f\"input/L1node_sequences.fasta\",'r'))\n",
    "\n",
    "# Read in table of all mutation events from SNPPar\n",
    "mutation_table = pd.read_csv(f\"input/L1all_mutation_events.tsv\", sep=\"\\t\")\n",
    "print(\"Number of unique mutations:\", len(mutation_table.Position.unique()))\n",
    "print(\"Number of mutations:\", len(mutation_table))\n",
    "\n",
    "# Count homoplasies on phylogeny relative to ancestral sequence\n",
    "lineage = 1\n",
    "ha = HomoplasyAnalysis(mutation_table, t1, lineage, ancestral_file, aln)\n",
    "ha.mutation_event_matrix = ha.create_mutation_matrix(mode=\"mutation\")\n",
    "ha.reversion_matrix = ha.create_mutation_matrix(mode=\"reversion\")\n",
    "ha.background_matrix = ha.create_background_matrix()\n",
    "\n",
    "# Save data\n",
    "ha.dump_sparse(prefix)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/tmp/ipykernel_30924/1369222226.py:1: DtypeWarning: Columns (7,8,9,16,17,19,20) have mixed types. Specify dtype option on import or set low_memory=False.\n",
      "  old_mutation_table = pd.read_csv(\"../../newRepo/DependentMutations/input/L1all_mutation_events_checked_against_seqs.tsv\", sep=\"\\t\")\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "204248"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "old_mutation_table = pd.read_csv(\"../../newRepo/DependentMutations/input/L1all_mutation_events_checked_against_seqs.tsv\", sep=\"\\t\")\n",
    "len(old_mutation_table)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Number of branches: 16180\n",
      "Number of tips: 8091\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/tmp/ipykernel_30924/509732897.py:14: DtypeWarning: Columns (7,8,9) have mixed types. Specify dtype option on import or set low_memory=False.\n",
      "  mutation_table = pd.read_csv(f\"input/L2all_mutation_events.tsv\", sep=\"\\t\")\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Number of unique mutations: 155809\n"
     ]
    }
   ],
   "source": [
    "prefix = \"output/00.mutation_data_by_lineage/Lineage2/L2\"\n",
    "\n",
    "ancestral_file=\"input/MTB_ancestor_reference_modified.fasta\"\n",
    "\n",
    "# Read in the phylogenetic tree\n",
    "t1 = ete3.Tree(f\"input/L2node_labelled_newick_cleaned.tre\", format=1)\n",
    "print(\"Number of branches:\", len(t1.get_edges()))\n",
    "print(\"Number of tips:\", len(t1.get_leaves()))\n",
    "\n",
    "# Read in concatenated SNP sequences used to do ASR on phylogeny\n",
    "aln = Alignment.from_file(open(f\"input/L2node_sequences.fasta\",'r'))\n",
    "\n",
    "# Read in table of all mutation events from SNPPar\n",
    "mutation_table = pd.read_csv(f\"input/L2all_mutation_events.tsv\", sep=\"\\t\")\n",
    "print(\"Number of unique mutations:\", len(mutation_table.Position.unique()))\n",
    "\n",
    "# Count homoplasies on phylogeny relative to ancestral sequence\n",
    "lineage = 2\n",
    "ha = HomoplasyAnalysis(mutation_table, t1, lineage, ancestral_file, aln)\n",
    "ha.mutation_event_matrix = ha.create_mutation_matrix(mode=\"mutation\")\n",
    "ha.reversion_matrix = ha.create_mutation_matrix(mode=\"reversion\")\n",
    "ha.background_matrix = ha.create_background_matrix()\n",
    "\n",
    "# Save data\n",
    "ha.dump_sparse(prefix)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Number of branches: 6796\n",
      "Number of tips: 3399\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/tmp/ipykernel_30924/3882852268.py:14: DtypeWarning: Columns (7,8,9,16,17,19,20) have mixed types. Specify dtype option on import or set low_memory=False.\n",
      "  mutation_table = pd.read_csv(f\"input/L3all_mutation_events.tsv\", sep=\"\\t\")\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Number of unique mutations: 128623\n"
     ]
    }
   ],
   "source": [
    "prefix = \"output/00.mutation_data_by_lineage/Lineage3/L3\"\n",
    "\n",
    "ancestral_file=\"input/MTB_ancestor_reference_modified.fasta\"\n",
    "\n",
    "# Read in the phylogenetic tree\n",
    "t1 = ete3.Tree(f\"input/L3node_labelled_newick_cleaned.tre\", format=1)\n",
    "print(\"Number of branches:\", len(t1.get_edges()))\n",
    "print(\"Number of tips:\", len(t1.get_leaves()))\n",
    "\n",
    "# Read in concatenated SNP sequences used to do ASR on phylogeny\n",
    "aln = Alignment.from_file(open(f\"input/L3node_sequences.fasta\",'r'))\n",
    "\n",
    "# Read in table of all mutation events from SNPPar\n",
    "mutation_table = pd.read_csv(f\"input/L3all_mutation_events.tsv\", sep=\"\\t\")\n",
    "print(\"Number of unique mutations:\", len(mutation_table.Position.unique()))\n",
    "\n",
    "# Count homoplasies on phylogeny relative to ancestral sequence\n",
    "lineage = 3\n",
    "ha = HomoplasyAnalysis(mutation_table, t1, lineage, ancestral_file, aln)\n",
    "ha.mutation_event_matrix = ha.create_mutation_matrix(mode=\"mutation\")\n",
    "ha.reversion_matrix = ha.create_mutation_matrix(mode=\"reversion\")\n",
    "ha.background_matrix = ha.create_background_matrix()\n",
    "\n",
    "# Save data\n",
    "ha.dump_sparse(prefix)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Number of branches: 11678\n",
      "Number of tips: 5840\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/tmp/ipykernel_30924/3416117973.py:14: DtypeWarning: Columns (7,8,9) have mixed types. Specify dtype option on import or set low_memory=False.\n",
      "  mutation_table = pd.read_csv(f\"input/L4Aall_mutation_events.tsv\", sep=\"\\t\")\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Number of unique mutations: 140430\n"
     ]
    }
   ],
   "source": [
    "prefix = \"output/00.mutation_data_by_lineage/Lineage4A/L4A\"\n",
    "\n",
    "ancestral_file=\"input/MTB_ancestor_reference_modified.fasta\"\n",
    "\n",
    "# Read in the phylogenetic tree\n",
    "t1 = ete3.Tree(f\"input/L4Anode_labelled_newick_cleaned.tre\", format=1)\n",
    "print(\"Number of branches:\", len(t1.get_edges()))\n",
    "print(\"Number of tips:\", len(t1.get_leaves()))\n",
    "\n",
    "# Read in concatenated SNP sequences used to do ASR on phylogeny\n",
    "aln = Alignment.from_file(open(f\"input/L4Anode_sequences.fasta\",'r'))\n",
    "\n",
    "# Read in table of all mutation events from SNPPar\n",
    "mutation_table = pd.read_csv(f\"input/L4Aall_mutation_events.tsv\", sep=\"\\t\")\n",
    "print(\"Number of unique mutations:\", len(mutation_table.Position.unique()))\n",
    "\n",
    "# Count homoplasies on phylogeny relative to ancestral sequence\n",
    "lineage = '4A'\n",
    "ha = HomoplasyAnalysis(mutation_table, t1, lineage, ancestral_file, aln)\n",
    "ha.mutation_event_matrix = ha.create_mutation_matrix(mode=\"mutation\")\n",
    "ha.reversion_matrix = ha.create_mutation_matrix(mode=\"reversion\")\n",
    "ha.background_matrix = ha.create_background_matrix()\n",
    "\n",
    "# Save data\n",
    "ha.dump_sparse(prefix)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Number of branches: 13914\n",
      "Number of tips: 6958\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/tmp/ipykernel_30924/1034920307.py:14: DtypeWarning: Columns (7,8,9,16,17,19,20) have mixed types. Specify dtype option on import or set low_memory=False.\n",
      "  mutation_table = pd.read_csv(f\"input/L4Ball_mutation_events.tsv\", sep=\"\\t\")\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Number of unique mutations: 126238\n"
     ]
    }
   ],
   "source": [
    "prefix = \"output/00.mutation_data_by_lineage/Lineage4B/L4B\"\n",
    "\n",
    "ancestral_file=\"input/MTB_ancestor_reference_modified.fasta\"\n",
    "\n",
    "# Read in the phylogenetic tree\n",
    "t1 = ete3.Tree(f\"input/L4Bnode_labelled_newick_cleaned.tre\", format=1)\n",
    "print(\"Number of branches:\", len(t1.get_edges()))\n",
    "print(\"Number of tips:\", len(t1.get_leaves()))\n",
    "\n",
    "# Read in concatenated SNP sequences used to do ASR on phylogeny\n",
    "aln = Alignment.from_file(open(f\"input/L4Bnode_sequences.fasta\",'r'))\n",
    "\n",
    "# Read in table of all mutation events from SNPPar\n",
    "mutation_table = pd.read_csv(f\"input/L4Ball_mutation_events.tsv\", sep=\"\\t\")\n",
    "print(\"Number of unique mutations:\", len(mutation_table.Position.unique()))\n",
    "\n",
    "# Count homoplasies on phylogeny relative to ancestral sequence\n",
    "lineage = '4B'\n",
    "ha = HomoplasyAnalysis(mutation_table, t1, lineage, ancestral_file, aln)\n",
    "ha.mutation_event_matrix = ha.create_mutation_matrix(mode=\"mutation\")\n",
    "ha.reversion_matrix = ha.create_mutation_matrix(mode=\"reversion\")\n",
    "ha.background_matrix = ha.create_background_matrix()\n",
    "\n",
    "# Save data\n",
    "ha.dump_sparse(prefix)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Number of branches: 8268\n",
      "Number of tips: 4135\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/tmp/ipykernel_30924/2790187550.py:14: DtypeWarning: Columns (7,8,9) have mixed types. Specify dtype option on import or set low_memory=False.\n",
      "  mutation_table = pd.read_csv(f\"input/L4Call_mutation_events.tsv\", sep=\"\\t\")\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Number of unique mutations: 142165\n"
     ]
    }
   ],
   "source": [
    "prefix = \"output/00.mutation_data_by_lineage/Lineage4C/L4C\"\n",
    "\n",
    "ancestral_file=\"input/MTB_ancestor_reference_modified.fasta\"\n",
    "\n",
    "# Read in the phylogenetic tree\n",
    "t1 = ete3.Tree(f\"input/L4Cnode_labelled_newick_cleaned.tre\", format=1)\n",
    "print(\"Number of branches:\", len(t1.get_edges()))\n",
    "print(\"Number of tips:\", len(t1.get_leaves()))\n",
    "\n",
    "# Read in concatenated SNP sequences used to do ASR on phylogeny\n",
    "aln = Alignment.from_file(open(f\"input/L4Cnode_sequences.fasta\",'r'))\n",
    "\n",
    "# Read in table of all mutation events from SNPPar\n",
    "mutation_table = pd.read_csv(f\"input/L4Call_mutation_events.tsv\", sep=\"\\t\")\n",
    "print(\"Number of unique mutations:\", len(mutation_table.Position.unique()))\n",
    "\n",
    "# Count homoplasies on phylogeny relative to ancestral sequence\n",
    "lineage = '4C'\n",
    "ha = HomoplasyAnalysis(mutation_table, t1, lineage, ancestral_file, aln)\n",
    "ha.mutation_event_matrix = ha.create_mutation_matrix(mode=\"mutation\")\n",
    "ha.reversion_matrix = ha.create_mutation_matrix(mode=\"reversion\")\n",
    "ha.background_matrix = ha.create_background_matrix()\n",
    "\n",
    "# Save data\n",
    "ha.dump_sparse(prefix)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Number of branches: 196\n",
      "Number of tips: 99\n",
      "Number of unique mutations: 8692\n"
     ]
    }
   ],
   "source": [
    "prefix = \"output/00.mutation_data_by_lineage/Lineage5/L5\"\n",
    "\n",
    "ancestral_file=\"input/MTB_ancestor_reference_modified.fasta\"\n",
    "\n",
    "# Read in the phylogenetic tree\n",
    "t1 = ete3.Tree(f\"input/L5node_labelled_newick_cleaned.tre\", format=1)\n",
    "print(\"Number of branches:\", len(t1.get_edges()))\n",
    "print(\"Number of tips:\", len(t1.get_leaves()))\n",
    "\n",
    "# Read in concatenated SNP sequences used to do ASR on phylogeny\n",
    "aln = Alignment.from_file(open(f\"input/L5node_sequences.fasta\",'r'))\n",
    "\n",
    "# Read in table of all mutation events from SNPPar\n",
    "mutation_table = pd.read_csv(f\"input/L5all_mutation_events.tsv\", sep=\"\\t\")\n",
    "print(\"Number of unique mutations:\", len(mutation_table.Position.unique()))\n",
    "\n",
    "# Count homoplasies on phylogeny relative to ancestral sequence\n",
    "lineage = '5'\n",
    "ha = HomoplasyAnalysis(mutation_table, t1, lineage, ancestral_file, aln)\n",
    "ha.mutation_event_matrix = ha.create_mutation_matrix(mode=\"mutation\")\n",
    "ha.reversion_matrix = ha.create_mutation_matrix(mode=\"reversion\")\n",
    "ha.background_matrix = ha.create_background_matrix()\n",
    "\n",
    "# Save data\n",
    "ha.dump_sparse(prefix)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Number of branches: 192\n",
      "Number of tips: 97\n",
      "Number of unique mutations: 11271\n"
     ]
    }
   ],
   "source": [
    "prefix = \"output/00.mutation_data_by_lineage/Lineage6/L6\"\n",
    "\n",
    "ancestral_file=\"input/MTB_ancestor_reference_modified.fasta\"\n",
    "\n",
    "# Read in the phylogenetic tree\n",
    "t1 = ete3.Tree(f\"input/L6node_labelled_newick_cleaned.tre\", format=1)\n",
    "print(\"Number of branches:\", len(t1.get_edges()))\n",
    "print(\"Number of tips:\", len(t1.get_leaves()))\n",
    "\n",
    "# Read in concatenated SNP sequences used to do ASR on phylogeny\n",
    "aln = Alignment.from_file(open(f\"input/L6node_sequences.fasta\",'r'))\n",
    "\n",
    "# Read in table of all mutation events from SNPPar\n",
    "mutation_table = pd.read_csv(f\"input/L6all_mutation_events.tsv\", sep=\"\\t\")\n",
    "print(\"Number of unique mutations:\", len(mutation_table.Position.unique()))\n",
    "\n",
    "# Count homoplasies on phylogeny relative to ancestral sequence\n",
    "lineage = '6'\n",
    "ha = HomoplasyAnalysis(mutation_table, t1, lineage, ancestral_file, aln)\n",
    "ha.mutation_event_matrix = ha.create_mutation_matrix(mode=\"mutation\")\n",
    "ha.reversion_matrix = ha.create_mutation_matrix(mode=\"reversion\")\n",
    "ha.background_matrix = ha.create_background_matrix()\n",
    "\n",
    "# Save data\n",
    "ha.dump_sparse(prefix)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Determine mutations to examine based on number of evolutionary events"
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {},
   "source": [
    "Construct a table with the frequency of all evolutionary events, and filter for 5 or more events in a particular site.\n",
    "Note that this doesn't distinguish between mutations and reversions. Sites with fewer than 5 mutations, excluding reversions, are filtered\n",
    "at a later step. Initial filtering done here to reduce the computational needs of next step"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/tmp/ipykernel_30924/3305117556.py:1: DtypeWarning: Columns (7,8,9,16,17,19,20) have mixed types. Specify dtype option on import or set low_memory=False.\n",
      "  mutation_table = pd.read_csv(f\"input/L1all_mutation_events.tsv\", sep=\"\\t\")\n",
      "/tmp/ipykernel_30924/3305117556.py:4: DtypeWarning: Columns (7,8,9) have mixed types. Specify dtype option on import or set low_memory=False.\n",
      "  mutation_table = pd.read_csv(f\"input/L2all_mutation_events.tsv\", sep=\"\\t\")\n",
      "/tmp/ipykernel_30924/3305117556.py:7: DtypeWarning: Columns (7,8,9,16,17,19,20) have mixed types. Specify dtype option on import or set low_memory=False.\n",
      "  mutation_table = pd.read_csv(f\"input/L3all_mutation_events.tsv\", sep=\"\\t\")\n",
      "/tmp/ipykernel_30924/3305117556.py:10: DtypeWarning: Columns (7,8,9) have mixed types. Specify dtype option on import or set low_memory=False.\n",
      "  mutation_table = pd.read_csv(f\"input/L4Aall_mutation_events.tsv\", sep=\"\\t\")\n",
      "/tmp/ipykernel_30924/3305117556.py:13: DtypeWarning: Columns (7,8,9,16,17,19,20) have mixed types. Specify dtype option on import or set low_memory=False.\n",
      "  mutation_table = pd.read_csv(f\"input/L4Ball_mutation_events.tsv\", sep=\"\\t\")\n",
      "/tmp/ipykernel_30924/3305117556.py:16: DtypeWarning: Columns (7,8,9) have mixed types. Specify dtype option on import or set low_memory=False.\n",
      "  mutation_table = pd.read_csv(f\"input/L4Call_mutation_events.tsv\", sep=\"\\t\")\n"
     ]
    }
   ],
   "source": [
    "mutation_table = pd.read_csv(f\"input/L1all_mutation_events.tsv\", sep=\"\\t\")\n",
    "t1 = mutation_table.groupby(\"Position\").size()\n",
    "\n",
    "mutation_table = pd.read_csv(f\"input/L2all_mutation_events.tsv\", sep=\"\\t\")\n",
    "t2 = mutation_table.groupby(\"Position\").size()\n",
    "\n",
    "mutation_table = pd.read_csv(f\"input/L3all_mutation_events.tsv\", sep=\"\\t\")\n",
    "t3 = mutation_table.groupby(\"Position\").size()\n",
    "\n",
    "mutation_table = pd.read_csv(f\"input/L4Aall_mutation_events.tsv\", sep=\"\\t\")\n",
    "t4a = mutation_table.groupby(\"Position\").size()\n",
    "\n",
    "mutation_table = pd.read_csv(f\"input/L4Ball_mutation_events.tsv\", sep=\"\\t\")\n",
    "t4b = mutation_table.groupby(\"Position\").size()\n",
    "\n",
    "mutation_table = pd.read_csv(f\"input/L4Call_mutation_events.tsv\", sep=\"\\t\")\n",
    "t4c = mutation_table.groupby(\"Position\").size()\n",
    "\n",
    "mutation_table = pd.read_csv(f\"input/L5all_mutation_events.tsv\", sep=\"\\t\")\n",
    "t5 = mutation_table.groupby(\"Position\").size()\n",
    "\n",
    "mutation_table = pd.read_csv(f\"input/L6all_mutation_events.tsv\", sep=\"\\t\")\n",
    "t6 = mutation_table.groupby(\"Position\").size()\n",
    "\n",
    "\n",
    "t = pd.DataFrame({\n",
    "    \"L1\": t1,\n",
    "    \"L2\": t2,\n",
    "    \"L3\": t3,\n",
    "    \"L4A\": t4a,\n",
    "    \"L4B\": t4b,\n",
    "    \"L4C\": t4c,\n",
    "    \"L5\": t5,\n",
    "    \"L6\": t6,\n",
    "})\n",
    "t = t.fillna(0)\n",
    "t[\"summed\"] = t.sum(axis=1)\n",
    "\n",
    "t.to_csv(\"output/00.mutation_data_by_lineage/event_frequency_table.csv\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Combine mutation matrices for each lineage into final matrix"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "def create_matrix(prefixes, matrix_name, all_branches, positions):\n",
    "    \"\"\"\n",
    "    Creates an N_branches x N_positions matrix tracking evolution on all the trees\n",
    "    \n",
    "    Parameters\n",
    "    ----------\n",
    "    prefixes: list of str\n",
    "        list of prefixes containing the lineage-specific data\n",
    "    matrix_name: str\n",
    "        prefix indicating which matrix to load\n",
    "        current options are mutation_matrix, background_matrix, and reversion_matrix\n",
    "    all_branches: list of str\n",
    "    \n",
    "    positions: list of int\n",
    "    \"\"\"\n",
    "\n",
    "    full_positions_to_index = {x:idx for idx,x in enumerate(positions)}\n",
    "    print(\"Number of positions after 5 event filtered\", len(full_positions_to_index))\n",
    "    \n",
    "    # Initialize the N_branches x N_positions matrix\n",
    "    full_mutation_matrix = np.zeros(shape=(len(all_branches), len(filtered_positions)),dtype='int8')\n",
    "    print(\"initialized matrix with shape:\", full_mutation_matrix.shape)\n",
    "\n",
    "    # For each lineage sub-tree\n",
    "    for prefix in prefixes:\n",
    "\n",
    "        #Load in the branches from current lineage and get their indices in the full matrix\n",
    "        branches = np.load(f\"{prefix}_branches.npy\")\n",
    "        branches_start_index = all_branches.index(branches[0])\n",
    "        branches_end_index= all_branches.index(branches[-1])\n",
    "\n",
    "        #Load in the positions found to be mutated in current lineage\n",
    "        positions_in_lineage = list(np.load(f\"{prefix}_positions.npy\"))\n",
    "        positions_in_lineage_to_index = {x:idx for idx,x in enumerate(positions_in_lineage)}\n",
    "\n",
    "        #Load in the matrix and condense\n",
    "        matrix = sparse.load_npz(f\"{prefix}_{matrix_name}.npz\").toarray()\n",
    "        matrix = matrix.astype(\"int8\")\n",
    "        print(f\"{prefix} matrix has shape {matrix.shape}\")\n",
    "\n",
    "        # Find intersection between positions found in current lineage and full matrix (which has been pre-filtered for frequency)\n",
    "        test_set = set(positions)\n",
    "        filtered_positions_in_lineage = [x for x in positions_in_lineage if x in test_set]\n",
    "        print(f\"{len(filtered_positions_in_lineage)} positions in the {prefix} matrix are found in the full matrix\")\n",
    "\n",
    "        # Find the indices\n",
    "        mutation_matrix_index = np.array([positions_in_lineage_to_index[pos] for pos in filtered_positions_in_lineage])\n",
    "        full_mutation_matrix_index = np.array([full_positions_to_index[pos] for pos in filtered_positions_in_lineage])\n",
    "\n",
    "        # Assert that indices are both in ascending order -> otherwise, results will be wrong\n",
    "        assert(list(mutation_matrix_index) == sorted(list(mutation_matrix_index)))\n",
    "        assert(list(full_mutation_matrix_index) == sorted(list(full_mutation_matrix_index)))\n",
    "\n",
    "        # fill in the branch, positions in the full matrix corresponding to current lineage matrix\n",
    "        full_mutation_matrix[branches_start_index:branches_end_index + 1,:][:,full_mutation_matrix_index] = matrix[:,mutation_matrix_index]\n",
    "\n",
    "    return full_mutation_matrix\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "mkdir: cannot create directory ‘output/00.mutation_data_by_lineage/00.combined_data’: File exists\n"
     ]
    }
   ],
   "source": [
    "!mkdir output/00.mutation_data_by_lineage/00.combined_data\n",
    "\n",
    "prefixes =  [\n",
    "    \"output/00.mutation_data_by_lineage/Lineage1/L1\", \n",
    "    \"output/00.mutation_data_by_lineage/Lineage2/L2\", \n",
    "    \"output/00.mutation_data_by_lineage/Lineage3/L3\", \n",
    "    \"output/00.mutation_data_by_lineage/Lineage4A/L4A\",\n",
    "    \"output/00.mutation_data_by_lineage/Lineage4B/L4B\",\n",
    "    \"output/00.mutation_data_by_lineage/Lineage4C/L4C\", \n",
    "    \"output/00.mutation_data_by_lineage/Lineage5/L5\", \n",
    "    \"output/00.mutation_data_by_lineage/Lineage6/L6\"]\n",
    "\n",
    "# Get a list of all the branches across all the trees\n",
    "all_branches = []\n",
    "\n",
    "for prefix in prefixes:\n",
    "    branches = np.load(f\"{prefix}_branches.npy\")\n",
    "    \n",
    "    all_branches = all_branches + list(branches)\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Number of positions before EBR filtering 5589\n",
      "Number of positions after EBR filtering 5543\n",
      "Number of positions after blindspot filtering 5489\n"
     ]
    }
   ],
   "source": [
    "\n",
    "# Get a list of all positions to be used\n",
    "filtered_positions = pd.read_csv(\"output/00.mutation_data_by_lineage/event_frequency_table.csv\").query(\"summed > 4\")\n",
    "positions = list(filtered_positions.Position)\n",
    "\n",
    "# Ensure that we've accurately filtered positions from Marin et al\n",
    "ebr = pd.read_csv(\n",
    "    \"input/EBR.csv\", index_col=0\n",
    ")\n",
    "filtered_positions = filtered_positions.merge(ebr, how=\"left\", left_on=\"Position\", right_on=\"position\")\n",
    "print(\"Number of positions before EBR filtering\", len(filtered_positions))\n",
    "filtered_positions = filtered_positions.query(\"EBR >= 0.9\")\n",
    "print(\"Number of positions after EBR filtering\", len(filtered_positions))\n",
    "\n",
    "# Accurately filter the positions from Modlin et al\n",
    "modlin = pd.read_csv(\"input/S7-seq-attributes.csv\")\n",
    "blindspots = modlin.query(\"blindspot==1\")\n",
    "blindspot_list = list(blindspots.rv_position)\n",
    "filtered_positions = filtered_positions.query(\"not Position in @blindspot_list\")\n",
    "print(\"Number of positions after blindspot filtering\", len(filtered_positions))\n",
    "\n",
    "positions = list(filtered_positions.Position)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Number of positions after 5 event filtered 5489\n",
      "initialized matrix with shape: (62846, 5489)\n",
      "output/00.mutation_data_by_lineage/Lineage1/L1 matrix has shape (5629, 189040)\n",
      "2937 positions in the output/00.mutation_data_by_lineage/Lineage1/L1 matrix are found in the full matrix\n",
      "output/00.mutation_data_by_lineage/Lineage2/L2 matrix has shape (16179, 155809)\n",
      "4069 positions in the output/00.mutation_data_by_lineage/Lineage2/L2 matrix are found in the full matrix\n",
      "output/00.mutation_data_by_lineage/Lineage3/L3 matrix has shape (6795, 128623)\n",
      "2779 positions in the output/00.mutation_data_by_lineage/Lineage3/L3 matrix are found in the full matrix\n",
      "output/00.mutation_data_by_lineage/Lineage4A/L4A matrix has shape (11677, 140430)\n",
      "3306 positions in the output/00.mutation_data_by_lineage/Lineage4A/L4A matrix are found in the full matrix\n",
      "output/00.mutation_data_by_lineage/Lineage4B/L4B matrix has shape (13913, 126238)\n",
      "3109 positions in the output/00.mutation_data_by_lineage/Lineage4B/L4B matrix are found in the full matrix\n",
      "output/00.mutation_data_by_lineage/Lineage4C/L4C matrix has shape (8267, 142165)\n",
      "3341 positions in the output/00.mutation_data_by_lineage/Lineage4C/L4C matrix are found in the full matrix\n",
      "output/00.mutation_data_by_lineage/Lineage5/L5 matrix has shape (195, 8692)\n",
      "356 positions in the output/00.mutation_data_by_lineage/Lineage5/L5 matrix are found in the full matrix\n",
      "output/00.mutation_data_by_lineage/Lineage6/L6 matrix has shape (191, 11271)\n",
      "325 positions in the output/00.mutation_data_by_lineage/Lineage6/L6 matrix are found in the full matrix\n",
      "Number of positions after 5 event filtered 5489\n",
      "initialized matrix with shape: (62846, 5489)\n",
      "output/00.mutation_data_by_lineage/Lineage1/L1 matrix has shape (5629, 189040)\n",
      "2937 positions in the output/00.mutation_data_by_lineage/Lineage1/L1 matrix are found in the full matrix\n",
      "output/00.mutation_data_by_lineage/Lineage2/L2 matrix has shape (16179, 155809)\n",
      "4069 positions in the output/00.mutation_data_by_lineage/Lineage2/L2 matrix are found in the full matrix\n",
      "output/00.mutation_data_by_lineage/Lineage3/L3 matrix has shape (6795, 128623)\n",
      "2779 positions in the output/00.mutation_data_by_lineage/Lineage3/L3 matrix are found in the full matrix\n",
      "output/00.mutation_data_by_lineage/Lineage4A/L4A matrix has shape (11677, 140430)\n",
      "3306 positions in the output/00.mutation_data_by_lineage/Lineage4A/L4A matrix are found in the full matrix\n",
      "output/00.mutation_data_by_lineage/Lineage4B/L4B matrix has shape (13913, 126238)\n",
      "3109 positions in the output/00.mutation_data_by_lineage/Lineage4B/L4B matrix are found in the full matrix\n",
      "output/00.mutation_data_by_lineage/Lineage4C/L4C matrix has shape (8267, 142165)\n",
      "3341 positions in the output/00.mutation_data_by_lineage/Lineage4C/L4C matrix are found in the full matrix\n",
      "output/00.mutation_data_by_lineage/Lineage5/L5 matrix has shape (195, 8692)\n",
      "356 positions in the output/00.mutation_data_by_lineage/Lineage5/L5 matrix are found in the full matrix\n",
      "output/00.mutation_data_by_lineage/Lineage6/L6 matrix has shape (191, 11271)\n",
      "325 positions in the output/00.mutation_data_by_lineage/Lineage6/L6 matrix are found in the full matrix\n",
      "Number of positions after 5 event filtered 5489\n",
      "initialized matrix with shape: (62846, 5489)\n",
      "output/00.mutation_data_by_lineage/Lineage1/L1 matrix has shape (5629, 189040)\n",
      "2937 positions in the output/00.mutation_data_by_lineage/Lineage1/L1 matrix are found in the full matrix\n",
      "output/00.mutation_data_by_lineage/Lineage2/L2 matrix has shape (16179, 155809)\n",
      "4069 positions in the output/00.mutation_data_by_lineage/Lineage2/L2 matrix are found in the full matrix\n",
      "output/00.mutation_data_by_lineage/Lineage3/L3 matrix has shape (6795, 128623)\n",
      "2779 positions in the output/00.mutation_data_by_lineage/Lineage3/L3 matrix are found in the full matrix\n",
      "output/00.mutation_data_by_lineage/Lineage4A/L4A matrix has shape (11677, 140430)\n",
      "3306 positions in the output/00.mutation_data_by_lineage/Lineage4A/L4A matrix are found in the full matrix\n",
      "output/00.mutation_data_by_lineage/Lineage4B/L4B matrix has shape (13913, 126238)\n",
      "3109 positions in the output/00.mutation_data_by_lineage/Lineage4B/L4B matrix are found in the full matrix\n",
      "output/00.mutation_data_by_lineage/Lineage4C/L4C matrix has shape (8267, 142165)\n",
      "3341 positions in the output/00.mutation_data_by_lineage/Lineage4C/L4C matrix are found in the full matrix\n",
      "output/00.mutation_data_by_lineage/Lineage5/L5 matrix has shape (195, 8692)\n",
      "356 positions in the output/00.mutation_data_by_lineage/Lineage5/L5 matrix are found in the full matrix\n",
      "output/00.mutation_data_by_lineage/Lineage6/L6 matrix has shape (191, 11271)\n",
      "325 positions in the output/00.mutation_data_by_lineage/Lineage6/L6 matrix are found in the full matrix\n"
     ]
    }
   ],
   "source": [
    "\n",
    "np.save(\"output/00.mutation_data_by_lineage/00.combined_data/all_positions_5strain_filtering.npy\", np.array(positions))\n",
    "np.save(\"output/00.mutation_data_by_lineage/00.combined_data/all_branches.npy\", np.array(all_branches))\n",
    "\n",
    "mutation_matrix = create_matrix(prefixes, \"mutation_matrix\", all_branches, positions)\n",
    "np.save(\"output/00.mutation_data_by_lineage/00.combined_data/all_mutation_matrix_5strain_filtering.npy\", mutation_matrix)\n",
    "\n",
    "background_matrix =  create_matrix(prefixes, \"background_matrix\", all_branches, positions)\n",
    "np.save(\"output/00.mutation_data_by_lineage/00.combined_data/all_background_matrix_5strain_filtering.npy\", background_matrix)\n",
    "\n",
    "reversion_matrix =  create_matrix(prefixes, \"reversion_matrix\", all_branches, positions)\n",
    "np.save(\"output/00.mutation_data_by_lineage/00.combined_data/all_reversion_matrix_5strain_filtering.npy\", reversion_matrix)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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.9.13"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}