OrthoSLC_Python_jupyter_interface.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "e19ff2b5",
   "metadata": {},
   "source": [
    "# Jupyter interface of OrthoSLC (0.2.4)\n",
    "\n",
    "<b>OrthoSLC</b> is a pipline that perfomrs Reciprocal Best Blast Hit (RBBH) Single Linkage Clustering to obtain Orthologous Genes which assist core genome construction. <br>\n",
    "\n",
    "**This interface** is python jupyter based to facilitate customized analysis. However, for datasets with 500 or more genomes, we still recommend using the binary files, which is mainly written in C++, for optimal performance.<br>\n",
    "\n",
    "It is: <br>\n",
    "* <b>independent</b> of relational database management systems (e.g., MySQL)\n",
    "* recommend to handle less genomes.\n",
    "\n",
    "It is recommended for sub-species level single copy core genome construction since RBBH may not work well for missions like Human-Microbe core genome construction. \n",
    "\n",
    "<b>Caveat:</b><br>\n",
    "The pipeline is currently only tested on Ubuntu 20.04 and 18.04 and Debian strech.\n",
    "\n",
    "<b>Requirement:</b><br>\n",
    "Python3 (suggest newest stable release or higher),<br>\n",
    "NCBI Blast+ (suggest 2.12 or higher) <br>\n",
    "\n",
    "<b>Note:</b><br>\n",
    "For all steps, users do not need to make the ourput directory manually, program will do that for you.\n",
    "\n",
    "Bug report: \n",
    "* <jingjie.chencharly@gmail.com>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6ac28ea7",
   "metadata": {},
   "outputs": [],
   "source": [
    "import os,sys\n",
    "import random\n",
    "import copy\n",
    "# import numpy as np\n",
    "# import pandas as pd\n",
    "\n",
    "from multiprocessing import Lock, Process, Manager\n",
    "import subprocess\n",
    "import shutil\n",
    "\n",
    "# pip install biopython\n",
    "from Bio import SeqIO\n",
    "from Bio.SeqRecord import SeqRecord\n",
    "\n",
    "from collections import defaultdict\n",
    "\n",
    "import time\n",
    "from typing import List\n",
    "\n",
    "# def mission spliter\n",
    "def mission_spliter(lst: List, \n",
    "                    num_splits: int) -> List[List]:\n",
    "    n = len(lst)\n",
    "    random.shuffle(lst)\n",
    "    \n",
    "    split_size = n // num_splits\n",
    "    remainder = n % num_splits\n",
    "    start = 0\n",
    "    result = []\n",
    "    \n",
    "    for i in range(num_splits):\n",
    "        end = start + split_size\n",
    "        if i < remainder:\n",
    "            end += 1\n",
    "        if lst[start:end] == []:\n",
    "                break\n",
    "        result.append(lst[start:end])\n",
    "        start = end\n",
    "        \n",
    "    return result"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fc3f6a5d",
   "metadata": {},
   "source": [
    "## Step 1 Genome information preparation\n",
    "The pipeline starts with annotated genomes in fasta format. FASTA file name require strain name and extension (e.g., `strain_A.ffn`, `strain_B.fna`, `strain_C.fasta` etc.).<br> \n",
    "Step 1 needs the **path to directory of annotated FASTA files** as input, to genereate a header less, tab separated table, in which the \n",
    "* first column is a short ID, \n",
    "* second column is the strain name, \n",
    "* third column as the absolute path. \n",
    "\n",
    "The short ID of each genome is generated to save computational resources and storage space. Since reciprocal BLAST generates a large volume of files (millions to billions of rows if large number of genomes participated), each row contains the names of the query and subject. If the user provides input FASTA file names like:\n",
    "\n",
    "* `GCA_900627445.1_PFR31F05_genomic.fasta`\n",
    "* `GCA_021980615.1_PDT001237823.1_genomic.fasta`\n",
    "\n",
    "and such file names become part of gene identifier instead of the short ID used in this program, large size of additional storage will be consumed for intermediate files and even more pressure on computing memory for analysis of 1000~ genomes."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "08f870c1",
   "metadata": {},
   "source": [
    "### Set path\n",
    "* Set variable `raw_annotated_dir_path` to the path of directory of annoated FASTA files\n",
    "* Set variable `preparation_output_path` to the path of output tab separated files.<br>\n",
    "Then run the codes below [Step 1 execution](###step-1-execution)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a9b35ccd",
   "metadata": {},
   "outputs": [],
   "source": [
    "# set your own\n",
    "# directory path of input fastas\n",
    "raw_annotated_dir_path = 'test_inputs/'\n",
    "# path of output tsv\n",
    "preparation_output_path = 'test_op/Step1_op.txt'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a18153c9",
   "metadata": {},
   "outputs": [],
   "source": [
    "# os.mkdir(\"test_op\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7c777e17",
   "metadata": {},
   "source": [
    "### Step 1 execution"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c6d38771",
   "metadata": {},
   "outputs": [],
   "source": [
    "start_c = 0\n",
    "\n",
    "with open(preparation_output_path, 'w') as pre_op:\n",
    "    for file_naam in os.listdir(raw_annotated_dir_path):\n",
    "        abs_path = os.path.join(raw_annotated_dir_path,\n",
    "                                file_naam\n",
    "                               )\n",
    "        abs_path = os.path.abspath(abs_path)\n",
    "        \n",
    "        strain_naam = os.path.basename(abs_path)\n",
    "        strain_naam = os.path.splitext(strain_naam)[0]\n",
    "        \n",
    "        start_cs = str(start_c)\n",
    "        if len(start_cs) == 4:\n",
    "            start_cs = '0' + start_cs\n",
    "        \n",
    "        pre_op.write(start_cs\n",
    "                     + '\\t'\n",
    "                     + strain_naam\n",
    "                     + '\\t'\n",
    "                     + abs_path\n",
    "                     + '\\n'\n",
    "                    )\n",
    "        start_c += 1\n",
    "        \n",
    "pre_op.close()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "da74db3a",
   "metadata": {},
   "source": [
    "## Step2 FASTA dereplication"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2aa0c565",
   "metadata": {},
   "source": [
    "Step 2 is to remove potential sequence duplication (e.g., copies of tRNA, some cds) within each genome. This dereplication is equivalent to 100% clustering, to obtain single copy.<br>\n",
    "Step 2 **requires the tab separated table output by Step 1 as input**, and specifying a directory for dereplicated files."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6f27499a",
   "metadata": {},
   "source": [
    "### define functions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3ca4a81a",
   "metadata": {
    "code_folding": []
   },
   "outputs": [],
   "source": [
    "def redundancy_rm_a_file(path_n_id_tup, # abs path of an input fasta\n",
    "                         rred_folder, # path to 100% clustered result\n",
    "                        ):\n",
    "    \n",
    "    in_fasta = SeqIO.parse(path_n_id_tup[0],\n",
    "                           'fasta'\n",
    "                          )\n",
    "    \n",
    "    strain_naam = os.path.basename(path_n_id_tup[0])\n",
    "    strain_naam = os.path.splitext(strain_naam)[0]\n",
    "    \n",
    "    short_id = path_n_id_tup[1]\n",
    "    \n",
    "    added_seq = set()\n",
    "    to_save = []\n",
    "    gene_id_c = 0\n",
    "    \n",
    "    for SEQ in in_fasta:\n",
    "        if str(SEQ.seq) in added_seq:\n",
    "            continue\n",
    "        else:\n",
    "            to_save.append(SeqRecord(seq = SEQ.seq,\n",
    "                                     id = short_id + '-' + str(gene_id_c),\n",
    "                                     description = SEQ.description\n",
    "                                    )\n",
    "                          )\n",
    "            gene_id_c += 1\n",
    "            added_seq = added_seq.union(set([str(SEQ.seq)\n",
    "                                            ] # set a list\n",
    "                                           ) # add to set\n",
    "                                       )\n",
    "    \n",
    "    SeqIO.write(to_save, \n",
    "                os.path.join(os.path.abspath(rred_folder), \n",
    "                             strain_naam + '.fasta'\n",
    "                            ),\n",
    "                'fasta')\n",
    "\n",
    "def redundancy_rm_files(m_lst,\n",
    "                        rred_folder_\n",
    "                       ):\n",
    "    \n",
    "    for ab_fp in m_lst:\n",
    "        redundancy_rm_a_file(ab_fp, # to provide abs path\n",
    "                             rred_folder_\n",
    "                            )"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "aa8fbca2",
   "metadata": {},
   "source": [
    "### Set path\n",
    "* Set variable `preparation_output_path` to the path of output tab separated files.\n",
    "* Set variable `dereped_dir_path` to the path to directory of dereplicated files.\n",
    "* Set variable `process_number` as thread number.<br>\n",
    "Then run the codes below [Step 2 execution](###step-2-execution)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "086bb622",
   "metadata": {},
   "outputs": [],
   "source": [
    "# set your own\n",
    "# path of output tsv from Step1\n",
    "preparation_output_path = 'test_op/Step1_op.txt'\n",
    "# directory path of output\n",
    "dereped_dir_path = 'test_op/S2_op_dereped'\n",
    "# how many jobs to parallel\n",
    "process_number = 16"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bf71addc",
   "metadata": {},
   "source": [
    "### Step 2 execution"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3b0bf28c",
   "metadata": {},
   "outputs": [],
   "source": [
    "print(time.asctime())\n",
    "\n",
    "if __name__ == \"__main__\":\n",
    "    \n",
    "    # mkdir or not ----\n",
    "    if os.path.exists(dereped_dir_path):\n",
    "        pass\n",
    "    else:\n",
    "        os.mkdir(dereped_dir_path)\n",
    "        \n",
    "    # read short id ----\n",
    "    strain_short_pth_n_id = []\n",
    "    \n",
    "    with open(preparation_output_path, 'r') as pre_op:\n",
    "        rls = pre_op.readlines()\n",
    "        \n",
    "        for l in rls:\n",
    "            l_ = l.split('\\t')\n",
    "            strain_short_pth_n_id.append((l_[2][0: -1], # key remove \\n\n",
    "                                          l_[0]\n",
    "                                         )\n",
    "                                        )\n",
    "        \n",
    "    pre_op.close()\n",
    "    \n",
    "    mission_lst = mission_spliter(strain_short_pth_n_id, \n",
    "                                  process_number\n",
    "                                 )\n",
    "    # mp ----\n",
    "    jobs = []\n",
    "    \n",
    "    for sub_mission_lst in mission_lst:\n",
    "    \n",
    "        p = Process(target = redundancy_rm_files,\n",
    "                    args = (sub_mission_lst, \n",
    "                            dereped_dir_path\n",
    "                           )\n",
    "                   )\n",
    "        p.start()\n",
    "        jobs.append(p)\n",
    "\n",
    "\n",
    "    for z in jobs:\n",
    "        z.join()\n",
    "        \n",
    "print(time.asctime())"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "00f0fcab",
   "metadata": {},
   "source": [
    "### <font color=\"red\">Note before next step</font>\n",
    "After dereplication, users should give a careful check of <b>size</b> of dereplicated fasta files. It is worth noting that if a fasta file with a very low amount of sequences, getting processed together with all the rest, the final \"core\" clusters will be heavily affected and may bias your analysis.\n",
    "    \n",
    "Since the core genome construction is similar with intersection construction. <font color=\"red\"><b>It is recommend to remove some very small dereplicated fasta files BEFORE NEXT STEP</b>, e.g., remove all dereplicated <i>E.coli </i> genomes with file size lower than 2.5MB as most should be larger than 4.0MB.</font>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d4a27804",
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "!grep -c \">\" test_op/S2_op_dereped/* | sort -t: -k2rn"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c39aa986",
   "metadata": {},
   "source": [
    "## Step 3 Pre-clustering of using all dereplicated FASTAs and non-redundant genome generation"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1cc02cf4",
   "metadata": {},
   "source": [
    "Step 3 is the new feature since version `0.1` onward comparing with `0.1Beta`.It performs 100% clustering on all dereplicated FASTAs made in Step 2.<br>\n",
    "\n",
    "The program of Step 3 will take the **dereplicated fasta files made in step 3 as input**, and produce:<br>\n",
    "* dereplicated concatenated FASTA prepared for reciprocal blast\n",
    "* length of each non-redundant sequence\n",
    "* pre-clustered gene id\n",
    "* each genome that is redundancy-removed\n",
    "* gene id information: a tab delimited table, first column as short id generated by software, second colums as description line of original fasta \n",
    "\n",
    "In previous versions, program would BLAST the concatenated FASTA against each dereplicated genome sequentially, as direct all-vs-all BLAST using concatenated FASTA would be too memory intensive to run (1~ G FASTA direct all-vs-all BLAST cost roughly 120 GB memory using `-mt_mode 1` with 36 threads). <br>\n",
    "\n",
    "As tested, BLAST the concatenated FASTA against each genome (`-mt_mode 1` and 36 threads) could also be very time consuming:\n",
    "* 500 <i>Listeria monocytogenes</i> costs > 7 hours.<br>\n",
    "* 1150 <i>E. coli</i> costs > 2.7 days.<br>\n",
    "\n",
    "However, when running the program for phylogneticlly close genmoes, there would be a high duplication level in the concatenated FASTA.<br>\n",
    "\n",
    "As tested, the size of concatenated FASTA could be significantly reduced after dereplication:\n",
    "* 500 <i>Listeria monocytogenes</i> before dereplication -> ~1.33GB, after -> ~187MB, \n",
    "* 1150 <i>E. coli</i> before dereplication -> ~5.2GB, after -> ~986MB\n",
    "\n",
    "Besides dereplication on concatenated FASTA, this version would also use this dereplicated concatenaion to re-generate each genome without overall redundancy (non-redundant genome). Which would significantly reduce task labor.<br>\n",
    "\n",
    "As tested, the BLAST time usage after dereplication could be significantly reduced:\n",
    "* 500 <i>Listeria monocytogenes</i> before dereplication -> ~7 hours, after -> ~22 mins\n",
    "* 1150 <i>E. coli</i> before dereplication -> ~2.7 days, after -> less than 3.5 hours"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1b54bb12",
   "metadata": {},
   "source": [
    "### define functions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "21e7f318",
   "metadata": {},
   "outputs": [],
   "source": [
    "def save_nr_each(in_genome_lst, fasta_dict, op):\n",
    "    \n",
    "    for nr_genome in in_genome_lst:\n",
    "        saver = []\n",
    "        \n",
    "        for id_ in nr_genome: \n",
    "            saver.append(fasta_dict[id_])\n",
    "        \n",
    "        SeqIO.write(saver, \n",
    "                    os.path.join(op,\n",
    "                                 nr_genome[0].partition(\"-\")[0] + '.fasta'\n",
    "                                ), \n",
    "                    'fasta'\n",
    "                   )"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ed537598",
   "metadata": {},
   "source": [
    "### Set path\n",
    "* Set variable `dir_dereped_fasta_path` as path to directory dereplicated FASTAs file made in Step 2.\n",
    "* Set variable `cated_dereped_fasta_path` as path to dereplicated concatenated FASTA file.\n",
    "* Set variable `nr_genomes_pth` as path to directory of each non-redundant genome file.\n",
    "* Set variable `len_infor_path` as path to file of sequence length information.\n",
    "* Set variable `pre_cluster_path` as path to file pre-clustered gene id.\n",
    "* Set variable `id_info_path` as path information gene id.\n",
    "* Set variable `process_number` as thread number.<br>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1f0397af",
   "metadata": {},
   "outputs": [],
   "source": [
    "# path to concatenated fasta made in Step 3\n",
    "dir_dereped_fasta_path = 'test_op/S2_op_dereped/'\n",
    "# path to the dereplicated concatenated fasta\n",
    "cated_dereped_fasta_path = 'test_op/S3_op_dereped_cated.fasta'\n",
    "# path to the dereplicated of non-redundant genomes\n",
    "nr_genomes_pth = 'test_op/S3_op_nr_genomes/'\n",
    "# path to save sequence lenght of each sequence\n",
    "len_infor_path = 'test_op/S3_op_seq_len.txt'\n",
    "# path to save pre-cluster\n",
    "pre_cluster_path = 'test_op/S3_op_pre_cluster.txt'\n",
    "# path to save gene id info\n",
    "id_info_path = 'test_op/S3_op_id_info.txt'\n",
    "# how many jobs to parallel\n",
    "process_number = 16"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "89cf3e80",
   "metadata": {},
   "outputs": [],
   "source": [
    "print(time.asctime())\n",
    "\n",
    "length_info = open(len_infor_path, 'w')\n",
    "\n",
    "pre_cluster = defaultdict(list)\n",
    "each_d = defaultdict(list)\n",
    "\n",
    "dereped_cated_dict = {}\n",
    "\n",
    "dereped_cated = []\n",
    "\n",
    "id_info = open(id_info_path, 'w')\n",
    "# derep of all ----\n",
    "for f_naam in os.listdir(dir_dereped_fasta_path):\n",
    "    a_fasta = SeqIO.parse(os.path.join(dir_dereped_fasta_path, \n",
    "                                       f_naam),\n",
    "                          'fasta')\n",
    "    for rec in a_fasta:\n",
    "        \n",
    "        id_info.write(rec.id + '\\t' + rec.description.partition(\" \")[2] + '\\n')\n",
    "        SEQ = str(rec.seq)\n",
    "\n",
    "        if SEQ in pre_cluster.keys():# if already presence\n",
    "            pre_cluster[SEQ].append(rec.id) # must after if\n",
    "\n",
    "        else:\n",
    "            length_info.write(rec.id \n",
    "                              + '\\t' \n",
    "                              +  str(len(SEQ)) \n",
    "                              + '\\n')\n",
    "            pre_cluster[SEQ].append(rec.id)\n",
    "\n",
    "            dereped_cated.append(rec)\n",
    "            dereped_cated_dict[rec.id] = rec\n",
    "            \n",
    "            each_d[rec.id.partition(\"-\")[0]].append(rec.id)\n",
    "        \n",
    "length_info.close()\n",
    "id_info.close()\n",
    "\n",
    "# save pre cluster ----\n",
    "pre_cluster_file = open(pre_cluster_path, 'w')\n",
    "\n",
    "for _100_clus in pre_cluster.values():\n",
    "    LEN = len(_100_clus)\n",
    "    LEN_1 = LEN - 1\n",
    "    for ele in range(LEN):\n",
    "        if ele == LEN_1:\n",
    "            pre_cluster_file.write(_100_clus[ele] + '\\n')\n",
    "        else:\n",
    "            pre_cluster_file.write(_100_clus[ele] + '\\t')\n",
    "\n",
    "pre_cluster_file.close()\n",
    "\n",
    "# save dereped_cated ----\n",
    "SeqIO.write(dereped_cated, \n",
    "            cated_dereped_fasta_path,\n",
    "            'fasta')\n",
    "\n",
    "# nr genomes ----\n",
    "if os.path.exists(nr_genomes_pth):\n",
    "    pass\n",
    "else:\n",
    "    os.mkdir(nr_genomes_pth)\n",
    "    \n",
    "if __name__ == \"__main__\":\n",
    "                                       \n",
    "        \n",
    "    mission_lst = list(each_d.values())\n",
    "    \n",
    "    mission_lst = mission_spliter(mission_lst, \n",
    "                                  process_number\n",
    "                                 )\n",
    "    \n",
    "    jobs = []\n",
    "    \n",
    "    for sub_mission_lst in mission_lst:\n",
    "        \n",
    "        p = Process(target = save_nr_each,\n",
    "                    args = (sub_mission_lst,\n",
    "                            dereped_cated_dict,\n",
    "                            nr_genomes_pth\n",
    "                           )\n",
    "                   )\n",
    "            \n",
    "        p.start()\n",
    "        jobs.append(p)\n",
    "\n",
    "\n",
    "    for z in jobs:\n",
    "        z.join()\n",
    "print(time.asctime())"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9b41d4f2",
   "metadata": {},
   "source": [
    "## Step 4 Reciprocal Blast"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8127a084",
   "metadata": {},
   "source": [
    "Step 4 will carry out the Reciprocal Blast using NCBI Blast. You can get it from [NCBI official](https://ftp.ncbi.nlm.nih.gov/blast/executables/blast+/LATEST/)\n",
    " \n",
    "The pipeline will assist you to:<br>\n",
    "1. Create databases for each of all dereplicated genomes using `makeblastdb`.\n",
    "2. Using `blastn` or `blastp` to align the `dereplicated concatenated FASTA` against each of the database just made and get tabular output.\n",
    "\n",
    "### makeblastdb\n",
    "To create database to BLAST, users should provide **path to directory where all non-redundant FASTA made in Step 3 is**, and a **path to output directory where BLAST database is to store**.<br>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f6742aea",
   "metadata": {},
   "source": [
    "#### define functions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "00b9513a",
   "metadata": {},
   "outputs": [],
   "source": [
    "def make_a_db(in_abs_pth, \n",
    "              db_dir_abs_pth, \n",
    "              dbtype,\n",
    "              makeblastdb_bin):\n",
    "#     if db folder exist\n",
    "    strain_naam = os.path.basename(in_abs_pth)\n",
    "    strain_naam = os.path.splitext(strain_naam)[0]\n",
    "    \n",
    "    save_dest = os.path.join(db_dir_abs_pth, strain_naam)\n",
    "    \n",
    "    if os.path.exists(save_dest):\n",
    "        pass\n",
    "    else:\n",
    "        os.mkdir(save_dest)\n",
    "        \n",
    "    save_naam = os.path.join(save_dest, strain_naam)\n",
    "        \n",
    "    subprocess.run([makeblastdb_bin, \n",
    "                    \"-dbtype\", dbtype, \n",
    "                    \"-in\", in_abs_pth,\n",
    "                    \"-out\", save_naam,\n",
    "                    \"-parse_seqids\"\n",
    "                   ],\n",
    "                   stdout = subprocess.DEVNULL\n",
    "                  )\n",
    "    \n",
    "def make_dbs(in_abs_pth_lst,\n",
    "             db_dir_abs_pth_,\n",
    "             dbtype_,\n",
    "             makeblastdb_bin_):\n",
    "    for s in in_abs_pth_lst:\n",
    "        make_a_db(s, \n",
    "                  db_dir_abs_pth_, \n",
    "                  dbtype_,\n",
    "                  makeblastdb_bin_)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "de5a6954",
   "metadata": {},
   "source": [
    "#### Set path\n",
    "* Set variable `nr_genomes_pth` as path to directory of each non-redundant genome file.\n",
    "* Set variable `blastdb_dir_path` as path to directory to store generated databases.\n",
    "* Set variable `process_number` thread number.\n",
    "* Set variable `dbt` as the tpye of database to make.\n",
    "* Set variable `mbdb` to the path to makeblastdb command.\n",
    "\n",
    "Then run the codes below [makeblastdb](###makeblastdb)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e81b8236",
   "metadata": {},
   "source": [
    "In case you have installed your blast but not exported to `$PATH`, you can simply input `whereis blastn` or `whereis makeblastdb` to get the full path to your blast binary file.<br>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b80fa9b6",
   "metadata": {},
   "outputs": [],
   "source": [
    "! whereis makeblastdb"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "72bba58c",
   "metadata": {},
   "outputs": [],
   "source": [
    "# set your own\n",
    "# directory path of output\n",
    "nr_genomes_pth = 'test_op/S3_op_nr_genomes/'\n",
    "# directory path of blastdb\n",
    "blastdb_dir_path = 'test_op/S4_op_dbs'\n",
    "# how many jobs to parallel\n",
    "process_number = 16\n",
    "# set dbtype: nucl or prot\n",
    "dbt = 'nucl'\n",
    "# set path to makblastdb bin file\n",
    "mbdb = 'makeblastdb'"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6a97c380",
   "metadata": {},
   "source": [
    "#### makeblastdb"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "940f8246",
   "metadata": {},
   "outputs": [],
   "source": [
    "print(time.asctime())\n",
    "\n",
    "if __name__ == \"__main__\":\n",
    "    \n",
    "    # mkdir or not\n",
    "    if os.path.exists(blastdb_dir_path):\n",
    "        pass\n",
    "    else:\n",
    "        os.mkdir(blastdb_dir_path)\n",
    "    \n",
    "    mission_lst = os.listdir(nr_genomes_pth)\n",
    "    \n",
    "    mission_lst = [os.path.join(nr_genomes_pth, x) for x in mission_lst]\n",
    "    \n",
    "    mission_lst = mission_spliter(mission_lst, \n",
    "                                  process_number\n",
    "                                 )\n",
    "    \n",
    "    # mp\n",
    "    jobs = []\n",
    "    \n",
    "    for sub_mission_lst in mission_lst:\n",
    "    \n",
    "        p = Process(target = make_dbs,\n",
    "                    args = (sub_mission_lst, \n",
    "                            blastdb_dir_path, \n",
    "                            dbt,\n",
    "                            mbdb\n",
    "                           )\n",
    "                   )\n",
    "        p.start()\n",
    "        jobs.append(p)\n",
    "\n",
    "\n",
    "    for z in jobs:\n",
    "        z.join()\n",
    "        \n",
    "print(time.asctime())"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "aa3263e0",
   "metadata": {},
   "source": [
    "### Reciprocal Blast"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c6f81848",
   "metadata": {},
   "source": [
    "To perform reciprocal BLAST, users should provide **path to dereplicated concatenated FASTA producd by step 3**, **path to directory where databases made by `makeblastdb`**, and a **path to output directory where BLAST tabular output** is to store.<br>\n",
    "\n",
    "The reason to BLAST against each database sequentially rather than directly using an all-vs-all approach is to reduce computational overhead. This can be very useful if the task involves many genomes. For example, if you have 1000 dereplicated genomes to analyze, the total size of concatenated FASTA may reach 5-10 GB. A multi-threaded BLAST job using the `-mt_mode 1` by all-vs-all style could be too memory-intensive to run for such a large dataset.\n",
    "\n",
    "In addition, sequentially running BLAST will produce one tabular output per database. This will be a better adaptation for the job parallelization of finding reciprocal best hits in later steps, which will apply the hash binning method.\n",
    "\n",
    "Since version `0.1` onward, program allow users to choose recoprocal BLAST executed under memory efficient mode or not. Under memory efficient mode, the real-time memory usage will be much lower but **much more** time consuming.\n",
    "\n",
    "Since version `0.2.1` onward, users are allowed to set customized blast output format for other analysis requirement. If a follow up OthoSLC analysis is needed, leave the parameter as unspecified as default - `\"6 qseqid sseqid score\"` so that follow up command could parse the blast output."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "866732aa",
   "metadata": {},
   "source": [
    "#### define function"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "87a80e5f",
   "metadata": {},
   "outputs": [],
   "source": [
    "def RBB(cmd_lst_,\n",
    "        db_pth, # alread abs\n",
    "        blast_res_dir, # already abs\n",
    "       ):\n",
    "    strain_naam = db_pth.split('/')[-1]\n",
    "    \n",
    "    db_name = os.path.join(db_pth, strain_naam)\n",
    "    save_dest = os.path.join(blast_res_dir, strain_naam + '.tab')\n",
    "    \n",
    "    cmd_lst_l = copy.deepcopy(cmd_lst_)\n",
    "\n",
    "    cmd_lst_l.extend([\n",
    "        \"-db\", db_name,\n",
    "        \"-out\", save_dest\n",
    "    ])\n",
    "    \n",
    "    subprocess.run(cmd_lst_l)\n",
    "    \n",
    "def RBB_high(cmd_lst_,\n",
    "             db_pth_lst, # alread abs\n",
    "             blast_res_dir_, # already abs\n",
    "             pro_id_, \n",
    "             s_t_d,\n",
    "             st_lock_,\n",
    "             available_threads):\n",
    "    for db in db_pth_lst:\n",
    "        st_lock_.acquire()\n",
    "        if sum(s_t_d.values()) >= available_threads:\n",
    "            useable_th = s_t_d[pro_id_]\n",
    "#             print(str(pro_id_) + ' uses ' + str(useable_th))\n",
    "        else:\n",
    "            s_t_d[pro_id_] = s_t_d[pro_id_] + (available_threads - sum(s_t_d.values()\n",
    "                                                                       )\n",
    "                                                )\n",
    "            useable_th = s_t_d[pro_id_]\n",
    "#             print(str(pro_id_) + ' uses ' + str(useable_th))\n",
    "        st_lock_.release()\n",
    "    \n",
    "        cmd_lst_l = copy.deepcopy(cmd_lst_)\n",
    "    \n",
    "        cmd_lst_l.extend([\"-num_threads\", str(useable_th)])\n",
    "        RBB(\n",
    "            cmd_lst_l,\n",
    "            db, # alread abs\n",
    "            blast_res_dir_, # already abs\n",
    "        )\n",
    "\n",
    "    st_lock_.acquire()\n",
    "    s_t_d[pro_id_] = 0\n",
    "#     print(str(pro_id_) + ' finished')\n",
    "    st_lock_.release()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "dd18271c",
   "metadata": {},
   "source": [
    "#### Set path\n",
    "* Set variable `cated_dereped_fasta_path` as path to dereplicated concatenated FASTA file.\n",
    "* Set variable `blastdb_dir_path` as path to directory to store generated databases.\n",
    "* Set variable `process_number` thread number.\n",
    "* Set variable `E_value` as e value parameter during BLAST.\n",
    "* Set variable `blast_bin_` to the path to BLAST command.\n",
    "* Set variable `blast_op_dir` to the path of directory to save BLAST output.\n",
    "* Set variable `mem_eff` as Ture or False to run under memory efficient mode or not.\n",
    "* Set variable `outfmt` to specify output format if other format needed.\n",
    "\n",
    "Then run the codes below [BLAST](###BLAST)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b4dd0748",
   "metadata": {},
   "outputs": [],
   "source": [
    "! whereis blastn"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "03f1f256",
   "metadata": {},
   "outputs": [],
   "source": [
    "# set your own\n",
    "# path to concatenated fasta\n",
    "cated_dereped_fasta_path = 'test_op/S3_op_dereped_cated.fasta'\n",
    "# directory path of blastdb\n",
    "blastdb_dir_path = 'test_op/S4_op_dbs/'\n",
    "# how many jobs to parallel\n",
    "process_number = 16\n",
    "# set e_value\n",
    "E_value = 1e-5\n",
    "# set path to blast bin file\n",
    "blast_bin_ = 'blastn'\n",
    "# set blast output directory\n",
    "blast_op_dir = 'test_op/S4_op_blast_op'\n",
    "# memory efficient mode \n",
    "mem_eff_mode = False\n",
    "# blast output format\n",
    "outfmt = str(\"6 qseqid sseqid score\")\n",
    "# blastp task, select from <'blastp' 'blastp-fast' 'blastp-short'>, unspecified means `'blastp'` as default\n",
    "blastp_task = ''\n",
    "# blastp task <blastn' 'blastn-short' 'dc-megablast' 'megablast' 'rmblastn' > , unspecified means `'megablast'` as default\n",
    "blastn_task = ''"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "de7aba86",
   "metadata": {},
   "source": [
    "#### BLAST"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "920faff1",
   "metadata": {},
   "outputs": [],
   "source": [
    "print(time.asctime())\n",
    "# mkdir or not\n",
    "if os.path.exists(blast_op_dir):\n",
    "    pass\n",
    "else:\n",
    "    os.mkdir(blast_op_dir)\n",
    "    \n",
    "    \n",
    "mission_lst = os.listdir(blastdb_dir_path)\n",
    "mission_lst = [os.path.join(blastdb_dir_path, x) for x in mission_lst]\n",
    "\n",
    "cmd_lst = [blast_bin_, \n",
    "            \"-query\", cated_dereped_fasta_path,\n",
    "            # \"-db\", db_name,\n",
    "            # \"-out\", save_dest,\n",
    "            \"-evalue\", str(E_value),\n",
    "            \"-max_hsps\", \"1\",\n",
    "            # \"-dust\", \"no\",\n",
    "            \"-outfmt\", outfmt,\n",
    "            \"-mt_mode\", \"1\"\n",
    "            # \"-num_threads\", str(t_num)\n",
    "            ]\n",
    "\n",
    "if 'blastn' in blast_bin_:\n",
    "    cmd_lst.extend([\"-dust\", \"no\"])\n",
    "    if blastn_task == '':\n",
    "        pass\n",
    "    else:\n",
    "        cmd_lst.extend([\"-task\", blastn_task])\n",
    "elif 'blastp' in blast_bin_ and blastp_task == \"\":\n",
    "    if blastp_task == '':\n",
    "        pass\n",
    "    else:\n",
    "        cmd_lst.extend([\"-task\", blastp_task])\n",
    "\n",
    "if mem_eff_mode:\n",
    "    blast_low = RBB\n",
    "\n",
    "    cmd_lst.extend([\"-num_threads\", str(process_number)])\n",
    "\n",
    "    for db in mission_lst:\n",
    "        \n",
    "        blast_low(cmd_lst,\n",
    "                  db, \n",
    "                  blast_op_dir\n",
    "                  )\n",
    "        \n",
    "else:\n",
    "    mission_lst = mission_spliter(mission_lst, \n",
    "                                  process_number\n",
    "                                 )\n",
    "    \n",
    "    blast_high = RBB_high\n",
    "\n",
    "    if __name__ == \"__main__\":\n",
    "        with Manager() as manager:\n",
    "            smart_threads_d = manager.dict()\n",
    "            \n",
    "            st_lock = Lock()\n",
    "            \n",
    "            for x in range(len(mission_lst)):\n",
    "                smart_threads_d[x] = 1\n",
    "        \n",
    "            # mp\n",
    "            jobs = []\n",
    "    \n",
    "            pro_id = 0\n",
    "            for db_lst in mission_lst:\n",
    "\n",
    "                p = Process(target = blast_high,\n",
    "                            args = (cmd_lst,\n",
    "                                    db_lst,\n",
    "                                    blast_op_dir,\n",
    "                                    pro_id,\n",
    "                                    smart_threads_d, \n",
    "                                    st_lock,\n",
    "                                    process_number,\n",
    "                                   )\n",
    "                           )\n",
    "                p.start()\n",
    "                jobs.append(p)\n",
    "                \n",
    "                pro_id += 1\n",
    "\n",
    "\n",
    "            for z in jobs:\n",
    "                z.join()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b6aadde9",
   "metadata": {},
   "source": [
    "## Step 5 query binning"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "06625672",
   "metadata": {},
   "source": [
    "This is the new feature in since version `0.1`. This step is to apply hash binning to bin all presence of a query into same file to facilitate next step filtering.\n",
    "\n",
    "<font color=\"red\"><b>Set bin level:</b></font><br>\n",
    "According to the amount of genomes to analysze, user should provide binning level, which is to set how many bins should be used. Level $L$ should be interger within range $0 < L \\le 9999$, and will generate $L$ bins. \n",
    "\n",
    "Suggestion is that do not set the bin level too high, especially when less than 200 genomes participated. If such amount of genomes participated analysis, bin level from 10 to 100 should work as most efficient way. \n",
    "\n",
    "As tested, an analysis of 30 genomes, has 30 BLAST output after step 4. \n",
    "* A bin level of 10, takes 7 seconds to finish, \n",
    "* a bin level of 100, takes 10 seconds to finish,\n",
    "* a bin level of 1000, takes 24 seconds to finish,\n",
    "\n",
    "<font color=\"red\"><b>When to set a high bin level:</b></font><br>\n",
    "Simply speaking, when you have really larger amount of genomes and not enough memory (e.g., more than 1000 genomes and less than 100 GB memory) <br>\n",
    "\n",
    "For example, if the output of BLAST for 1000 genomes reach 150 GB in size, and if the bin level is set to 10, there will be 10 bins to evenly distribute the data. On average, each bin will contain 1.5 GB of data, which may be too memory-intensive to process in step 5 (where requires approximately 1.5 GB of memory per bin). However, if the number of bins is increased to 1000, the size of each bin will be reduced to between 100-200 MB, it will then facilitate step 6 parallelization.\n",
    "\n",
    "<font color=\"red\"><b>No lock mode:</b></font><br>\n",
    "we provide **no lock mode** in all steps that apply hash binning to speed up the process. We allow users to turn off mutex lock which is to safely write into files when multi-threading. In ours tests, program can generate files without data corruption when multi-threading with no lock (data corruption were rarely observed, the possiblity of data corruption may vary between computation platform)."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bf897311",
   "metadata": {},
   "source": [
    "### define functions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "235f0ca3",
   "metadata": {},
   "outputs": [],
   "source": [
    "def q_binning(in_f_lst, bin_level, bin_op_dir, no_lock_mode, LOCKS):\n",
    "    for in_f in in_f_lst:\n",
    "        f = open(in_f, 'r')\n",
    "\n",
    "        save_dict = defaultdict(list)\n",
    "\n",
    "        while True:\n",
    "            l_ = f.readline()\n",
    "            if l_ == '':\n",
    "                break\n",
    "\n",
    "            query = l_.split('\\t')[0]\n",
    "\n",
    "            bin_ = hash(query) % bin_level\n",
    "\n",
    "            save_dict[bin_].append(l_)\n",
    "        f.close()\n",
    "        \n",
    "        # saving\n",
    "        for b in save_dict.keys():\n",
    "            if not no_lock_mode:\n",
    "                LOCKS[int(b)].acquire()\n",
    "\n",
    "            f = open(os.path.join(bin_op_dir, str(b) + '.txt'), 'a')\n",
    "\n",
    "            f.writelines(save_dict[b])\n",
    "\n",
    "            f.close()\n",
    "            if not no_lock_mode:\n",
    "                LOCKS[int(b)].release()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "533201e2",
   "metadata": {},
   "source": [
    "### Set path\n",
    "* Set variable `blast_op_dir` to the path of directory to save BLAST output.\n",
    "* Set variable `process_number` thread number.\n",
    "* Set variable `b_level` as how many bins to generate.\n",
    "* Set len_infor_path `query_bin_op_dir_` to the path of output directory.\n",
    "\n",
    "Then run the codes below [Step 5 execution](###Step-6-exeution)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1cb5dc3f",
   "metadata": {},
   "outputs": [],
   "source": [
    "# set your own\n",
    "# set blast output directory\n",
    "blast_op_dir = 'test_op/S4_op_blast_op/'\n",
    "# how many jobs to parallel\n",
    "process_number = 16\n",
    "# set bin level\n",
    "b_level = 64\n",
    "# set blast output directory\n",
    "query_bin_op_dir_ = 'test_op/S5_op'\n",
    "# set whether turn off lock\n",
    "no_lock_mode = True"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "43ae5996",
   "metadata": {},
   "source": [
    "### Step 5 execution"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "566de097",
   "metadata": {},
   "outputs": [],
   "source": [
    "if b_level < 0 or b_level > 9999 or not isinstance(b_level, int):\n",
    "    print('b_level must be integer within range 0 < L <= 9999')\n",
    "    sys.exit()\n",
    "    \n",
    "# mkdir or not\n",
    "\n",
    "if os.path.exists(query_bin_op_dir_):\n",
    "    pass\n",
    "else:\n",
    "    os.mkdir(query_bin_op_dir_)\n",
    "    \n",
    "print(time.asctime())\n",
    "\n",
    "if __name__ == \"__main__\":\n",
    "    \n",
    "    mission_lst = os.listdir(blast_op_dir)\n",
    "    \n",
    "    mission_lst = [os.path.join(blast_op_dir, x) for x in mission_lst]\n",
    "    \n",
    "    mission_lst = mission_spliter(mission_lst, \n",
    "                                  process_number\n",
    "                                 )\n",
    "    \n",
    "    # mp\n",
    "    locks = [Lock() for i in range(b_level)]\n",
    "    \n",
    "    jobs = []\n",
    "    \n",
    "    for sub_mission_lst in mission_lst:\n",
    "    \n",
    "        p = Process(target = q_binning,\n",
    "                    args = (sub_mission_lst, \n",
    "                            b_level, \n",
    "                            query_bin_op_dir_,\n",
    "                            no_lock_mode,\n",
    "                            locks\n",
    "                           )\n",
    "                   )\n",
    "        p.start()\n",
    "        jobs.append(p)\n",
    "\n",
    "\n",
    "    for z in jobs:\n",
    "        z.join()\n",
    "        \n",
    "print(time.asctime())"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c407698f",
   "metadata": {},
   "source": [
    "## Step 6 Filtering and binning"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a8b83db5",
   "metadata": {},
   "source": [
    "This step is to filter the blast output and to apply hash binning, in order to provide best preparation for reciprocal best find.<br>\n",
    "\n",
    "Step 6 requires **path to directory of query binning output (Step 5)**, **sequence length information, pre-cluster information output by Step 3** as input.\n",
    "\n",
    "The pipeline will carry out following treatment to BLAST output:\n",
    "1. Paralog removal: <br>\n",
    "If query and subject is from same strain, the hit will be skipped, as to remove paralog.\n",
    "2. Length ratio filtering:<br>\n",
    "Within a hit, query length $Q$ and subject length $S$, the ratio $v$ of this 2 length\n",
    "\n",
    "$$v = \\frac{Q}{S}$$\n",
    "\n",
    "should be within a range, according to [L. Salichos et al](https://journals.plos.org/plosone/article?id=10.1371/journal.pone.0018755), $r (0 < r \\le 1)$ is recommended to be higher than 0.3 which means the shorter sequence should not be shorter than 30% of the longer sequence:\n",
    "\n",
    "$$r \\le v \\le \\frac{1}{r}$$<br>\n",
    "\n",
    "If above condition not met, the hit will be removed from analysis.\n",
    "\n",
    "3. Non-best-hit removal: <br>\n",
    "* Identical sequences are always regarded as best hit.\n",
    "* If a query has more than 1 subject hits, only the query-subject pair with highest score will then be kept.\n",
    "* if pairs are of same score, the pair whose query and subject are of more similar length will be kept.\n",
    "4. Sorting and binning:<br>\n",
    "For every kept hit, its query and subject will be sorted using Python or C++ built in sort algorithm. This is because in a sequential blast output file, only \"<b>single direction best hit</b>\" can be obtained, its \"<b>reciprocal best hit</b>\" only exist in other files, which poses difficulty doing \"<b>repriprocal finding</b>\". <br>\n",
    "However, if a query $a$ and its best suject hit $b$, passed filter above, and form $(a, b)$, and in the mean time we sort its rericprocal hit $(b, a)$ from another file into $(a, b)$, then both $(a, b)$ will generate same hash value. This hashed value with last several digits will allow us to bin them into same new file. Therefore, after this binning, \"<b>reciprocal finding</b>\" will be turned into \"<b>duplication finding</b>\" within one same file.<br>\n",
    "\n",
    "<font color=\"red\"><b>Set bin level:</b></font><br>\n",
    "According to the amount of genomes to analysze, user should provide binning level, which is to set how many bins should be used. Level $L$ should be interger within range $0 < L \\le 9999$, and will generate $L$ bins. \n",
    "\n",
    "Suggestion is that do not set the bin level too high, especially when less than 200 genomes participated. If such amount of genomes participated analysis, bin level from 10 to 100 should work as most efficient way. \n",
    "\n",
    "<font color=\"red\"><b>When to set a high bin level:</b></font><br>\n",
    "Simply speaking, when you have really larger amount of genomes and not enough memory (e.g., more than 1000 genomes and less than 100 GB memory) <br>\n",
    "\n",
    "for example, if the output of BLAST for 1000 genomes reach 200 GB in size, and if the bin level is set to 100, there will be 100 bins to evenly distribute the data. On average, each bin will contain 1.7 GB of data, which may be too memory-intensive to process in step 7, where reciprocal find is performed (which requires approximately 1.7 GB of memory per bin). However, if the number of bins is increased to 1000, the size of each bin will be reduced to between 100-200 MB, which then facilitate step 7 parallelization.\n",
    "\n",
    "<font color=\"red\"><b>Note:</b></font><br>\n",
    "This is the one of the most computation and I/O intensive step, use the C++ based binary file to process for better efficiency."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "72027951",
   "metadata": {},
   "outputs": [],
   "source": [
    "def fileter_n_bin(blast_op_pt_lst, # to avoid redundant copy of  seqlen info, abs pth\n",
    "                  r, \n",
    "                  bin_level, \n",
    "                  bin_op_dir, # already abs\n",
    "                  seq_len_info,\n",
    "                  bin_spe, # pre_cluster\n",
    "                  no_lock_mode,\n",
    "                  LOCKS):\n",
    "    for fs in blast_op_pt_lst:\n",
    "        \n",
    "        f = open(fs, 'r')\n",
    "            \n",
    "        q_s_d = defaultdict(dict)\n",
    "        \n",
    "        sspe_m = defaultdict(lambda: defaultdict(list\n",
    "                                                )\n",
    "                            )\n",
    "        \n",
    "        while True:\n",
    "            l_ = f.readline()\n",
    "            if l_ == '':\n",
    "                break\n",
    "                \n",
    "#             l_ = l_.rstrip('\\n') '\\n' do not affect number  \n",
    "            l_ = l_.split('\\t')\n",
    "            \n",
    "            query = l_[0]\n",
    "            subject = l_[1]\n",
    "            \n",
    "            if query.partition(\"-\")[0] == subject.partition(\"-\")[0]: # if self hit\n",
    "                continue\n",
    "            else:\n",
    "                score = int(l_[2])\n",
    "                \n",
    "                q_s_d[query][subject] = score\n",
    "#                 {q1: {s1: score,\n",
    "#                       s2: score\n",
    "#                      },\n",
    "#                  q2\n",
    "#                 }\n",
    "                \n",
    "                # add all spe\n",
    "                for a_spe_of_clus in bin_spe[subject]:\n",
    "                    sspe_m[query][a_spe_of_clus].append(subject)\n",
    "#                     {\n",
    "#                         q1: {\n",
    "#                             a1_spe_in_cluster_of_sub: [subject1 id,\n",
    "#                                                     subject2 id\n",
    "#                                                    ],\n",
    "#                             a2_spe_in_cluster_of_sub: [subject3 id,\n",
    "#                                                     subject4 id\n",
    "#                                                    ],\n",
    "#                         },\n",
    "#                         q2\n",
    "#                     }\n",
    "        f.close()\n",
    "        # how to save\n",
    "        saver = defaultdict(list)\n",
    "        \n",
    "        for qs in q_s_d.keys():\n",
    "            q_len = seq_len_info[qs]\n",
    "            \n",
    "            for sub in q_s_d[qs].keys():\n",
    "                s_spe = sub.partition(\"-\")[0]\n",
    "                \n",
    "                if s_spe in bin_spe[qs]:# if qs pre cluster has that spe\n",
    "                    continue\n",
    "                    \n",
    "                else:\n",
    "                    best_hit = True\n",
    "                    \n",
    "                    same_score = []\n",
    "                    \n",
    "                    if len(sspe_m[qs][s_spe]) > 1: #if multiple hit in a spe\n",
    "                        for id_ in sspe_m[qs][s_spe]:\n",
    "                            if sub == id_: # skip self\n",
    "                                continue\n",
    "                                \n",
    "                            if q_s_d[qs][sub] < q_s_d[qs][id_]: # if there is a better hit\n",
    "                                best_hit = False\n",
    "                                break\n",
    "                            elif q_s_d[qs][sub] == q_s_d[qs][id_]: # same score \n",
    "                                same_score.append(id_)\n",
    "                    \n",
    "                    if len(same_score) > 0 and best_hit: # need more similar length # print(qs + \"<>\" + sub)\n",
    "                        \n",
    "                        qs_len_diff = abs(q_len - seq_len_info[sub])\n",
    "                        \n",
    "                        for id_ in same_score:\n",
    "                            q_id_len_diff = abs(q_len - seq_len_info[id_])\n",
    "                            \n",
    "                            if q_id_len_diff < qs_len_diff:\n",
    "#                                 print(qs + \"< >\" + sub + \"< >\" + id_)\n",
    "                                best_hit = False\n",
    "                                break\n",
    "                                \n",
    "                    if best_hit: \n",
    "                        \n",
    "                        len_ratio = q_len/seq_len_info[sub]\n",
    "                        \n",
    "                        c1 = len_ratio <= r\n",
    "                        c2 = len_ratio >= (1/r)\n",
    "                        \n",
    "                        if c1 or c2: # len ratio filter\n",
    "#                             print(qs + \"< >\" + sub)\n",
    "                            continue\n",
    "                        else:\n",
    "                            sorting = sorted((qs, sub)) \n",
    "                            \n",
    "                            bin_ = hash(tuple(sorting\n",
    "                                             )\n",
    "                                       ) % bin_level\n",
    "                            \n",
    "                            saver[str(bin_)].append(sorting[0] + '\\t' + sorting[1] + '\\n')\n",
    "                            \n",
    "        for b in saver.keys():\n",
    "            if not no_lock_mode:\n",
    "                LOCKS[int(b)].acquire()\n",
    "\n",
    "            f = open(os.path.join(bin_op_dir, b + '.txt'), 'a')\n",
    "\n",
    "            f.writelines(saver[b])\n",
    "\n",
    "            f.close()\n",
    "            \n",
    "            if not no_lock_mode:\n",
    "                LOCKS[int(b)].release()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0c0379bb",
   "metadata": {},
   "source": [
    "### Set path\n",
    "* Set variable `process_number` thread number.\n",
    "* Set variable `b_level` as how many bins to generate.\n",
    "* Set variable `r_` as the limit of length difference ratio.\n",
    "* Set variable `blast_op_dir` to the path of directory to save BLAST output.\n",
    "* Set variable `len_infor_path` as path to file of sequence length information.\n",
    "* Set variable `pre_cluster_path` as path to file pre-clustered gene id.\n",
    "\n",
    "Then run the codes below [Step 6 execution](###Step-6-execution)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d251d4c6",
   "metadata": {},
   "outputs": [],
   "source": [
    "# set your own\n",
    "# set blast output directory\n",
    "qbin_op_dir = 'test_op/S5_op/'\n",
    "# how many jobs to parallel\n",
    "process_number = 16\n",
    "# set bin level\n",
    "b_level = 48\n",
    "# set length ration limit r value \n",
    "r_ = 0.3\n",
    "# set blast output directory\n",
    "bin_op_dir_ = 'test_op/S6_op'\n",
    "# path save sequence lenght of each sequence\n",
    "len_infor_path = 'test_op/S3_op_seq_len.txt'\n",
    "# path to save pre-cluster\n",
    "pre_cluster_path = 'test_op/S3_op_pre_cluster.txt'\n",
    "# set whether turn off lock\n",
    "no_lock_mode = True"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "125b682c",
   "metadata": {},
   "source": [
    "### Step 6 execution"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ce5220dc",
   "metadata": {},
   "outputs": [],
   "source": [
    "if b_level < 0 or b_level > 9999 or not isinstance(b_level, int):\n",
    "    print('b_level must be integer within range 0 < L <= 9999')\n",
    "    sys.exit()\n",
    "\n",
    "# mkdir or not\n",
    "\n",
    "if os.path.exists(bin_op_dir_):\n",
    "    pass\n",
    "else:\n",
    "    os.mkdir(bin_op_dir_)\n",
    "\n",
    "# seq len\n",
    "seq_len_d = {}\n",
    "seq_len = open(len_infor_path, 'r')\n",
    "seq_len_rls = seq_len.readlines()\n",
    "seq_len.close()\n",
    "\n",
    "for ros in seq_len_rls:\n",
    "    l_ = ros.rstrip('\\n').split('\\t')\n",
    "    seq_len_d[l_[0]] = float(l_[1])\n",
    "\n",
    "\n",
    "# bin spe\n",
    "bin_spe = {}\n",
    "\n",
    "f = open(pre_cluster_path, 'r')\n",
    "\n",
    "while True:\n",
    "    l_ = f.readline()\n",
    "    if l_ == '':\n",
    "        break\n",
    "    \n",
    "    \n",
    "    l_ = l_.rstrip('\\n').split('\\t')\n",
    "    \n",
    "    bin_ = l_[0]\n",
    "    \n",
    "    bin_spe[bin_] = [x.partition(\"-\")[0] for x in l_]\n",
    "    \n",
    "\n",
    "f.close()\n",
    "\n",
    "print(time.asctime())\n",
    "\n",
    "if __name__ == \"__main__\":\n",
    "    \n",
    "    mission_lst = os.listdir(qbin_op_dir)\n",
    "    \n",
    "    mission_lst = [os.path.join(qbin_op_dir, x) for x in mission_lst]\n",
    "    \n",
    "    mission_lst = mission_spliter(mission_lst, \n",
    "                                  process_number\n",
    "                                 )\n",
    "    \n",
    "    # mp\n",
    "    locks = [Lock() for i in range(b_level)]\n",
    "    \n",
    "    jobs = []\n",
    "    \n",
    "    for sub_mission_lst in mission_lst:\n",
    "    \n",
    "        p = Process(target = fileter_n_bin,\n",
    "                    args = (sub_mission_lst, \n",
    "                            r_, \n",
    "                            b_level,\n",
    "                            bin_op_dir_,\n",
    "                            seq_len_d,\n",
    "                            bin_spe,\n",
    "                            no_lock_mode,\n",
    "                            locks\n",
    "                           )\n",
    "                   )\n",
    "        p.start()\n",
    "        jobs.append(p)\n",
    "\n",
    "\n",
    "    for z in jobs:\n",
    "        z.join()\n",
    "        \n",
    "print(time.asctime())"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c88d89fe",
   "metadata": {},
   "source": [
    "## Step 7 Reciprocal Best find"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4b3e6332",
   "metadata": {},
   "source": [
    "This Step is to find reciprocal best hits. In `Step 6`, query-subject pairs had been binned into different files according to their hash value, therefore, pair $(a, b)$ and its reciprocal pair $(b, a)$ (which was sorted into $(a, b)$), will be in the same bin. Thus, a pair found twice in a bin will be reported as a reciprocal best blast pair.\n",
    "\n",
    "In addition, Step 7 also does hash binning after a reciprocal best hit is comfirmed. Query-subject pairs will be binned by the hash value of query ID, which then put pairs with common elements into same bin to assist faster clustering in next step.\n",
    "\n",
    "Step 7 requires **path to directory of bins output by Step 6**, and path to output directory.\n",
    "\n",
    "\n",
    "<font color=\"red\"><b>Set bin level:</b></font><br>\n",
    "According to the amount of genomes to analysze, user should provide binning level, which is to set how many bins should be used. Level $0 < L \\le 9999$, and will generate $L$ bins. \n",
    "\n",
    "Suggestion is that do not set the bin level too high, especially when less than 200 genomes participated. If such amount of genomes participated analysis, bin level from 10 to 100 should work as most efficient way. \n",
    "\n",
    "As tested, 30 genomes, if 10 bins generated by Step 7: \n",
    "* A bin level of 10, takes 1 seconds to finish, \n",
    "* a bin level of 100, takes 2 seconds to finish,\n",
    "* a bin level of 1000, takes 7 seconds to finish,\n",
    "\n",
    "<font color=\"red\"><b>When to set a high bin level:</b></font><br>\n",
    "Simply speaking, when you have really larger amount of genomes and not enough memory (e.g., more than 1000 genomes and less than 100 GB memory) <br>\n",
    "\n",
    "Less bins could make step 7 faster, but step 8 more memory intensive.\n",
    "\n",
    "<font color=\"red\"><b>Note:</b></font><br>\n",
    "This is the one of the most computation and I/O intensive step, use the C++ based binary file to process for better efficiency."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0f993312",
   "metadata": {},
   "source": [
    "### define functions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f3684157",
   "metadata": {
    "code_folding": []
   },
   "outputs": [],
   "source": [
    "def dup_find_mem(in_file_lst, RBB_op_dir, bin_level, no_lock_mode, LOCKS):\n",
    "    for in_file in in_file_lst:\n",
    "        naam = os.path.basename(in_file)\n",
    "\n",
    "        f = open(in_file, 'r')\n",
    "\n",
    "        rec = set()\n",
    "        d_bin = defaultdict(list)\n",
    "\n",
    "        while True:\n",
    "            l_ = f.readline()\n",
    "            if l_ == '':\n",
    "                break\n",
    "\n",
    "            if l_ in rec:\n",
    "                hash_res = str(hash(l_.split('\\t')[0]\n",
    "                                   ) % bin_level\n",
    "                              )\n",
    "                d_bin[hash_res].append(l_)\n",
    "                rec.remove(l_)\n",
    "            else:\n",
    "                rec |= set([l_])\n",
    "\n",
    "        f.close()\n",
    "\n",
    "        for b in d_bin.keys():\n",
    "            if not no_lock_mode:\n",
    "                LOCKS[int(b)].acquire()\n",
    "\n",
    "            f2 = open(os.path.join(RBB_op_dir, b + '.txt'), 'a')\n",
    "\n",
    "            f2.writelines(d_bin[b])\n",
    "\n",
    "            f2.close()\n",
    "            if not no_lock_mode:\n",
    "                LOCKS[int(b)].release()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c9b583b2",
   "metadata": {},
   "source": [
    "#### Set path\n",
    "* Set variable `bin_op_dir_` as path to directory to bins generated by Step 6.\n",
    "* Set variable `RBB_op_dir_pth` as path to output directory.\n",
    "* Set variable `process_number` thread number.\n",
    "* Set variable `b_level` as how many bins to generate\n",
    "\n",
    "Then run the codes below [Step 7 execution](###Step-7-execution)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "706fe1aa",
   "metadata": {},
   "outputs": [],
   "source": [
    "# set blast output directory\n",
    "bin_op_dir_ = 'test_op/S6_op/'\n",
    "# set reciprocal best output directory\n",
    "RBB_op_dir_pth = 'test_op/S7_op'\n",
    "# how many jobs to parallel\n",
    "process_number = 16\n",
    "# set bin level\n",
    "b_level = 36\n",
    "# set whether turn off lock\n",
    "no_lock_mode = True"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "911cf269",
   "metadata": {},
   "source": [
    "### Step 7 execution "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "472d946b",
   "metadata": {},
   "outputs": [],
   "source": [
    "if b_level < 0 or b_level > 9999 or not isinstance(b_level, int):\n",
    "    print('b_level must be integer within range 0 < L <= 9999')\n",
    "    sys.exit()\n",
    "\n",
    "# mkdir or not\n",
    "if os.path.exists(RBB_op_dir_pth):\n",
    "    pass\n",
    "else:\n",
    "    os.mkdir(RBB_op_dir_pth)\n",
    "\n",
    "print(time.asctime())\n",
    "\n",
    "if __name__ == \"__main__\":\n",
    "    \n",
    "    mission_lst = os.listdir(bin_op_dir_)\n",
    "    \n",
    "    mission_lst = [os.path.join(bin_op_dir_, x) for x in mission_lst]\n",
    "    \n",
    "    mission_lst = mission_spliter(mission_lst, \n",
    "                                  process_number\n",
    "                                 )\n",
    "    \n",
    "    # mp\n",
    "    locks = [Lock() for i in range(b_level)\n",
    "            ]\n",
    "             \n",
    "    jobs = []\n",
    "    \n",
    "    for sub_mission_lst in mission_lst:\n",
    "        \n",
    "        p = Process(target = dup_find_mem,\n",
    "                    args = (sub_mission_lst,\n",
    "                            RBB_op_dir_pth,\n",
    "                            b_level,\n",
    "                            no_lock_mode,\n",
    "                            locks\n",
    "                           )\n",
    "                   )\n",
    "            \n",
    "        p.start()\n",
    "        jobs.append(p)\n",
    "\n",
    "\n",
    "    for z in jobs:\n",
    "        z.join()\n",
    "        \n",
    "print(time.asctime())"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ad7800e9",
   "metadata": {},
   "source": [
    "## Step 8 Single Linkage Clustering"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8d738370",
   "metadata": {},
   "source": [
    "This step will carry out single linkage clustering on output from step 7. Users may perform \"<b>multi-step-to-final</b>\" or \"<b>one-step-to-final</b>\" clustering by adjusting the `compression_size` parameter. In the output files, each row is a cluster (stopped by \"\\n\") and each gene ID is separated by \"\\t\".\n",
    "\n",
    "In case that large amount genomes participated analysis, it could be memory intensive to reach final cluster in a single step. The pipeline provide ability to extenuate such pressure by reaching final cluster with multiple steps. For example, if `compression_size` = 5 is provided, program will perform clustering using 5 files at a time and shrink the output file number by a factor of 5.\n",
    "\n",
    "<font color=\"red\"><b>Note Before Start</b></font><br>\n",
    "User <font color=\"red\"><b>must</b></font> specify the path to `pre-cluster file` produced in `Step 3`, when running the <font color=\"red\"><b>LAST step of multi step to final</b></font>, or when running direct <font color=\"red\"><b>one step to final</b></font>."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "eb3e5b87",
   "metadata": {},
   "source": [
    "### define functions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2c44acc2",
   "metadata": {},
   "outputs": [],
   "source": [
    "def connect_bin(bin_ind_lst, ind):\n",
    "    while bin_ind_lst[ind] != ind:\n",
    "        ind = bin_ind_lst[ind]\n",
    "    return ind\n",
    "\n",
    "def ptr_like_merger(in_p_lst, op_dir):\n",
    "    \n",
    "    ptr_coll = dict()\n",
    "    bin_inds_connector = []\n",
    "    \n",
    "    to_save_lst = []\n",
    "    \n",
    "    const_ind_every = 0\n",
    "    \n",
    "    for in_p in in_p_lst[1]:\n",
    "        f = open(in_p, 'r')\n",
    "    \n",
    "        while True:\n",
    "            l_ = f.readline()\n",
    "\n",
    "            if l_ == '':\n",
    "                break\n",
    "\n",
    "            l_ = set(l_.rstrip('\\n').split('\\t'))\n",
    "\n",
    "            ind_to_go = const_ind_every\n",
    "            to_save_lst.append(l_)\n",
    "            bin_inds_connector.append(ind_to_go)\n",
    "\n",
    "            for ele in l_:\n",
    "                if ele not in ptr_coll.keys():\n",
    "                    ptr_coll[ele] = ind_to_go\n",
    "\n",
    "                else:\n",
    "                    its_lowest_bin_ind = connect_bin(bin_inds_connector, \n",
    "                                                     ptr_coll[ele]\n",
    "                                                    )\n",
    "\n",
    "                    if its_lowest_bin_ind == ind_to_go:\n",
    "                        continue\n",
    "\n",
    "                    if its_lowest_bin_ind > ind_to_go:\n",
    "                        its_lowest_bin_ind, ind_to_go = ind_to_go, its_lowest_bin_ind\n",
    "\n",
    "                    to_save_lst[its_lowest_bin_ind] |= to_save_lst[ind_to_go]\n",
    "                    to_save_lst[ind_to_go] = None\n",
    "\n",
    "                    bin_inds_connector[ind_to_go] = its_lowest_bin_ind\n",
    "                    ind_to_go = its_lowest_bin_ind\n",
    "\n",
    "\n",
    "            const_ind_every += 1\n",
    "\n",
    "        f.close()\n",
    "        \n",
    "    f2 = open(os.path.join(op_dir, str(in_p_lst[0]) + '.txt'), 'w')\n",
    "    for i in to_save_lst:\n",
    "        if i:\n",
    "            i = tuple(i)\n",
    "            Len = len(i)\n",
    "            for j in list(range(Len)):\n",
    "                if j == Len - 1: \n",
    "                    f2.write(i[j] + '\\n')\n",
    "                else:\n",
    "                    f2.write(i[j] + '\\t')\n",
    "        else:\n",
    "            pass\n",
    "    f2.close()\n",
    "                \n",
    "def files_clustering_s(in_p_lst_lst, op_dir_):\n",
    "    \n",
    "    for inp_lst in in_p_lst_lst:\n",
    "        ptr_like_merger(inp_lst, op_dir_)\n",
    "        \n",
    "def SLC(ip, op, proc_num, c_s, pre_clus_pth = ''):\n",
    "    if os.path.exists(op):\n",
    "        pass\n",
    "    else:\n",
    "        os.mkdir(op)\n",
    "\n",
    "    print(time.asctime())\n",
    "\n",
    "    if __name__ == \"__main__\":\n",
    "\n",
    "        mission_lst = os.listdir(ip)\n",
    "\n",
    "        if c_s > len(mission_lst):\n",
    "            print('Error: Compression size can only be smaller than ' + str(len(mission_lst)))\n",
    "        else:\n",
    "\n",
    "            mission_lst = [os.path.join(ip, x) for x in mission_lst]\n",
    "            \n",
    "            if pre_clus_pth != '':\n",
    "                mission_lst.append(pre_clus_pth)\n",
    "            \n",
    "            per = len(mission_lst)//c_s\n",
    "\n",
    "            mission_lst = mission_spliter(mission_lst, \n",
    "                                          per\n",
    "                                         )\n",
    "\n",
    "            x = 0\n",
    "            mission_lst_ = []\n",
    "            for i in mission_lst:\n",
    "                mission_lst_.append([x, i])\n",
    "                x += 1\n",
    "\n",
    "            if len(mission_lst_) < proc_num:\n",
    "                proc_num = len(mission_lst_)\n",
    "\n",
    "\n",
    "            mission_lst_ = mission_spliter(mission_lst_, \n",
    "                                           proc_num\n",
    "                                          )\n",
    "            # mp\n",
    "            jobs = []\n",
    "            \n",
    "\n",
    "            for sub_mission_lst in mission_lst_:\n",
    "\n",
    "                p = Process(target = files_clustering_s,\n",
    "                            args = (sub_mission_lst,\n",
    "                                    op\n",
    "                                   )\n",
    "                           )\n",
    "                p.start()\n",
    "                jobs.append(p)\n",
    "\n",
    "\n",
    "            for z in jobs:\n",
    "                z.join()\n",
    "\n",
    "        print(time.asctime())"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "84acc2a9",
   "metadata": {},
   "source": [
    "### Example to reach final cluster by multi step, starting with 1000 input files."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b9782ded",
   "metadata": {},
   "source": [
    "#### Set path\n",
    "* Set variable `RBB_op_dir_pth` as path to output directory.\n",
    "* Set variable `SLC_inter_op_dir_pth` as path to output directory.\n",
    "* Set variable `process_number` thread number.\n",
    "* Set variable `compression_size` as how many files to cluster into one file.\n",
    "* Set variable `pre_cluster_path` as path to file pre-clustered gene id."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ffeb35a6",
   "metadata": {},
   "outputs": [],
   "source": [
    "# set reciprocal best output directory\n",
    "RBB_op_dir_pth = 'test_op/S7_op/'\n",
    "# set path to directory of clustering output\n",
    "SLC_inter_op_dir_pth = 'test_op/SLC_1/'\n",
    "# how many jobs to parallel\n",
    "process_number = 16\n",
    "# Compression size\n",
    "compression_size = 1"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "42455e3c",
   "metadata": {},
   "source": [
    "Then run the codes below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5f381a6e",
   "metadata": {},
   "outputs": [],
   "source": [
    "SLC(RBB_op_dir_pth, SLC_inter_op_dir_pth, process_number, compression_size)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "29dedd17",
   "metadata": {},
   "source": [
    "#### Set path\n",
    "* Set variable `ip_dir_pth` as path to output directory of last round clustering.\n",
    "* Set variable `SLC_inter_op_dir_pth` as path to output directory.\n",
    "* Set variable `process_number` thread number.\n",
    "* Set variable `compression_size` as how many files to cluster into one file."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "932f7639",
   "metadata": {},
   "outputs": [],
   "source": [
    "# set path to last step clustering output directory\n",
    "ip_dir_pth = 'test_op/SLC_1/'\n",
    "# set path to directory of clustering output\n",
    "SLC_inter_op_dir_pth = 'test_op/SLC_2/'\n",
    "# how many jobs to parallel\n",
    "process_number = 16\n",
    "# Compression size\n",
    "compression_size = 2"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "32ff543b",
   "metadata": {},
   "source": [
    "Then run the codes below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "40f032b5",
   "metadata": {},
   "outputs": [],
   "source": [
    "SLC(ip_dir_pth, SLC_inter_op_dir_pth, process_number, compression_size)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0bc9f47f",
   "metadata": {},
   "source": [
    "#### Set path\n",
    "* Set variable `ip_dir_pth` as path to output directory of last round clustering.\n",
    "* Set variable `SLC_inter_op_dir_pth` as path to output directory.\n",
    "* Set variable `process_number` thread number.\n",
    "* Set variable `compression_size` as how many files to cluster into one file.\n",
    "* Set variable `pre_cluster_path` path to pre-clustered id file made in **step 3**."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ff1d2dee",
   "metadata": {},
   "outputs": [],
   "source": [
    "# set path to last step clustering output directory\n",
    "ip_dir_pth = 'test_op/SLC_2/'\n",
    "# set path to directory of each clustering output\n",
    "SLC_inter_op_dir_pth = 'test_op/SLC_final/'\n",
    "# how many jobs to parallel\n",
    "process_number = 16\n",
    "# Compression size\n",
    "compression_size = len(os.listdir(ip_dir_pth))\n",
    "# path to pre-cluster\n",
    "pre_cluster_path = 'test_op/S3_op_pre_cluster.txt'"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1dc3d822",
   "metadata": {},
   "source": [
    "Then run the codes below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "41851593",
   "metadata": {},
   "outputs": [],
   "source": [
    "SLC(ip_dir_pth, SLC_inter_op_dir_pth, process_number, compression_size, pre_cluster_path)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f4c02193",
   "metadata": {},
   "source": [
    "## Step 9 Write clusters into FASTA\n",
    "In Step 9, program will help user to generate FASTA file for each cluster. By providing the final one cluster file generated by Step 8 as input, program produces 3 types of clusters into 3 directories separately.<br>\n",
    "\n",
    "Noteably, those genomes \n",
    "* depreplicated in **Step 2**,  \n",
    "* not removed because of too low genome size\n",
    "* participated processes up to this step, \n",
    "\n",
    "are used to separate 3 types of clusters.<br>\n",
    "\n",
    "1. In drectory `accessory_cluster`(a cluster not shared by all genomes), FASTA files of clusters, which **do not have genes from all genomes** participated analysis, will be output in this drectory. For example, there are 100 genomes in analysis, a cluster with less than 100 genes will have its FASTA output here. Also, if a cluster has >= 100 genes, but all these genes are from less than 100 genomes, its FASTA will be in this directory.\n",
    "2. In drectory `strict_core`, each cluster has **exactly 1 gene from every genome** to analyze. Such clusters will have their FASTA files here.\n",
    "3. In drectory `surplus_core`, each cluster has **at least 1 gene from every genome** to analyze, and **some genomes has more than 1 genes** in this cluster. Such clusters will have their FASTA files here.\n",
    "\n",
    "This step also requires the concatenated FASTA made in Step 3 as input."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b1789811",
   "metadata": {},
   "outputs": [],
   "source": [
    "def write_a_cluster(a_cluster_of_genes, \n",
    "                    total_amount, \n",
    "                    which_to_write, \n",
    "                    op_path,\n",
    "                    fasta_dict,\n",
    "                    id_cluster_d,\n",
    "                    id_info_d\n",
    "                    ):\n",
    "    spe_count = len(set([x.partition(\"-\")[0] for x in a_cluster_of_genes]))\n",
    "    \n",
    "    if spe_count < total_amount:\n",
    "        if which_to_write[\"accessory_cluster\"]:\n",
    "            file_naam = os.path.join(op_path, 'accessory_cluster/' \n",
    "                                     + a_cluster_of_genes[0] \n",
    "                                     + '.fasta'\n",
    "                                    )\n",
    "        else:\n",
    "            return\n",
    "    else:\n",
    "        if len(a_cluster_of_genes) == total_amount:\n",
    "            if which_to_write[\"strict_core\"]:\n",
    "                file_naam = os.path.join(op_path, 'strict_core/' \n",
    "                                         + a_cluster_of_genes[0] \n",
    "                                         + '.fasta'\n",
    "                                        )\n",
    "            else:\n",
    "                return\n",
    "        else:\n",
    "            if which_to_write[\"surplus_core\"]:\n",
    "                file_naam = os.path.join(op_path, 'surplus_core/' \n",
    "                                         + a_cluster_of_genes[0] \n",
    "                                         + '.fasta'\n",
    "                                        )\n",
    "            else:\n",
    "                return\n",
    "    \n",
    "    to_save = []\n",
    "    for  recs in a_cluster_of_genes:\n",
    "        to_save.append(SeqRecord(id = recs, \n",
    "                                 description = id_info_d[recs],\n",
    "                                 seq = fasta_dict[id_cluster_d[recs]].seq\n",
    "                                )\n",
    "                      )\n",
    "        \n",
    "    SeqIO.write(to_save, \n",
    "                file_naam, \n",
    "                'fasta')\n",
    "\n",
    "def write_cluster_s(lst_of_a_clusters, \n",
    "                    total_amount_, \n",
    "                    which_to_write_, \n",
    "                    op_path_, \n",
    "                    fasta_dict_,\n",
    "                    id_cluster_d_,\n",
    "                    id_info_d_):\n",
    "    \n",
    "    for a_cl in lst_of_a_clusters:\n",
    "        write_a_cluster(a_cl, \n",
    "                        total_amount_,\n",
    "                        which_to_write_,\n",
    "                        op_path_,\n",
    "                        fasta_dict_,\n",
    "                        id_cluster_d_,\n",
    "                        id_info_d_)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d4959d6e",
   "metadata": {},
   "source": [
    "#### Set path\n",
    "* Set variable `final_cluster_path` as path to the final cluster file made by Step 8.\n",
    "* Set variable `output_path` as path to output directory.\n",
    "* Set variable `process_number` thread number.\n",
    "* Set variable `cluster_type` select from `accessory/strict/surplus`, separate by comma (`,`).\n",
    "\n",
    "\n",
    "* Set variable `dereped_fasta_path` as path to globally dereplicated FASTAs made in Step 3.\n",
    "* Set variable `id_info_path` path to save gene id info made in Step3.\n",
    "* Set variable `pre_cluster_path` path to pre-clustered gene id info made in Step3."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c4c3da1c",
   "metadata": {},
   "outputs": [],
   "source": [
    "# set path to the final cluster file made by Step 7\n",
    "final_cluster_path = 'test_op/SLC_final/0.txt'\n",
    "# set path to output directory\n",
    "output_path = 'test_op/S9_write_fasta'\n",
    "# how many jobs to parallel\n",
    "process_number = 16\n",
    "\n",
    "# select at least 1 from `accessory/strict/surplus`, separate by comma (`,`), \n",
    "# specifying types of clusters to write\n",
    "cluster_type = 'accessory,strict,surplus'\n",
    "\n",
    "# set path to concatenated dereped FASTA file made in Step 3.\n",
    "dereped_fasta_path = 'test_op/S3_op_dereped_cated.fasta'\n",
    "# set path to save gene id info made in Step3\n",
    "id_info_path = 'test_op/S3_op_id_info.txt'\n",
    "# set path to save pre-cluster\n",
    "pre_cluster_path = 'test_op/S3_op_pre_cluster.txt'"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e1100006",
   "metadata": {},
   "source": [
    "* Set variable `total_amount` of the how many genomes participated analysis."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8a19fac1",
   "metadata": {},
   "outputs": [],
   "source": [
    "# directory path of output\n",
    "print(len(os.listdir(dir_dereped_fasta_path)\n",
    "         )\n",
    "     )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b3791d82",
   "metadata": {},
   "outputs": [],
   "source": [
    "# how many genomes participated analysis.\n",
    "total_amount = len(os.listdir(dir_dereped_fasta_path)\n",
    "                      )"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e1cd860a",
   "metadata": {},
   "source": [
    "#### run the codes below get final FASTA files"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4c2aea60",
   "metadata": {
    "code_folding": [],
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "print(time.asctime())\n",
    "# mkdir or not\n",
    "if os.path.exists(output_path):\n",
    "    pass\n",
    "else:\n",
    "    os.mkdir(output_path)\n",
    "    \n",
    "\n",
    "# readin ----\n",
    "cated_fasta = SeqIO.to_dict(SeqIO.parse(dereped_fasta_path, \n",
    "                                        'fasta')\n",
    "                           )\n",
    "\n",
    "# id: cluster_id dict of 100% same seqs ----\n",
    "f = open(pre_cluster_path, 'r')\n",
    "id_cluster = f.readlines()\n",
    "f.close()\n",
    "\n",
    "id_cluster_dict = {}\n",
    "for cls in id_cluster:\n",
    "    l_ = cls.rstrip('\\n').split('\\t')\n",
    "    \n",
    "    for a_id in l_:\n",
    "        id_cluster_dict[a_id] = l_[0]\n",
    "\n",
    "# id: sqe_info ----\n",
    "f = open(id_info_path, 'r')\n",
    "id_info = f.readlines()\n",
    "f.close()\n",
    "\n",
    "id_info_dict = {}\n",
    "for pair in id_info:\n",
    "    l_ = pair.rstrip('\\n').split('\\t')\n",
    "    id_info_dict[l_[0]] = l_[1]\n",
    "        \n",
    "# which to write\n",
    "types_to_write = {\n",
    "    'accessory_cluster': False,\n",
    "    'strict_core': False,\n",
    "    'surplus_core': False\n",
    "} \n",
    "for x in cluster_type.split(','):\n",
    "    if x == 'accessory':\n",
    "        types_to_write['accessory_cluster'] = True\n",
    "        if os.path.exists(os.path.join(output_path, 'accessory_cluster')):\n",
    "            pass\n",
    "        else:\n",
    "            os.mkdir(os.path.join(output_path, 'accessory_cluster'))\n",
    "    \n",
    "    elif x == 'strict':\n",
    "        types_to_write['strict_core'] = True\n",
    "        if os.path.exists(os.path.join(output_path, 'strict_core')):\n",
    "            pass\n",
    "        else:\n",
    "            os.mkdir(os.path.join(output_path, 'strict_core'))\n",
    "            \n",
    "    elif x == 'surplus':\n",
    "        types_to_write['surplus_core'] = True\n",
    "        if os.path.exists(os.path.join(output_path, 'surplus_core')):\n",
    "            pass\n",
    "        else:\n",
    "            os.mkdir(os.path.join(output_path, 'surplus_core'))\n",
    "        \n",
    "    \n",
    "# read cluster\n",
    "f = open(final_cluster_path, 'r')\n",
    "CLUSTERS = f.readlines()\n",
    "f.close()\n",
    "\n",
    "if __name__ == \"__main__\":\n",
    "    \n",
    "    mission_lst = [x.rstrip('\\n').split('\\t') for x in CLUSTERS]\n",
    "    \n",
    "    mission_lst = mission_spliter(mission_lst, \n",
    "                                  process_number\n",
    "                                 )\n",
    "    \n",
    "    # mp\n",
    "             \n",
    "    jobs = []\n",
    "    \n",
    "    for sub_mission_lst in mission_lst:\n",
    "        \n",
    "        p = Process(target = write_cluster_s,\n",
    "                    args = (sub_mission_lst,\n",
    "                            total_amount,\n",
    "                            types_to_write,\n",
    "                            output_path,\n",
    "                            cated_fasta,\n",
    "                            id_cluster_dict,\n",
    "                            id_info_dict\n",
    "                           )\n",
    "                   )\n",
    "            \n",
    "        p.start()\n",
    "        jobs.append(p)\n",
    "\n",
    "\n",
    "    for z in jobs:\n",
    "        z.join()\n",
    "        \n",
    "print(time.asctime())"
   ]
  }
 ],
 "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.8.15"
  },
  "toc": {
   "base_numbering": 1,
   "nav_menu": {},
   "number_sections": true,
   "sideBar": true,
   "skip_h1_title": false,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": false,
   "toc_position": {
    "height": "calc(100% - 180px)",
    "left": "10px",
    "top": "150px",
    "width": "181px"
   },
   "toc_section_display": true,
   "toc_window_display": true
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}