Jack_SSNPipeline.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This notebook will get a take a protein family of choice from UniProt, download the sequences with annotations such as domains, PDB codes and taxonomy and then cluster by similarity using mmSEQs. \n",
    "\n",
    "The sequence space will be plotted and the most representative sequences will be taken from across the sequence space \n",
    "\n",
    "Finally the selected sequences will be analysed for solubility and expressebility in e coli \n",
    "\n",
    "Eventually this will be combined with automated analysis of the proteins in the lab to train a reinforcement learning network "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First you will need to go to uniprot, search for the interpro or pfam ID you are interested in. We will download 2 files, a Fasta file and a TSV file\n",
    "for the fasta file just download the default option. for TSV follow below.\n",
    "Select download, format: TSV and select the annotations you are interested in. I generally add taxnomic lineage, PDB, and domain [FT], ensure you also. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Instead of using the API we can just download the TSV and move it to the directory you are running this notebook from. Then we will use mmseqs to clean the data create an all by all comparison, you will need to install mmseqs2 first. https://github.com/soedinglab/mmseqs2/wiki"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "easy-cluster GH18_input.fasta clusterRes tmp --min-seq-id 0.8 \n",
      "\n",
      "MMseqs Version:                     \t15-6f452\n",
      "Substitution matrix                 \taa:blosum62.out,nucl:nucleotide.out\n",
      "Seed substitution matrix            \taa:VTML80.out,nucl:nucleotide.out\n",
      "Sensitivity                         \t4\n",
      "k-mer length                        \t0\n",
      "Target search mode                  \t0\n",
      "k-score                             \tseq:2147483647,prof:2147483647\n",
      "Alphabet size                       \taa:21,nucl:5\n",
      "Max sequence length                 \t65535\n",
      "Max results per query               \t20\n",
      "Split database                      \t0\n",
      "Split mode                          \t2\n",
      "Split memory limit                  \t0\n",
      "Coverage threshold                  \t0.8\n",
      "Coverage mode                       \t0\n",
      "Compositional bias                  \t1\n",
      "Compositional bias                  \t1\n",
      "Diagonal scoring                    \ttrue\n",
      "Exact k-mer matching                \t0\n",
      "Mask residues                       \t1\n",
      "Mask residues probability           \t0.9\n",
      "Mask lower case residues            \t0\n",
      "Minimum diagonal score              \t15\n",
      "Selected taxa                       \t\n",
      "Include identical seq. id.          \tfalse\n",
      "Spaced k-mers                       \t1\n",
      "Preload mode                        \t0\n",
      "Pseudo count a                      \tsubstitution:1.100,context:1.400\n",
      "Pseudo count b                      \tsubstitution:4.100,context:5.800\n",
      "Spaced k-mer pattern                \t\n",
      "Local temporary path                \t\n",
      "Threads                             \t10\n",
      "Compressed                          \t0\n",
      "Verbosity                           \t3\n",
      "Add backtrace                       \tfalse\n",
      "Alignment mode                      \t3\n",
      "Alignment mode                      \t0\n",
      "Allow wrapped scoring               \tfalse\n",
      "E-value threshold                   \t0.001\n",
      "Seq. id. threshold                  \t0.8\n",
      "Min alignment length                \t0\n",
      "Seq. id. mode                       \t0\n",
      "Alternative alignments              \t0\n",
      "Max reject                          \t2147483647\n",
      "Max accept                          \t2147483647\n",
      "Score bias                          \t0\n",
      "Realign hits                        \tfalse\n",
      "Realign score bias                  \t-0.2\n",
      "Realign max seqs                    \t2147483647\n",
      "Correlation score weight            \t0\n",
      "Gap open cost                       \taa:11,nucl:5\n",
      "Gap extension cost                  \taa:1,nucl:2\n",
      "Zdrop                               \t40\n",
      "Rescore mode                        \t0\n",
      "Remove hits by seq. id. and coverage\tfalse\n",
      "Sort results                        \t0\n",
      "Cluster mode                        \t0\n",
      "Max connected component depth       \t1000\n",
      "Similarity type                     \t2\n",
      "Weight file name                    \t\n",
      "Cluster Weight threshold            \t0.9\n",
      "Single step clustering              \tfalse\n",
      "Cascaded clustering steps           \t3\n",
      "Cluster reassign                    \tfalse\n",
      "Remove temporary files              \ttrue\n",
      "Force restart with latest tmp       \tfalse\n",
      "MPI runner                          \t\n",
      "k-mers per sequence                 \t21\n",
      "Scale k-mers per sequence           \taa:0.000,nucl:0.200\n",
      "Adjust k-mer length                 \tfalse\n",
      "Shift hash                          \t67\n",
      "Include only extendable             \tfalse\n",
      "Skip repeating k-mers               \tfalse\n",
      "Database type                       \t0\n",
      "Shuffle input database              \ttrue\n",
      "Createdb mode                       \t1\n",
      "Write lookup file                   \t0\n",
      "Offset of numeric ids               \t0\n",
      "\n",
      "createdb GH18_input.fasta tmp/15846222258004601135/input --dbtype 0 --shuffle 1 --createdb-mode 1 --write-lookup 0 --id-offset 0 --compressed 0 -v 3 \n",
      "\n",
      "\u001b[33mShuffle database cannot be combined with --createdb-mode 0\n",
      "\u001b[39m\u001b[33mWe recompute with --shuffle 0\n",
      "\u001b[39mConverting sequences\n",
      "[41108] 0s 12ms\n",
      "Time for merging to input_h: 0h 0m 0s 0ms\n",
      "Time for merging to input: 0h 0m 0s 0ms\n",
      "Database type: Aminoacid\n",
      "Time for processing: 0h 0m 0s 14ms\n",
      "Create directory tmp/15846222258004601135/clu_tmp\n",
      "cluster tmp/15846222258004601135/input tmp/15846222258004601135/clu tmp/15846222258004601135/clu_tmp --max-seqs 20 -c 0.8 --spaced-kmer-mode 1 --alignment-mode 3 -e 0.001 --min-seq-id 0.8 --remove-tmp-files 1 \n",
      "\n",
      "Set cluster sensitivity to -s 1.000000\n",
      "Set cluster mode SET COVER\n",
      "Set cluster iterations to 1\n",
      "linclust tmp/15846222258004601135/input tmp/15846222258004601135/clu_tmp/6091375675523680812/clu_redundancy tmp/15846222258004601135/clu_tmp/6091375675523680812/linclust --cluster-mode 0 --max-iterations 1000 --similarity-type 2 --threads 10 --compressed 0 -v 3 --cluster-weight-threshold 0.9 --sub-mat 'aa:blosum62.out,nucl:nucleotide.out' -a 0 --alignment-mode 3 --alignment-output-mode 0 --wrapped-scoring 0 -e 0.001 --min-seq-id 0.8 --min-aln-len 0 --seq-id-mode 0 --alt-ali 0 -c 0.8 --cov-mode 0 --max-seq-len 65535 --comp-bias-corr 0 --comp-bias-corr-scale 1 --max-rejected 2147483647 --max-accept 2147483647 --add-self-matches 0 --db-load-mode 0 --pca substitution:1.100,context:1.400 --pcb substitution:4.100,context:5.800 --score-bias 0 --realign 0 --realign-score-bias -0.2 --realign-max-seqs 2147483647 --corr-score-weight 0 --gap-open aa:11,nucl:5 --gap-extend aa:1,nucl:2 --zdrop 40 --alph-size aa:13,nucl:5 --kmer-per-seq 21 --spaced-kmer-mode 1 --kmer-per-seq-scale aa:0.000,nucl:0.200 --adjust-kmer-len 0 --mask 0 --mask-prob 0.9 --mask-lower-case 0 -k 0 --hash-shift 67 --split-memory-limit 0 --include-only-extendable 0 --ignore-multi-kmer 0 --rescore-mode 0 --filter-hits 0 --sort-results 0 --remove-tmp-files 1 --force-reuse 0 \n",
      "\n",
      "kmermatcher tmp/15846222258004601135/input tmp/15846222258004601135/clu_tmp/6091375675523680812/linclust/5499051141078875745/pref --sub-mat 'aa:blosum62.out,nucl:nucleotide.out' --alph-size aa:13,nucl:5 --min-seq-id 0.8 --kmer-per-seq 21 --spaced-kmer-mode 1 --kmer-per-seq-scale aa:0.000,nucl:0.200 --adjust-kmer-len 0 --mask 0 --mask-prob 0.9 --mask-lower-case 0 --cov-mode 0 -k 0 -c 0.8 --max-seq-len 65535 --hash-shift 67 --split-memory-limit 0 --include-only-extendable 0 --ignore-multi-kmer 0 --threads 10 --compressed 0 -v 3 --cluster-weight-threshold 0.9 \n",
      "\n",
      "kmermatcher tmp/15846222258004601135/input tmp/15846222258004601135/clu_tmp/6091375675523680812/linclust/5499051141078875745/pref --sub-mat 'aa:blosum62.out,nucl:nucleotide.out' --alph-size aa:13,nucl:5 --min-seq-id 0.8 --kmer-per-seq 21 --spaced-kmer-mode 1 --kmer-per-seq-scale aa:0.000,nucl:0.200 --adjust-kmer-len 0 --mask 0 --mask-prob 0.9 --mask-lower-case 0 --cov-mode 0 -k 0 -c 0.8 --max-seq-len 65535 --hash-shift 67 --split-memory-limit 0 --include-only-extendable 0 --ignore-multi-kmer 0 --threads 10 --compressed 0 -v 3 --cluster-weight-threshold 0.9 \n",
      "\n",
      "Database size: 41191 type: Aminoacid\n",
      "Reduced amino acid alphabet: (A S T) (C) (D B N) (E Q Z) (F Y) (G) (H) (I V) (K R) (L J M) (P) (W) (X) \n",
      "\n",
      "Generate k-mers list for 1 split\n",
      "[=================================================================] 100.00% 41.19K 0s 80ms     \n",
      "Sort kmer 0h 0m 0s 19ms\n",
      "Sort by rep. sequence 0h 0m 0s 9ms\n",
      "Time for fill: 0h 0m 0s 8ms\n",
      "Time for merging to pref: 0h 0m 0s 0ms\n",
      "Time for processing: 0h 0m 0s 137ms\n",
      "rescorediagonal tmp/15846222258004601135/input tmp/15846222258004601135/input tmp/15846222258004601135/clu_tmp/6091375675523680812/linclust/5499051141078875745/pref tmp/15846222258004601135/clu_tmp/6091375675523680812/linclust/5499051141078875745/pref_rescore1 --sub-mat 'aa:blosum62.out,nucl:nucleotide.out' --rescore-mode 0 --wrapped-scoring 0 --filter-hits 0 -e 0.001 -c 0.8 -a 0 --cov-mode 0 --min-seq-id 0.8 --min-aln-len 0 --seq-id-mode 0 --add-self-matches 0 --sort-results 0 --db-load-mode 0 --threads 10 --compressed 0 -v 3 \n",
      "\n",
      "[=================================================================] 100.00% 41.19K 0s 22ms       ] 39.82% 16.40K eta 0s       \n",
      "Time for merging to pref_rescore1: 0h 0m 0s 11ms=================>] 99.94% 41.17K eta 0s       \n",
      "Time for processing: 0h 0m 0s 44ms\n",
      "clust tmp/15846222258004601135/input tmp/15846222258004601135/clu_tmp/6091375675523680812/linclust/5499051141078875745/pref_rescore1 tmp/15846222258004601135/clu_tmp/6091375675523680812/linclust/5499051141078875745/pre_clust --cluster-mode 0 --max-iterations 1000 --similarity-type 2 --threads 10 --compressed 0 -v 3 --cluster-weight-threshold 0.9 \n",
      "\n",
      "Clustering mode: Set Cover\n",
      "[=================================================================] 100.00% 41.19K 0s 14ms     ==>                               ] 51.51% 21.22K eta 0s       \n",
      "Sort entries\n",
      "Find missing connections\n",
      "Found 70586 new connections.\n",
      "Reconstruct initial order\n",
      "[=================================================================] 100.00% 41.19K 0s 13ms     \n",
      "Add missing connections\n",
      "[=================================================================] 100.00% 41.19K 0s 1ms      \n",
      "\n",
      "Time for read in: 0h 0m 0s 32ms\n",
      "Total time: 0h 0m 0s 35ms\n",
      "\n",
      "Size of the sequence database: 41191\n",
      "Size of the alignment database: 41191\n",
      "Number of clusters: 15115\n",
      "\n",
      "Writing results 0h 0m 0s 1ms\n",
      "Time for merging to pre_clust: 0h 0m 0s 0ms\n",
      "Time for processing: 0h 0m 0s 43ms\n",
      "createsubdb tmp/15846222258004601135/clu_tmp/6091375675523680812/linclust/5499051141078875745/order_redundancy tmp/15846222258004601135/input tmp/15846222258004601135/clu_tmp/6091375675523680812/linclust/5499051141078875745/input_step_redundancy -v 3 --subdb-mode 1 \n",
      "\n",
      "Time for merging to input_step_redundancy: 0h 0m 0s 0ms\n",
      "Time for processing: 0h 0m 0s 5ms\n",
      "createsubdb tmp/15846222258004601135/clu_tmp/6091375675523680812/linclust/5499051141078875745/order_redundancy tmp/15846222258004601135/clu_tmp/6091375675523680812/linclust/5499051141078875745/pref tmp/15846222258004601135/clu_tmp/6091375675523680812/linclust/5499051141078875745/pref_filter1 -v 3 --subdb-mode 1 \n",
      "\n",
      "Time for merging to pref_filter1: 0h 0m 0s 0ms\n",
      "Time for processing: 0h 0m 0s 7ms\n",
      "filterdb tmp/15846222258004601135/clu_tmp/6091375675523680812/linclust/5499051141078875745/pref_filter1 tmp/15846222258004601135/clu_tmp/6091375675523680812/linclust/5499051141078875745/pref_filter2 --filter-file tmp/15846222258004601135/clu_tmp/6091375675523680812/linclust/5499051141078875745/order_redundancy --threads 10 --compressed 0 -v 3 \n",
      "\n",
      "Filtering using file(s)\n",
      "[=================================================================] 100.00% 15.11K 0s 14ms     \n",
      "Time for merging to pref_filter2: 0h 0m 0s 6ms\n",
      "Time for processing: 0h 0m 0s 29ms\n",
      "rescorediagonal tmp/15846222258004601135/clu_tmp/6091375675523680812/linclust/5499051141078875745/input_step_redundancy tmp/15846222258004601135/clu_tmp/6091375675523680812/linclust/5499051141078875745/input_step_redundancy tmp/15846222258004601135/clu_tmp/6091375675523680812/linclust/5499051141078875745/pref_filter2 tmp/15846222258004601135/clu_tmp/6091375675523680812/linclust/5499051141078875745/pref_rescore2 --sub-mat 'aa:blosum62.out,nucl:nucleotide.out' --rescore-mode 1 --wrapped-scoring 0 --filter-hits 1 -e 0.001 -c 0.8 -a 0 --cov-mode 0 --min-seq-id 0.8 --min-aln-len 0 --seq-id-mode 0 --add-self-matches 0 --sort-results 0 --db-load-mode 0 --threads 10 --compressed 0 -v 3 \n",
      "\n",
      "[=================================================================] 100.00% 15.11K 0s 14ms     \n",
      "Time for merging to pref_rescore2: 0h 0m 0s 7ms==================>] 98.94% 14.95K eta 0s       \n",
      "Time for processing: 0h 0m 0s 31ms\n",
      "align tmp/15846222258004601135/clu_tmp/6091375675523680812/linclust/5499051141078875745/input_step_redundancy tmp/15846222258004601135/clu_tmp/6091375675523680812/linclust/5499051141078875745/input_step_redundancy tmp/15846222258004601135/clu_tmp/6091375675523680812/linclust/5499051141078875745/pref_rescore2 tmp/15846222258004601135/clu_tmp/6091375675523680812/linclust/5499051141078875745/aln --sub-mat 'aa:blosum62.out,nucl:nucleotide.out' -a 0 --alignment-mode 3 --alignment-output-mode 0 --wrapped-scoring 0 -e 0.001 --min-seq-id 0.8 --min-aln-len 0 --seq-id-mode 0 --alt-ali 0 -c 0.8 --cov-mode 0 --max-seq-len 65535 --comp-bias-corr 0 --comp-bias-corr-scale 1 --max-rejected 2147483647 --max-accept 2147483647 --add-self-matches 0 --db-load-mode 0 --pca substitution:1.100,context:1.400 --pcb substitution:4.100,context:5.800 --score-bias 0 --realign 0 --realign-score-bias -0.2 --realign-max-seqs 2147483647 --corr-score-weight 0 --gap-open aa:11,nucl:5 --gap-extend aa:1,nucl:2 --zdrop 40 --threads 10 --compressed 0 -v 3 \n",
      "\n",
      "Compute score, coverage and sequence identity\n",
      "Query database size: 15115 type: Aminoacid\n",
      "Target database size: 15115 type: Aminoacid\n",
      "Calculation of alignments\n",
      "[=================================================================] 100.00% 15.11K 1s 976ms    \n",
      "Time for merging to aln: 0h 0m 0s 4ms\n",
      "38164 alignments calculated\n",
      "16938 sequence pairs passed the thresholds (0.443821 of overall calculated)\n",
      "1.120609 hits per query sequence\n",
      "Time for processing: 0h 0m 1s 988ms\n",
      "clust tmp/15846222258004601135/clu_tmp/6091375675523680812/linclust/5499051141078875745/input_step_redundancy tmp/15846222258004601135/clu_tmp/6091375675523680812/linclust/5499051141078875745/aln tmp/15846222258004601135/clu_tmp/6091375675523680812/linclust/5499051141078875745/clust --cluster-mode 0 --max-iterations 1000 --similarity-type 2 --threads 10 --compressed 0 -v 3 --cluster-weight-threshold 0.9 \n",
      "\n",
      "Clustering mode: Set Cover\n",
      "[=================================================================] 100.00% 15.11K 0s 7ms      \n",
      "Sort entries\n",
      "Find missing connections\n",
      "Found 1823 new connections.\n",
      "Reconstruct initial order\n",
      "[=================================================================] 100.00% 15.11K 0s 9ms      \n",
      "Add missing connections\n",
      "[=================================================================] 100.00% 15.11K 0s 1ms      \n",
      "\n",
      "Time for read in: 0h 0m 0s 20ms\n",
      "Total time: 0h 0m 0s 23ms\n",
      "\n",
      "Size of the sequence database: 15115\n",
      "Size of the alignment database: 15115\n",
      "Number of clusters: 13825\n",
      "\n",
      "Writing results 0h 0m 0s 2ms\n",
      "Time for merging to clust: 0h 0m 0s 0ms\n",
      "Time for processing: 0h 0m 0s 28ms\n",
      "mergeclusters tmp/15846222258004601135/input tmp/15846222258004601135/clu_tmp/6091375675523680812/clu_redundancy tmp/15846222258004601135/clu_tmp/6091375675523680812/linclust/5499051141078875745/pre_clust tmp/15846222258004601135/clu_tmp/6091375675523680812/linclust/5499051141078875745/clust --threads 10 --compressed 0 -v 3 \n",
      "\n",
      "Clustering step 1\n",
      "[=================================================================] 100.00% 15.11K 0s 8ms      \n",
      "Clustering step 2\n",
      "[=================================================================] 100.00% 13.82K 0s 14ms     \n",
      "Write merged clustering\n",
      "[=================================================================] 100.00% 41.19K 0s 24ms     \n",
      "Time for merging to clu_redundancy: 0h 0m 0s 3ms\n",
      "Time for processing: 0h 0m 0s 32ms\n",
      "rmdb tmp/15846222258004601135/clu_tmp/6091375675523680812/linclust/5499051141078875745/pref_filter1 -v 3 \n",
      "\n",
      "Time for processing: 0h 0m 0s 0ms\n",
      "rmdb tmp/15846222258004601135/clu_tmp/6091375675523680812/linclust/5499051141078875745/pref -v 3 \n",
      "\n",
      "Time for processing: 0h 0m 0s 2ms\n",
      "rmdb tmp/15846222258004601135/clu_tmp/6091375675523680812/linclust/5499051141078875745/pref_rescore1 -v 3 \n",
      "\n",
      "Time for processing: 0h 0m 0s 0ms\n",
      "rmdb tmp/15846222258004601135/clu_tmp/6091375675523680812/linclust/5499051141078875745/pre_clust -v 3 \n",
      "\n",
      "Time for processing: 0h 0m 0s 0ms\n",
      "rmdb tmp/15846222258004601135/clu_tmp/6091375675523680812/linclust/5499051141078875745/input_step_redundancy -v 3 \n",
      "\n",
      "Time for processing: 0h 0m 0s 0ms\n",
      "rmdb tmp/15846222258004601135/clu_tmp/6091375675523680812/linclust/5499051141078875745/input_step_redundancy_h -v 3 \n",
      "\n",
      "Time for processing: 0h 0m 0s 0ms\n",
      "rmdb tmp/15846222258004601135/clu_tmp/6091375675523680812/linclust/5499051141078875745/pref_filter2 -v 3 \n",
      "\n",
      "Time for processing: 0h 0m 0s 0ms\n",
      "rmdb tmp/15846222258004601135/clu_tmp/6091375675523680812/linclust/5499051141078875745/pref_rescore2 -v 3 \n",
      "\n",
      "Time for processing: 0h 0m 0s 0ms\n",
      "rmdb tmp/15846222258004601135/clu_tmp/6091375675523680812/linclust/5499051141078875745/aln -v 3 \n",
      "\n",
      "Time for processing: 0h 0m 0s 0ms\n",
      "rmdb tmp/15846222258004601135/clu_tmp/6091375675523680812/linclust/5499051141078875745/clust -v 3 \n",
      "\n",
      "Time for processing: 0h 0m 0s 0ms\n",
      "createsubdb tmp/15846222258004601135/clu_tmp/6091375675523680812/clu_redundancy tmp/15846222258004601135/input tmp/15846222258004601135/clu_tmp/6091375675523680812/input_step_redundancy -v 3 --subdb-mode 1 \n",
      "\n",
      "Time for merging to input_step_redundancy: 0h 0m 0s 0ms\n",
      "Time for processing: 0h 0m 0s 5ms\n",
      "prefilter tmp/15846222258004601135/clu_tmp/6091375675523680812/input_step_redundancy tmp/15846222258004601135/clu_tmp/6091375675523680812/input_step_redundancy tmp/15846222258004601135/clu_tmp/6091375675523680812/pref_step0 --sub-mat 'aa:blosum62.out,nucl:nucleotide.out' --seed-sub-mat 'aa:VTML80.out,nucl:nucleotide.out' -s 1 -k 0 --target-search-mode 0 --k-score seq:2147483647,prof:2147483647 --alph-size aa:21,nucl:5 --max-seq-len 65535 --max-seqs 20 --split 0 --split-mode 2 --split-memory-limit 0 -c 0.8 --cov-mode 0 --comp-bias-corr 0 --comp-bias-corr-scale 1 --diag-score 0 --exact-kmer-matching 0 --mask 1 --mask-prob 0.9 --mask-lower-case 0 --min-ungapped-score 0 --add-self-matches 0 --spaced-kmer-mode 1 --db-load-mode 0 --pca substitution:1.100,context:1.400 --pcb substitution:4.100,context:5.800 --threads 10 --compressed 0 -v 3 \n",
      "\n",
      "Query database size: 13825 type: Aminoacid\n",
      "Estimated memory consumption: 1012M\n",
      "Target database size: 13825 type: Aminoacid\n",
      "Index table k-mer threshold: 154 at k-mer size 6 \n",
      "Index table: counting k-mers\n",
      "[=================================================================] 100.00% 13.82K 0s 131ms    \n",
      "Index table: Masked residues: 1922\n",
      "Index table: fill\n",
      "[=================================================================] 100.00% 13.82K 0s 56ms                                                       ] 5.75% 796 eta 0s       \n",
      "Index statistics\n",
      "Entries:          2085886\n",
      "DB size:          500 MB\n",
      "Avg k-mer size:   0.032592\n",
      "Top 10 k-mers\n",
      "    DGDDYP\t2253\n",
      "    LKLSGW\t847\n",
      "    FDIIEY\t564\n",
      "    EKVGFY\t436\n",
      "    MGPYDW\t409\n",
      "    THNANV\t393\n",
      "    LGGYGW\t391\n",
      "    GIFGWT\t341\n",
      "    DWYNGL\t322\n",
      "    LGPYGW\t318\n",
      "Time for index table init: 0h 0m 0s 434ms\n",
      "Process prefiltering step 1 of 1\n",
      "\n",
      "k-mer similarity threshold: 154\n",
      "Starting prefiltering scores calculation (step 1 of 1)\n",
      "Query db start 1 to 13825\n",
      "Target db start 1 to 13825\n",
      "[=================================================================] 100.00% 13.82K 0s 250ms    \n",
      "\n",
      "2.317447 k-mers per position\n",
      "4585 DB matches per sequence\n",
      "0 overflows\n",
      "19 sequences passed prefiltering per query sequence\n",
      "20 median result list length\n",
      "0 sequences with 0 size result lists\n",
      "Time for merging to pref_step0: 0h 0m 0s 4ms\n",
      "Time for processing: 0h 0m 1s 168ms\n",
      "align tmp/15846222258004601135/clu_tmp/6091375675523680812/input_step_redundancy tmp/15846222258004601135/clu_tmp/6091375675523680812/input_step_redundancy tmp/15846222258004601135/clu_tmp/6091375675523680812/pref_step0 tmp/15846222258004601135/clu_tmp/6091375675523680812/aln_step0 --sub-mat 'aa:blosum62.out,nucl:nucleotide.out' -a 0 --alignment-mode 3 --alignment-output-mode 0 --wrapped-scoring 0 -e 0.001 --min-seq-id 0.8 --min-aln-len 0 --seq-id-mode 0 --alt-ali 0 -c 0.8 --cov-mode 0 --max-seq-len 65535 --comp-bias-corr 0 --comp-bias-corr-scale 1 --max-rejected 2147483647 --max-accept 2147483647 --add-self-matches 0 --db-load-mode 0 --pca substitution:1.100,context:1.400 --pcb substitution:4.100,context:5.800 --score-bias 0 --realign 0 --realign-score-bias -0.2 --realign-max-seqs 2147483647 --corr-score-weight 0 --gap-open aa:11,nucl:5 --gap-extend aa:1,nucl:2 --zdrop 40 --threads 10 --compressed 0 -v 3 \n",
      "\n",
      "Compute score, coverage and sequence identity\n",
      "Query database size: 13825 type: Aminoacid\n",
      "Target database size: 13825 type: Aminoacid\n",
      "Calculation of alignments\n",
      "[=================================================================] 100.00% 13.82K 14s 580ms   \n",
      "Time for merging to aln_step0: 0h 0m 0s 3ms\n",
      "219210 alignments calculated\n",
      "22759 sequence pairs passed the thresholds (0.103823 of overall calculated)\n",
      "1.646221 hits per query sequence\n",
      "Time for processing: 0h 0m 14s 595ms\n",
      "clust tmp/15846222258004601135/clu_tmp/6091375675523680812/input_step_redundancy tmp/15846222258004601135/clu_tmp/6091375675523680812/aln_step0 tmp/15846222258004601135/clu_tmp/6091375675523680812/clu_step0 --cluster-mode 0 --max-iterations 1000 --similarity-type 2 --threads 10 --compressed 0 -v 3 --cluster-weight-threshold 0.9 \n",
      "\n",
      "Clustering mode: Set Cover\n",
      "[=================================================================] 100.00% 13.82K 0s 8ms      \n",
      "Sort entries\n",
      "Find missing connections\n",
      "Found 1756 new connections.\n",
      "Reconstruct initial order\n",
      "[=================================================================] 100.00% 13.82K 0s 9ms      \n",
      "Add missing connections\n",
      "[=================================================================] 100.00% 13.82K 0s 1ms      \n",
      "\n",
      "Time for read in: 0h 0m 0s 20ms\n",
      "Total time: 0h 0m 0s 22ms\n",
      "\n",
      "Size of the sequence database: 13825\n",
      "Size of the alignment database: 13825\n",
      "Number of clusters: 12516\n",
      "\n",
      "Writing results 0h 0m 0s 1ms\n",
      "Time for merging to clu_step0: 0h 0m 0s 0ms\n",
      "Time for processing: 0h 0m 0s 27ms\n",
      "mergeclusters tmp/15846222258004601135/input tmp/15846222258004601135/clu tmp/15846222258004601135/clu_tmp/6091375675523680812/clu_redundancy tmp/15846222258004601135/clu_tmp/6091375675523680812/clu_step0 --threads 10 --compressed 0 -v 3 \n",
      "\n",
      "Clustering step 1\n",
      "[=================================================================] 100.00% 13.82K 0s 15ms     ============>                    ] 68.32% 9.44K eta 0s       \n",
      "Clustering step 2\n",
      "[=================================================================] 100.00% 12.52K 0s 24ms     \n",
      "Write merged clustering\n",
      "[=================================================================] 100.00% 41.19K 0s 36ms     \n",
      "Time for merging to clu: 0h 0m 0s 9ms\n",
      "Time for processing: 0h 0m 0s 60ms\n",
      "rmdb tmp/15846222258004601135/clu_tmp/6091375675523680812/clu_redundancy -v 3 \n",
      "\n",
      "Time for processing: 0h 0m 0s 0ms\n",
      "rmdb tmp/15846222258004601135/clu_tmp/6091375675523680812/input_step_redundancy -v 3 \n",
      "\n",
      "Time for processing: 0h 0m 0s 0ms\n",
      "rmdb tmp/15846222258004601135/clu_tmp/6091375675523680812/input_step_redundancy_h -v 3 \n",
      "\n",
      "Time for processing: 0h 0m 0s 0ms\n",
      "rmdb tmp/15846222258004601135/clu_tmp/6091375675523680812/pref_step0 -v 3 \n",
      "\n",
      "Time for processing: 0h 0m 0s 1ms\n",
      "rmdb tmp/15846222258004601135/clu_tmp/6091375675523680812/aln_step0 -v 3 \n",
      "\n",
      "Time for processing: 0h 0m 0s 1ms\n",
      "rmdb tmp/15846222258004601135/clu_tmp/6091375675523680812/clu_step0 -v 3 \n",
      "\n",
      "Time for processing: 0h 0m 0s 0ms\n",
      "createtsv tmp/15846222258004601135/input tmp/15846222258004601135/input tmp/15846222258004601135/clu tmp/15846222258004601135/cluster.tsv --threads 10 -v 3 \n",
      "\n",
      "Time for merging to cluster.tsv: 0h 0m 0s 5ms\n",
      "Time for processing: 0h 0m 0s 22ms\n",
      "result2repseq tmp/15846222258004601135/input tmp/15846222258004601135/clu tmp/15846222258004601135/clu_rep --db-load-mode 0 --compressed 0 --threads 10 -v 3 \n",
      "\n",
      "[=================================================================] 100.00% 12.52K 0s 7ms      \n",
      "Time for merging to clu_rep: 0h 0m 0s 23ms\n",
      "Time for processing: 0h 0m 0s 41ms\n",
      "result2flat tmp/15846222258004601135/input tmp/15846222258004601135/input tmp/15846222258004601135/clu_rep tmp/15846222258004601135/rep_seq.fasta --use-fasta-header -v 3 \n",
      "\n",
      "Time for processing: 0h 0m 0s 23ms\n",
      "createseqfiledb tmp/15846222258004601135/input tmp/15846222258004601135/clu tmp/15846222258004601135/clu_seqs --threads 10 -v 3 \n",
      "\n",
      "[=================================================================] 100.00% 12.52K 0s 11ms     \n",
      "Time for merging to clu_seqs: 0h 0m 0s 3ms\n",
      "Time for processing: 0h 0m 0s 37ms\n",
      "result2flat tmp/15846222258004601135/input tmp/15846222258004601135/input tmp/15846222258004601135/clu_seqs tmp/15846222258004601135/all_seqs.fasta -v 3 \n",
      "\n",
      "Time for processing: 0h 0m 0s 103ms\n",
      "rmdb tmp/15846222258004601135/input -v 3 \n",
      "\n",
      "Time for processing: 0h 0m 0s 0ms\n",
      "rmdb tmp/15846222258004601135/input_h -v 3 \n",
      "\n",
      "Time for processing: 0h 0m 0s 5ms\n",
      "rmdb tmp/15846222258004601135/clu_seqs -v 3 \n",
      "\n",
      "Time for processing: 0h 0m 0s 1ms\n",
      "rmdb tmp/15846222258004601135/clu_rep -v 3 \n",
      "\n",
      "Time for processing: 0h 0m 0s 0ms\n",
      "rmdb tmp/15846222258004601135/clu -v 3 \n",
      "\n",
      "Time for processing: 0h 0m 0s 1ms\n"
     ]
    }
   ],
   "source": [
    "# remove sequences with greater than 0.8 sequence identity\n",
    "!mmseqs easy-cluster GH18_input.fasta clusterRes tmp --min-seq-id 0.8 #clusters sequences with 80% similarity"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "After cleaning the dataset we need to do an all by all comparison to cluster the sequences by similarity "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\u001b[33mall_by_all.tsv exists and will be overwritten\n",
      "\u001b[39measy-search clusterRes_rep_seq.fasta clusterRes_rep_seq.fasta all_by_all.tsv /tmp --threads 8 \n",
      "\n",
      "MMseqs Version:                        \t15-6f452\n",
      "Substitution matrix                    \taa:blosum62.out,nucl:nucleotide.out\n",
      "Add backtrace                          \tfalse\n",
      "Alignment mode                         \t3\n",
      "Alignment mode                         \t0\n",
      "Allow wrapped scoring                  \tfalse\n",
      "E-value threshold                      \t0.001\n",
      "Seq. id. threshold                     \t0\n",
      "Min alignment length                   \t0\n",
      "Seq. id. mode                          \t0\n",
      "Alternative alignments                 \t0\n",
      "Coverage threshold                     \t0\n",
      "Coverage mode                          \t0\n",
      "Max sequence length                    \t65535\n",
      "Compositional bias                     \t1\n",
      "Compositional bias                     \t1\n",
      "Max reject                             \t2147483647\n",
      "Max accept                             \t2147483647\n",
      "Include identical seq. id.             \tfalse\n",
      "Preload mode                           \t0\n",
      "Pseudo count a                         \tsubstitution:1.100,context:1.400\n",
      "Pseudo count b                         \tsubstitution:4.100,context:5.800\n",
      "Score bias                             \t0\n",
      "Realign hits                           \tfalse\n",
      "Realign score bias                     \t-0.2\n",
      "Realign max seqs                       \t2147483647\n",
      "Correlation score weight               \t0\n",
      "Gap open cost                          \taa:11,nucl:5\n",
      "Gap extension cost                     \taa:1,nucl:2\n",
      "Zdrop                                  \t40\n",
      "Threads                                \t8\n",
      "Compressed                             \t0\n",
      "Verbosity                              \t3\n",
      "Seed substitution matrix               \taa:VTML80.out,nucl:nucleotide.out\n",
      "Sensitivity                            \t5.7\n",
      "k-mer length                           \t0\n",
      "Target search mode                     \t0\n",
      "k-score                                \tseq:2147483647,prof:2147483647\n",
      "Alphabet size                          \taa:21,nucl:5\n",
      "Max results per query                  \t300\n",
      "Split database                         \t0\n",
      "Split mode                             \t2\n",
      "Split memory limit                     \t0\n",
      "Diagonal scoring                       \ttrue\n",
      "Exact k-mer matching                   \t0\n",
      "Mask residues                          \t1\n",
      "Mask residues probability              \t0.9\n",
      "Mask lower case residues               \t0\n",
      "Minimum diagonal score                 \t15\n",
      "Selected taxa                          \t\n",
      "Spaced k-mers                          \t1\n",
      "Spaced k-mer pattern                   \t\n",
      "Local temporary path                   \t\n",
      "Rescore mode                           \t0\n",
      "Remove hits by seq. id. and coverage   \tfalse\n",
      "Sort results                           \t0\n",
      "Mask profile                           \t1\n",
      "Profile E-value threshold              \t0.001\n",
      "Global sequence weighting              \tfalse\n",
      "Allow deletions                        \tfalse\n",
      "Filter MSA                             \t1\n",
      "Use filter only at N seqs              \t0\n",
      "Maximum seq. id. threshold             \t0.9\n",
      "Minimum seq. id.                       \t0.0\n",
      "Minimum score per column               \t-20\n",
      "Minimum coverage                       \t0\n",
      "Select N most diverse seqs             \t1000\n",
      "Pseudo count mode                      \t0\n",
      "Min codons in orf                      \t30\n",
      "Max codons in length                   \t32734\n",
      "Max orf gaps                           \t2147483647\n",
      "Contig start mode                      \t2\n",
      "Contig end mode                        \t2\n",
      "Orf start mode                         \t1\n",
      "Forward frames                         \t1,2,3\n",
      "Reverse frames                         \t1,2,3\n",
      "Translation table                      \t1\n",
      "Translate orf                          \t0\n",
      "Use all table starts                   \tfalse\n",
      "Offset of numeric ids                  \t0\n",
      "Create lookup                          \t0\n",
      "Add orf stop                           \tfalse\n",
      "Overlap between sequences              \t0\n",
      "Sequence split mode                    \t1\n",
      "Header split mode                      \t0\n",
      "Chain overlapping alignments           \t0\n",
      "Merge query                            \t1\n",
      "Search type                            \t0\n",
      "Search iterations                      \t1\n",
      "Start sensitivity                      \t4\n",
      "Search steps                           \t1\n",
      "Prefilter mode                         \t0\n",
      "Exhaustive search mode                 \tfalse\n",
      "Filter results during exhaustive search\t0\n",
      "Strand selection                       \t1\n",
      "LCA search mode                        \tfalse\n",
      "Disk space limit                       \t0\n",
      "MPI runner                             \t\n",
      "Force restart with latest tmp          \tfalse\n",
      "Remove temporary files                 \ttrue\n",
      "Alignment format                       \t0\n",
      "Format alignment output                \tquery,target,fident,alnlen,mismatch,gapopen,qstart,qend,tstart,tend,evalue,bits\n",
      "Database output                        \tfalse\n",
      "Overlap threshold                      \t0\n",
      "Database type                          \t0\n",
      "Shuffle input database                 \ttrue\n",
      "Createdb mode                          \t0\n",
      "Write lookup file                      \t0\n",
      "Greedy best hits                       \tfalse\n",
      "\n",
      "createdb clusterRes_rep_seq.fasta /tmp/2199896455259100104/query --dbtype 0 --shuffle 1 --createdb-mode 0 --write-lookup 0 --id-offset 0 --compressed 0 -v 3 \n",
      "\n",
      "Converting sequences\n",
      "[12424] 0s 16ms\n",
      "Time for merging to query_h: 0h 0m 0s 6ms\n",
      "Time for merging to query: 0h 0m 0s 6ms\n",
      "Database type: Aminoacid\n",
      "Time for processing: 0h 0m 0s 46ms\n",
      "createdb clusterRes_rep_seq.fasta /tmp/2199896455259100104/target --dbtype 0 --shuffle 1 --createdb-mode 0 --write-lookup 0 --id-offset 0 --compressed 0 -v 3 \n",
      "\n",
      "Converting sequences\n",
      "[12424] 0s 14ms\n",
      "Time for merging to target_h: 0h 0m 0s 6ms\n",
      "Time for merging to target: 0h 0m 0s 6ms\n",
      "Database type: Aminoacid\n",
      "Time for processing: 0h 0m 0s 47ms\n",
      "Create directory /tmp/2199896455259100104/search_tmp\n",
      "search /tmp/2199896455259100104/query /tmp/2199896455259100104/target /tmp/2199896455259100104/result /tmp/2199896455259100104/search_tmp --alignment-mode 3 --threads 8 -s 5.7 --remove-tmp-files 1 \n",
      "\n",
      "prefilter /tmp/2199896455259100104/query /tmp/2199896455259100104/target /tmp/2199896455259100104/search_tmp/1160442705650774285/pref_0 --sub-mat 'aa:blosum62.out,nucl:nucleotide.out' --seed-sub-mat 'aa:VTML80.out,nucl:nucleotide.out' -k 0 --target-search-mode 0 --k-score seq:2147483647,prof:2147483647 --alph-size aa:21,nucl:5 --max-seq-len 65535 --max-seqs 300 --split 0 --split-mode 2 --split-memory-limit 0 -c 0 --cov-mode 0 --comp-bias-corr 1 --comp-bias-corr-scale 1 --diag-score 1 --exact-kmer-matching 0 --mask 1 --mask-prob 0.9 --mask-lower-case 0 --min-ungapped-score 15 --add-self-matches 0 --spaced-kmer-mode 1 --db-load-mode 0 --pca substitution:1.100,context:1.400 --pcb substitution:4.100,context:5.800 --threads 8 --compressed 0 -v 3 -s 5.7 \n",
      "\n",
      "Query database size: 12516 type: Aminoacid\n",
      "Estimated memory consumption: 1008M\n",
      "Target database size: 12516 type: Aminoacid\n",
      "Index table k-mer threshold: 112 at k-mer size 6 \n",
      "Index table: counting k-mers\n",
      "[=================================================================] 100.00% 12.52K 0s 164ms    \n",
      "Index table: Masked residues: 1772\n",
      "Index table: fill\n",
      "[=================================================================] 100.00% 12.52K 0s 83ms     \n",
      "Index statistics\n",
      "Entries:          3627912\n",
      "DB size:          509 MB\n",
      "Avg k-mer size:   0.056686\n",
      "Top 10 k-mers\n",
      "    DGDDYP\t1858\n",
      "    LKLSGW\t647\n",
      "    HEWERS\t503\n",
      "    EKVGFY\t414\n",
      "    MGPYDW\t388\n",
      "    FDIIEY\t322\n",
      "    LGPYGW\t281\n",
      "    GTEGDY\t279\n",
      "    DGNDNV\t271\n",
      "    GGSSFS\t247\n",
      "Time for index table init: 0h 0m 0s 495ms\n",
      "Process prefiltering step 1 of 1\n",
      "\n",
      "k-mer similarity threshold: 112\n",
      "Starting prefiltering scores calculation (step 1 of 1)\n",
      "Query db start 1 to 12516\n",
      "Target db start 1 to 12516\n",
      "[=================================================================] 100.00% 12.52K 9s 840ms    \n",
      "\n",
      "397.895709 k-mers per position\n",
      "44372 DB matches per sequence\n",
      "0 overflows\n",
      "299 sequences passed prefiltering per query sequence\n",
      "300 median result list length\n",
      "0 sequences with 0 size result lists\n",
      "Time for merging to pref_0: 0h 0m 0s 4ms\n",
      "Time for processing: 0h 0m 10s 789ms\n",
      "align /tmp/2199896455259100104/query /tmp/2199896455259100104/target /tmp/2199896455259100104/search_tmp/1160442705650774285/pref_0 /tmp/2199896455259100104/result --sub-mat 'aa:blosum62.out,nucl:nucleotide.out' -a 0 --alignment-mode 3 --alignment-output-mode 0 --wrapped-scoring 0 -e 0.001 --min-seq-id 0 --min-aln-len 0 --seq-id-mode 0 --alt-ali 0 -c 0 --cov-mode 0 --max-seq-len 65535 --comp-bias-corr 1 --comp-bias-corr-scale 1 --max-rejected 2147483647 --max-accept 2147483647 --add-self-matches 0 --db-load-mode 0 --pca substitution:1.100,context:1.400 --pcb substitution:4.100,context:5.800 --score-bias 0 --realign 0 --realign-score-bias -0.2 --realign-max-seqs 2147483647 --corr-score-weight 0 --gap-open aa:11,nucl:5 --gap-extend aa:1,nucl:2 --zdrop 40 --threads 8 --compressed 0 -v 3 \n",
      "\n",
      "Compute score, coverage and sequence identity\n",
      "Query database size: 12516 type: Aminoacid\n",
      "Target database size: 12516 type: Aminoacid\n",
      "Calculation of alignments\n",
      "[=================================================================] 100.00% 12.52K 3m 14s 127ms   \n",
      "Time for merging to result: 0h 0m 0s 6ms\n",
      "3742986 alignments calculated\n",
      "3579753 sequence pairs passed the thresholds (0.956390 of overall calculated)\n",
      "286.014130 hits per query sequence\n",
      "Time for processing: 0h 3m 14s 488ms\n",
      "rmdb /tmp/2199896455259100104/search_tmp/1160442705650774285/pref_0 -v 3 \n",
      "\n",
      "Time for processing: 0h 0m 0s 12ms\n",
      "rmdb /tmp/2199896455259100104/search_tmp/1160442705650774285/aln_0 -v 3 \n",
      "\n",
      "Time for processing: 0h 0m 0s 0ms\n",
      "rmdb /tmp/2199896455259100104/search_tmp/1160442705650774285/input_0 -v 3 \n",
      "\n",
      "Time for processing: 0h 0m 0s 0ms\n",
      "rmdb /tmp/2199896455259100104/search_tmp/1160442705650774285/aln_merge -v 3 \n",
      "\n",
      "Time for processing: 0h 0m 0s 0ms\n",
      "\u001b[33mall_by_all.tsv exists and will be overwritten\n",
      "\u001b[39mconvertalis /tmp/2199896455259100104/query /tmp/2199896455259100104/target /tmp/2199896455259100104/result all_by_all.tsv --sub-mat 'aa:blosum62.out,nucl:nucleotide.out' --format-mode 0 --format-output query,target,fident,alnlen,mismatch,gapopen,qstart,qend,tstart,tend,evalue,bits --translation-table 1 --gap-open aa:11,nucl:5 --gap-extend aa:1,nucl:2 --db-output 0 --db-load-mode 0 --search-type 0 --threads 8 --compressed 0 -v 3 \n",
      "\n",
      "[=================================================================] 100.00% 12.52K 2s 577ms    \n",
      "Time for merging to all_by_all.tsv: 0h 0m 0s 294ms\n",
      "Time for processing: 0h 0m 2s 903ms\n",
      "rmdb /tmp/2199896455259100104/result -v 3 \n",
      "\n",
      "Time for processing: 0h 0m 0s 3ms\n",
      "rmdb /tmp/2199896455259100104/target -v 3 \n",
      "\n",
      "Time for processing: 0h 0m 0s 1ms\n",
      "rmdb /tmp/2199896455259100104/target_h -v 3 \n",
      "\n",
      "Time for processing: 0h 0m 0s 0ms\n",
      "rmdb /tmp/2199896455259100104/query -v 3 \n",
      "\n",
      "Time for processing: 0h 0m 0s 1ms\n",
      "rmdb /tmp/2199896455259100104/query_h -v 3 \n",
      "\n",
      "Time for processing: 0h 0m 0s 1ms\n"
     ]
    }
   ],
   "source": [
    "# perform all by all search on cleaned sequences\n",
    "!mmseqs easy-search clusterRes_rep_seq.fasta clusterRes_rep_seq.fasta all_by_all.tsv /tmp --threads 8 "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we will use John's Cluster tools package to analyse the results of the all by all search"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "#first import the modules we need\n",
    "\n",
    "from ClusterTools.cluster_stream import cluster_stream\n",
    "from ClusterTools.cluster_progression import cluster_progression\n",
    "#from ClusterTools.circle_plot import PlotCirclesBox2D\n",
    "from ClusterTools.cluster_size_table import MakeTable\n",
    "from ClusterTools.cluster_consensus import ParseClusters\n",
    "import sys\n",
    "sys.path.append(\"./ClusterTools/\")  # replace with the actual path\n",
    "import pandas as pd\n",
    "import re\n",
    "import os\n",
    "import zipfile\n",
    "from pathlib import Path\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The code below analyses how the cluster sizes vary by bitscore cutoff, other scoring metrics can be used by changing the score setting to a different column in the all_by_all.tsv file"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "---cutoff: 100\n",
      "1 clusters with 12167 sequences\n",
      "349 single nodes added\n",
      "0 unconnected nodes added\n",
      "---cutoff: 200\n",
      "1 clusters with 11004 sequences\n",
      "1512 single nodes added\n",
      "0 unconnected nodes added\n",
      "---cutoff: 300\n",
      "6 clusters with 9797 sequences\n",
      "2719 single nodes added\n",
      "0 unconnected nodes added\n",
      "---cutoff: 400\n",
      "8 clusters with 3071 sequences\n",
      "9445 single nodes added\n",
      "0 unconnected nodes added\n",
      "---cutoff: 500\n",
      "3 clusters with 602 sequences\n",
      "11914 single nodes added\n",
      "0 unconnected nodes added\n",
      "---cutoff: 600\n",
      "1 clusters with 2 sequences\n",
      "12514 single nodes added\n",
      "0 unconnected nodes added\n",
      "---cutoff: 700\n",
      "1 clusters with 2 sequences\n",
      "12514 single nodes added\n",
      "0 unconnected nodes added\n",
      "3579753 \n",
      "27.50s\n",
      "\n",
      "['cluster_400.csv', 'cluster_600.csv', 'cluster_200.csv', 'cluster_700.csv', 'cluster_500.csv', 'cluster_100.csv', 'cluster_300.csv']\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjcAAAGwCAYAAABVdURTAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8pXeV/AAAACXBIWXMAAA9hAAAPYQGoP6dpAABXHElEQVR4nO3deXxM5/4H8M9km+yTRVYJCWLLZkkRRak1gipVXVXXn5aqWoveKl1i6aLubakWrev26iWonShCrUU2kYQsJCKLJJLJOlnm/P4IU2mCDDM5s3zer9e82jnnzMxnTkm/Od/neY5EEAQBRERERAbCROwARERERJrE4oaIiIgMCosbIiIiMigsboiIiMigsLghIiIig8LihoiIiAwKixsiIiIyKEZX3AiCALlcDi7vQ0REZJiMrrgpLS2FTCZDaWmp2FGIiIhIC4yuuCEiIiLDxuKGiIiIDAqLGyIiIjIoLG6IiIjIoLC4ISIiIoPC4oaIiIgMCosbIiIiMigsboiIiMigsLghIiIig8LihoiIiAwKixsiIiIyKCxuiIiIyKCIWtysXr0aQUFBsLe3h729PUJDQ7Fv3777viY6Oho9e/aEpaUl2rVrhzVr1rRQWiIiItIHohY3Xl5eWLp0Kc6dO4dz587hySefxFNPPYXExMQmj8/IyMDIkSPRv39/xMTEYMGCBZg+fToiIyNbODkRERHpKokgCILYIe7m5OSEFStW4PXXX2+0b968edi5cyeSkpJU26ZMmYK4uDicOnWqWe8vl8shk8lQUlICe3t7jeUmIt1QWV0HKwtTsWMQkYh0ZsxNXV0dNm/ejPLycoSGhjZ5zKlTpzBs2LAG24YPH45z586hpqamydcoFArI5fIGDyIyTMv3J6PLR/txNCVf7ChEJCLRi5uEhATY2tpCKpViypQp2L59O7p27drksbm5uXBzc2uwzc3NDbW1tSgoKGjyNREREZDJZKqHt7e3xr8DEYnvdHohvjuaBgDYdPqayGmISEyiFzedOnVCbGwsTp8+jbfffhuvvPIKLl26dM/jJRJJg+d3ump/337H/PnzUVJSonpkZWVpLjwR6YRyRS3mbI1TPT92uQDyqqav5hKR4TMTO4CFhQU6dOgAAAgJCcGff/6Jb775Bt9//32jY93d3ZGbm9tgW35+PszMzODs7Nzk+0ulUkilUs0HJyKdsXRfMrKKKtHawQoWZibIKCjHoUt5GNfDS+xoRCQC0a/c/J0gCFAoFE3uCw0NRVRUVINtBw8eREhICMzNzVsiHhHpmJOpBfj37TbUsvFBeKqbJwBgT3yOmLGISESiFjcLFizA8ePHcfXqVSQkJGDhwoU4evQoXnzxRQD1LaVJkyapjp8yZQquXbuGmTNnIikpCevXr8e6deswe/Zssb4CEYmoTFGLOVvjAQAv9WmDfn6tEB7oAQA4duUmSirZmiIyRqIWN3l5eXj55ZfRqVMnDB48GGfOnMH+/fsxdOhQAEBOTg4yMzNVx/v6+mLv3r04evQounXrhk8++QSrVq3C+PHjxfoKRCSiz/cmIbu4El6OVpgf1gUA4Odmh05udqipExB1KU/khEQkBp1b50bbuM4NkWE4dvkmJq0/CwD475t9ENr+r3F33xy6gq8PXcagTi7Y8GovsSISkUh0bswNEdGDyKtq8EFkfTvqldC2DQobAAgPcgcAHL9SgJIKtqaIjA2LGyLSO5/vScKNkiq0cbLGvLDOjfZ3cLVDZ3c71CoFHLiU28Q7EJEhY3FDRHrlaEo+Nv+ZBYkE+GJCMKwtml7R4s7AYs6aIjI+LG6ISG+UVNbgg8gEAMCrfX3Ry9fpnseODKovbk6kFqC4orpF8hGRbmBxQ0R645Pdl5Arr4JvKxvMGd7pvse2d7FFFw971CoFHEzkrCkiY8Lihoj0wuHkPGw9fx0SCbDimaBm3fl71O2rN7sT2JoiMiYsbohI55VU/NWOeqOfL0J87t2OutvIwL9aU7fK2ZoiMhYsbohI5y3elYj8UgXaudhg1rD7t6Pu5tvKBv6e9qhTCjiQyFlTRMaCxQ0R6bSDibnYFpMNk9uzoyzNH9yOutudqzd72JoiMhosbohIZ90qr8aC7RcBAG8OaIcebRzVfo87U8JPphWisKzpm/ISkWFhcUNEOuvjXYkoKFOgg6st3h/S8aHew6eVDQJa32lNcdYUkTFgcUNEOmn/xRz8FnsDpiYSfPkQ7ai7hQd6AgD2JNzQVDwi0mEsbohI5xSWKbDwdjtqyhPtEOzt8Ejvd6c1dSqtEAVsTREZPBY3RKRzPtqZiMLyanRys8P0wX6P/H5tnK0R5CWDUgBnTREZARY3RKRT9sTnYE98DkxNJPhiQjCkZg/fjrob7zVFZDxY3BCRzigoU+Afv9W3o6YObI9AL5nG3vvOlPDT6YW4WcrWFJEhY3FDRDpBEAT8Y8dFFJVXo7O7HaY9+ejtqLt5O1kj2NsBSgHYz9YUkUFjcUNEOmFXfA72XcyFmYkEXz4bDAszzf94Cg90BwDsieesKSJDxuKGiESXX1qFj263o6Y92QH+npprR93tTmvqTEYR8kurtPIZRCQ+FjdEJCpBELBw+0UUV9Sgq4c9pg7qoLXP8nK0RjdvBwgCsP8iW1NEhorFDRGJ6rfYG4i6lAdz0/p2lLmpdn8sjQqqv3qzm7OmiAwWixsiEk2evAqLdiYCAN4b7IcuHvZa/8yw262pP68WIV/O1hSRIWJxQ0SiEAQBC7YloKSyBoGtZZjyRPsW+dzWDlbo0aa+NbWPrSkig8TihohEEXkhG78n58PC1ARfTAiGmZbbUXcLD7p9rym2pogMEosbImpxuSVVWLyrvh01Y6gfOrnbtejnj7w9JfzPa0XILWFrisjQsLghohYlCAI+2BaP0qpaBHs74K3+7Vo8g4fMCj3bOt5uTfHqDZGhYXFDRC1qy7nrOJpyExZmJvhyQlCLtqPuxntNERkuFjdE1GJuFFfik92XAACzhnZEB9eWbUfd7c6Cfueu3UJOSaVoOYhI81jcEFGLEAQB8yLjUaqoRfc2DnhDhHbU3dxllnjMxxEAsDeBs6aIDAmLGyJqEZv/zMLxKwWQmtXPjjI1kYgd6a7WFO81RWRIWNwQkdZdv1WBT2+3o+YM74T2LrYiJ6oXFugBiQS4kFmMG8VsTREZChY3RKRVSqWAuVvjUV5dh5C2jnj1cV+xI6m42VviMR8nAMDeBA4sJjIULG6ISKv+czYTJ9MKYWlughU60o662517Te1hcUNkMFjcEJHWZBVVIGJvEgBg3ojO8G1lI3KixkYEuEMiAWIyi3H9VoXYcYhIA1jcEJFWKJUC5myNQ0V1HXr5OuGVUB+xIzXJ1c4SvW63pvZx1hSRQWBxQ0Ra8e/T13A6vQjWFqb44plgmOhYO+pud1pTu9maIjIILG6ISOOuFZZj6b5kAMAHYZ3Rxtla5ET3NzzAHSYSIC6rGFlFbE0R6TsWN0SkUUqlgDlb4lFZU4fQds54qXdbsSM9kKudJXr7OgPgrCkiQ8Dihog06qeTV3H2ahFsLEyx/JkgnW5H3S2cs6aIDAaLGyLSmPSbZVh+oL4dtSC8C7yddLsddbcRt1tT8ddL2Joi0nMsbohII+qUAuZsjUdVjRL9OrTCC73aiB1JLa1spQhtX9+a4tUbIv3G4oaINGL9Hxk4f+0WbKVmWDo+EBKJfrSj7hYe6AkA2BPP4oZIn7G4IaJHlppfhi8OpgAAPgzvAi9H/WlH3W24vxtMTSRIyC7BtcJyseMQ0UNicUNEj6ROKWD2ljgoapUY0NEFEx/zFjvSQ3O2lSK0HVtTRPqOxQ0RPZIfjqcjNqsYdpZmWKan7ai7qWZNsTVFpLdELW4iIiLw2GOPwc7ODq6urhg7dixSUlLu+5qjR49CIpE0eiQnJ7dQaiK640peKb46eBkA8I9RXeEhsxI50aMb7u8OUxMJEm/IkVHA1hSRPhK1uImOjsbUqVNx+vRpREVFoba2FsOGDUN5+YN/oKSkpCAnJ0f18PPza4HERHRHbZ0Ss7fEobpOiUGdXDChp5fYkTTCycYCfdtzQT8ifWYm5ofv37+/wfMNGzbA1dUV58+fx4ABA+77WldXVzg4OGgxHRHdz/fH0hF3vQT2lmaIGBek9+2ou40K8sDxKwXYHZ+DqYM6iB2HiNSkU2NuSkpKAABOTk4PPLZ79+7w8PDA4MGDceTIkXsep1AoIJfLGzyI6NEk58qx8lB9O+rjMf5wl1mKnEizhnV1h5mJBEk5cqTfLBM7DhGpSWeKG0EQMHPmTPTr1w8BAQH3PM7DwwNr165FZGQktm3bhk6dOmHw4ME4duxYk8dHRERAJpOpHt7e+juTg0gX1NxuR9XUCRjSxQ1Pd28tdiSNc7SxwOMdWgFga4pIH0kEQRDEDgEAU6dOxZ49e/DHH3/Ay0u93v3o0aMhkUiwc+fORvsUCgUUCoXquVwuh7e3N0pKSmBvb//IuYmMzarfr+CrqMuQWZkj6v0BcLU3rKs2d/zvXBbmbo1HZ3c77J9x/zY5EekWnbhy8+6772Lnzp04cuSI2oUNAPTp0wdXrlxpcp9UKoW9vX2DBxE9nEs35Pjn4fq/a0ue8jfYwgYAhnd1h7mpBMm5pUjNZ2uKSJ+IWtwIgoBp06Zh27ZtOHz4MHx9fR/qfWJiYuDh4aHhdER0t+rav9pRw/3dMCbYU+xIWiWzNmdrikhPiTpbaurUqfjll1/w22+/wc7ODrm5uQAAmUwGK6v69TLmz5+P7OxsbNy4EQCwcuVK+Pj4wN/fH9XV1di0aRMiIyMRGRkp2vcgMgbfHknFpRw5HK3N8elY/V+srznCAz1wNOUm9sTnYPpgLjdBpC9ELW5Wr14NABg4cGCD7Rs2bMDkyZMBADk5OcjMzFTtq66uxuzZs5GdnQ0rKyv4+/tjz549GDlyZEvFJjI6F7NL8O2RVADAkqcC4GInFTlRyxjW1R0LTBOQkleKK3ml8HOzEzsSETWDzgwobilyuRwymYwDiomaqbpWiTH/+gPJuaUYGeiOb1/oYRRXbe547ac/cTg5HzOG+GHGkI5ixyGiZtCJAcVEpLv+efgKknNL4WxjgU+eCjCqwgaob00BvNcUkT5hcUNE9xR/vRjfHU0DAHw6NgDOtsbRjrrbkK5usDA1wZX8MlzOKxU7DhE1A4sbImqSorYOs7fEoU4pYHSwJ8ICjXNGoszKHAM61s+a4tUbIv3A4oaImvTNoSu4nFeGVrYWWDzGX+w4ogoPut2aSsiBkQ1TJNJLGiluiouLNfE2RKQjYrOKsSb6TjsqEE42FiInEtfgLvWtqdT8MlzO44J+RLpO7eJm2bJl+PXXX1XPn332WTg7O6N169aIi4vTaDgianlVNXWY9b9YKAVgbDdPjAhwFzuS6OwtzTGgowsAYE/8DZHTENGDqF3cfP/996qbT0ZFRSEqKgr79u1DWFgY5syZo/GARNSyvo66jLSb5XCxk+JjI29H3W3U7dbUbramiHSe2ov45eTkqIqb3bt349lnn8WwYcPg4+OD3r17azwgEbWc89duYe3xdABAxNOBcLA27nbU3QZ3cYWFmQnSb5YjObcUXTy4ThaRrlL7yo2joyOysrIAAPv378eQIUMA1N8nqq6uTrPpiKjFVNXUYc6WOAgCMK5Hawzp6iZ2JJ1iZ2mOgarWFGdNEekytYubcePG4YUXXsDQoUNRWFiIsLAwAEBsbCw6dOig8YBE1DK+OJCC9IJyuNlLsWgU21FN4awpIv2gdlvq66+/ho+PD7KysrB8+XLY2toCqG9XvfPOOxoPSETa9+fVIqw7kQEAWDouCDJrc5ET6abBXdwgNTNBRkE5knJK0dWTrSkiXcR7SxEZucrqOoR9cwxXCyvwbIgXlj8TLHYknTbl3+exPzEXUwe1x5zhncWOQ0RNeKh1blJSUjBt2jQMHjwYQ4YMwbRp05CSkqLpbETUApYfSMbVwgp4yCzx4aiuYsfRearWVDxbU0S6Su3iZuvWrQgICMD58+cRHByMoKAgXLhwAQEBAdiyZYs2MhKRlpxOL8SGE1cBAEvHB8Heku2oB3mysyukZia4WliBxBtyseMQURPUHnMzd+5czJ8/H0uWLGmwfdGiRZg3bx4mTJigsXBEpD3lilrM3RoPAHi+lzeeuD0TiO7PRmqGJzu7Yt/FXOxJyEFAa5nYkYjob9S+cpObm4tJkyY12v7SSy8hNzdXI6GISPuW7U9GZlEFWjtYYcHILmLH0StsTRHpNrWLm4EDB+L48eONtv/xxx/o37+/RkIRkXadTC3AxlPXAADLxgfBju0otTzZ2RWW5ibILKrAxWy2poh0jdptqTFjxmDevHk4f/48+vTpAwA4ffo0tmzZgsWLF2Pnzp0NjiUi3VKmqMXcyPp21Iu926CfXyuRE+kfawszDO7shj0JOdidcAOBXmxNEekStaeCm5g072KPRCLRyRWLORWcjN3C7Qn4z5lMeDlaYf+MAbCVqv07DgHYm5CDd/5zAd5OVjg2ZxAkEonYkYjoNrXbUkqlslkPXSxsiIzd8Ss38Z8zmQCA5c8EsbB5BIM6ucLK3BRZRZVIyC4ROw4R3eWh1rm5o6qqSlM5iEjLSqtqMO/27KhXQtuib3u2ox6FlYUpBndxBcB7TRHpGrWLm7q6OnzyySdo3bo1bG1tkZ5efwfhf/zjH1i3bp3GAxKRZny+Nwk3SqrQxska88K4sq4mjLo9a2o3Z00R6RS1i5vPPvsMP/30E5YvXw4LCwvV9sDAQPz4448aDUdEmhF9+Sb+ezYLALDimSBYW7AdpQkDO7nC2sIU2cWViLvO1hSRrlC7uNm4cSPWrl2LF198EaampqrtQUFBSE5O1mg4Inp0JZV/taNefdwHvds5i5zIcFiam2JwFzcAwJ74GyKnIaI71C5usrOz0aFDh0bblUolampqNBKKiDTn092XkCuvgo+zNebyRo8aFx7IBf2IdI3axY2/v3+Ti/ht2bIF3bt310goItKMw8l52HL+OiQS4IsJwbCyMH3wi0gtAzu5wMbCFDdKqhCTVSx2HCLCQyzit2jRIrz88svIzs6GUqnEtm3bkJKSgo0bN2L37t3ayEhED6GkogbztyUAAF5/3BchPk4iJzJMluamGNLVDb/F3sCe+Bz0aOModiQio6f2lZvRo0fj119/xd69eyGRSPDRRx8hKSkJu3btwtChQ7WRkYgewuLdiciTK9CulQ1mD+8kdhyDdqc1tS8hB0olW1NEYlPryk1tbS0+++wzvPbaa4iOjtZWJiJ6RFGX8rDtQjZMJMAXzwbD0pztKG0a0NEFtlIzVWuqZ1tevSESk1pXbszMzLBixQquPkykw4orqrFge3076s0B7dgmaQGW5qYY2vXOrCku6EckNrXbUkOGDMHRo0e1EIWINOHjnYm4WapAB1dbvD+ko9hxjMad1tRetqaIRKf2gOKwsDDMnz8fFy9eRM+ePWFjY9NgP+8ETiSe/RdzsSP2Rn07agLbUS2pf8dWsJOaIVdehQuZtziAm0hEGr0ruK7eCfxuvCs4Gaqi8moM+zoaBWXVeGdge8wdwTVtWtrMX2OxLSYbk/v64OMx/mLHITJaGr0ruK4XNkSG7KPfLqKgrBod3Wzx3hA/seMYpfAgtqaIdMEj3RWciHTD3oQc7I7PgamJBF9O6AapGdtRYujn1wp2lmbIL1Xg3LVbYschMlosboj0XEGZAh/uuAgAeGdgewR6yUROZLykZqYY1tUdAO81RSQmFjdEekwQBPxjx0UUlVejs7sd3n2S7Sixjbrdmtp3MRd1bE0RiYLFDZEe2x2fg30Xc2FmIsEXE4JhYca/0mJ7vEMr2N9pTV0tEjsOkVHiT0IiPXWzVIGPfqtvR017sgMCWrMdpQsszEww3P92ayqBC/oRiUHtdW7Kyspw/vx55ObmQiKRwM3NDT179oStra028hFREwRBwIc7EnCrogZdPewxdVAHsSPRXUYGeWDL+evYm5CLRaP9YWoiETsSkVFpdnFTW1uLWbNm4YcffkBVVRUsLCwgCAJqampgaWmJt956CytWrIC5ubk28xIRgJ1xN3AgMQ/mpvXtKHNTXoTVJY+3bwWZlTkKyhQ4m1GE0PbOYkciMirN/ok4a9YsREZGYsOGDSgqKkJVVRUUCgWKioqwYcMGbNu2DXPmzNFmViICkC+vwke/JQIApj/ph66eXIxS19S3pm7fayqBs6aIWlqzi5tffvkFGzduxMSJE+Hg4KDa7uDggIkTJ2LDhg34z3/+o42MRHSbIAhYsD0BJZU1CGwtw5SB7cWORPcQHuQJoP6WGLV1SpHTEBmXZhc3lZWVaNWq1T33Ozs7o7KyUiOhiKhp2y5k41BSPixMTdiO0nF92zvDwdocBWXVOJvBWVNELanZPxkHDRqEmTNnIi8vr9G+vLw8zJ07F08++aRGwxHRX3JLqrB4V3076r0hfujkbidyIrofc1MTjLg9a2o3Z00RtahmFzffffcd8vLy4OXlhe7du2PEiBEICwtD9+7d4eXlhby8PHz33XdqfXhERAQee+wx2NnZwdXVFWPHjkVKSsoDXxcdHY2ePXvC0tIS7dq1w5o1a9T6XCJ9IwgC5m+Lh7yqFsFeMvzfgHZiR6JmuHOvqQNsTRG1qGYXN97e3oiLi8POnTsxZswYtG3bFm3atMGYMWOwa9cuxMTEwMvLS60Pj46OxtSpU3H69GlERUWhtrYWw4YNQ3l5+T1fk5GRgZEjR6J///6IiYnBggULMH36dERGRqr12UT6ZMv56ziSchMWZvXtKDO2o/RCaDtnOFqbo7C8GmfYmiJqMRJBEHRmffCbN2/C1dUV0dHRGDBgQJPHzJs3Dzt37kRSUpJq25QpUxAXF4dTp0498DPkcjlkMhlKSkpgb89ZJqT7bhRXYvjXx1CqqMX8sM74vyc4iFifzN+WgP+ezcTzvdogYlyg2HGIjILGfv0rLy/HsWPHHuk9SkpKAABOTk73PObUqVMYNmxYg23Dhw/HuXPnUFNT0+h4hUIBuVze4EGkLwRBwAfbElCqqEX3Ng54oz/bUfomPLC+NbX/Yg5bU0QtRGPFTWpqKgYNGvTQrxcEATNnzkS/fv0QEBBwz+Nyc3Ph5ubWYJubmxtqa2tRUFDQ6PiIiAjIZDLVw9vb+6EzErW0X//MwrHLNyG93Y7iSrf6p087JzjZWOBWRQ1OpReKHYfIKOhM437atGmIj4/Hf//73wceK5E0/AF/p7P29+0AMH/+fJSUlKgeWVlZmglMpGXXb1Xg0z317dc5wzuhvQtvcaKPzExNMCLg9r2m4jlriqglNPv2C/drFQFAXV3dQ4d49913sXPnThw7duyBg5Ld3d2Rm5vbYFt+fj7MzMzg7Nx4iXOpVAqpVPrQ2YjEIAgC5kXGo0xRi5C2jnj1cV+xI9EjGBXogV/OZGJ/Yi4+GRvA9YmItKzZxY1CocDbb7+NwMCmB8Rdu3YNixcvVuvDBUHAu+++i+3bt+Po0aPw9X3wD/DQ0FDs2rWrwbaDBw8iJCSE97Uig/GfM5k4kVoIS3MTrGA7Su/18nVCK1sLFJRV42RaIZ7o6CJ2JCKD1uziplu3bvD29sYrr7zS5P64uDi1i5upU6fil19+wW+//QY7OzvVFRmZTAYrKysA9W2l7OxsbNy4EUD9zKh//etfmDlzJt58802cOnUK69ata1Y7i0gfZBVV4PO99e2oucM7w7eVjciJ6FHdaU1tOp2JvfE5LG6ItKzZ10bDw8NRXFx8z/1OTk6YNGmSWh++evVqlJSUYODAgfDw8FA9fv31V9UxOTk5yMzMVD339fXF3r17cfToUXTr1g2ffPIJVq1ahfHjx6v12US6SKkUMHdrPCqq69DLxwmT+/qIHYk0JDzw9r2mEnNRw1lTRFqlU+vctASuc0O6bOOpq/jot0RYmZti/4z+aOvMqzaGok4poPfnv6OgTIGfXn0MAzu5ih2JyGBxVBuRjrhWWI6IvckAgPkjO7OwMTCmJhKMDOSsKaKWwOKGSAcolQLmbI1HZU0d+rRzwku924odibRg5O0F/Q4k5qK6lq0pIm1hcUOkA34+dRVnM4pgbWGKFc8Ew4SzowzSYz5OcLGTQl5VixOpjRcdJSLNYHFDJLKMgnIs21/fjlowsgu8naxFTkTaYmoiwcjbC/rtZmuKSGtY3BCJqE4pYM6WOFTVKNGvQyu82LuN2JFIy8KD6mdNHbyUC0Xtwy9+SkT3xuKGSEQbTmTg3LVbsJWaYen4wCZvIUKGJaStI1ztpCitqsUfV9iaItKGhy5uEhMT0aNHD5w7d06TeYiMRtrNMqw4kAIAWBjeBV6ObEcZAxMTiWpg8Z4EtqaItOGhi5uffvoJcXFxWL9+vSbzEBmFOqWA2VvioKhVor9fKzz3GO9Wb0xGBdUXN1GJeWxNEWnBQxU3dXV1+OWXX/D+++/j119/RXV1taZzERm0H4+nIyazGHZSMywbH8R2lJHp0cYR7vaWKFXU4vhltqaINO2hipsDBw6gtrYWERERsLe3x2+//abpXEQGKzW/FF9GXQYA/GN0V3g6WImciFoaW1NE2vVQxc3PP/+M5557Dubm5njxxRfx008/aTgWkWGqrVNi1pZ4VNcqMaiTCyb09BI7EokkPKh+SnjUpTxU1bA1RaRJahc3xcXF2LVrl+ommS+//DKioqKQl5en8XBEhmbt8XTEZRXDztIMEePYjjJm3b0d4SGzRJmiFscu3xQ7DpFBUbu42bx5M3x9fdGzZ08AQKdOndCjRw9s2rRJ4+GIDElKbilWRl0BAHw82h/uMkuRE5GY2Joi0h61i5uff/4ZL7/8coNtL730EltTRPdRU6fE7C1xqK5TYkgXV4zr0VrsSKQDwm/PmjrE1hSRRqlV3GRlZSEvL69RcfP888+jsrISly9f1mg4IkOx5mgaErJLILMyx+dPc7E+qtfd2wGtHaxQXl2HoylsTRFpilrFjbe3N9LT09G6dcPfOp2dnZGamoqOHTtqNByRIUjKkWPV4fp21OIx/nC1ZzuK6kkkEowMrB9YvJetKSKN4e0XiLSopk6JWf+LQ02dgGFd3fBUN0+xI5GOuXOvqUNJbE0RaQqLGyIt+vZIKi7lyOFobY7P2I6iJgR7ydDawQoV1XU4mpIvdhwig8DihkhLEm+U4F+HUwEAS54KgIudVOREpIskEonqdgy749maItIEFjdEWlBdW9+OqlUKCAtwV/3Pi6gpd6aE/56Uj8pqtqaIHhWLGyIt+NfhK0jOLYWTjQU+GRvAdhTdV5CXDF6OVqisqcMRtqaIHpnaxU1WVhauX7+uen727FnMmDEDa9eu1WgwIn2VcL0E3x5NAwB88lQAWtmyHUX3J5FIVGve7GFriuiRqV3cvPDCCzhy5AgAIDc3F0OHDsXZs2exYMECLFmyROMBifSJorYOs7bEok4pYFSQh+p/WEQPMiqwftbU78l5qKiuFTkNkX5Tu7i5ePEievXqBQD43//+h4CAAJw8eRK//PILVykmo/fNoSu4nFeGVrYWWPJUgNhxSI8EtLZHGydrVNUocTiZrSmiR6F2cVNTUwOptP4y+6FDhzBmzBgAQOfOnZGTw8upZLzisoqxJrq+HfXp2EA42ViInIj0yd2tKS7oR/Ro1C5u/P39sWbNGhw/fhxRUVEYMWIEAODGjRtwdnbWeEAifVBVU4dZW+KgFICnunliRIC72JFID4XfnjV1ODkf5Qq2pogeltrFzbJly/D9999j4MCBeP755xEcHAwA2Llzp6pdRWRsvj50Gan5ZXCxk+Lj0f5ixyE95e9pDx9ntqaIHpWZOgcLggBfX19cu3YNdXV1cHR0VO176623YG1trfGARLruQuYt/HAsHQDw+dOBcGQ7ih5S/b2mPPDd0TTsic/B6GDeroPoYah15UYQBPj5+SEvL69BYQMAPj4+cHV11Wg4Il1XVVOH2bfbUeO6t8bQrm5iRyI9d2fczZGUfJSxNUX0UNQqbkxMTODn54fCwkJt5SHSK18eTEH6zXK42kmxiO0o0oCuHvbwbWUDRa0SvyfliR2HSC+pPeZm+fLlmDNnDi5evKiNPER649zVIvz4RwYAYOn4QMiszUVORIZAIpGoBhZzQT+ih6PWmBsAeOmll1BRUYHg4GBYWFjAysqqwf6ioiKNhSPSVZXV9e0oQQAm9PTCk53ZjiLNCQ/ywL+OpOLo5ZsoraqBnSULZyJ1qF3crFy5UgsxiPTL8gPJuFpYAXd7S3w4qqvYccjAdHa3QzsXG6TfLMfh5Hw81a212JGI9Iraxc0rr7yijRxEeuNMeiF+OnkVwO12lBV/qybNkkgkGBXogVWHU7E7PofFDZGaHuqu4Glpafjwww/x/PPPIz+/fi2G/fv3IzExUaPhiHRNRXUt5myNhyAAzz3mjYGdOEOQtCM8qH4aeHRKfWuKiJpP7eImOjoagYGBOHPmDLZt24aysjIAQHx8PBYtWqTxgES6ZNm+ZGQWVcBTZomF4V3EjkMGrKObLTq42qK6TolDnDVFpBa1i5sPPvgAn376KaKiomBh8ddiZYMGDcKpU6c0Go5Il5xMK8DPp64BAJY/E8xBnqRVdxb0Azhrikhdahc3CQkJePrppxttd3Fx4fo3ZLDKFbWYuzUeAPBC7zbo59dK5ERkDEbdXtDv2OUClFSyNUXUXGoXNw4ODk3e/TsmJgatW3PQGxmmiH1JuH6rEq0drLBgJNtR1DI6utnB705r6hJbU0TNpXZx88ILL2DevHnIzc2FRCKBUqnEiRMnMHv2bEyaNEkbGYlE9ceVAmw6nQkAWPFMEGylak8yJHpod27HsCeBrSmi5lK7uPnss8/Qpk0btG7dGmVlZejatSsGDBiAvn374sMPP9RGRiLRlFbVYF5kfTtqUmhb9O3AdhS1rDurFR+/cpOtKaJmkgiCIDzMC9PT03HhwgUolUp0794dfn5+ms6mFXK5HDKZDCUlJbC3txc7Dum4+dvi8d+zWfB2ssL+9wbAhldtSATDvz6GlLxSfDEhGM/09BI7DpHOU/vKzZIlS1BRUYF27drhmWeewbPPPgs/Pz9UVlZiyZIl2shIJIpjl2/iv2ezAAArnglmYUOiUbWm4m+InIRIP6hd3CxevFi1ts3dKioqsHjxYo2EIhKb/K521OS+PujTzlnkRGTMRqpaUwUoqWBriuhB1C5uBEGARCJptD0uLg5OTk4aCUUktk93X0JOSRV8nK0xd0QnseOQkevgaovO7naoVQo4cClX7DhEOq/ZxY2joyOcnJwgkUjQsWNHODk5qR4ymQxDhw7Fs88+q9aHHzt2DKNHj4anpyckEgl27Nhx3+OPHj0KiUTS6JGcnKzW5xLdz5HkfPzv3HVIJMCKCcGwtmA7isQXzgX9iJqt2T+1V65cCUEQ8Nprr2Hx4sWQyWSqfRYWFvDx8UFoaKhaH15eXo7g4GC8+uqrGD9+fLNfl5KS0mAwsIuLi1qfS3QvJRU1+GBbfTvqtcd98ZgPr0aSbhgZ5IEvoy7jRGoBbpVXw9HG4sEvIjJSzS5u7twN3NfXF48//jjMzB79t9mwsDCEhYWp/TpXV1c4ODg88ucT/d2S3ZeQJ1egXSsbzB7GdhTpjvYutujiYY+kHDkOXsrFxMfaiB2JSGepPebGzs4OSUlJque//fYbxo4diwULFqC6ulqj4e6le/fu8PDwwODBg3HkyJH7HqtQKCCXyxs8iJpy6FIeIi9ch8ntdpSVhanYkYgauHM7ht1sTRHdl9rFzf/93//h8uXLAOrXupk4cSKsra2xZcsWzJ07V+MB7+bh4YG1a9ciMjIS27ZtQ6dOnTB48GAcO3bsnq+JiIiATCZTPby9vbWakfRTcUU15m9PAAC82b8derZ1FDkRUWN3Zk2dTCtEUXnL/DJJpI/UXsRPJpPhwoULaN++PZYtW4bDhw/jwIEDOHHiBJ577jlkZWU9XBCJBNu3b8fYsWPVet3o0aMhkUiwc+fOJvcrFAooFArVc7lcDm9vby7iRw3M2ByDHbE30N7FBnum94elOa/akG4KX3UciTfkWDouEM/1YmuKqCkPNRVcqVQCAA4dOoSRI0cCALy9vVFQUKDZdM3Qp08fXLly5Z77pVIp7O3tGzyI7nYgMRc7Ym/ARAJ8MSGYhQ3pNN5riujB1C5uQkJC8Omnn+Lf//43oqOjER4eDgDIyMiAm5ubxgM+SExMDDw8PFr8c8kwFJVXY+HtdtT/PdEe3duwHUW6Lfyu1lRhmeIBRxMZJ7WnPK1cuRIvvvgiduzYgYULF6JDhw4AgK1bt6Jv375qvVdZWRlSU1NVzzMyMhAbGwsnJye0adMG8+fPR3Z2NjZu3Kj6bB8fH/j7+6O6uhqbNm1CZGQkIiMj1f0aRACARTsTUVBWjY5utpgxRD/uj0bGra2zDQJby5CQXYIDiXl4oTdbU0R/p3ZxExQUhISEhEbbV6xYAVNT9S7nnzt3DoMGDVI9nzlzJoD6aec//fQTcnJykJmZqdpfXV2N2bNnIzs7G1ZWVvD398eePXtUrTEidexLyMGuuBswNZHgiwnBkJqxHUX6YWSgBxKyS7An4QaLG6ImPPRdwfUV7wpOAFBYpsCwr4+hsLwa0wZ1wOzhXNOG9EdmYQUGrDgCEwlwduEQtLKVih2JSKeoPebGxMQEpqam93wQ6YOPfktEYXk1Orvb4d3BHcSOQ6SWNs7WCPKSQSkA+y/yXlNEf6d2W2r79u0NntfU1CAmJgY///wz7wpOemF3/A3sSciBGdtRpMfCAz0Qf70Ee+Jz8FKftmLHIdIpGmtL/fLLL/j111/x22+/aeLttIZtKeN2s1SBYV9H41ZFDaYP9sPMoR3FjkT0ULKKKtB/eX1r6syCIXCxY2uK6A6121L30rt3bxw6dEhTb0ekcYIg4MMdCbhVUYMuHvaYNojtKNJf3k7WCPZ2qG9NJbI1RXQ3jRQ3lZWV+Oc//wkvLy9NvB2RVuyMu4EDiXkwM5HgywnBsDDTWG1PJIpRt9e82RN/Q+QkRLpF7TE3jo6OkEgkqueCIKC0tBTW1tbYtGmTRsMRaUq+vAof/ZYIAJg+2A9dPdmSJP0XFuiOz/Ym4UxGEfJLq+BqZyl2JCKdoHZx8/XXXzcobkxMTODi4oLevXvD0ZGru5LuEQQBC7YnoKSyBgGt7fH2wPZiRyLSCC9Ha3TzdkBsVjH2X8zFpFAfsSMR6QS1i5vJkydrIQaR9myPycahpHyYm9bPjjI3ZTuKDMeoIA/EZhVjd3wOixui25pV3MTHxzf7DYOCgh46DJGm5cmr8PHO+nbUjCEd0dmd7SgyLGGBHvh0TxL+vFqEPHkV3OzZmiJqVnHTrVs3SCQSPGjWuEQiQV1dnUaCET0qQRAwf1sC5FW1CPKS4f8GtBM7EpHGtXawQo82DriQWYx9CTmY/Liv2JGIRNes4iYjI0PbOYg0buv56zicnA8LUxN8OSEYZmxHkYEKD/LEhcxi7GFxQwSgmcVN27Zc/ZL0S05JJZbsugQAeH9oR/i52YmciEh7Rga645Pdl3Du2i3kllTBXcbWFBk3tX+VjYiIwPr16xttX79+PZYtW6aRUESPQhAEfBCZgFJFLbp5O+DN/vxNlgybh8wKIW0dIQjAvos5YschEp3axc3333+Pzp07N9ru7++PNWvWaCQU0aP437ksRF++CQszE3zBdhQZifCgOwv6sbghUvunfm5uLjw8PBptd3FxQU4O/1KRuLKLK/HJ7iQAwJxhndDB1VbkREQtIyzAAxIJcO7aLeSUVIodh0hUahc33t7eOHHiRKPtJ06cgKenp0ZCET0MQRAwb2s8yhS16NnWEa/1YzuKjIe7zBIhbesXUt2bwHtNkXFTexG/N954AzNmzEBNTQ2efPJJAMDvv/+OuXPnYtasWRoPSNRcv5zNxB+pBZCamWDFM0EwNZE8+EVEBiQ80AN/Xr2FPfE38DqLezJiahc3c+fORVFREd555x1UV1cDACwtLTFv3jzMnz9f4wGJmiOrqAKf76lvR80d0RntXNiOIuMTFuiBxbsv4UJmMbKLK9HawUrsSESikAgPWpnvHsrKypCUlAQrKyv4+flBKpVqOptWyOVyyGQylJSUwN6eq9UaAqVSwEvrzuBkWiF6+Thh81t9YMKrNmSknv3+FM5mFOHD8C54oz8XriTj9NDTSGxtbfHYY48hICBAbwobMkz/OXMNJ9MKYWVuiuXPBLGwIaM26vasqd2cNUVGjHNkSa9lFlbg873JAIB5IzrBp5WNyImIxDUiwB0SCRCbVYzrtyrEjkMkChY3pLeUSgFztsahsqYOvX2deEdkIgCudpbo7esEANjHWVNkpFjckN7aeOoqzmQUwdrCFCueCWY7iui28KD6ZTl2J7A1RcapWcVNjx49cOvWLQDAkiVLUFHBS50krqsF5Vi6v74dNX9kF7RxthY5EZHuGOHvDhMJEJdVjKwi/rwm49Os4iYpKQnl5eUAgMWLF6OsrEyroYjup04pYPaWOFTVKNG3vTNe7NVG7EhEOsXFTorevs4AgL28ekNGqFnr3HTr1g2vvvoq+vXrB0EQ8MUXX8DWtul1RD766CONBiT6uw0nMnDu2i3YWJhi2XjOjiJqSniQB06lF2JPQg7+74n2YschalHNWucmJSUFixYtQlpaGi5cuICuXbvCzKxxXSSRSHDhwgWtBNUUrnOj39JvliHsm+NQ1Crx+dOBeKE3r9oQNaWgTIFenx2CUgCOzRnE1i0ZlWZduenUqRM2b94MADAxMcHvv/8OV1dXrQYj+rs77ShFrRL9/Vrh+V7eYkci0lmtbKUIbe+ME6n1V2/eHsirN2Q81J4tpVQqWdiQKNb9kY4LmcWwk5ph2fggSCRsRxHdT3hg/aypPQk3RE5C1LIeaip4Wloa3n33XQwZMgRDhw7F9OnTkZaWpulsRCqp+aX44uBlAMCHo7rAk/fMIXqg4f5uMDWR4GK2HNcKy8WOQ9Ri1C5uDhw4gK5du+Ls2bMICgpCQEAAzpw5A39/f0RFRWkjIxm52jolZm2JR3WtEk90dMGzIWxHETWHs60UfdvXz5raw1lTZETUvnFm9+7dMXz4cCxdurTB9g8++AAHDx7kgGLSuNVH07BsfzLsLM1w8P0B8JDxqg1Rc20+m4kPtiXA39Mee6b3FzsOUYtQ+8pNUlISXn/99UbbX3vtNVy6dEkjoYjuuJxXiq+j6ttRi0b7s7AhUtNwf3eYmkiQeEOOjAK2psg4qF3cuLi4IDY2ttH22NhYDjQmjaqpU2LW/+JQXafEk51dMb5Ha7EjEekdRxsLVWuKC/qRsWjWVPC7vfnmm3jrrbeQnp6Ovn37QiKR4I8//sCyZcswa9YsbWQkI/V9dBoSsktgb2mGiHGBnB1F9JBGBXng+JUC7I7PwdRBHcSOQ6R1ao+5EQQBK1euxJdffokbN+qnF3p6emLOnDmYPn26zv8PiGNu9ENyrhyj//kHauoEfD0xGE939xI7EpHeulVejcc+O4RapYDfZz2B9i5NrzBPZCjULm7uVlpaCgCws7PTWCBtY3Gj+2rqlBj77Qkk3pBjaFc3rH25p84XzUS67pX1ZxF9+SZmDe2Idwf7iR2HSKseap2bO+zs7PSqsCH98N2RNCTekMPB2hyfPR3AwoZIA8KDPABwSjgZh0cqbog0LfFGCf55+AoAYPEYf7jaWYqciMgwDO/qDnNTCZJzS5GaXyZ2HCKtYnFDOqO6VonZW+JRqxQwwt8dY4I9xY5EZDBk1ubo16EVAM6aIsPH4oZ0xr+OpCIpRw4nGwt8ynYUkcaFB92+11Q8ixsybGoVNzU1NRg0aBAuX76srTxkpC5ml+DbI6kAgE+eCkArW6nIiYgMz9CubjA3lSAlrxRX8krFjkOkNWoVN+bm5rh48SJ/oyaNUtTWYdb/4lCnFBAe6KEa+EhEmiWzMkd/PxcAHFhMhk3tttSkSZOwbt06bWQhI7Xq9ytIySuFs40FljzlL3YcIoMWHnh71hRbU2TA1F6huLq6Gj/++COioqIQEhICGxubBvu/+uorjYUjwxeXVYw10ekAgE/HBsCZ7SgirRrS1Q0Wpia4kl+Gy3ml6OjG5TzI8Kh95ebixYvo0aMH7O3tcfnyZcTExKgeTd1z6n6OHTuG0aNHw9PTExKJBDt27Hjga6Kjo9GzZ09YWlqiXbt2WLNmjbpfgXREVU0dZm+pb0eNCfZEWCDbUUTaJrMyx4CO9bOmdvPqDRkota/cHDlyRGMfXl5ejuDgYLz66qsYP378A4/PyMjAyJEj8eabb2LTpk04ceIE3nnnHbi4uDTr9aRbVh66giv5ZWhlK8XiMWxHEbWU8CAPHErKx96EHLw/xI/jKMngqF3c3JGamoq0tDQMGDAAVlZWEARB7b8gYWFhCAsLa/bxa9asQZs2bbBy5UoAQJcuXXDu3Dl88cUXohc3dUoBOSWVombQJ2k3y7H2WBoA4POnA+BoYyFyIiLjMaSLGyzMTJCaX4bLeWXo5M7WFBkWtYubwsJCPPvsszhy5AgkEgmuXLmCdu3a4Y033oCDgwO+/PJLbeQEAJw6dQrDhg1rsG348OFYt24dampqYG5u3ug1CoUCCoVC9Vwul2slW2G5Av2Wae6qlrF4untrDPN3FzsGkVGxszTHEx1dEHUpD5v/zMSi0bxySoZF7TE377//PszNzZGZmQlra2vV9okTJ2L//v0aDfd3ubm5cHNza7DNzc0NtbW1KCgoaPI1ERERkMlkqoe3t7fW8knNTPhoxsPSvP6fXT3ssWh0V6399yCie3uhVxsAwM8nryI2q1jcMEQapvaVm4MHD+LAgQPw8vJqsN3Pzw/Xrl3TWLB7+Xvr685Nze/VEps/fz5mzpypei6Xy7VS4LjaWSLl0+a32IiIxDSosyvGBHtiZ9wNzPpfLPZM7w9Lc1OxYxFphNpXbsrLyxtcsbmjoKAAUql2p/G6u7sjNze3wbb8/HyYmZnB2dm5yddIpVLY29s3eBARUf3NaVvZSpF2sxxfH+LK82Q41C5uBgwYgI0bN6qeSyQSKJVKrFixAoMGDdJouL8LDQ1FVFRUg20HDx5ESEhIk+NtiIjo3hxtLPD50wEAgB+OpeP8tVsiJyLSDLWLmxUrVuD7779HWFgYqqurMXfuXAQEBODYsWNYtmyZWu9VVlaG2NhY1fo4GRkZiI2NRWZmJoD6ltKkSZNUx0+ZMgXXrl3DzJkzkZSUhPXr12PdunWYPXu2ul+DiIgADPN3x7juraEUgDlb4lBVUyd2JKJHpnZx07VrV8THx6NXr14YOnQoysvLMW7cOMTExKB9+/Zqvde5c+fQvXt3dO/eHQAwc+ZMdO/eHR999BEAICcnR1XoAICvry/27t2Lo0ePolu3bvjkk0+watUq0aeBExHps0Wj/eFqJ0V6QTm+OJAidhyiRyYR7ozINRJyuRwymQwlJSUcf0NEdNvh5Dy89tM5SCTAlv8LRYiPk9iRiB7aQxU3t27dwrp165CUlASJRIIuXbrg1VdfhZOT7v9lYHFDRNS02VvisPX8dfg4W2PfewNgZcHZU6Sf1G5LRUdHw9fXF6tWrcKtW7dQVFSEVatWwdfXF9HR0drISERELeAfo7rC3d4SVwsrsPxAsthxiB6a2lduAgIC0LdvX6xevRqmpvVVfV1dHd555x2cOHECFy9e1EpQTeGVm5bzMLfkICJxHU3Jx+QNfwIANr/VB33aNb3MBpEuU/vKTVpaGmbNmqUqbADA1NQUM2fORFpamkbDkf5iYUOknwZ2csXEkPqFTudujUdFda3IiYjUp3Zx06NHDyQlJTXanpSUhG7dumkiExkAFjZE+mvhqC7wlFkis6gCy/axPUX6p1m3X4iPj1f9+/Tp0/Hee+8hNTUVffr0AQCcPn0a3377LZYuXaqdlERE1GLsLc2x7JkgvLzuLH4+dQ3DA9zRt30rsWMRNVuzxtyYmJhAIpHgQYdKJBLU1en2AlAcc0NE1DwLtifglzOZ8HK0wv4ZA2ArVft2hESiaNaf1IyMDG3nICIiHbNgZBdEp9zE9VuViNibhM+eDhQ7ElGzcBE/IiK6p5OpBXjhxzMAgE2v90Y/P7anSPc9VHGTnZ2NEydOID8/H0qlssG+6dOnayycNrC4ISJSzz92XMS/T19Dawcr7J/RH3aWvFEx6Ta1i5sNGzZgypQpsLCwgLOzc4NZMRKJBOnp6RoPqUksboiI1FOuqMWIb44hq6gSz/fyRsS4ILEjEd2X2sWNt7c3pkyZgvnz58PERO2Z5KJjcaN5XNOGyPCdTi/Ec2tPAwB+fq0XnujoInIiontTuzqpqKjAc889p5eFDWkHCxsiw9ennTMm9/UBAHwQGQ95VY24gYjuQ+0K5fXXX8eWLVu0kYWIiHTY3BGd0NbZGjklVfh09yWx4xDdk9ptqbq6OowaNQqVlZUIDAyEuXnDgWVfffWVRgNqGttSREQP72xGESauPQVBADZMfgyDOruKHYmoEbVXZPr8889x4MABdOrUCQAaDSgmIiLD1cvXCa897ot1f2Tgg23xODjjCcisOXuKdIvaV24cHR3x9ddfY/LkyVqKpF28ckNE9Ggqq+sQvuo40gvKMa5Ha3z1bDexIxE1oPaYG6lUiscff1wbWYiISA9YWZhixYQgSCTAtgvZOHQpT+xIRA2oXdy89957+Oc//6mNLEREpCd6tnXCm/3bAQDmb09AcUW1yImI/qJ2W+rpp5/G4cOH4ezsDH9//0YDirdt26bRgJrGthQRkWZU1dS3p9JulmNsN0+sfK672JGIADzEgGIHBweMGzdOG1lIx3GxPiK6m6W5Kb6YEIzxq09iR+wNjAjwwIgAd7FjEfHGmURE9GiW7kvGmug0tLK1wMH3n4CTjYXYkcjIcZlhIiJ6JDOG+MHP1RYFZdVYtDNR7DhE6l+58fX1vW9rgjfOJCIyPnFZxRi3+iTqlAK+e7EHRgZ6iB2JjJjaY25mzJjR4HlNTQ1iYmKwf/9+zJkzR1O5SEdwnA0RNUewtwPefqI9/nUkFR/uuIhevk5oZSsVOxYZKY2Nufn2229x7tw5bNiwQRNvpzW8ckNEpB2K2jo89a8TSM4txchAd3z3Yk+xI5GR0tiYm7CwMERGRmrq7YiISM9IzepnT5maSLA3IRe742+IHYmMlMaKm61bt8LJyUlTb0dERHoooLUMUwd1AAD8Y8dF3CxViJyIjJHaY266d+/eYAyGIAjIzc3FzZs38d1332k0HBER6Z9pgzog6lIeknLk+HBHAta81JNj96hFqV3cjB07tsFzExMTuLi4YODAgejcubOmchERkZ6yMDPBlxOCMeZff+BAYh52xt3AU91aix2LjAgX8SMiIq345tAVfH3oMmRW5oh6fwBc7S3FjkRGgov4ERGRVrwzqD38Pe1RUlmDBdsTYGS/S5OIml3cmJiYwNTU9L4PMzO1u1xERGSgzE1N8OWzwTA3leBQUj62x2SLHYmMRLOrke3bt99z38mTJ/HPf/6TVbke42J9RKQNnd3tMWNIR6w4kIKPdyaib/tWcJexPUXa9UhjbpKTkzF//nzs2rULL774Ij755BO0adNGk/k0jmNuiIhaVm2dEuNWn0T89RIM6uSC9ZMf4y9TpFUPNebmxo0bePPNNxEUFITa2lrExsbi559/1vnChoiIWp6Zaf3sKQtTExxJuYkt56+LHYkMnFrFTUlJCebNm4cOHTogMTERv//+O3bt2oWAgABt5SMiIgPg52aH94d2BAB8susSckoqRU5EhqzZxc3y5cvRrl077N69G//9739x8uRJ9O/fX5vZiIjIgLzZ3xfdvB1QqqjFvEjOniLtafaYGxMTE1hZWWHIkCEwNTW953Hbtm3TWDht4JgbIiLxpOaXYeSq46iuVWLpuEA814vDGUjzmj1batKkSRwARkREj6SDqy1mD+uIz/cm49M9Sejf0QWtHazEjkUGhisUExFRi6pTCpiw5iQuZBajX4dW+PfrvfjLM2kUVygmIqIWZWoiwRcTgiE1M8EfqQX45Wym2JHIwLC4MTJGdqGOiHRUOxdbzB1Rf7Plz/YkIauoQuREZEhY3BgZXvolIl3xal8f9PJxQkV1HeZujYdSyV++SDNY3BARkShMTCRY/kwQLM1NcCq9EJvOXBM7EhkI0Yub7777Dr6+vrC0tETPnj1x/Pjxex579OhRSCSSRo/k5OQWTExERJri08oGH9xuT0XsTUZmIdtT9OhELW5+/fVXzJgxAwsXLkRMTAz69++PsLAwZGbef3BZSkoKcnJyVA8/P78WSkxERJo2KdQHvX2dUFlTh9lb49ieokcmanHz1Vdf4fXXX8cbb7yBLl26YOXKlfD29sbq1avv+zpXV1e4u7urHvdbVJCIiHSbiYkEK54JhrWFKc5mFOHnU1fFjkR6TrTiprq6GufPn8ewYcMabB82bBhOnjx539d2794dHh4eGDx4MI4cOXLfYxUKBeRyeYMHERHpljbO1pgfVt+eWrY/GRkF5SInIn0mWnFTUFCAuro6uLm5Ndju5uaG3NzcJl/j4eGBtWvXIjIyEtu2bUOnTp0wePBgHDt27J6fExERAZlMpnp4e3tr9HsQEZFmvNi7Lfq2d0ZVjRJztsShju0pekiirVB848YNtG7dGidPnkRoaKhq+2effYZ///vfzR4kPHr0aEgkEuzcubPJ/QqFAgqFQvVcLpfD29ubKxQTEemgrKIKjFh5DOXVdfgwvAve6N9O7Eikh0S7ctOqVSuYmpo2ukqTn5/f6GrO/fTp0wdXrly5536pVAp7e/sGDyIi0k3eTtZYGN4VALDiQArSbpaJnIj0kWjFjYWFBXr27ImoqKgG26OiotC3b99mv09MTAw8PDw0HY+IiETyfC9v9PdrBUWtErPZnqKH0Oy7gmvDzJkz8fLLLyMkJAShoaFYu3YtMjMzMWXKFADA/PnzkZ2djY0bNwIAVq5cCR8fH/j7+6O6uhqbNm1CZGQkIiMjxfwaRESkQRKJBEvHB2H418cQk1mMH4+n4/+eaC92LNIjohY3EydORGFhIZYsWYKcnBwEBARg7969aNu2LQAgJyenwZo31dXVmD17NrKzs2FlZQV/f3/s2bMHI0eOFOsrEBGRFrR2sMI/RnXBvMgEfBl1GYO7uKKDq53YsUhPiDagWCxyuRwymYwDiomIdJwgCJi84U9EX76JYG8HRE4JhZmp6Avrkx7gnxIiItJJ9e2pQNhZmiEuqxhrj6eLHYn0BIsbIiLSWR4yK3w0qn721MqoK0jJLRU5EekDFjdERKTTnunphSc7u6K6rn72VE2dUuxIpONY3BARkU6TSCSIGBcIe0szJGSX4PvoNLEjkY5jcUNERDrPzd4Si5/yBwB88/sVJOXwPoF0byxuiIhIL4zt1hpDu7qhpk7ArP+xPUX3xuKGiIj0gkQiwWdPB8DB2hyXcuT49kiq2JFIR7G4ISIiveFqZ4nFY+rbU/86nIrEGyUiJyJdxOKGiIj0yphgT4zwd0etsr49VV3L9hQ1xOKGiIj0ikQiwadPB8DJxgLJuaX41+ErYkciHcPihoiI9E4rWymW3J499e3RNCRcZ3uK/sLihoiI9NKoIE+EB3qgTilg9pY4KGrrxI5EOoLFDRER6a0lT/nD2cYCKXmlWPU721NUj8UNERHpLWdbKT4dGwAAWH00DXFZxeIGIp3A4oaIiPRaWKAHxgR7QikAs7bEoaqG7Sljx+KGiIj03uIx/mhlK0Vqfhm+PnRZ7DgkMhY3RESk9xxtLPD50/XtqR+OpeNC5i2RE5GYWNwQEZFBGObvjqe7t4ZSAGazPWXUWNwQEZHBWDS6K1ztpEi/WY4vD6aIHYdEwuKGiIgMhoO1BSLGBQIAfvwjA+euFomciMTA4oaIiAzK4C5uGN/DC4IAzNkaj8pqtqeMDYsbIiIyOB+N7go3eykyCsqx4gDbU8aGxQ0RERkcmZU5lo4PAgBsOJmBM+mFIieilsTihoiIDNKgTq6YGOKtak9VVNeKHYlaCIsbIiIyWAtHdYGHzBKZRRVYti9Z7DjUQljcEBGRwbK3NMey2+2pn09dw6k0tqeMAYsbIiIyaAM6uuD5Xm0AAHO2xqFcwfaUoWNxQ0REBm9heBe0drDC9VuViNiXJHYc0jIWN0REZPBspWZY/kx9e2rT6UycSC0QORFpE4sbIiIyCo93aIWX+tS3p+ZujUdpVY3IiUhbWNwQEZHRmB/WBV6OVsgursTnezl7ylCxuCEiIqNhIzXDimeCAQD/PZuJY5dvipyItIHFDRERGZXQ9s6Y3NcHADAvMh5ytqcMDosbIiIyOnNHdEJbZ2vklFTh092XxI5DGsbihoiIjI61RX17SiIB/nfuOo6k5IsdiTSIxQ0RERmlXr5OeLWvLwDgg8h4lFSwPWUoWNwQEZHRmjO8E3xb2SBPrsAStqcMBosbIiIyWlYWpljxTBAkEiDywnUcupQndiTSABY3RERk1EJ8nPBGv/r21PztCSiuqBY5ET0qFjdERGT0Zg3rhHYuNrhZqsDiXWxP6TsWN0REZPQszU3xxYRgmEiA7THZOJCYK3YkegQsboiIiAD0aOOItwa0BwAs3J6AonK2p/QVixsiIqLbZgzxg5+rLQrKqrFoZ6LYceghsbghIiK67U57ytREgl1xN7AvIUfsSPQQWNwQERHdJdjbAVOeaAcA+HDHRRSWKUROROoSvbj57rvv4OvrC0tLS/Ts2RPHjx+/7/HR0dHo2bMnLC0t0a5dO6xZs6aFkhIRkbGYPtgPndzsUFhejY9+Y3tK34ha3Pz666+YMWMGFi5ciJiYGPTv3x9hYWHIzMxs8viMjAyMHDkS/fv3R0xMDBYsWIDp06cjMjKyhZMTEZEhk5r91Z7ak5CD3fE3xI5EapAIgiCI9eG9e/dGjx49sHr1atW2Ll26YOzYsYiIiGh0/Lx587Bz504kJSWptk2ZMgVxcXE4depUsz5TLpdDJpOhpKQE9vb2j/4liIjIYH11MAWrDqfC0doc/32rD2ylZmJH0gumJhJ4yKxE+3zR/itVV1fj/Pnz+OCDDxpsHzZsGE6ePNnka06dOoVhw4Y12DZ8+HCsW7cONTU1MDc3b/QahUIBheKvfqlcLtdAeiIiMgbTnvTDwUt5SM4txYiV9x82QX9xtZPi7MIhon2+aMVNQUEB6urq4Obm1mC7m5sbcnObXjwpNze3yeNra2tRUFAADw+PRq+JiIjA4sWLNReciIiMhoWZCb6e2A2v/fQn171Rg9Rc3CG9ol9fk0gkDZ4LgtBo24OOb2r7HfPnz8fMmTNVz+VyOby9vR82LhERGZkuHvY4NX+w2DFIDaIVN61atYKpqWmjqzT5+fmNrs7c4e7u3uTxZmZmcHZ2bvI1UqkUUqlUM6GJiIhI54l23cjCwgI9e/ZEVFRUg+1RUVHo27dvk68JDQ1tdPzBgwcREhLS5HgbIiIiMj6iNsVmzpyJH3/8EevXr0dSUhLef/99ZGZmYsqUKQDqW0qTJk1SHT9lyhRcu3YNM2fORFJSEtavX49169Zh9uzZYn0FIiIi0jGijrmZOHEiCgsLsWTJEuTk5CAgIAB79+5F27ZtAQA5OTkN1rzx9fXF3r178f777+Pbb7+Fp6cnVq1ahfHjx4v1FYiIiEjHiLrOjRi4zg0REZFhE/32C0RERESaxOKGiIiIDAqLGyIiIjIoLG6IiIjIoLC4ISIiIoPC4oaIiIgMCosbIiIiMigsboiIiMigsLghIiIigyLq7RfEcGdBZrlcLnISIiIiUpednR0kEsl9jzG64qa0tBQA4O3tLXISIiIiUldzbp9kdPeWUiqVuHHjRrMqP3XJ5XJ4e3sjKyuL9616AJ6r5uO5aj6eK/XwfDUfz1Xzaftc8cpNE0xMTODl5aXVz7C3t+cf/mbiuWo+nqvm47lSD89X8/FcNZ+Y54oDiomIiMigsLghIiIig8LiRoOkUikWLVoEqVQqdhSdx3PVfDxXzcdzpR6er+bjuWo+XThXRjegmIiIiAwbr9wQERGRQWFxQ0RERAaFxQ0REREZFBY3REREZFBY3DzAsWPHMHr0aHh6ekIikWDHjh0N9guCgI8//hienp6wsrLCwIEDkZiY2OAYhUKBd999F61atYKNjQ3GjBmD69evt+C3aBkRERF47LHHYGdnB1dXV4wdOxYpKSkNjuH5qrd69WoEBQWpFrkKDQ3Fvn37VPt5nu4tIiICEokEM2bMUG3j+ar38ccfQyKRNHi4u7ur9vM8NZSdnY2XXnoJzs7OsLa2Rrdu3XD+/HnVfp6vv/j4+DT6syWRSDB16lQAOniuBLqvvXv3CgsXLhQiIyMFAML27dsb7F+6dKlgZ2cnREZGCgkJCcLEiRMFDw8PQS6Xq46ZMmWK0Lp1ayEqKkq4cOGCMGjQICE4OFiora1t4W+jXcOHDxc2bNggXLx4UYiNjRXCw8OFNm3aCGVlZapjeL7q7dy5U9izZ4+QkpIipKSkCAsWLBDMzc2FixcvCoLA83QvZ8+eFXx8fISgoCDhvffeU23n+aq3aNEiwd/fX8jJyVE98vPzVft5nv5SVFQktG3bVpg8ebJw5swZISMjQzh06JCQmpqqOobn6y/5+fkN/lxFRUUJAIQjR44IgqB754rFjRr+XtwolUrB3d1dWLp0qWpbVVWVIJPJhDVr1giCIAjFxcWCubm5sHnzZtUx2dnZgomJibB///4Wyy6G/Px8AYAQHR0tCALP14M4OjoKP/74I8/TPZSWlgp+fn5CVFSU8MQTT6iKG56vvyxatEgIDg5uch/PU0Pz5s0T+vXrd8/9PF/399577wnt27cXlEqlTp4rtqUeQUZGBnJzczFs2DDVNqlUiieeeAInT54EAJw/fx41NTUNjvH09ERAQIDqGENVUlICAHBycgLA83UvdXV12Lx5M8rLyxEaGsrzdA9Tp05FeHg4hgwZ0mA7z1dDV65cgaenJ3x9ffHcc88hPT0dAM/T3+3cuRMhISGYMGECXF1d0b17d/zwww+q/Txf91ZdXY1Nmzbhtddeg0Qi0clzxeLmEeTm5gIA3NzcGmx3c3NT7cvNzYWFhQUcHR3veYwhEgQBM2fORL9+/RAQEACA5+vvEhISYGtrC6lUiilTpmD79u3o2rUrz1MTNm/ejAsXLiAiIqLRPp6vv/Tu3RsbN27EgQMH8MMPPyA3Nxd9+/ZFYWEhz9PfpKenY/Xq1fDz88OBAwcwZcoUTJ8+HRs3bgTAP1f3s2PHDhQXF2Py5MkAdPNcGd1dwbXh77deFwThgbdjb84x+mzatGmIj4/HH3/80Wgfz1e9Tp06ITY2FsXFxYiMjMQrr7yC6Oho1X6ep3pZWVl47733cPDgQVhaWt7zOJ4vICwsTPXvgYGBCA0NRfv27fHzzz+jT58+AHie7lAqlQgJCcHnn38OAOjevTsSExOxevVqTJo0SXUcz1dj69atQ1hYGDw9PRts16VzxSs3j+DOLIS/V535+fmqCtbd3R3V1dW4devWPY8xNO+++y527tyJI0eOwMvLS7Wd56shCwsLdOjQASEhIYiIiEBwcDC++eYbnqe/OX/+PPLz89GzZ0+YmZnBzMwM0dHRWLVqFczMzFTfl+erMRsbGwQGBuLKlSv8c/U3Hh4e6Nq1a4NtXbp0QWZmJgD+vLqXa9eu4dChQ3jjjTdU23TxXLG4eQS+vr5wd3dHVFSUalt1dTWio6PRt29fAEDPnj1hbm7e4JicnBxcvHhRdYyhEAQB06ZNw7Zt23D48GH4+vo22M/zdX+CIEChUPA8/c3gwYORkJCA2NhY1SMkJAQvvvgiYmNj0a5dO56ve1AoFEhKSoKHhwf/XP3N448/3mipisuXL6Nt27YA+PPqXjZs2ABXV1eEh4ertunkudL4EGUDU1paKsTExAgxMTECAOGrr74SYmJihGvXrgmCUD/9TSaTCdu2bRMSEhKE559/vsnpb15eXsKhQ4eECxcuCE8++aRBThV8++23BZlMJhw9erTBlMGKigrVMTxf9ebPny8cO3ZMyMjIEOLj44UFCxYIJiYmwsGDBwVB4Hl6kLtnSwkCz9cds2bNEo4ePSqkp6cLp0+fFkaNGiXY2dkJV69eFQSB5+luZ8+eFczMzITPPvtMuHLlivCf//xHsLa2FjZt2qQ6huerobq6OqFNmzbCvHnzGu3TtXPF4uYBjhw5IgBo9HjllVcEQaifLrho0SLB3d1dkEqlwoABA4SEhIQG71FZWSlMmzZNcHJyEqysrIRRo0YJmZmZInwb7WrqPAEQNmzYoDqG56vea6+9JrRt21awsLAQXFxchMGDB6sKG0HgeXqQvxc3PF/17qwtYm5uLnh6egrjxo0TEhMTVft5nhratWuXEBAQIEilUqFz587C2rVrG+zn+WrowIEDAgAhJSWl0T5dO1cSQRAEzV8PIiIiIhIHx9wQERGRQWFxQ0RERAaFxQ0REREZFBY3REREZFBY3BAREZFBYXFDREREBoXFDRERERkUFjdERERkUFjcEJHRO3HiBAIDA2Fubo6xY8fecxsR6QcWN0SkE3Jzc/Huu++iXbt2kEql8Pb2xujRo/H777836/VHjx6FRCJBcXGx2p89c+ZMdOvWDRkZGfjpp5/uuY2I9IOZ2AGIiK5evYrHH38cDg4OWL58OYKCglBTU4MDBw5g6tSpSE5O1urnp6WlYcqUKfDy8rrvNiLSD7y3FBGJbuTIkYiPj0dKSgpsbGwa7CsuLkZxcTF8fX0RExODbt26qbY7OjriyJEj8PHxga+vb4PXvfLKK/jpp5+gUCgwZ84cbN68GXK5HCEhIfj666/x2GOP4erVq41et2HDBrz66quNtk2ePFnj35uItINtKSISVVFREfbv34+pU6c2KmwAwMHB4YHv4e3tjcjISABASkoKcnJy8M033wAA5s6di8jISPz888+4cOECOnTogOHDh6OoqAje3t7IycmBvb09Vq5ciZycHEyYMKHRtokTJ2r0OxORdrG4ISJRpaamQhAEdO7c+aHfw9TUFE5OTgAAV1dXuLu7QyaToby8HKtXr8aKFSsQFhaGrl274ocffoCVlRXWrVsHU1NTuLu7QyKRQCaTwd3dHTY2No22WVlZaerrElELYHFDRKK60xmXSCQaf++0tDTU1NTg8ccfV20zNzdHr169kJSUpPHPIyLdwOKGiETl5+cHiURy32LDxKT+R9XdQwRramoe+N73KpwEQdBKMUVEuoHFDRGJysnJCcOHD8e3336L8vLyRvuLi4vh4uICAMjJyVFtj42NbXCchYUFAKCurk61rUOHDrCwsMAff/yh2lZTU4Nz586hS5cumvwaRKRDWNwQkei+++471NXVoVevXoiMjMSVK1eQlJSEVatWITQ0FFZWVujTpw+WLl2KS5cu4dixY/jwww8bvEfbtm0hkUiwe/du3Lx5E2VlZbCxscHbb7+NOXPmYP/+/bh06RLefPNNVFRU4PXXXxfp2xKRtrG4ISLR+fr64sKFCxg0aBBmzZqFgIAADB06FL///jtWr14NAFi/fj1qamoQEhKC9957D59++mmD92jdujUWL16MDz74AG5ubpg2bRoAYOnSpRg/fjxefvll9OjRA6mpqThw4AAcHR1b/HsSUcvgOjdERERkUHjlhoiIiAwKixsiIiIyKCxuiIiIyKCwuCEiIiKDwuKGiIiIDAqLGyIiIjIoLG6IiIjIoLC4ISIiIoPC4oaIiIgMCosbIiIiMigsboiIiMig/D8ISDVnx5b/bAAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "cluster_stream('all_by_all.tsv',cutoffs=[100,200,300,400,500,600,700], score = 3, directory='cluster_stream_res')\n",
    "cluster_progression('cluster_stream_res', filter = 10) #creates the plot"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "the next function creates a table with the cluster IDs and the numbers of members. list files containing uniprot IDs of interest can be added as additional arguments to locate which cluster these proteins are in."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "      cluster  rep_count\n",
      "1           1       2468\n",
      "0           0        443\n",
      "2           2        149\n",
      "7           7          3\n",
      "3           3          2\n",
      "...       ...        ...\n",
      "3155     3155          1\n",
      "3156     3156          1\n",
      "3157     3157          1\n",
      "3158     3158          1\n",
      "9452     9452          1\n",
      "\n",
      "[9453 rows x 2 columns]\n"
     ]
    }
   ],
   "source": [
    "clusters = MakeTable('cluster_stream_res/cluster_400.csv', )\n",
    "print(clusters)\n",
    "clusters.to_csv('cluster_size_table.csv', index=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The next function, ParseClusters will create a zip file containing an Hmm and an alginment file for each cluster with the most representative sequence at the top"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "['0', '2', '1', '4', '3', '5']\n"
     ]
    },
    {
     "ename": "IndexError",
     "evalue": "list index out of range",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mIndexError\u001b[0m                                Traceback (most recent call last)",
      "Cell \u001b[0;32mIn[8], line 10\u001b[0m\n\u001b[1;32m      7\u001b[0m \u001b[38;5;28mprint\u001b[39m(cluster_list)\n\u001b[1;32m      9\u001b[0m \u001b[38;5;66;03m#This function will parse the clusters and write the most representative sequences to a file\u001b[39;00m\n\u001b[0;32m---> 10\u001b[0m \u001b[43mParseClusters\u001b[49m\u001b[43m(\u001b[49m\u001b[43mcluster_file\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[38;5;124;43m'\u001b[39;49m\u001b[38;5;124;43mcluster_stream_res/cluster_400.csv\u001b[39;49m\u001b[38;5;124;43m'\u001b[39;49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mfasta_file\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43m \u001b[49m\u001b[38;5;124;43m'\u001b[39;49m\u001b[38;5;124;43mGH18_input.fasta\u001b[39;49m\u001b[38;5;124;43m'\u001b[39;49m\u001b[43m,\u001b[49m\u001b[43mtarget_clusters\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43mcluster_list\u001b[49m\u001b[43m,\u001b[49m\u001b[43mzip_path\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[38;5;124;43m'\u001b[39;49m\u001b[38;5;124;43mconsensus_cutoff400.zip\u001b[39;49m\u001b[38;5;124;43m'\u001b[39;49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mid_index\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43m \u001b[49m\u001b[38;5;241;43m1\u001b[39;49m\u001b[43m)\u001b[49m\n\u001b[1;32m     11\u001b[0m \u001b[38;5;28mprint\u001b[39m(\u001b[38;5;124mf\u001b[39m\u001b[38;5;124m'\u001b[39m\u001b[38;5;124mClusters collected: \u001b[39m\u001b[38;5;132;01m{\u001b[39;00mclusters\u001b[38;5;132;01m}\u001b[39;00m\u001b[38;5;124m'\u001b[39m)\n",
      "File \u001b[0;32m~/Library/CloudStorage/OneDrive-Nexus365/GitHub/honours-chitinase/ClusterTools/cluster_consensus.py:143\u001b[0m, in \u001b[0;36mParseClusters\u001b[0;34m(cluster_file, fasta_file, target_clusters, zip_path, id_index, num_processes)\u001b[0m\n\u001b[1;32m    141\u001b[0m \u001b[38;5;28;01mdef\u001b[39;00m \u001b[38;5;21mParseClusters\u001b[39m(cluster_file, fasta_file, target_clusters, zip_path, id_index\u001b[38;5;241m=\u001b[39m\u001b[38;5;28;01mNone\u001b[39;00m, num_processes\u001b[38;5;241m=\u001b[39m\u001b[38;5;28;01mNone\u001b[39;00m):\n\u001b[1;32m    142\u001b[0m     clusters \u001b[38;5;241m=\u001b[39m defaultdict(\u001b[38;5;28mlist\u001b[39m)\n\u001b[0;32m--> 143\u001b[0m     fa \u001b[38;5;241m=\u001b[39m \u001b[43mReadFasta\u001b[49m\u001b[43m(\u001b[49m\u001b[43mfasta_file\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mid_index\u001b[49m\u001b[43m)\u001b[49m\n\u001b[1;32m    145\u001b[0m     \u001b[38;5;66;03m# Collect clusters\u001b[39;00m\n\u001b[1;32m    146\u001b[0m     \u001b[38;5;28;01mwith\u001b[39;00m \u001b[38;5;28mopen\u001b[39m(cluster_file) \u001b[38;5;28;01mas\u001b[39;00m f:\n",
      "File \u001b[0;32m~/Library/CloudStorage/OneDrive-Nexus365/GitHub/honours-chitinase/ClusterTools/cluster_consensus.py:25\u001b[0m, in \u001b[0;36mReadFasta\u001b[0;34m(file, id_index)\u001b[0m\n\u001b[1;32m     23\u001b[0m \u001b[38;5;28mid\u001b[39m\u001b[38;5;241m=\u001b[39mrecord\u001b[38;5;241m.\u001b[39mid\n\u001b[1;32m     24\u001b[0m \u001b[38;5;28;01mif\u001b[39;00m id_index \u001b[38;5;129;01mis\u001b[39;00m \u001b[38;5;129;01mnot\u001b[39;00m \u001b[38;5;28;01mNone\u001b[39;00m:\n\u001b[0;32m---> 25\u001b[0m     \u001b[38;5;28mid\u001b[39m\u001b[38;5;241m=\u001b[39m\u001b[38;5;28;43mid\u001b[39;49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43msplit\u001b[49m\u001b[43m(\u001b[49m\u001b[38;5;124;43m'\u001b[39;49m\u001b[38;5;124;43m|\u001b[39;49m\u001b[38;5;124;43m'\u001b[39;49m\u001b[43m)\u001b[49m\u001b[43m[\u001b[49m\u001b[43mid_index\u001b[49m\u001b[43m]\u001b[49m\n\u001b[1;32m     26\u001b[0m      \u001b[38;5;66;03m# Print the extracted ID to check what is being read\u001b[39;00m\n\u001b[1;32m     27\u001b[0m \u001b[38;5;66;03m#print(f'Extracted ID: {id}')\u001b[39;00m\n\u001b[1;32m     28\u001b[0m \u001b[38;5;28;01mif\u001b[39;00m _non_aa_re\u001b[38;5;241m.\u001b[39msearch(\u001b[38;5;28mstr\u001b[39m(record\u001b[38;5;241m.\u001b[39mseq)):\n",
      "\u001b[0;31mIndexError\u001b[0m: list index out of range"
     ]
    }
   ],
   "source": [
    "def clusters_to_string_list(clusters):\n",
    "    # Filter clusters with at least 2 members\n",
    "    filtered_clusters = clusters[clusters['rep_count'] >= 2]  # Assuming 'rep_count' contains the member count\n",
    "    return filtered_clusters['cluster'].astype(str).tolist()\n",
    "\n",
    "cluster_list = clusters_to_string_list(clusters)\n",
    "print(cluster_list)\n",
    "\n",
    "#This function will parse the clusters and write the most representative sequences to a file\n",
    "ParseClusters(cluster_file='cluster_stream_res/cluster_400.csv', fasta_file= 'GH18_input.fasta',target_clusters=cluster_list,zip_path='consensus_cutoff400.zip', id_index= 1)\n",
    "print(f'Clusters collected: {clusters}') #debugging\n",
    "\n",
    "# unzip the consensus to seperate folder\n",
    "#import zipfile\n",
    "#with zipfile.ZipFile('consensus_cutoff500.zip', 'r') as zip_ref:\n",
    "    #zip_ref.extractall('consensus_cutoff500')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The below cell takes clusters of varying bit score trhesholds from JOhn's SSN workflow and compiles representatives into an 'allcluster.fasta' master file. This can be used as the basis of a phylogeny or for experimental anlaysis of the sequnence space"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "---cutoff: 100\n",
      "1 clusters with 12167 sequences\n",
      "349 single nodes added\n",
      "0 unconnected nodes added\n",
      "---cutoff: 200\n",
      "1 clusters with 11004 sequences\n",
      "1512 single nodes added\n",
      "0 unconnected nodes added\n",
      "---cutoff: 300\n",
      "6 clusters with 9797 sequences\n",
      "2719 single nodes added\n",
      "0 unconnected nodes added\n",
      "---cutoff: 400\n",
      "8 clusters with 3071 sequences\n",
      "9445 single nodes added\n",
      "0 unconnected nodes added\n",
      "---cutoff: 500\n",
      "3 clusters with 602 sequences\n",
      "11914 single nodes added\n",
      "0 unconnected nodes added\n",
      "---cutoff: 600\n",
      "1 clusters with 2 sequences\n",
      "12514 single nodes added\n",
      "0 unconnected nodes added\n",
      "3579753 \n",
      "21.51s\n",
      "\n",
      "Warning: Cluster 1 is empty. Skipping.\n",
      "Warning: Cluster 5 is empty. Skipping.\n",
      "Warning: Cluster 0 is empty. Skipping.\n",
      "Warning: Cluster 2 is empty. Skipping.\n",
      "Warning: Cluster 6 is empty. Skipping.\n",
      "Warning: Cluster 4 is empty. Skipping.Warning: Cluster 3 is empty. Skipping.\n",
      "\n",
      "Warning: Cluster 7 is empty. Skipping.\n",
      "Clusters collected:       cluster  rep_count\n",
      "1           1       2468\n",
      "0           0        443\n",
      "2           2        149\n",
      "7           7          3\n",
      "3           3          2\n",
      "...       ...        ...\n",
      "3155     3155          1\n",
      "3156     3156          1\n",
      "3157     3157          1\n",
      "3158     3158          1\n",
      "9452     9452          1\n",
      "\n",
      "[9453 rows x 2 columns]\n",
      "Warning: Cluster 0 is empty. Skipping.\n",
      "Clusters collected:        cluster  rep_count\n",
      "0            0          2\n",
      "8359      8359          1\n",
      "8337      8337          1\n",
      "8338      8338          1\n",
      "8339      8339          1\n",
      "...        ...        ...\n",
      "4174      4174          1\n",
      "4175      4175          1\n",
      "4176      4176          1\n",
      "4177      4177          1\n",
      "12514    12514          1\n",
      "\n",
      "[12515 rows x 2 columns]\n",
      "Warning: Cluster 0 is empty. Skipping.\n",
      "Clusters collected:       cluster  rep_count\n",
      "0           0      11004\n",
      "1005     1005          1\n",
      "1014     1014          1\n",
      "1013     1013          1\n",
      "1012     1012          1\n",
      "...       ...        ...\n",
      "503       503          1\n",
      "502       502          1\n",
      "501       501          1\n",
      "500       500          1\n",
      "1512     1512          1\n",
      "\n",
      "[1513 rows x 2 columns]\n",
      "Warning: Cluster 0 is empty. Skipping.\n",
      "Clusters collected:        cluster  rep_count\n",
      "0            0          2\n",
      "8359      8359          1\n",
      "8337      8337          1\n",
      "8338      8338          1\n",
      "8339      8339          1\n",
      "...        ...        ...\n",
      "4174      4174          1\n",
      "4175      4175          1\n",
      "4176      4176          1\n",
      "4177      4177          1\n",
      "12514    12514          1\n",
      "\n",
      "[12515 rows x 2 columns]\n",
      "Warning: Cluster 2 is empty. Skipping.\n",
      "Warning: Cluster 1 is empty. Skipping.\n",
      "Warning: Cluster 0 is empty. Skipping.\n",
      "Clusters collected:        cluster  rep_count\n",
      "0            0        597\n",
      "2            2          3\n",
      "1            1          2\n",
      "11915    11915          1\n",
      "7947      7947          1\n",
      "...        ...        ...\n",
      "3975      3975          1\n",
      "3976      3976          1\n",
      "3977      3977          1\n",
      "3978      3978          1\n",
      "11916    11916          1\n",
      "\n",
      "[11917 rows x 2 columns]\n",
      "Warning: Cluster 0 is empty. Skipping.\n",
      "Clusters collected:      cluster  rep_count\n",
      "0          0      12167\n",
      "230      230          1\n",
      "238      238          1\n",
      "237      237          1\n",
      "236      236          1\n",
      "..       ...        ...\n",
      "115      115          1\n",
      "114      114          1\n",
      "113      113          1\n",
      "112      112          1\n",
      "349      349          1\n",
      "\n",
      "[350 rows x 2 columns]\n",
      "Warning: Cluster 3 is empty. Skipping.\n",
      "Warning: Cluster 1 is empty. Skipping.\n",
      "Warning: Cluster 5 is empty. Skipping.\n",
      "Warning: Cluster 4 is empty. Skipping.\n",
      "Warning: Cluster 2 is empty. Skipping.\n",
      "Warning: Cluster 0 is empty. Skipping.\n",
      "Clusters collected:       cluster  rep_count\n",
      "0           0       9771\n",
      "2           2         15\n",
      "1           1          4\n",
      "4           4          3\n",
      "3           3          2\n",
      "...       ...        ...\n",
      "912       912          1\n",
      "913       913          1\n",
      "914       914          1\n",
      "915       915          1\n",
      "2724     2724          1\n",
      "\n",
      "[2725 rows x 2 columns]\n",
      "cutoff500\n",
      "Extracted ./extracted_fasta_files/cluster_cutoff500_top_hits.fasta\n",
      "cutoff400\n",
      "Extracted ./extracted_fasta_files/cluster_cutoff400_top_hits.fasta\n",
      "200\n",
      "Extracted ./extracted_fasta_files/cluster_200_top_hits.fasta\n",
      "300\n",
      "Extracted ./extracted_fasta_files/cluster_300_top_hits.fasta\n",
      "100\n",
      "Extracted ./extracted_fasta_files/cluster_100_top_hits.fasta\n",
      "700\n",
      "Extracted ./extracted_fasta_files/cluster_700_top_hits.fasta\n",
      "600\n",
      "Extracted ./extracted_fasta_files/cluster_600_top_hits.fasta\n",
      "500\n",
      "Extracted ./extracted_fasta_files/cluster_500_top_hits.fasta\n",
      "400\n",
      "Extracted ./extracted_fasta_files/cluster_400_top_hits.fasta\n"
     ]
    }
   ],
   "source": [
    "\n",
    "\n",
    "start_cluster = 100\n",
    "end_cluster = 700\n",
    "step = 100\n",
    "cutoffs = range(start_cluster, end_cluster, step)\n",
    "\n",
    "\n",
    "cluster_stream('all_by_all.tsv',cutoffs=(cutoffs), score = 3, directory='cluster_stream_res')\n",
    "\n",
    "directory_path = 'cluster_stream_res/'\n",
    "#create hmms and representative alignments for each cluster at each cutoff\n",
    "for filename in os.listdir(directory_path):\n",
    "    file_path = os.path.join(directory_path, filename)\n",
    "    clusters = MakeTable(file_path, ) \n",
    "    cluster_list = clusters_to_string_list(clusters)\n",
    "    ParseClusters(cluster_file=file_path, fasta_file= 'input_clean.fasta',target_clusters=cluster_list,zip_path=f'{filename}_consensus.zip', id_index= 1)\n",
    "    print(f'Clusters collected: {clusters}') #debugging\n",
    "\n",
    "\n",
    "#extract the fasta file containing the top hits at each cutoff and combine in allcluster.fasta output file \n",
    "# Paths\n",
    "output_directory = './extracted_fasta_files'\n",
    "combined_fasta_path = 'allcluster.fasta'\n",
    "\n",
    "# Ensure the output directory exists\n",
    "os.makedirs(output_directory, exist_ok=True)\n",
    "\n",
    "# Step 1: Extract and rename each fasta file\n",
    "for filename in os.listdir():\n",
    "    if filename.endswith('.zip'):\n",
    "        # Extract the number from the zip filename, e.g., '100' from 'cluster_100.csv_consensus.zip'\n",
    "        cluster_number = filename.split('_')[1].split('.')[0]\n",
    "        print(cluster_number)\n",
    "        with zipfile.ZipFile(filename, 'r') as zip_ref:\n",
    "            for zip_info in zip_ref.infolist():\n",
    "                if zip_info.filename == 'cluster_top_hits.fasta':\n",
    "                    # Save each extracted fasta file with the cluster number in its filename\n",
    "                    extracted_fasta_path = os.path.join(output_directory, f'cluster_{cluster_number}_top_hits.fasta')\n",
    "                    with zip_ref.open(zip_info) as fasta_file, open(extracted_fasta_path, 'wb') as out_fasta:\n",
    "                        out_fasta.write(fasta_file.read())\n",
    "                    print(f\"Extracted {extracted_fasta_path}\")\n",
    "\n",
    "\n",
    "with open(combined_fasta_path, 'w') as combined_fasta:\n",
    "    # Loop through each .fasta file in the output directory\n",
    "    for fasta_file in Path(output_directory).glob('*.fasta'):\n",
    "        with open(fasta_file, 'r') as f:\n",
    "            combined_fasta.write(f.read())  # Write each file's content to `allcluster.fasta`\n",
    "\n",
    " "
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "bio",
   "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.18"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}