dev.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "63e057ac-51aa-4a96-9ede-b76fe7ba7407",
   "metadata": {},
   "source": [
    "# Dependencies"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "251cef2a-37b5-46e6-9fda-4546d1c14315",
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "import argparse\n",
    "import pathlib\n",
    "import sys\n",
    "import subprocess\n",
    "import math\n",
    "import warnings\n",
    "from collections import defaultdict\n",
    "from multiprocessing import Pool, Manager\n",
    "from itertools import repeat\n",
    "warnings.filterwarnings(\"ignore\")\n",
    "\n",
    "import pandas as pd\n",
    "from ete3 import NCBITaxa\n",
    "from Bio import SeqIO"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "3911e5c3-e6b4-41f6-a06b-0fc2337e0e7b",
   "metadata": {},
   "outputs": [],
   "source": [
    "from utils.reference_finder import prepare_reference_genomes\n",
    "from utils.alignment import run_minimap2, run_bwa, sort_samfile, samtools_calculate_coverage\n",
    "from utils.summary import alignment_summary, merge_reference_fasta, call_present_absent\n",
    "from utils.ani import samtools_merged_consensus, ani_summary\n",
    "from utils.input_parsing import parsing_input_f, filter_input_df, get_seq2assembly_dict"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "6941614a-c1d2-46c1-875a-d2e4ba163d61",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "73ec3bc9-0ad1-48c9-bf76-82a2203da268",
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.cluster import AgglomerativeClustering"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "dea0e3d3-1fb4-4fc7-9d0e-7aed5b7b05f5",
   "metadata": {},
   "outputs": [],
   "source": [
    "from ast import literal_eval"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "f7bacebe-e16c-4fd6-858a-8d1677bf830d",
   "metadata": {},
   "outputs": [],
   "source": [
    "def samtools_calculate_coverage(output_dir, include_supp=False):\n",
    "    coverage_files = os.path.join(output_dir, \"coverage_files\")\n",
    "    bam_files = os.path.join(output_dir, \"bam_files\")\n",
    "\n",
    "    command = [\"samtools\",\n",
    "               \"coverage\",\n",
    "               os.path.join(bam_files, f\"merged.sorted.bam\")]\n",
    "    \n",
    "    if include_supp:\n",
    "        # samtools coverage --ff UNMAP,QCFAIL,DUP -q 0 merged.sorted.bam > secondary_coverage.tsv\n",
    "        coverage_file = os.path.join(coverage_files, f\"secondary_coverage.tsv\")\n",
    "        command += ['--ff', 'UNMAP,QCFAIL,DUP', '-q', str(0)]\n",
    "    else:\n",
    "        # samtools coverage merged.sorted.bam -q 20 > primary_coverage.tsv\n",
    "        coverage_file = os.path.join(coverage_files, f\"primary_coverage.tsv\")\n",
    "        command += ['-q', str(20)]\n",
    "        # command += ['--ff', 'UNMAP,QCFAIL,DUP', '-q', str(20)]\n",
    "\n",
    "    subprocess.run(command,\n",
    "                    check=True,\n",
    "                    stdout=open(coverage_file, \"w\"))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "c7fc4610-9595-4d62-aa9b-9a7a9b8e9b06",
   "metadata": {},
   "outputs": [],
   "source": [
    "def parse_fastani_line(line):\n",
    "    ret_list = []\n",
    "    for item in line.strip().split('\\t')[1:]:\n",
    "        if item == 'NA':\n",
    "            ret_list.append(float(0))\n",
    "        else:\n",
    "            ret_list.append(float(item))\n",
    "    return ret_list"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "77f2abbd-c444-4a3b-92c2-890189448434",
   "metadata": {},
   "outputs": [],
   "source": [
    "def find_representative_genome(fastani_path, fastani_assemblies, downloaded_assemblies):\n",
    "    fastani_result = os.path.join(fastani_path, f'pairwise_ani.matrix')\n",
    "    #fastani_assemblies = downloaded_assemblies[downloaded_assemblies['Genus Taxid'] == genus_taxid]['Assembly Accession ID'].values\n",
    "    \n",
    "\n",
    "    ani_matrix = []\n",
    "    with open(fastani_result, 'r') as fastani_out_f:\n",
    "        lines = fastani_out_f.readlines()\n",
    "        num_seqs = int(lines[0].strip())\n",
    "        \n",
    "        for idx, line in enumerate(lines[1:]):\n",
    "            ani_matrix.append(parse_fastani_line(line)+list(np.ones(num_seqs-idx)*100))\n",
    "\n",
    "    dist_nparray = 1-np.array(ani_matrix)/100\n",
    "    dist_df = pd.DataFrame(dist_nparray.T + dist_nparray,\n",
    "                 columns=fastani_assemblies,\n",
    "                 index=fastani_assemblies)\n",
    "\n",
    "    model = AgglomerativeClustering(affinity='precomputed', n_clusters=None, compute_full_tree=True,\n",
    "                                    linkage='complete', \n",
    "                                    distance_threshold=0.05).fit(dist_df)\n",
    "    cluster_df = dist_df.copy()\n",
    "    #print(cluster_df)\n",
    "    cluster_df['Cluster Label'] = model.labels_\n",
    "    \n",
    "    representative_genomes = defaultdict(list)\n",
    "    member2representative = dict()\n",
    "    for cluster_idx in cluster_df['Cluster Label'].unique():\n",
    "        cluster_members = cluster_df[cluster_df['Cluster Label'] == cluster_idx].index.to_list()\n",
    "        selected_df = downloaded_assemblies[downloaded_assemblies['Assembly Accession ID'].isin(cluster_members)].copy()\n",
    "        if selected_df[selected_df['Assembly Level'] == 'Complete Genome'].shape[0] == 1:\n",
    "            representative_genome = selected_df[selected_df['Assembly Level'] == 'Complete Genome']['Assembly Accession ID'].values[0]\n",
    "        elif selected_df[selected_df['Assembly Level'] == 'Complete Genome'].shape[0] > 1:\n",
    "            complete_genomes = list(selected_df[selected_df['Assembly Level'] == 'Complete Genome']['Assembly Accession ID'].values)\n",
    "            representative_genome = dist_df.loc[complete_genomes].sum().idxmin()\n",
    "        else:\n",
    "            representative_genome = dist_df.loc[cluster_members].sum().idxmin()\n",
    "        \n",
    "        representative_genomes[representative_genome] = cluster_members\n",
    "        for member in cluster_members:\n",
    "            member2representative[member] = representative_genome\n",
    "            \n",
    "    return representative_genomes, member2representative"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "158caf83-1fe5-417e-9b09-c93373755ed3",
   "metadata": {},
   "outputs": [],
   "source": [
    "def samtools_merged_consensus(output_directory, threads):\n",
    "    merged_bam = os.path.join(output_directory, 'bam_files', 'merged.sorted.bam')\n",
    "    subprocess.run(['samtools', 'consensus', \n",
    "                    '--show-ins', 'no', \n",
    "                    '--show-del', 'yes',\n",
    "                    '--min-MQ', str(20),\n",
    "                    '-a',\n",
    "                    '--mode', \"simple\",\n",
    "                    '--threads', str(threads),\n",
    "                    merged_bam, \n",
    "                    '-o', os.path.join(output_directory, 'merged_consensus.fasta')],\n",
    "                  check=True)\n",
    "    \n",
    "    consensus_record_dict = SeqIO.to_dict(SeqIO.parse(os.path.join(output_directory, 'merged_consensus.fasta'), \"fasta\"))\n",
    "    return consensus_record_dict"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "33836dc5-d693-4a16-99a4-179ebd87805f",
   "metadata": {},
   "source": [
    "# Full Run"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "d400aad9-4a57-425d-b97b-af99b8217362",
   "metadata": {},
   "outputs": [],
   "source": [
    "def run_magnet(cmd_args):\n",
    "    args = parser.parse_args(cmd_args)\n",
    "\n",
    "    input_tsv = args.classification\n",
    "    input_fastq = args.fastq\n",
    "    input_fastq2 = args.fastq2\n",
    "    mode = args.mode\n",
    "    working_directory = args.output\n",
    "\n",
    "    taxid_col_idx = args.taxid_idx\n",
    "    abundance_col_idx = args.abundance_idx\n",
    "    min_abundance = args.min_abundance\n",
    "    min_mapq = args.min_mapq\n",
    "    min_coverage_score = args.min_covscore\n",
    "    threads = args.threads\n",
    "    valid_kingdom_str = args.kingdom\n",
    "\n",
    "    valid_kingdom = set()\n",
    "    for i in valid_kingdom_str.split(','):\n",
    "        valid_kingdom.add(int(i))\n",
    "\n",
    "    if args.include_mag:\n",
    "        mag_flag = 'all'\n",
    "    else:\n",
    "        mag_flag = 'exclude'\n",
    "\n",
    "    accession_flag = args.accession\n",
    "    call_subspecies = args.subspecies\n",
    "\n",
    "    sep = '\\t'\n",
    "    if str(input_tsv)[-3:] == 'csv':\n",
    "        sep = ','\n",
    "\n",
    "    ncbi_taxa_db = NCBITaxa()\n",
    "    if not os.path.exists(working_directory):\n",
    "        os.mkdir(working_directory)\n",
    "\n",
    "    if accession_flag:\n",
    "        abundance_col_idx = None\n",
    "        min_abundance = 0\n",
    "        input_df, min_abundance = parsing_input_f(input_tsv, sep, taxid_col_idx, abundance_col_idx, min_abundance)\n",
    "        valid_taxids = list(input_df['tax_id'].values)\n",
    "    else:\n",
    "        input_df, min_abundance = parsing_input_f(input_tsv, sep, taxid_col_idx, abundance_col_idx, min_abundance)\n",
    "        # make valid_kingdom a variable?\n",
    "        valid_taxids = filter_input_df(input_df, min_abundance, ncbi_taxa_db, valid_kingdom=valid_kingdom, ret_subspecies=call_subspecies)\n",
    "\n",
    "    reference_metadata = prepare_reference_genomes(valid_taxids, working_directory, ncbi_taxa_db, accession_flag=accession_flag, mag_flag=mag_flag)\n",
    "    downloaded_assemblies = reference_metadata[reference_metadata['Downloaded']]\n",
    "\n",
    "    reference_genome_path = os.path.join(working_directory, 'reference_genomes')\n",
    "\n",
    "    fastani_path = os.path.join(working_directory, 'fastANI')\n",
    "    if not os.path.exists(fastani_path):\n",
    "        os.mkdir(fastani_path)\n",
    "\n",
    "    fastani_assemblies = downloaded_assemblies['Assembly Accession ID'].values\n",
    "    assemblie_list_f = os.path.join(fastani_path, f\"assemblie_list.txt\")\n",
    "    with open(assemblie_list_f, \"w\") as rl_f:\n",
    "        for accession in fastani_assemblies:\n",
    "            reference_genome = os.path.join(reference_genome_path, f'{accession}.fasta')\n",
    "            rl_f.write(f\"{reference_genome}\\n\")\n",
    "\n",
    "    subprocess.run(['fastANI',\n",
    "                    '--rl', assemblie_list_f,\n",
    "                    '--ql', assemblie_list_f,\n",
    "                    \"--matrix\",\n",
    "                    '--threads', str(threads),\n",
    "                    '-o', os.path.join(fastani_path, f'pairwise_ani')],\n",
    "                   check=True,\n",
    "                   stdout=open(os.path.join(fastani_path, \"fastani.log\"), \"a\"),\n",
    "                   stderr=open(os.path.join(fastani_path, \"fastani.err\"), \"a\"))\n",
    "\n",
    "    representative_genomes, member2representative = find_representative_genome(fastani_path, fastani_assemblies, downloaded_assemblies)\n",
    "\n",
    "    representative_labels = []\n",
    "    cluster_members = []\n",
    "    for idx, row in downloaded_assemblies.iterrows():\n",
    "        accession = row['Assembly Accession ID']\n",
    "        if accession in representative_genomes.keys():\n",
    "            representative_labels.append(True)\n",
    "            cluster_members.append(','.join(representative_genomes[accession]))\n",
    "        else:\n",
    "            representative_labels.append(False)\n",
    "            cluster_members.append(','.join(representative_genomes[member2representative[accession]]))\n",
    "    downloaded_assemblies['Cluster Representative'] = representative_labels\n",
    "    downloaded_assemblies['Cluster Members'] = cluster_members\n",
    "\n",
    "    representative_df = downloaded_assemblies[downloaded_assemblies['Cluster Representative']]\n",
    "\n",
    "    seq2assembly_dict = get_seq2assembly_dict(working_directory, representative_df)\n",
    "    reference_fasta = merge_reference_fasta(list(representative_df['Assembly Accession ID']), working_directory)\n",
    "\n",
    "    if mode == 'ont':\n",
    "        aligner_output = run_minimap2(input_fastq, reference_fasta, 'merged', working_directory, threads=threads)\n",
    "    if mode == 'illumina':\n",
    "        aligner_output = run_bowtie2(input_fastq, input_fastq2, reference_fasta, 'merged', working_directory, threads=threads)\n",
    "    sort_samfile('merged', aligner_output, working_directory, min_mapq=0, threads=threads)\n",
    "\n",
    "    coverage_files = os.path.join(working_directory, \"coverage_files\")\n",
    "    if not os.path.exists(coverage_files):\n",
    "        os.mkdir(coverage_files)\n",
    "\n",
    "    pool = Pool(processes=threads)\n",
    "    pool.starmap(samtools_calculate_coverage, zip(repeat(working_directory), [True, False]))\n",
    "    pool.close()\n",
    "    pool.join()\n",
    "\n",
    "\n",
    "    representative_df = alignment_summary(representative_df,\n",
    "                                          working_directory,\n",
    "                                          seq2assembly_dict,\n",
    "                                          include_supp=True)\n",
    "\n",
    "    representative_df = alignment_summary(representative_df,\n",
    "                                          working_directory,\n",
    "                                          seq2assembly_dict,\n",
    "                                          include_supp=False)\n",
    "\n",
    "    consensus_record_dict = samtools_merged_consensus(working_directory, threads)\n",
    "    representative_df = ani_summary(representative_df, consensus_record_dict, working_directory, threads)\n",
    "    representative_df = call_present_absent(representative_df, min_coverage_score)\n",
    "    representative_df.sort_values(['Primary Score'], ascending=False).to_csv(os.path.join(working_directory, \n",
    "                                                                                          'cluster_representative.csv'), \n",
    "                                                                             index=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "83f474b5-c9ae-495e-a99e-728f4310c7b0",
   "metadata": {},
   "source": [
    "# Argparse"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "f12642b9-1fab-4399-88c2-251f0f237b64",
   "metadata": {},
   "outputs": [],
   "source": [
    "parser = argparse.ArgumentParser(description=\"Universal Taxonomic Classification Verifier.\")\n",
    "\n",
    "parser.add_argument(\"-c\", \"--classification\", type=pathlib.Path, required=True, help=\"Path to the Taxonomic Classification Report. Accepting csv/tsv file format, other text formats are treated as tsv.\")\n",
    "parser.add_argument(\"-i\", \"--fastq\", type=pathlib.Path, required=True, help=\"Path to the first fastq file.\")\n",
    "parser.add_argument(\"-I\", \"--fastq2\", type=pathlib.Path, required=False, help=\"Path to the second fastq file for paired-end reads.\")\n",
    "parser.add_argument(\"-m\", \"--mode\", type=str, required=False, choices=['ont', 'illumina'], help=\"Modes for different sequencing platforms [ont, illumina]. Default:[ont]\",  default='ont')\n",
    "parser.add_argument(\"-o\", \"--output\", type=pathlib.Path, required=True, help=\"Path to the output directory.\")\n",
    "parser.add_argument(\"-t\", \"--taxid-idx\", type=int, required=False, help=\"The column index (0-based) of the taxids. Default:[0]\", default=0)\n",
    "parser.add_argument(\"-a\", \"--abundance-idx\", type=int, required=False, help=\"The column index (0-based) of the abundance. Default:[None]\")\n",
    "parser.add_argument(\"--min-abundance\", type=float, required=False, help=\"Minimum abundance (0-1) for pre-filtering, exclude taxa below the threshold.\", default=0)\n",
    "parser.add_argument(\"--min-mapq\", type=int, required=False, help=\"Minimum MAPQ for primary alignments. Default:[20]\", default=20)\n",
    "parser.add_argument(\"--min-covscore\", type=float, required=False, help=\"Minimum Coverage Score for supplementary alignments. Default:[0.7]\", default=0.7)\n",
    "parser.add_argument(\"--threads\", type=int, required=False, help=\"Number of threads for Multi-threading. Default:[1]\", default=1)\n",
    "parser.add_argument(\"--kingdom\", type=str, help=\"A comma separated list of taxids of valid kingdoms. Default:[2,4751,2157,10239]\", default='2,4751,2157,10239')\n",
    "parser.add_argument(\"--include-mag\", action='store_true', required=False, help=\"Include metagenomic assemble genomes. Default:[off]\")\n",
    "parser.set_defaults(include_mag=False)\n",
    "parser.add_argument(\"--subspecies\", action='store_true', required=False, help=\"Verify taxonomic classification at subspecies rank. Default:[off]\")\n",
    "parser.set_defaults(subspecies=False)\n",
    "parser.add_argument(\"--accession\", action='store_true', required=False, help=\"Take accession ids as taxids. Does not work with min-abundance. Default:[off]\")\n",
    "parser.set_defaults(accession=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a80e5d8d-c952-4b5c-ad33-2a56f4d645a8",
   "metadata": {},
   "source": [
    "## Metamap-sim"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "f51c4556-ac9e-4ed5-9790-860262d17c4a",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "-c /home/Users/ns58/Mob-experiments/Lemur-outputs/MetaMaps-sim/nof-i100-relative_abundance.tsv -i /home/Users/ns58/Mob-experiments/Data/Metamaps-sim/simulations_i100_specifiedFrequencies/0/reads.fastq -m ont -t 4 -a 0 --min-abundance 0.01 -o /home/Users/yl181/memu/magnet_v222_MetaMapsSim/nof-i100/nof-i100 --threads 40\n"
     ]
    }
   ],
   "source": [
    "cmd_records = []\n",
    "\n",
    "magnet_output = '/home/Users/yl181/memu/magnet_v222_MetaMapsSim'\n",
    "if not os.path.exists(magnet_output):\n",
    "    os.mkdir(magnet_output)\n",
    "\n",
    "classification_output = '/home/Users/ns58/Mob-experiments/Lemur-outputs/MetaMaps-sim/nof-i100-relative_abundance.tsv'\n",
    "if not os.path.exists(classification_output):\n",
    "    print(\"classification_output missing:\", classification_output)\n",
    "input_fastq = '/home/Users/ns58/Mob-experiments/Data/Metamaps-sim/simulations_i100_specifiedFrequencies/0/reads.fastq'\n",
    "if not os.path.exists(input_fastq):\n",
    "    print(\"input_fastq missing:\", input_fastq)\n",
    "    \n",
    "magnet_output_prefix = 'nof-i100'\n",
    "if not os.path.exists(os.path.join(magnet_output, magnet_output_prefix)):\n",
    "    os.mkdir(os.path.join(magnet_output, magnet_output_prefix))\n",
    "magnet_sample_output = os.path.join(magnet_output, magnet_output_prefix, magnet_output_prefix)\n",
    "\n",
    "cmd_args = ['-c', classification_output,\n",
    "'-i', input_fastq,\n",
    "'-m', 'ont',\n",
    "'-t', str(4),\n",
    "'-a', str(0),\n",
    "'--min-abundance', str(0.01),\n",
    "'-o', magnet_sample_output,\n",
    "'--threads', '40']\n",
    "\n",
    "print(\" \".join(cmd_args))\n",
    "cmd_records.append(cmd_args)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4f6cfdcd-7ea8-4aca-8a8a-81c0b9850601",
   "metadata": {},
   "source": [
    "## Kraken"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "f603d325-7d70-43d9-8526-8dd3879fad9a",
   "metadata": {},
   "outputs": [],
   "source": [
    "mob_output_paths = ['/home/Users/ns58/Mob-experiments/Kraken2-outputs/Zymo-Log-0.10']\n",
    "input_fastq_paths = ['/home/Users/ns58/Mob-experiments/Data-grouped/Zymo-Log-p0.10']\n",
    "\n",
    "magnet_output = '/home/Users/yl181/memu/magnet_v222_ZYMO_Log'\n",
    "if not os.path.exists(magnet_output):\n",
    "    os.mkdir(magnet_output)\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "id": "844377c0-23f3-4c70-b0af-4721c23dbd39",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "-c /home/Users/ns58/Mob-experiments/Kraken2-outputs/Zymo-Log-0.10/Zymo-Log-1.report.txt -i /home/Users/ns58/Mob-experiments/Data-grouped/Zymo-Log-p0.10/Zymo-Log-1.fastq -m ont -t 4 -a 0 --min-abundance 0.01 -o /home/Users/yl181/memu/magnet_v222_ZYMO_Log/Zymo-Log-p0.10-kraken/Zymo-Log-1 --threads 40\n",
      "-c /home/Users/ns58/Mob-experiments/Kraken2-outputs/Zymo-Log-0.10/Zymo-Log-5.report.txt -i /home/Users/ns58/Mob-experiments/Data-grouped/Zymo-Log-p0.10/Zymo-Log-5.fastq -m ont -t 4 -a 0 --min-abundance 0.01 -o /home/Users/yl181/memu/magnet_v222_ZYMO_Log/Zymo-Log-p0.10-kraken/Zymo-Log-5 --threads 40\n",
      "-c /home/Users/ns58/Mob-experiments/Kraken2-outputs/Zymo-Log-0.10/Zymo-Log-2.report.txt -i /home/Users/ns58/Mob-experiments/Data-grouped/Zymo-Log-p0.10/Zymo-Log-2.fastq -m ont -t 4 -a 0 --min-abundance 0.01 -o /home/Users/yl181/memu/magnet_v222_ZYMO_Log/Zymo-Log-p0.10-kraken/Zymo-Log-2 --threads 40\n",
      "-c /home/Users/ns58/Mob-experiments/Kraken2-outputs/Zymo-Log-0.10/Zymo-Log-3.report.txt -i /home/Users/ns58/Mob-experiments/Data-grouped/Zymo-Log-p0.10/Zymo-Log-3.fastq -m ont -t 4 -a 0 --min-abundance 0.01 -o /home/Users/yl181/memu/magnet_v222_ZYMO_Log/Zymo-Log-p0.10-kraken/Zymo-Log-3 --threads 40\n",
      "-c /home/Users/ns58/Mob-experiments/Kraken2-outputs/Zymo-Log-0.10/Zymo-Log-4.report.txt -i /home/Users/ns58/Mob-experiments/Data-grouped/Zymo-Log-p0.10/Zymo-Log-4.fastq -m ont -t 4 -a 0 --min-abundance 0.01 -o /home/Users/yl181/memu/magnet_v222_ZYMO_Log/Zymo-Log-p0.10-kraken/Zymo-Log-4 --threads 40\n"
     ]
    }
   ],
   "source": [
    "cmd_records = []\n",
    "\n",
    "for idx, mob_output_path in enumerate(mob_output_paths):\n",
    "    for f in os.listdir(mob_output_path):\n",
    "        if f.startswith('Zymo-Log-') and 'report' in f:\n",
    "            sample_id = f.split('.')[0]\n",
    "            input_fastq = os.path.join(input_fastq_paths[idx], f'{sample_id}.fastq')\n",
    "            if not os.path.exists(input_fastq):\n",
    "                print(\"input_fastq missing:\", mob_output_path, sample_id)\n",
    "            classification_output = os.path.join(mob_output_path, f)\n",
    "            if not os.path.exists(classification_output):\n",
    "                print(\"classification_output missing:\", input_fastq_paths[idx], sample_id)\n",
    "                \n",
    "            magnet_output_prefix = input_fastq_paths[idx].split('/')[-1]+'-kraken'\n",
    "            if not os.path.exists(os.path.join(magnet_output, magnet_output_prefix)):\n",
    "                os.mkdir(os.path.join(magnet_output, magnet_output_prefix))\n",
    "            magnet_sample_output = os.path.join(magnet_output, magnet_output_prefix, sample_id)\n",
    "            \n",
    "            # print(input_fastq)\n",
    "            # print(classification_output)\n",
    "            # print(magnet_sample_output)\n",
    "            \n",
    "            cmd_args = ['-c', classification_output,\n",
    "            '-i', input_fastq,\n",
    "            '-m', 'ont',\n",
    "            '-t', str(4),\n",
    "            '-a', str(0),\n",
    "            '--min-abundance', str(0.01),\n",
    "            '-o', magnet_sample_output,\n",
    "            '--threads', '40']\n",
    "            \n",
    "            print(\" \".join(cmd_args))\n",
    "            cmd_records.append(cmd_args)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a252fde1-9c00-46d5-ad21-fc3348c40222",
   "metadata": {},
   "source": [
    "## SRR17687125"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "5e1f732e-6810-4ea5-8792-63271c55d994",
   "metadata": {},
   "outputs": [],
   "source": [
    "mob_output_paths = ['/home/Users/ns58/Mob-experiments/Mob-outputs/SRR17687125']\n",
    "input_fastq_paths = ['/home/Users/ns58/Mob-experiments/Data/Gut']\n",
    "\n",
    "magnet_output = '/home/Users/yl181/memu/magnet_v222_SRR17687125'\n",
    "if not os.path.exists(magnet_output):\n",
    "    os.mkdir(magnet_output)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "id": "ce030bfb-a3a7-4d29-8894-1dcc6c8621aa",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "-c /home/Users/ns58/Mob-experiments/Mob-outputs/SRR17687125/SRR17687125-relative_abundance.tsv -i /home/Users/ns58/Mob-experiments/Data/Gut/SRR17687125.fastq -m ont -a 12 --min-abundance 0.0 -o /home/Users/yl181/memu/magnet_v222_SRR17687125/Gut/SRR17687125 --threads 40\n",
      "-c /home/Users/ns58/Mob-experiments/Mob-outputs/SRR17687125/nof-SRR17687125-relative_abundance.tsv -i /home/Users/ns58/Mob-experiments/Data/Gut/SRR17687125.fastq -m ont -a 12 --min-abundance 0.0 -o /home/Users/yl181/memu/magnet_v222_SRR17687125/nof-SRR17687125/SRR17687125 --threads 40\n"
     ]
    }
   ],
   "source": [
    "cmd_records = []\n",
    "\n",
    "for idx, mob_output_path in enumerate(mob_output_paths):\n",
    "    for f in os.listdir(mob_output_path):\n",
    "        if 'relative_abundance.tsv' in f and not 'nof' in f:\n",
    "            sample_id = f.split('-')[0]\n",
    "            input_fastq = os.path.join(input_fastq_paths[idx], f'{sample_id}.fastq')\n",
    "            if not os.path.exists(input_fastq):\n",
    "                print(\"input_fastq missing:\", mob_output_path, sample_id)\n",
    "            classification_output = os.path.join(mob_output_path, f)\n",
    "            if not os.path.exists(classification_output):\n",
    "                print(\"classification_output missing:\", input_fastq_paths[idx], sample_id)\n",
    "                \n",
    "            magnet_output_prefix = input_fastq_paths[idx].split('/')[-1]\n",
    "            if not os.path.exists(os.path.join(magnet_output, magnet_output_prefix)):\n",
    "                os.mkdir(os.path.join(magnet_output, magnet_output_prefix))\n",
    "            magnet_sample_output = os.path.join(magnet_output, magnet_output_prefix, sample_id)\n",
    "            \n",
    "            cmd_args = ['-c', classification_output,\n",
    "            '-i', input_fastq,\n",
    "            '-m', 'ont',\n",
    "            '-a', str(12),\n",
    "            '--min-abundance', str(0.000),\n",
    "            '-o', magnet_sample_output,\n",
    "            '--threads', '40']\n",
    "            \n",
    "            #print(sample_id)\n",
    "            print(\" \".join(cmd_args))\n",
    "            cmd_records.append(cmd_args)\n",
    "            \n",
    "for idx, mob_output_path in enumerate(mob_output_paths):\n",
    "    for f in os.listdir(mob_output_path):\n",
    "        if 'relative_abundance.tsv' in f and 'nof' in f:\n",
    "            sample_id = f.split('-')[1]\n",
    "            input_fastq = os.path.join(input_fastq_paths[idx], f'{sample_id}.fastq')\n",
    "            if not os.path.exists(input_fastq):\n",
    "                print(\"input_fastq missing:\", mob_output_path, sample_id)\n",
    "            classification_output = os.path.join(mob_output_path, f)\n",
    "            if not os.path.exists(classification_output):\n",
    "                print(\"classification_output missing:\", input_fastq_paths[idx], sample_id)\n",
    "                \n",
    "            magnet_output_prefix = '-'.join(f.split('-')[0:2])\n",
    "            if not os.path.exists(os.path.join(magnet_output, magnet_output_prefix)):\n",
    "                os.mkdir(os.path.join(magnet_output, magnet_output_prefix))\n",
    "            magnet_sample_output = os.path.join(magnet_output, magnet_output_prefix, sample_id)\n",
    "            \n",
    "            cmd_args = ['-c', classification_output,\n",
    "            '-i', input_fastq,\n",
    "            '-m', 'ont',\n",
    "            '-a', str(12),\n",
    "            '--min-abundance', str(0.000),\n",
    "            '-o', magnet_sample_output,\n",
    "            '--threads', '40']\n",
    "            \n",
    "            #print(sample_id)\n",
    "            print(\" \".join(cmd_args))\n",
    "            cmd_records.append(cmd_args)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5cec7a03-36d6-44ab-92df-efbde100f6ac",
   "metadata": {},
   "source": [
    "## Zymo TM"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "dc88addd-0c02-4cad-8521-b86d799aa230",
   "metadata": {},
   "outputs": [],
   "source": [
    "mob_output_paths = ['/home/Users/ns58/Mob-experiments/Mob-outputs/Zymo-TM']\n",
    "input_fastq_paths = ['/home/Users/ns58/Mob-experiments/Data-grouped/Zymo-TM']\n",
    "\n",
    "magnet_output = '/home/Users/yl181/memu/magnet_v222_ZYMO_TM'\n",
    "if not os.path.exists(magnet_output):\n",
    "    os.mkdir(magnet_output)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "f78a1b95-04f7-4988-a71a-2a219e6a481b",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "-c /home/Users/ns58/Mob-experiments/Mob-outputs/Zymo-TM/Zymo-TM-3-relative_abundance.tsv -i /home/Users/ns58/Mob-experiments/Data-grouped/Zymo-TM/Zymo-TM-3.fastq -m ont -a 12 --min-abundance 0.0 -o /home/Users/yl181/memu/magnet_v222_ZYMO_TM/Zymo-TM/Zymo-TM-3 --threads 40\n",
      "-c /home/Users/ns58/Mob-experiments/Mob-outputs/Zymo-TM/Zymo-TM-1-relative_abundance.tsv -i /home/Users/ns58/Mob-experiments/Data-grouped/Zymo-TM/Zymo-TM-1.fastq -m ont -a 12 --min-abundance 0.0 -o /home/Users/yl181/memu/magnet_v222_ZYMO_TM/Zymo-TM/Zymo-TM-1 --threads 40\n",
      "-c /home/Users/ns58/Mob-experiments/Mob-outputs/Zymo-TM/Zymo-TM-2-relative_abundance.tsv -i /home/Users/ns58/Mob-experiments/Data-grouped/Zymo-TM/Zymo-TM-2.fastq -m ont -a 12 --min-abundance 0.0 -o /home/Users/yl181/memu/magnet_v222_ZYMO_TM/Zymo-TM/Zymo-TM-2 --threads 40\n",
      "-c /home/Users/ns58/Mob-experiments/Mob-outputs/Zymo-TM/Zymo-TM-5-relative_abundance.tsv -i /home/Users/ns58/Mob-experiments/Data-grouped/Zymo-TM/Zymo-TM-5.fastq -m ont -a 12 --min-abundance 0.0 -o /home/Users/yl181/memu/magnet_v222_ZYMO_TM/Zymo-TM/Zymo-TM-5 --threads 40\n",
      "-c /home/Users/ns58/Mob-experiments/Mob-outputs/Zymo-TM/Zymo-TM-6-relative_abundance.tsv -i /home/Users/ns58/Mob-experiments/Data-grouped/Zymo-TM/Zymo-TM-6.fastq -m ont -a 12 --min-abundance 0.0 -o /home/Users/yl181/memu/magnet_v222_ZYMO_TM/Zymo-TM/Zymo-TM-6 --threads 40\n",
      "-c /home/Users/ns58/Mob-experiments/Mob-outputs/Zymo-TM/Zymo-TM-4-relative_abundance.tsv -i /home/Users/ns58/Mob-experiments/Data-grouped/Zymo-TM/Zymo-TM-4.fastq -m ont -a 12 --min-abundance 0.0 -o /home/Users/yl181/memu/magnet_v222_ZYMO_TM/Zymo-TM/Zymo-TM-4 --threads 40\n"
     ]
    }
   ],
   "source": [
    "cmd_records = []\n",
    "\n",
    "for idx, mob_output_path in enumerate(mob_output_paths):\n",
    "    for f in os.listdir(mob_output_path):\n",
    "        if 'relative_abundance.tsv' in f and not 'nof' in f:\n",
    "            sample_id = '-'.join(f.split('-')[0:3])\n",
    "            input_fastq = os.path.join(input_fastq_paths[idx], f'{sample_id}.fastq')\n",
    "            if not os.path.exists(input_fastq):\n",
    "                print(\"input_fastq missing:\", mob_output_path, sample_id)\n",
    "            classification_output = os.path.join(mob_output_path, f)\n",
    "            if not os.path.exists(classification_output):\n",
    "                print(\"classification_output missing:\", input_fastq_paths[idx], sample_id)\n",
    "                \n",
    "            magnet_output_prefix = input_fastq_paths[idx].split('/')[-1]\n",
    "            if not os.path.exists(os.path.join(magnet_output, magnet_output_prefix)):\n",
    "                os.mkdir(os.path.join(magnet_output, magnet_output_prefix))\n",
    "            magnet_sample_output = os.path.join(magnet_output, magnet_output_prefix, sample_id)\n",
    "            \n",
    "            cmd_args = ['-c', classification_output,\n",
    "            '-i', input_fastq,\n",
    "            '-m', 'ont',\n",
    "            '-a', str(12),\n",
    "            '--min-abundance', str(0.000),\n",
    "            '-o', magnet_sample_output,\n",
    "            '--threads', '40']\n",
    "            \n",
    "            #print(sample_id)\n",
    "            print(\" \".join(cmd_args))\n",
    "            cmd_records.append(cmd_args)\n",
    "            "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "id": "f5c178be-cfbb-4746-822a-664fdf2d1843",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "-c /home/Users/ns58/Mob-experiments/Mob-outputs/Zymo-TM/nof-Zymo-TM-1-relative_abundance.tsv -i /home/Users/ns58/Mob-experiments/Data-grouped/Zymo-TM/Zymo-TM-1.fastq -m ont -a 12 --min-abundance 0.0 -o /home/Users/yl181/memu/magnet_v222_ZYMO_TM/nof-Zymo-TM/Zymo-TM-1 --threads 40\n",
      "-c /home/Users/ns58/Mob-experiments/Mob-outputs/Zymo-TM/nof-Zymo-TM-2-relative_abundance.tsv -i /home/Users/ns58/Mob-experiments/Data-grouped/Zymo-TM/Zymo-TM-2.fastq -m ont -a 12 --min-abundance 0.0 -o /home/Users/yl181/memu/magnet_v222_ZYMO_TM/nof-Zymo-TM/Zymo-TM-2 --threads 40\n",
      "-c /home/Users/ns58/Mob-experiments/Mob-outputs/Zymo-TM/nof-Zymo-TM-4-relative_abundance.tsv -i /home/Users/ns58/Mob-experiments/Data-grouped/Zymo-TM/Zymo-TM-4.fastq -m ont -a 12 --min-abundance 0.0 -o /home/Users/yl181/memu/magnet_v222_ZYMO_TM/nof-Zymo-TM/Zymo-TM-4 --threads 40\n",
      "-c /home/Users/ns58/Mob-experiments/Mob-outputs/Zymo-TM/nof-Zymo-TM-5-relative_abundance.tsv -i /home/Users/ns58/Mob-experiments/Data-grouped/Zymo-TM/Zymo-TM-5.fastq -m ont -a 12 --min-abundance 0.0 -o /home/Users/yl181/memu/magnet_v222_ZYMO_TM/nof-Zymo-TM/Zymo-TM-5 --threads 40\n",
      "-c /home/Users/ns58/Mob-experiments/Mob-outputs/Zymo-TM/nof-Zymo-TM-3-relative_abundance.tsv -i /home/Users/ns58/Mob-experiments/Data-grouped/Zymo-TM/Zymo-TM-3.fastq -m ont -a 12 --min-abundance 0.0 -o /home/Users/yl181/memu/magnet_v222_ZYMO_TM/nof-Zymo-TM/Zymo-TM-3 --threads 40\n",
      "-c /home/Users/ns58/Mob-experiments/Mob-outputs/Zymo-TM/nof-Zymo-TM-6-relative_abundance.tsv -i /home/Users/ns58/Mob-experiments/Data-grouped/Zymo-TM/Zymo-TM-6.fastq -m ont -a 12 --min-abundance 0.0 -o /home/Users/yl181/memu/magnet_v222_ZYMO_TM/nof-Zymo-TM/Zymo-TM-6 --threads 40\n"
     ]
    }
   ],
   "source": [
    "for idx, mob_output_path in enumerate(mob_output_paths):\n",
    "    for f in os.listdir(mob_output_path):\n",
    "        if 'relative_abundance.tsv' in f and 'nof' in f:\n",
    "            sample_id = '-'.join(f.split('-')[1:4])\n",
    "            input_fastq = os.path.join(input_fastq_paths[idx], f'{sample_id}.fastq')\n",
    "            if not os.path.exists(input_fastq):\n",
    "                print(\"input_fastq missing:\", mob_output_path, sample_id)\n",
    "            classification_output = os.path.join(mob_output_path, f)\n",
    "            if not os.path.exists(classification_output):\n",
    "                print(\"classification_output missing:\", input_fastq_paths[idx], sample_id)\n",
    "                \n",
    "            magnet_output_prefix = '-'.join(f.split('-')[0:3])\n",
    "            if not os.path.exists(os.path.join(magnet_output, magnet_output_prefix)):\n",
    "                os.mkdir(os.path.join(magnet_output, magnet_output_prefix))\n",
    "            magnet_sample_output = os.path.join(magnet_output, magnet_output_prefix, sample_id)\n",
    "            \n",
    "            cmd_args = ['-c', classification_output,\n",
    "            '-i', input_fastq,\n",
    "            '-m', 'ont',\n",
    "            '-a', str(12),\n",
    "            '--min-abundance', str(0.000),\n",
    "            '-o', magnet_sample_output,\n",
    "            '--threads', '40']\n",
    "            \n",
    "            #print(sample_id)\n",
    "            print(\" \".join(cmd_args))\n",
    "            cmd_records.append(cmd_args)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2fffed63-4279-476c-a46f-1e566fbde6a9",
   "metadata": {},
   "source": [
    "## Simulations"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "397b1b19-eae6-4edf-8306-93c5d1852c07",
   "metadata": {},
   "outputs": [],
   "source": [
    "mob_output_paths = ['/home/Users/ns58/Mob-experiments/Mob-outputs/Sim-Even-s20', '/home/Users/ns58/Mob-experiments/Mob-outputs/Sim-Log-s12']\n",
    "input_fastq_paths = ['/home/Users/ns58/Mob-experiments/Data-grouped/Sim-Even-s20', '/home/Users/ns58/Mob-experiments/Data-grouped/Sim-Log-s12']\n",
    "\n",
    "magnet_output = '/home/Users/yl181/memu/magnet_v222_simulations'\n",
    "if not os.path.exists(magnet_output):\n",
    "    os.mkdir(magnet_output)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "73677f03-cef2-46f7-a05e-c1599587dd52",
   "metadata": {},
   "outputs": [],
   "source": [
    "mob_output_paths = ['/home/Users/ns58/Mob-experiments/Mob-outputs/Sim-MAG-Even-s20']\n",
    "input_fastq_paths = ['/home/Users/ns58/Mob-experiments/Data-grouped/Sim-MAG-Even-s20']\n",
    "\n",
    "magnet_output = '/home/Users/yl181/memu/magnet_v222_simulations_mags'\n",
    "if not os.path.exists(magnet_output):\n",
    "    os.mkdir(magnet_output)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "0655c34b-0c8f-4b2c-8da5-603fab736dde",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "-c /home/Users/ns58/Mob-experiments/Mob-outputs/Sim-MAG-Even-s20/Sim-Even-4-relative_abundance.tsv -i /home/Users/ns58/Mob-experiments/Data-grouped/Sim-MAG-Even-s20/Sim-Even-4.fastq -m ont -a 12 --min-abundance 0.0 -o /home/Users/yl181/memu/magnet_v222_simulations_mags/Sim-MAG-Even-s20/Sim-Even-4 --threads 40\n",
      "-c /home/Users/ns58/Mob-experiments/Mob-outputs/Sim-MAG-Even-s20/Sim-Even-3-relative_abundance.tsv -i /home/Users/ns58/Mob-experiments/Data-grouped/Sim-MAG-Even-s20/Sim-Even-3.fastq -m ont -a 12 --min-abundance 0.0 -o /home/Users/yl181/memu/magnet_v222_simulations_mags/Sim-MAG-Even-s20/Sim-Even-3 --threads 40\n",
      "-c /home/Users/ns58/Mob-experiments/Mob-outputs/Sim-MAG-Even-s20/Sim-Even-1-relative_abundance.tsv -i /home/Users/ns58/Mob-experiments/Data-grouped/Sim-MAG-Even-s20/Sim-Even-1.fastq -m ont -a 12 --min-abundance 0.0 -o /home/Users/yl181/memu/magnet_v222_simulations_mags/Sim-MAG-Even-s20/Sim-Even-1 --threads 40\n",
      "-c /home/Users/ns58/Mob-experiments/Mob-outputs/Sim-MAG-Even-s20/Sim-Even-5-relative_abundance.tsv -i /home/Users/ns58/Mob-experiments/Data-grouped/Sim-MAG-Even-s20/Sim-Even-5.fastq -m ont -a 12 --min-abundance 0.0 -o /home/Users/yl181/memu/magnet_v222_simulations_mags/Sim-MAG-Even-s20/Sim-Even-5 --threads 40\n",
      "-c /home/Users/ns58/Mob-experiments/Mob-outputs/Sim-MAG-Even-s20/Sim-Even-2-relative_abundance.tsv -i /home/Users/ns58/Mob-experiments/Data-grouped/Sim-MAG-Even-s20/Sim-Even-2.fastq -m ont -a 12 --min-abundance 0.0 -o /home/Users/yl181/memu/magnet_v222_simulations_mags/Sim-MAG-Even-s20/Sim-Even-2 --threads 40\n"
     ]
    }
   ],
   "source": [
    "cmd_records = []\n",
    "\n",
    "for idx, mob_output_path in enumerate(mob_output_paths):\n",
    "    for f in os.listdir(mob_output_path):\n",
    "        if 'relative_abundance.tsv' in f:\n",
    "            sample_id = '-'.join(f.split('-')[0:3])\n",
    "            input_fastq = os.path.join(input_fastq_paths[idx], f'{sample_id}.fastq')\n",
    "            if not os.path.exists(input_fastq):\n",
    "                print(\"input_fastq missing:\", mob_output_path, sample_id)\n",
    "            classification_output = os.path.join(mob_output_path, f)\n",
    "            if not os.path.exists(classification_output):\n",
    "                print(\"classification_output missing:\", input_fastq_paths[idx], sample_id)\n",
    "                \n",
    "            magnet_output_prefix = input_fastq_paths[idx].split('/')[-1]\n",
    "            if not os.path.exists(os.path.join(magnet_output, magnet_output_prefix)):\n",
    "                os.mkdir(os.path.join(magnet_output, magnet_output_prefix))\n",
    "            magnet_sample_output = os.path.join(magnet_output, magnet_output_prefix, sample_id)\n",
    "            \n",
    "            cmd_args = ['-c', classification_output,\n",
    "            '-i', input_fastq,\n",
    "            '-m', 'ont',\n",
    "            '-a', str(12),\n",
    "            '--min-abundance', str(0.000),\n",
    "            '-o', magnet_sample_output,\n",
    "            '--threads', '40']\n",
    "            \n",
    "            #print(sample_id)\n",
    "            print(\" \".join(cmd_args))\n",
    "            cmd_records.append(cmd_args)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4c43c889-7d1f-4081-8d22-5b901b673bd8",
   "metadata": {},
   "source": [
    "## Zymo Log"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "3054ed3a-1a3d-4b3b-a50e-a35f86c6f6e9",
   "metadata": {},
   "outputs": [],
   "source": [
    "mob_output_paths = ['/home/Users/ns58/Mob-experiments/Mob-outputs/Zymo-Log-0.10', '/home/Users/ns58/Mob-experiments/Mob-outputs/Zymo-Log-0.25']\n",
    "input_fastq_paths = ['/home/Users/ns58/Mob-experiments/Data-grouped/Zymo-Log-p0.10', '/home/Users/ns58/Mob-experiments/Data-grouped/Zymo-Log-p0.25']\n",
    "\n",
    "magnet_output = '/home/Users/yl181/memu/magnet_v222_ZYMO_Log'\n",
    "if not os.path.exists(magnet_output):\n",
    "    os.mkdir(magnet_output)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "584da22f-bc63-4493-89bb-8cdc7e7978dc",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "-c /home/Users/ns58/Mob-experiments/Mob-outputs/Zymo-Log-0.10/Zymo-Log-2-relative_abundance.tsv -i /home/Users/ns58/Mob-experiments/Data-grouped/Zymo-Log-p0.10/Zymo-Log-2.fastq -m ont -a 12 --min-abundance 0.0 -o /home/Users/yl181/memu/magnet_v222_ZYMO_Log/Zymo-Log-p0.10/Zymo-Log-2 --threads 40\n",
      "-c /home/Users/ns58/Mob-experiments/Mob-outputs/Zymo-Log-0.10/Zymo-Log-5-relative_abundance.tsv -i /home/Users/ns58/Mob-experiments/Data-grouped/Zymo-Log-p0.10/Zymo-Log-5.fastq -m ont -a 12 --min-abundance 0.0 -o /home/Users/yl181/memu/magnet_v222_ZYMO_Log/Zymo-Log-p0.10/Zymo-Log-5 --threads 40\n",
      "-c /home/Users/ns58/Mob-experiments/Mob-outputs/Zymo-Log-0.10/Zymo-Log-1-relative_abundance.tsv -i /home/Users/ns58/Mob-experiments/Data-grouped/Zymo-Log-p0.10/Zymo-Log-1.fastq -m ont -a 12 --min-abundance 0.0 -o /home/Users/yl181/memu/magnet_v222_ZYMO_Log/Zymo-Log-p0.10/Zymo-Log-1 --threads 40\n",
      "-c /home/Users/ns58/Mob-experiments/Mob-outputs/Zymo-Log-0.10/Zymo-Log-4-relative_abundance.tsv -i /home/Users/ns58/Mob-experiments/Data-grouped/Zymo-Log-p0.10/Zymo-Log-4.fastq -m ont -a 12 --min-abundance 0.0 -o /home/Users/yl181/memu/magnet_v222_ZYMO_Log/Zymo-Log-p0.10/Zymo-Log-4 --threads 40\n",
      "-c /home/Users/ns58/Mob-experiments/Mob-outputs/Zymo-Log-0.10/Zymo-Log-3-relative_abundance.tsv -i /home/Users/ns58/Mob-experiments/Data-grouped/Zymo-Log-p0.10/Zymo-Log-3.fastq -m ont -a 12 --min-abundance 0.0 -o /home/Users/yl181/memu/magnet_v222_ZYMO_Log/Zymo-Log-p0.10/Zymo-Log-3 --threads 40\n",
      "-c /home/Users/ns58/Mob-experiments/Mob-outputs/Zymo-Log-0.25/Zymo-Log-2-relative_abundance.tsv -i /home/Users/ns58/Mob-experiments/Data-grouped/Zymo-Log-p0.25/Zymo-Log-2.fastq -m ont -a 12 --min-abundance 0.0 -o /home/Users/yl181/memu/magnet_v222_ZYMO_Log/Zymo-Log-p0.25/Zymo-Log-2 --threads 40\n",
      "-c /home/Users/ns58/Mob-experiments/Mob-outputs/Zymo-Log-0.25/Zymo-Log-5-relative_abundance.tsv -i /home/Users/ns58/Mob-experiments/Data-grouped/Zymo-Log-p0.25/Zymo-Log-5.fastq -m ont -a 12 --min-abundance 0.0 -o /home/Users/yl181/memu/magnet_v222_ZYMO_Log/Zymo-Log-p0.25/Zymo-Log-5 --threads 40\n",
      "-c /home/Users/ns58/Mob-experiments/Mob-outputs/Zymo-Log-0.25/Zymo-Log-1-relative_abundance.tsv -i /home/Users/ns58/Mob-experiments/Data-grouped/Zymo-Log-p0.25/Zymo-Log-1.fastq -m ont -a 12 --min-abundance 0.0 -o /home/Users/yl181/memu/magnet_v222_ZYMO_Log/Zymo-Log-p0.25/Zymo-Log-1 --threads 40\n",
      "-c /home/Users/ns58/Mob-experiments/Mob-outputs/Zymo-Log-0.25/Zymo-Log-4-relative_abundance.tsv -i /home/Users/ns58/Mob-experiments/Data-grouped/Zymo-Log-p0.25/Zymo-Log-4.fastq -m ont -a 12 --min-abundance 0.0 -o /home/Users/yl181/memu/magnet_v222_ZYMO_Log/Zymo-Log-p0.25/Zymo-Log-4 --threads 40\n",
      "-c /home/Users/ns58/Mob-experiments/Mob-outputs/Zymo-Log-0.25/Zymo-Log-3-relative_abundance.tsv -i /home/Users/ns58/Mob-experiments/Data-grouped/Zymo-Log-p0.25/Zymo-Log-3.fastq -m ont -a 12 --min-abundance 0.0 -o /home/Users/yl181/memu/magnet_v222_ZYMO_Log/Zymo-Log-p0.25/Zymo-Log-3 --threads 40\n"
     ]
    }
   ],
   "source": [
    "cmd_records = []\n",
    "\n",
    "for idx, mob_output_path in enumerate(mob_output_paths):\n",
    "    for f in os.listdir(mob_output_path):\n",
    "        if f.startswith('Zymo-Log-') and 'relative_abundance.tsv' in f:\n",
    "            sample_id = '-'.join(f.split('-')[0:3])\n",
    "            input_fastq = os.path.join(input_fastq_paths[idx], f'{sample_id}.fastq')\n",
    "            if not os.path.exists(input_fastq):\n",
    "                print(\"input_fastq missing:\", mob_output_path, sample_id)\n",
    "            classification_output = os.path.join(mob_output_path, f)\n",
    "            if not os.path.exists(classification_output):\n",
    "                print(\"classification_output missing:\", input_fastq_paths[idx], sample_id)\n",
    "                \n",
    "            magnet_output_prefix = input_fastq_paths[idx].split('/')[-1]\n",
    "            if not os.path.exists(os.path.join(magnet_output, magnet_output_prefix)):\n",
    "                os.mkdir(os.path.join(magnet_output, magnet_output_prefix))\n",
    "            magnet_sample_output = os.path.join(magnet_output, magnet_output_prefix, sample_id)\n",
    "            \n",
    "            # print(input_fastq)\n",
    "            # print(classification_output)\n",
    "            # print(magnet_sample_output)\n",
    "            \n",
    "            cmd_args = ['-c', classification_output,\n",
    "            '-i', input_fastq,\n",
    "            '-m', 'ont',\n",
    "            '-a', str(12),\n",
    "            '--min-abundance', str(0.000),\n",
    "            '-o', magnet_sample_output,\n",
    "            '--threads', '40']\n",
    "            \n",
    "            print(\" \".join(cmd_args))\n",
    "            cmd_records.append(cmd_args)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d2100ec0-64bc-4604-8dd7-df9a9df48adc",
   "metadata": {},
   "source": [
    "## Cheese"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "id": "08a784c5-fc3e-402c-9b34-295300dbc677",
   "metadata": {},
   "outputs": [],
   "source": [
    "mob_output_paths = ['/home/Users/ns58/Mob-experiments/Mob-outputs/Cheese']\n",
    "input_fastq_paths = ['/home/Users/ns58/Mob-experiments/Data-grouped/Cheese']\n",
    "\n",
    "magnet_output = '/home/Users/yl181/memu/magnet_v222_cheese'\n",
    "if not os.path.exists(magnet_output):\n",
    "    os.mkdir(magnet_output)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "id": "55bb62d9-bfad-49d1-80e3-169b1c1e891b",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "-c /home/Users/ns58/Mob-experiments/Mob-outputs/Cheese/Cheese-4-relative_abundance.tsv -i /home/Users/ns58/Mob-experiments/Data-grouped/Cheese/Cheese-4.fastq -m ont -a 12 --min-abundance 0.0 -o /home/Users/yl181/memu/magnet_v222_cheese/Cheese/Cheese-4 --threads 40\n",
      "-c /home/Users/ns58/Mob-experiments/Mob-outputs/Cheese/Cheese-2-relative_abundance.tsv -i /home/Users/ns58/Mob-experiments/Data-grouped/Cheese/Cheese-2.fastq -m ont -a 12 --min-abundance 0.0 -o /home/Users/yl181/memu/magnet_v222_cheese/Cheese/Cheese-2 --threads 40\n",
      "-c /home/Users/ns58/Mob-experiments/Mob-outputs/Cheese/Cheese-1-relative_abundance.tsv -i /home/Users/ns58/Mob-experiments/Data-grouped/Cheese/Cheese-1.fastq -m ont -a 12 --min-abundance 0.0 -o /home/Users/yl181/memu/magnet_v222_cheese/Cheese/Cheese-1 --threads 40\n",
      "-c /home/Users/ns58/Mob-experiments/Mob-outputs/Cheese/Cheese-8-relative_abundance.tsv -i /home/Users/ns58/Mob-experiments/Data-grouped/Cheese/Cheese-8.fastq -m ont -a 12 --min-abundance 0.0 -o /home/Users/yl181/memu/magnet_v222_cheese/Cheese/Cheese-8 --threads 40\n",
      "-c /home/Users/ns58/Mob-experiments/Mob-outputs/Cheese/Cheese-5-relative_abundance.tsv -i /home/Users/ns58/Mob-experiments/Data-grouped/Cheese/Cheese-5.fastq -m ont -a 12 --min-abundance 0.0 -o /home/Users/yl181/memu/magnet_v222_cheese/Cheese/Cheese-5 --threads 40\n",
      "-c /home/Users/ns58/Mob-experiments/Mob-outputs/Cheese/Cheese-14-relative_abundance.tsv -i /home/Users/ns58/Mob-experiments/Data-grouped/Cheese/Cheese-14.fastq -m ont -a 12 --min-abundance 0.0 -o /home/Users/yl181/memu/magnet_v222_cheese/Cheese/Cheese-14 --threads 40\n",
      "-c /home/Users/ns58/Mob-experiments/Mob-outputs/Cheese/Cheese-3-relative_abundance.tsv -i /home/Users/ns58/Mob-experiments/Data-grouped/Cheese/Cheese-3.fastq -m ont -a 12 --min-abundance 0.0 -o /home/Users/yl181/memu/magnet_v222_cheese/Cheese/Cheese-3 --threads 40\n",
      "-c /home/Users/ns58/Mob-experiments/Mob-outputs/Cheese/Cheese-7-relative_abundance.tsv -i /home/Users/ns58/Mob-experiments/Data-grouped/Cheese/Cheese-7.fastq -m ont -a 12 --min-abundance 0.0 -o /home/Users/yl181/memu/magnet_v222_cheese/Cheese/Cheese-7 --threads 40\n",
      "-c /home/Users/ns58/Mob-experiments/Mob-outputs/Cheese/Cheese-13-relative_abundance.tsv -i /home/Users/ns58/Mob-experiments/Data-grouped/Cheese/Cheese-13.fastq -m ont -a 12 --min-abundance 0.0 -o /home/Users/yl181/memu/magnet_v222_cheese/Cheese/Cheese-13 --threads 40\n",
      "-c /home/Users/ns58/Mob-experiments/Mob-outputs/Cheese/Cheese-9-relative_abundance.tsv -i /home/Users/ns58/Mob-experiments/Data-grouped/Cheese/Cheese-9.fastq -m ont -a 12 --min-abundance 0.0 -o /home/Users/yl181/memu/magnet_v222_cheese/Cheese/Cheese-9 --threads 40\n",
      "-c /home/Users/ns58/Mob-experiments/Mob-outputs/Cheese/Cheese-15-relative_abundance.tsv -i /home/Users/ns58/Mob-experiments/Data-grouped/Cheese/Cheese-15.fastq -m ont -a 12 --min-abundance 0.0 -o /home/Users/yl181/memu/magnet_v222_cheese/Cheese/Cheese-15 --threads 40\n",
      "-c /home/Users/ns58/Mob-experiments/Mob-outputs/Cheese/Cheese-6-relative_abundance.tsv -i /home/Users/ns58/Mob-experiments/Data-grouped/Cheese/Cheese-6.fastq -m ont -a 12 --min-abundance 0.0 -o /home/Users/yl181/memu/magnet_v222_cheese/Cheese/Cheese-6 --threads 40\n",
      "-c /home/Users/ns58/Mob-experiments/Mob-outputs/Cheese/Cheese-12-relative_abundance.tsv -i /home/Users/ns58/Mob-experiments/Data-grouped/Cheese/Cheese-12.fastq -m ont -a 12 --min-abundance 0.0 -o /home/Users/yl181/memu/magnet_v222_cheese/Cheese/Cheese-12 --threads 40\n",
      "-c /home/Users/ns58/Mob-experiments/Mob-outputs/Cheese/Cheese-10-relative_abundance.tsv -i /home/Users/ns58/Mob-experiments/Data-grouped/Cheese/Cheese-10.fastq -m ont -a 12 --min-abundance 0.0 -o /home/Users/yl181/memu/magnet_v222_cheese/Cheese/Cheese-10 --threads 40\n",
      "-c /home/Users/ns58/Mob-experiments/Mob-outputs/Cheese/Cheese-11-relative_abundance.tsv -i /home/Users/ns58/Mob-experiments/Data-grouped/Cheese/Cheese-11.fastq -m ont -a 12 --min-abundance 0.0 -o /home/Users/yl181/memu/magnet_v222_cheese/Cheese/Cheese-11 --threads 40\n"
     ]
    }
   ],
   "source": [
    "cmd_records = []\n",
    "\n",
    "for idx, mob_output_path in enumerate(mob_output_paths):\n",
    "    for f in os.listdir(mob_output_path):\n",
    "        if f.startswith('Cheese-') and 'relative_abundance.tsv' in f:\n",
    "            sample_id = '-'.join(f.split('-')[0:2])\n",
    "            input_fastq = os.path.join(input_fastq_paths[idx], f'{sample_id}.fastq')\n",
    "            if not os.path.exists(input_fastq):\n",
    "                print(\"input_fastq missing:\", mob_output_path, sample_id)\n",
    "            classification_output = os.path.join(mob_output_path, f)\n",
    "            if not os.path.exists(classification_output):\n",
    "                print(\"classification_output missing:\", input_fastq_paths[idx], sample_id)\n",
    "                \n",
    "            magnet_output_prefix = input_fastq_paths[idx].split('/')[-1]\n",
    "            if not os.path.exists(os.path.join(magnet_output, magnet_output_prefix)):\n",
    "                os.mkdir(os.path.join(magnet_output, magnet_output_prefix))\n",
    "            magnet_sample_output = os.path.join(magnet_output, magnet_output_prefix, sample_id)\n",
    "            \n",
    "            # print(input_fastq)\n",
    "            # print(classification_output)\n",
    "            # print(magnet_sample_output)\n",
    "            \n",
    "            cmd_args = ['-c', classification_output,\n",
    "            '-i', input_fastq,\n",
    "            '-m', 'ont',\n",
    "            '-a', str(12),\n",
    "            '--min-abundance', str(0.000),\n",
    "            '-o', magnet_sample_output,\n",
    "            '--threads', '40']\n",
    "            \n",
    "            print(\" \".join(cmd_args))\n",
    "            cmd_records.append(cmd_args)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "704786a6-c0a4-4d1c-8107-e7d1d14c248f",
   "metadata": {},
   "source": [
    "## iterate ZYMO"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "09489679-4055-49ae-93c2-070b945add08",
   "metadata": {},
   "outputs": [],
   "source": [
    "mob_output_paths = ['/home/Users/ns58/Mob-experiments/Mob-outputs/Zymo-Even',\n",
    "                    '/home/Users/ns58/Mob-experiments/Mob-outputs/Zymo-Even-0.25',\n",
    "                    '/home/Users/ns58/Mob-experiments/Mob-outputs/Zymo-Even-0.50',\n",
    "                    '/home/Users/ns58/Mob-experiments/Mob-outputs/Zymo-Even-0.75']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "7348bc91-eefb-445b-9aac-43f17218207b",
   "metadata": {},
   "outputs": [],
   "source": [
    "input_fastq_paths =['/home/Users/ns58/Mob-experiments/Data-grouped/Zymo-EVEN-p0.01',\n",
    "                    '/home/Users/ns58/Mob-experiments/Data-grouped/Zymo-EVEN-p0.25',\n",
    "                    '/home/Users/ns58/Mob-experiments/Data-grouped/Zymo-EVEN-p0.50',\n",
    "                    '/home/Users/ns58/Mob-experiments/Data-grouped/Zymo-EVEN-p0.75']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "08e26b08-cecc-4a84-9b05-6c542a67f332",
   "metadata": {},
   "outputs": [],
   "source": [
    "mob_output_paths = ['/home/Users/ns58/Mob-experiments/Mob-outputs/Zymo-Even-0.75']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "fa289727-91c2-4585-8772-8080805c9b5d",
   "metadata": {},
   "outputs": [],
   "source": [
    "input_fastq_paths =['/home/Users/ns58/Mob-experiments/Data-grouped/Zymo-EVEN-p0.75']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "fa4ca408-6c6b-4133-baf1-07c0d0381689",
   "metadata": {},
   "outputs": [],
   "source": [
    "magnet_output = '/home/Users/yl181/memu/magnet_v222_zymo_even_full'\n",
    "if not os.path.exists(magnet_output):\n",
    "    os.mkdir(magnet_output)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "id": "1f1a2350-96bd-40d0-895e-8673f1cb38a2",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "-c /home/Users/ns58/Mob-experiments/Mob-outputs/Zymo-Even-0.75/Zymo-Even-3-relative_abundance.tsv -i /home/Users/ns58/Mob-experiments/Data-grouped/Zymo-EVEN-p0.75/Zymo-Even-3.fastq -m ont -a 12 --min-abundance 0.0 -o /home/Users/yl181/memu/magnet_v222_zymo_even_full/Zymo-EVEN-p0.75/Zymo-Even-3 --threads 40\n",
      "-c /home/Users/ns58/Mob-experiments/Mob-outputs/Zymo-Even-0.75/Zymo-Even-1-relative_abundance.tsv -i /home/Users/ns58/Mob-experiments/Data-grouped/Zymo-EVEN-p0.75/Zymo-Even-1.fastq -m ont -a 12 --min-abundance 0.0 -o /home/Users/yl181/memu/magnet_v222_zymo_even_full/Zymo-EVEN-p0.75/Zymo-Even-1 --threads 40\n",
      "-c /home/Users/ns58/Mob-experiments/Mob-outputs/Zymo-Even-0.75/Zymo-Even-2-relative_abundance.tsv -i /home/Users/ns58/Mob-experiments/Data-grouped/Zymo-EVEN-p0.75/Zymo-Even-2.fastq -m ont -a 12 --min-abundance 0.0 -o /home/Users/yl181/memu/magnet_v222_zymo_even_full/Zymo-EVEN-p0.75/Zymo-Even-2 --threads 40\n",
      "-c /home/Users/ns58/Mob-experiments/Mob-outputs/Zymo-Even-0.75/Zymo-Even-5-relative_abundance.tsv -i /home/Users/ns58/Mob-experiments/Data-grouped/Zymo-EVEN-p0.75/Zymo-Even-5.fastq -m ont -a 12 --min-abundance 0.0 -o /home/Users/yl181/memu/magnet_v222_zymo_even_full/Zymo-EVEN-p0.75/Zymo-Even-5 --threads 40\n",
      "-c /home/Users/ns58/Mob-experiments/Mob-outputs/Zymo-Even-0.75/Zymo-Even-4-relative_abundance.tsv -i /home/Users/ns58/Mob-experiments/Data-grouped/Zymo-EVEN-p0.75/Zymo-Even-4.fastq -m ont -a 12 --min-abundance 0.0 -o /home/Users/yl181/memu/magnet_v222_zymo_even_full/Zymo-EVEN-p0.75/Zymo-Even-4 --threads 40\n"
     ]
    }
   ],
   "source": [
    "cmd_records = []\n",
    "\n",
    "for idx, mob_output_path in enumerate(mob_output_paths):\n",
    "    for f in os.listdir(mob_output_path):\n",
    "        if f.startswith('Zymo-Even-') and 'relative_abundance.tsv' in f:\n",
    "            sample_id = '-'.join(f.split('-')[0:3])\n",
    "            input_fastq = os.path.join(input_fastq_paths[idx], f'{sample_id}.fastq')\n",
    "            if not os.path.exists(input_fastq):\n",
    "                print(\"input_fastq missing:\", mob_output_path, sample_id)\n",
    "            classification_output = os.path.join(mob_output_path, f)\n",
    "            if not os.path.exists(classification_output):\n",
    "                print(\"classification_output missing:\", input_fastq_paths[idx], sample_id)\n",
    "                \n",
    "            magnet_output_prefix = input_fastq_paths[idx].split('/')[-1]\n",
    "            if not os.path.exists(os.path.join(magnet_output, magnet_output_prefix)):\n",
    "                os.mkdir(os.path.join(magnet_output, magnet_output_prefix))\n",
    "            magnet_sample_output = os.path.join(magnet_output, magnet_output_prefix, sample_id)\n",
    "            \n",
    "            # print(input_fastq)\n",
    "            # print(classification_output)\n",
    "            # print(magnet_sample_output)\n",
    "            \n",
    "            cmd_args = ['-c', classification_output,\n",
    "            '-i', input_fastq,\n",
    "            '-m', 'ont',\n",
    "            '-a', str(12),\n",
    "            '--min-abundance', str(0.000),\n",
    "            '-o', magnet_sample_output,\n",
    "            '--threads', '40']\n",
    "            \n",
    "            print(\" \".join(cmd_args))\n",
    "            cmd_records.append(cmd_args)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "id": "0507669c-b6ac-41c1-9d2e-631d51e07155",
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "min_abundance: 0.0\n",
      "num_species: 82\n",
      "2923522    \t GCF_031142765.1 \t Enterococcus sp. IsoGale005    \t IsoGale0.. \t Contig\n",
      "2555396    \t GCF_013350625.1 \t Salmonella sp. SG220           \t SG220      \t Contig\n",
      "1715212    \t GCF_001837115.1 \t Staphylococcus sp. HMSC063H1.. \t HMSC063H.. \t Scaffold\n",
      "562        \t GCF_000008865.2 \t Escherichia coli O157:H7 str.. \t Sakai su.. \t Complete Genome\n",
      "1739314    \t GCF_001811595.1 \t Enterococcus sp. HMSC078F03    \t HMSC078F.. \t Scaffold\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Process ForkPoolWorker-61:\n",
      "Process ForkPoolWorker-63:\n",
      "Process ForkPoolWorker-58:\n",
      "Process ForkPoolWorker-46:\n",
      "Process ForkPoolWorker-43:\n",
      "Exception ignored in: <Finalize object, dead>\n",
      "Traceback (most recent call last):\n",
      "  File \"/home/Users/yl181/miniconda3/envs/magnet/lib/python3.9/multiprocessing/util.py\", line 224, in __call__\n",
      "    res = self._callback(*self._args, **self._kwargs)\n",
      "  File \"/home/Users/yl181/miniconda3/envs/magnet/lib/python3.9/multiprocessing/pool.py\", line 692, in _terminate_pool\n",
      "    cls._help_stuff_finish(inqueue, task_handler, len(pool))\n",
      "  File \"/home/Users/yl181/miniconda3/envs/magnet/lib/python3.9/multiprocessing/pool.py\", line 672, in _help_stuff_finish\n",
      "    inqueue._rlock.acquire()\n",
      "KeyboardInterrupt: \n",
      "Process ForkPoolWorker-77:\n",
      "Process ForkPoolWorker-74:\n",
      "Process ForkPoolWorker-80:\n",
      "Process ForkPoolWorker-70:\n",
      "Process ForkPoolWorker-75:\n",
      "Process ForkPoolWorker-50:\n",
      "Process ForkPoolWorker-71:\n",
      "Process ForkPoolWorker-79:\n",
      "Process ForkPoolWorker-69:\n",
      "Process ForkPoolWorker-72:\n",
      "Process ForkPoolWorker-57:\n",
      "Process ForkPoolWorker-55:\n",
      "Process ForkPoolWorker-45:\n",
      "Process ForkPoolWorker-68:\n",
      "Process ForkPoolWorker-42:\n",
      "Process ForkPoolWorker-41:\n",
      "Process ForkPoolWorker-67:\n",
      "Process ForkPoolWorker-73:\n",
      "Process ForkPoolWorker-76:\n",
      "Process ForkPoolWorker-47:\n",
      "Process ForkPoolWorker-54:\n",
      "Process ForkPoolWorker-49:\n",
      "Process ForkPoolWorker-48:\n",
      "Process ForkPoolWorker-53:\n",
      "Process ForkPoolWorker-44:\n",
      "Process ForkPoolWorker-78:\n",
      "Process ForkPoolWorker-56:\n",
      "Process ForkPoolWorker-60:\n",
      "Process ForkPoolWorker-59:\n",
      "Process ForkPoolWorker-62:\n",
      "Process ForkPoolWorker-51:\n",
      "Process ForkPoolWorker-64:\n",
      "Process ForkPoolWorker-65:\n",
      "Process ForkPoolWorker-66:\n",
      "Traceback (most recent call last):\n",
      "Process ForkPoolWorker-52:\n",
      "Traceback (most recent call last):\n",
      "Traceback (most recent call last):\n",
      "Traceback (most recent call last):\n",
      "Traceback (most recent call last):\n",
      "Traceback (most recent call last):\n",
      "Traceback (most recent call last):\n",
      "Traceback (most recent call last):\n",
      "Traceback (most recent call last):\n",
      "Traceback (most recent call last):\n",
      "Traceback (most recent call last):\n",
      "Traceback (most recent call last):\n",
      "Traceback (most recent call last):\n",
      "Traceback (most recent call last):\n",
      "Traceback (most recent call last):\n",
      "Traceback (most recent call last):\n",
      "Traceback (most recent call last):\n",
      "Traceback (most recent call last):\n",
      "Traceback (most recent call last):\n",
      "Traceback (most recent call last):\n",
      "Traceback (most recent call last):\n",
      "Traceback (most recent call last):\n",
      "Traceback (most recent call last):\n",
      "Traceback (most recent call last):\n",
      "Traceback (most recent call last):\n",
      "Traceback (most recent call last):\n",
      "Traceback (most recent call last):\n",
      "  File \"/home/Users/yl181/miniconda3/envs/magnet/lib/python3.9/multiprocessing/process.py\", line 315, in _bootstrap\n",
      "    self.run()\n",
      "Traceback (most recent call last):\n",
      "Traceback (most recent call last):\n",
      "  File \"/home/Users/yl181/miniconda3/envs/magnet/lib/python3.9/multiprocessing/process.py\", line 315, in _bootstrap\n",
      "    self.run()\n",
      "  File \"/home/Users/yl181/miniconda3/envs/magnet/lib/python3.9/multiprocessing/process.py\", line 315, in _bootstrap\n",
      "    self.run()\n",
      "Traceback (most recent call last):\n",
      "Traceback (most recent call last):\n"
     ]
    },
    {
     "ename": "KeyboardInterrupt",
     "evalue": "",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mKeyboardInterrupt\u001b[0m                         Traceback (most recent call last)",
      "Cell \u001b[0;32mIn [32], line 2\u001b[0m\n\u001b[1;32m      1\u001b[0m \u001b[38;5;28;01mfor\u001b[39;00m cmd_args \u001b[38;5;129;01min\u001b[39;00m cmd_records:\n\u001b[0;32m----> 2\u001b[0m     \u001b[43mrun_magnet\u001b[49m\u001b[43m(\u001b[49m\u001b[43mcmd_args\u001b[49m\u001b[43m)\u001b[49m\n",
      "Cell \u001b[0;32mIn [16], line 48\u001b[0m, in \u001b[0;36mrun_magnet\u001b[0;34m(cmd_args)\u001b[0m\n\u001b[1;32m     45\u001b[0m     \u001b[38;5;66;03m# make valid_kingdom a variable?\u001b[39;00m\n\u001b[1;32m     46\u001b[0m     valid_taxids \u001b[38;5;241m=\u001b[39m filter_input_df(input_df, min_abundance, ncbi_taxa_db, valid_kingdom\u001b[38;5;241m=\u001b[39mvalid_kingdom, ret_subspecies\u001b[38;5;241m=\u001b[39mcall_subspecies)\n\u001b[0;32m---> 48\u001b[0m reference_metadata \u001b[38;5;241m=\u001b[39m \u001b[43mprepare_reference_genomes\u001b[49m\u001b[43m(\u001b[49m\u001b[43mvalid_taxids\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mworking_directory\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mncbi_taxa_db\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43maccession_flag\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43maccession_flag\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mmag_flag\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43mmag_flag\u001b[49m\u001b[43m)\u001b[49m\n\u001b[1;32m     49\u001b[0m downloaded_assemblies \u001b[38;5;241m=\u001b[39m reference_metadata[reference_metadata[\u001b[38;5;124m'\u001b[39m\u001b[38;5;124mDownloaded\u001b[39m\u001b[38;5;124m'\u001b[39m]]\n\u001b[1;32m     51\u001b[0m reference_genome_path \u001b[38;5;241m=\u001b[39m os\u001b[38;5;241m.\u001b[39mpath\u001b[38;5;241m.\u001b[39mjoin(working_directory, \u001b[38;5;124m'\u001b[39m\u001b[38;5;124mreference_genomes\u001b[39m\u001b[38;5;124m'\u001b[39m)\n",
      "File \u001b[0;32m~/memu/magnet/utils/reference_finder.py:173\u001b[0m, in \u001b[0;36mprepare_reference_genomes\u001b[0;34m(taxid_queries, output_directory, ncbi_taxa_db, accession_flag, mag_flag)\u001b[0m\n\u001b[1;32m    171\u001b[0m             download_result\u001b[38;5;241m.\u001b[39mappend(download_reference_accession(taxid, working_dir))\n\u001b[1;32m    172\u001b[0m         \u001b[38;5;28;01melse\u001b[39;00m:\n\u001b[0;32m--> 173\u001b[0m             download_result\u001b[38;5;241m.\u001b[39mappend(\u001b[43mdownload_reference_genome\u001b[49m\u001b[43m(\u001b[49m\u001b[43mtaxid\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mworking_dir\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mmag_flag\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43mmag_flag\u001b[49m\u001b[43m)\u001b[49m)\n\u001b[1;32m    175\u001b[0m     unpack(working_dir, output_directory)\n\u001b[1;32m    176\u001b[0m     reference_metadata \u001b[38;5;241m=\u001b[39m pd\u001b[38;5;241m.\u001b[39mDataFrame(download_result,\n\u001b[1;32m    177\u001b[0m                                       columns\u001b[38;5;241m=\u001b[39m[\u001b[38;5;124m'\u001b[39m\u001b[38;5;124mTaxonomy ID\u001b[39m\u001b[38;5;124m'\u001b[39m, \n\u001b[1;32m    178\u001b[0m                                                \u001b[38;5;124m'\u001b[39m\u001b[38;5;124mAssembly Accession ID\u001b[39m\u001b[38;5;124m'\u001b[39m, \n\u001b[0;32m   (...)\u001b[0m\n\u001b[1;32m    184\u001b[0m \t\t\t\t\t\t\t\t\t\t\t   \u001b[38;5;124m'\u001b[39m\u001b[38;5;124mTotal Length\u001b[39m\u001b[38;5;124m'\u001b[39m,\n\u001b[1;32m    185\u001b[0m                                                \u001b[38;5;124m'\u001b[39m\u001b[38;5;124mDownloaded\u001b[39m\u001b[38;5;124m'\u001b[39m])\n",
      "File \u001b[0;32m~/memu/magnet/utils/reference_finder.py:74\u001b[0m, in \u001b[0;36mdownload_reference_genome\u001b[0;34m(taxid, working_dir, mag_flag)\u001b[0m\n\u001b[1;32m     72\u001b[0m \u001b[38;5;66;03m# check complete assembly\u001b[39;00m\n\u001b[1;32m     73\u001b[0m \u001b[38;5;28;01mif\u001b[39;00m res[\u001b[38;5;124m'\u001b[39m\u001b[38;5;124mtotal_count\u001b[39m\u001b[38;5;124m'\u001b[39m] \u001b[38;5;241m==\u001b[39m \u001b[38;5;241m0\u001b[39m:\n\u001b[0;32m---> 74\u001b[0m     res \u001b[38;5;241m=\u001b[39m \u001b[43mrun_datasets_summary\u001b[49m\u001b[43m(\u001b[49m\u001b[43mtaxid\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43m[\u001b[49m\u001b[43m]\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mmag_flag\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43mmag_flag\u001b[49m\u001b[43m)\u001b[49m\n\u001b[1;32m     75\u001b[0m     \u001b[38;5;66;03m#source = 'Assembly'\u001b[39;00m\n\u001b[1;32m     76\u001b[0m \n\u001b[1;32m     77\u001b[0m \u001b[38;5;66;03m# check chromosome assembly\u001b[39;00m\n\u001b[1;32m     78\u001b[0m \u001b[38;5;28;01mif\u001b[39;00m res[\u001b[38;5;124m'\u001b[39m\u001b[38;5;124mtotal_count\u001b[39m\u001b[38;5;124m'\u001b[39m] \u001b[38;5;241m==\u001b[39m \u001b[38;5;241m0\u001b[39m:\n",
      "File \u001b[0;32m~/memu/magnet/utils/reference_finder.py:30\u001b[0m, in \u001b[0;36mrun_datasets_summary\u001b[0;34m(taxid, flags, assembly_level, mag_flag)\u001b[0m\n\u001b[1;32m     20\u001b[0m args \u001b[38;5;241m=\u001b[39m [\u001b[38;5;124m'\u001b[39m\u001b[38;5;124mdatasets\u001b[39m\u001b[38;5;124m'\u001b[39m, \u001b[38;5;124m'\u001b[39m\u001b[38;5;124msummary\u001b[39m\u001b[38;5;124m'\u001b[39m, \n\u001b[1;32m     21\u001b[0m         \u001b[38;5;124m'\u001b[39m\u001b[38;5;124mgenome\u001b[39m\u001b[38;5;124m'\u001b[39m,\n\u001b[1;32m     22\u001b[0m         \u001b[38;5;124m'\u001b[39m\u001b[38;5;124mtaxon\u001b[39m\u001b[38;5;124m'\u001b[39m, \u001b[38;5;28mstr\u001b[39m(taxid),\n\u001b[0;32m   (...)\u001b[0m\n\u001b[1;32m     25\u001b[0m         \u001b[38;5;124m'\u001b[39m\u001b[38;5;124m--assembly-level\u001b[39m\u001b[38;5;124m'\u001b[39m, assembly_level,\n\u001b[1;32m     26\u001b[0m         \u001b[38;5;124m'\u001b[39m\u001b[38;5;124m--limit\u001b[39m\u001b[38;5;124m'\u001b[39m, \u001b[38;5;124m'\u001b[39m\u001b[38;5;124m1\u001b[39m\u001b[38;5;124m'\u001b[39m]\n\u001b[1;32m     28\u001b[0m args \u001b[38;5;241m=\u001b[39m args \u001b[38;5;241m+\u001b[39m flags\n\u001b[0;32m---> 30\u001b[0m grepOut  \u001b[38;5;241m=\u001b[39m \u001b[43msubprocess\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mrun\u001b[49m\u001b[43m(\u001b[49m\u001b[43margs\u001b[49m\u001b[43m,\u001b[49m\n\u001b[1;32m     31\u001b[0m \u001b[43m                          \u001b[49m\u001b[43muniversal_newlines\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[38;5;28;43;01mTrue\u001b[39;49;00m\u001b[43m,\u001b[49m\n\u001b[1;32m     32\u001b[0m \u001b[43m                          \u001b[49m\u001b[43mstdout\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43msubprocess\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mPIPE\u001b[49m\u001b[43m,\u001b[49m\n\u001b[1;32m     33\u001b[0m \u001b[43m                          \u001b[49m\u001b[43mstderr\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43msubprocess\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mDEVNULL\u001b[49m\u001b[43m)\u001b[49m\n\u001b[1;32m     35\u001b[0m \u001b[38;5;28;01mreturn\u001b[39;00m json\u001b[38;5;241m.\u001b[39mloads(grepOut\u001b[38;5;241m.\u001b[39mstdout\u001b[38;5;241m.\u001b[39mstrip())\n",
      "File \u001b[0;32m~/miniconda3/envs/magnet/lib/python3.9/subprocess.py:507\u001b[0m, in \u001b[0;36mrun\u001b[0;34m(input, capture_output, timeout, check, *popenargs, **kwargs)\u001b[0m\n\u001b[1;32m    505\u001b[0m \u001b[38;5;28;01mwith\u001b[39;00m Popen(\u001b[38;5;241m*\u001b[39mpopenargs, \u001b[38;5;241m*\u001b[39m\u001b[38;5;241m*\u001b[39mkwargs) \u001b[38;5;28;01mas\u001b[39;00m process:\n\u001b[1;32m    506\u001b[0m     \u001b[38;5;28;01mtry\u001b[39;00m:\n\u001b[0;32m--> 507\u001b[0m         stdout, stderr \u001b[38;5;241m=\u001b[39m \u001b[43mprocess\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mcommunicate\u001b[49m\u001b[43m(\u001b[49m\u001b[38;5;28;43minput\u001b[39;49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mtimeout\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43mtimeout\u001b[49m\u001b[43m)\u001b[49m\n\u001b[1;32m    508\u001b[0m     \u001b[38;5;28;01mexcept\u001b[39;00m TimeoutExpired \u001b[38;5;28;01mas\u001b[39;00m exc:\n\u001b[1;32m    509\u001b[0m         process\u001b[38;5;241m.\u001b[39mkill()\n",
      "File \u001b[0;32m~/miniconda3/envs/magnet/lib/python3.9/subprocess.py:1121\u001b[0m, in \u001b[0;36mPopen.communicate\u001b[0;34m(self, input, timeout)\u001b[0m\n\u001b[1;32m   1119\u001b[0m     \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39m_stdin_write(\u001b[38;5;28minput\u001b[39m)\n\u001b[1;32m   1120\u001b[0m \u001b[38;5;28;01melif\u001b[39;00m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mstdout:\n\u001b[0;32m-> 1121\u001b[0m     stdout \u001b[38;5;241m=\u001b[39m \u001b[38;5;28;43mself\u001b[39;49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mstdout\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mread\u001b[49m\u001b[43m(\u001b[49m\u001b[43m)\u001b[49m\n\u001b[1;32m   1122\u001b[0m     \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mstdout\u001b[38;5;241m.\u001b[39mclose()\n\u001b[1;32m   1123\u001b[0m \u001b[38;5;28;01melif\u001b[39;00m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mstderr:\n",
      "\u001b[0;31mKeyboardInterrupt\u001b[0m: "
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "  File \"/home/Users/yl181/miniconda3/envs/magnet/lib/python3.9/multiprocessing/process.py\", line 315, in _bootstrap\n",
      "    self.run()\n",
      "  File \"/home/Users/yl181/miniconda3/envs/magnet/lib/python3.9/multiprocessing/process.py\", line 315, in _bootstrap\n",
      "    self.run()\n",
      "Traceback (most recent call last):\n",
      "  File \"/home/Users/yl181/miniconda3/envs/magnet/lib/python3.9/multiprocessing/process.py\", line 315, in _bootstrap\n",
      "    self.run()\n",
      "  File \"/home/Users/yl181/miniconda3/envs/magnet/lib/python3.9/multiprocessing/process.py\", line 315, in _bootstrap\n",
      "    self.run()\n",
      "Traceback (most recent call last):\n",
      "Traceback (most recent call last):\n",
      "  File \"/home/Users/yl181/miniconda3/envs/magnet/lib/python3.9/multiprocessing/process.py\", line 315, in _bootstrap\n",
      "    self.run()\n",
      "Traceback (most recent call last):\n",
      "Traceback (most recent call last):\n",
      "Traceback (most recent call last):\n",
      "  File \"/home/Users/yl181/miniconda3/envs/magnet/lib/python3.9/multiprocessing/process.py\", line 315, in _bootstrap\n",
      "    self.run()\n",
      "Traceback (most recent call last):\n",
      "  File \"/home/Users/yl181/miniconda3/envs/magnet/lib/python3.9/multiprocessing/process.py\", line 315, in _bootstrap\n",
      "    self.run()\n",
      "  File \"/home/Users/yl181/miniconda3/envs/magnet/lib/python3.9/multiprocessing/process.py\", line 315, in _bootstrap\n",
      "    self.run()\n",
      "Traceback (most recent call last):\n",
      "  File \"/home/Users/yl181/miniconda3/envs/magnet/lib/python3.9/multiprocessing/process.py\", line 315, in _bootstrap\n",
      "    self.run()\n",
      "  File \"/home/Users/yl181/miniconda3/envs/magnet/lib/python3.9/multiprocessing/process.py\", line 315, in _bootstrap\n",
      "    self.run()\n",
      "Traceback (most recent call last):\n",
      "  File \"/home/Users/yl181/miniconda3/envs/magnet/lib/python3.9/multiprocessing/process.py\", line 315, in _bootstrap\n",
      "    self.run()\n",
      "  File \"/home/Users/yl181/miniconda3/envs/magnet/lib/python3.9/multiprocessing/process.py\", line 315, in _bootstrap\n",
      "    self.run()\n",
      "  File \"/home/Users/yl181/miniconda3/envs/magnet/lib/python3.9/multiprocessing/process.py\", line 315, in _bootstrap\n",
      "    self.run()\n",
      "  File \"/home/Users/yl181/miniconda3/envs/magnet/lib/python3.9/multiprocessing/process.py\", line 315, in _bootstrap\n",
      "    self.run()\n",
      "  File \"/home/Users/yl181/miniconda3/envs/magnet/lib/python3.9/multiprocessing/process.py\", line 315, in _bootstrap\n",
      "    self.run()\n",
      "  File \"/home/Users/yl181/miniconda3/envs/magnet/lib/python3.9/multiprocessing/process.py\", line 315, in _bootstrap\n",
      "    self.run()\n",
      "  File \"/home/Users/yl181/miniconda3/envs/magnet/lib/python3.9/multiprocessing/process.py\", line 315, in _bootstrap\n",
      "    self.run()\n",
      "  File \"/home/Users/yl181/miniconda3/envs/magnet/lib/python3.9/multiprocessing/process.py\", line 315, in _bootstrap\n",
      "    self.run()\n",
      "  File \"/home/Users/yl181/miniconda3/envs/magnet/lib/python3.9/multiprocessing/process.py\", line 315, in _bootstrap\n",
      "    self.run()\n",
      "  File \"/home/Users/yl181/miniconda3/envs/magnet/lib/python3.9/multiprocessing/process.py\", line 315, in _bootstrap\n",
      "    self.run()\n",
      "  File \"/home/Users/yl181/miniconda3/envs/magnet/lib/python3.9/multiprocessing/process.py\", line 315, in _bootstrap\n",
      "    self.run()\n",
      "  File \"/home/Users/yl181/miniconda3/envs/magnet/lib/python3.9/multiprocessing/process.py\", line 315, in _bootstrap\n",
      "    self.run()\n",
      "  File \"/home/Users/yl181/miniconda3/envs/magnet/lib/python3.9/multiprocessing/process.py\", line 315, in _bootstrap\n",
      "    self.run()\n",
      "  File \"/home/Users/yl181/miniconda3/envs/magnet/lib/python3.9/multiprocessing/process.py\", line 108, in run\n",
      "    self._target(*self._args, **self._kwargs)\n",
      "  File \"/home/Users/yl181/miniconda3/envs/magnet/lib/python3.9/multiprocessing/process.py\", line 315, in _bootstrap\n",
      "    self.run()\n",
      "  File \"/home/Users/yl181/miniconda3/envs/magnet/lib/python3.9/multiprocessing/process.py\", line 315, in _bootstrap\n",
      "    self.run()\n",
      "  File \"/home/Users/yl181/miniconda3/envs/magnet/lib/python3.9/multiprocessing/process.py\", line 315, in _bootstrap\n",
      "    self.run()\n",
      "  File \"/home/Users/yl181/miniconda3/envs/magnet/lib/python3.9/multiprocessing/process.py\", line 108, in run\n",
      "    self._target(*self._args, **self._kwargs)\n",
      "  File \"/home/Users/yl181/miniconda3/envs/magnet/lib/python3.9/multiprocessing/process.py\", line 108, in run\n",
      "    self._target(*self._args, **self._kwargs)\n",
      "  File \"/home/Users/yl181/miniconda3/envs/magnet/lib/python3.9/multiprocessing/process.py\", line 315, in _bootstrap\n",
      "    self.run()\n",
      "  File \"/home/Users/yl181/miniconda3/envs/magnet/lib/python3.9/multiprocessing/process.py\", line 315, in _bootstrap\n",
      "    self.run()\n",
      "  File \"/home/Users/yl181/miniconda3/envs/magnet/lib/python3.9/multiprocessing/process.py\", line 108, in run\n",
      "    self._target(*self._args, **self._kwargs)\n",
      "  File \"/home/Users/yl181/miniconda3/envs/magnet/lib/python3.9/multiprocessing/process.py\", line 108, in run\n",
      "    self._target(*self._args, **self._kwargs)\n",
      "  File \"/home/Users/yl181/miniconda3/envs/magnet/lib/python3.9/multiprocessing/process.py\", line 315, in _bootstrap\n",
      "    self.run()\n",
      "  File \"/home/Users/yl181/miniconda3/envs/magnet/lib/python3.9/multiprocessing/process.py\", line 108, in run\n",
      "    self._target(*self._args, **self._kwargs)\n",
      "  File \"/home/Users/yl181/miniconda3/envs/magnet/lib/python3.9/multiprocessing/process.py\", line 315, in _bootstrap\n",
      "    self.run()\n",
      "  File \"/home/Users/yl181/miniconda3/envs/magnet/lib/python3.9/multiprocessing/process.py\", line 108, in run\n",
      "    self._target(*self._args, **self._kwargs)\n",
      "  File \"/home/Users/yl181/miniconda3/envs/magnet/lib/python3.9/multiprocessing/process.py\", line 108, in run\n",
      "    self._target(*self._args, **self._kwargs)\n",
      "  File \"/home/Users/yl181/miniconda3/envs/magnet/lib/python3.9/multiprocessing/process.py\", line 315, in _bootstrap\n",
      "    self.run()\n",
      "  File \"/home/Users/yl181/miniconda3/envs/magnet/lib/python3.9/multiprocessing/process.py\", line 315, in _bootstrap\n",
      "    self.run()\n",
      "  File \"/home/Users/yl181/miniconda3/envs/magnet/lib/python3.9/multiprocessing/process.py\", line 315, in _bootstrap\n",
      "    self.run()\n",
      "  File \"/home/Users/yl181/miniconda3/envs/magnet/lib/python3.9/multiprocessing/process.py\", line 108, in run\n",
      "    self._target(*self._args, **self._kwargs)\n",
      "  File \"/home/Users/yl181/miniconda3/envs/magnet/lib/python3.9/multiprocessing/process.py\", line 315, in _bootstrap\n",
      "    self.run()\n",
      "  File \"/home/Users/yl181/miniconda3/envs/magnet/lib/python3.9/multiprocessing/process.py\", line 108, in run\n",
      "    self._target(*self._args, **self._kwargs)\n",
      "  File \"/home/Users/yl181/miniconda3/envs/magnet/lib/python3.9/multiprocessing/process.py\", line 108, in run\n",
      "    self._target(*self._args, **self._kwargs)\n",
      "  File \"/home/Users/yl181/miniconda3/envs/magnet/lib/python3.9/multiprocessing/process.py\", line 315, in _bootstrap\n",
      "    self.run()\n",
      "  File \"/home/Users/yl181/miniconda3/envs/magnet/lib/python3.9/multiprocessing/process.py\", line 108, in run\n",
      "    self._target(*self._args, **self._kwargs)\n",
      "  File \"/home/Users/yl181/miniconda3/envs/magnet/lib/python3.9/multiprocessing/process.py\", line 315, in _bootstrap\n",
      "    self.run()\n",
      "  File \"/home/Users/yl181/miniconda3/envs/magnet/lib/python3.9/multiprocessing/process.py\", line 108, in run\n",
      "    self._target(*self._args, **self._kwargs)\n",
      "  File \"/home/Users/yl181/miniconda3/envs/magnet/lib/python3.9/multiprocessing/process.py\", line 108, in run\n",
      "    self._target(*self._args, **self._kwargs)\n",
      "  File \"/home/Users/yl181/miniconda3/envs/magnet/lib/python3.9/multiprocessing/process.py\", line 315, in _bootstrap\n",
      "    self.run()\n",
      "  File \"/home/Users/yl181/miniconda3/envs/magnet/lib/python3.9/multiprocessing/process.py\", line 108, in run\n",
      "    self._target(*self._args, **self._kwargs)\n",
      "  File \"/home/Users/yl181/miniconda3/envs/magnet/lib/python3.9/multiprocessing/process.py\", line 108, in run\n",
      "    self._target(*self._args, **self._kwargs)\n",
      "  File \"/home/Users/yl181/miniconda3/envs/magnet/lib/python3.9/multiprocessing/process.py\", line 108, in run\n",
      "    self._target(*self._args, **self._kwargs)\n",
      "  File \"/home/Users/yl181/miniconda3/envs/magnet/lib/python3.9/multiprocessing/process.py\", line 108, in run\n",
      "    self._target(*self._args, **self._kwargs)\n",
      "  File \"/home/Users/yl181/miniconda3/envs/magnet/lib/python3.9/multiprocessing/process.py\", line 108, in run\n",
      "    self._target(*self._args, **self._kwargs)\n",
      "  File \"/home/Users/yl181/miniconda3/envs/magnet/lib/python3.9/multiprocessing/process.py\", line 108, in run\n",
      "    self._target(*self._args, **self._kwargs)\n",
      "  File \"/home/Users/yl181/miniconda3/envs/magnet/lib/python3.9/multiprocessing/process.py\", line 108, in run\n",
      "    self._target(*self._args, **self._kwargs)\n",
      "  File \"/home/Users/yl181/miniconda3/envs/magnet/lib/python3.9/multiprocessing/process.py\", line 108, in run\n",
      "    self._target(*self._args, **self._kwargs)\n",
      "  File \"/home/Users/yl181/miniconda3/envs/magnet/lib/python3.9/multiprocessing/process.py\", line 108, in run\n",
      "    self._target(*self._args, **self._kwargs)\n",
      "  File \"/home/Users/yl181/miniconda3/envs/magnet/lib/python3.9/multiprocessing/process.py\", line 108, in run\n",
      "    self._target(*self._args, **self._kwargs)\n",
      "  File \"/home/Users/yl181/miniconda3/envs/magnet/lib/python3.9/multiprocessing/process.py\", line 108, in run\n",
      "    self._target(*self._args, **self._kwargs)\n",
      "  File \"/home/Users/yl181/miniconda3/envs/magnet/lib/python3.9/multiprocessing/pool.py\", line 114, in worker\n",
      "    task = get()\n",
      "  File \"/home/Users/yl181/miniconda3/envs/magnet/lib/python3.9/multiprocessing/process.py\", line 108, in run\n",
      "    self._target(*self._args, **self._kwargs)\n",
      "  File \"/home/Users/yl181/miniconda3/envs/magnet/lib/python3.9/multiprocessing/process.py\", line 108, in run\n",
      "    self._target(*self._args, **self._kwargs)\n",
      "  File \"/home/Users/yl181/miniconda3/envs/magnet/lib/python3.9/multiprocessing/process.py\", line 108, in run\n",
      "    self._target(*self._args, **self._kwargs)\n",
      "  File \"/home/Users/yl181/miniconda3/envs/magnet/lib/python3.9/multiprocessing/process.py\", line 108, in run\n",
      "    self._target(*self._args, **self._kwargs)\n",
      "  File \"/home/Users/yl181/miniconda3/envs/magnet/lib/python3.9/multiprocessing/pool.py\", line 114, in worker\n",
      "    task = get()\n",
      "  File \"/home/Users/yl181/miniconda3/envs/magnet/lib/python3.9/multiprocessing/pool.py\", line 114, in worker\n",
      "    task = get()\n",
      "  File \"/home/Users/yl181/miniconda3/envs/magnet/lib/python3.9/multiprocessing/process.py\", line 108, in run\n",
      "    self._target(*self._args, **self._kwargs)\n",
      "  File \"/home/Users/yl181/miniconda3/envs/magnet/lib/python3.9/multiprocessing/process.py\", line 108, in run\n",
      "    self._target(*self._args, **self._kwargs)\n",
      "  File \"/home/Users/yl181/miniconda3/envs/magnet/lib/python3.9/multiprocessing/pool.py\", line 114, in worker\n",
      "    task = get()\n",
      "  File \"/home/Users/yl181/miniconda3/envs/magnet/lib/python3.9/multiprocessing/pool.py\", line 114, in worker\n",
      "    task = get()\n",
      "  File \"/home/Users/yl181/miniconda3/envs/magnet/lib/python3.9/multiprocessing/pool.py\", line 114, in worker\n",
      "    task = get()\n",
      "  File \"/home/Users/yl181/miniconda3/envs/magnet/lib/python3.9/multiprocessing/process.py\", line 108, in run\n",
      "    self._target(*self._args, **self._kwargs)\n",
      "  File \"/home/Users/yl181/miniconda3/envs/magnet/lib/python3.9/multiprocessing/process.py\", line 108, in run\n",
      "    self._target(*self._args, **self._kwargs)\n",
      "  File \"/home/Users/yl181/miniconda3/envs/magnet/lib/python3.9/multiprocessing/pool.py\", line 114, in worker\n",
      "    task = get()\n",
      "  File \"/home/Users/yl181/miniconda3/envs/magnet/lib/python3.9/multiprocessing/pool.py\", line 114, in worker\n",
      "    task = get()\n",
      "  File \"/home/Users/yl181/miniconda3/envs/magnet/lib/python3.9/multiprocessing/process.py\", line 108, in run\n",
      "    self._target(*self._args, **self._kwargs)\n",
      "  File \"/home/Users/yl181/miniconda3/envs/magnet/lib/python3.9/multiprocessing/process.py\", line 108, in run\n",
      "    self._target(*self._args, **self._kwargs)\n",
      "  File \"/home/Users/yl181/miniconda3/envs/magnet/lib/python3.9/multiprocessing/process.py\", line 108, in run\n",
      "    self._target(*self._args, **self._kwargs)\n",
      "  File \"/home/Users/yl181/miniconda3/envs/magnet/lib/python3.9/multiprocessing/pool.py\", line 114, in worker\n",
      "    task = get()\n",
      "  File \"/home/Users/yl181/miniconda3/envs/magnet/lib/python3.9/multiprocessing/process.py\", line 108, in run\n",
      "    self._target(*self._args, **self._kwargs)\n",
      "  File \"/home/Users/yl181/miniconda3/envs/magnet/lib/python3.9/multiprocessing/pool.py\", line 114, in worker\n",
      "    task = get()\n",
      "  File \"/home/Users/yl181/miniconda3/envs/magnet/lib/python3.9/multiprocessing/pool.py\", line 114, in worker\n",
      "    task = get()\n",
      "  File \"/home/Users/yl181/miniconda3/envs/magnet/lib/python3.9/multiprocessing/pool.py\", line 114, in worker\n",
      "    task = get()\n",
      "  File \"/home/Users/yl181/miniconda3/envs/magnet/lib/python3.9/multiprocessing/process.py\", line 108, in run\n",
      "    self._target(*self._args, **self._kwargs)\n",
      "  File \"/home/Users/yl181/miniconda3/envs/magnet/lib/python3.9/multiprocessing/process.py\", line 108, in run\n",
      "    self._target(*self._args, **self._kwargs)\n",
      "  File \"/home/Users/yl181/miniconda3/envs/magnet/lib/python3.9/multiprocessing/pool.py\", line 114, in worker\n",
      "    task = get()\n",
      "  File \"/home/Users/yl181/miniconda3/envs/magnet/lib/python3.9/multiprocessing/pool.py\", line 114, in worker\n",
      "    task = get()\n",
      "  File \"/home/Users/yl181/miniconda3/envs/magnet/lib/python3.9/multiprocessing/process.py\", line 108, in run\n",
      "    self._target(*self._args, **self._kwargs)\n",
      "  File \"/home/Users/yl181/miniconda3/envs/magnet/lib/python3.9/multiprocessing/pool.py\", line 114, in worker\n",
      "    task = get()\n",
      "  File \"/home/Users/yl181/miniconda3/envs/magnet/lib/python3.9/multiprocessing/pool.py\", line 114, in worker\n",
      "    task = get()\n",
      "  File \"/home/Users/yl181/miniconda3/envs/magnet/lib/python3.9/multiprocessing/pool.py\", line 114, in worker\n",
      "    task = get()\n",
      "  File \"/home/Users/yl181/miniconda3/envs/magnet/lib/python3.9/multiprocessing/pool.py\", line 114, in worker\n",
      "    task = get()\n",
      "  File \"/home/Users/yl181/miniconda3/envs/magnet/lib/python3.9/multiprocessing/pool.py\", line 114, in worker\n",
      "    task = get()\n",
      "  File \"/home/Users/yl181/miniconda3/envs/magnet/lib/python3.9/multiprocessing/pool.py\", line 114, in worker\n",
      "    task = get()\n",
      "  File \"/home/Users/yl181/miniconda3/envs/magnet/lib/python3.9/multiprocessing/pool.py\", line 114, in worker\n",
      "    task = get()\n",
      "  File \"/home/Users/yl181/miniconda3/envs/magnet/lib/python3.9/multiprocessing/pool.py\", line 114, in worker\n",
      "    task = get()\n",
      "  File \"/home/Users/yl181/miniconda3/envs/magnet/lib/python3.9/multiprocessing/pool.py\", line 114, in worker\n",
      "    task = get()\n",
      "  File \"/home/Users/yl181/miniconda3/envs/magnet/lib/python3.9/multiprocessing/pool.py\", line 114, in worker\n",
      "    task = get()\n",
      "  File \"/home/Users/yl181/miniconda3/envs/magnet/lib/python3.9/multiprocessing/pool.py\", line 114, in worker\n",
      "    task = get()\n",
      "  File \"/home/Users/yl181/miniconda3/envs/magnet/lib/python3.9/multiprocessing/queues.py\", line 364, in get\n",
      "    with self._rlock:\n",
      "  File \"/home/Users/yl181/miniconda3/envs/magnet/lib/python3.9/multiprocessing/pool.py\", line 114, in worker\n",
      "    task = get()\n",
      "  File \"/home/Users/yl181/miniconda3/envs/magnet/lib/python3.9/multiprocessing/pool.py\", line 114, in worker\n",
      "    task = get()\n",
      "  File \"/home/Users/yl181/miniconda3/envs/magnet/lib/python3.9/multiprocessing/pool.py\", line 114, in worker\n",
      "    task = get()\n",
      "  File \"/home/Users/yl181/miniconda3/envs/magnet/lib/python3.9/multiprocessing/pool.py\", line 114, in worker\n",
      "    task = get()\n",
      "  File \"/home/Users/yl181/miniconda3/envs/magnet/lib/python3.9/multiprocessing/queues.py\", line 364, in get\n",
      "    with self._rlock:\n",
      "  File \"/home/Users/yl181/miniconda3/envs/magnet/lib/python3.9/multiprocessing/queues.py\", line 364, in get\n",
      "    with self._rlock:\n",
      "  File \"/home/Users/yl181/miniconda3/envs/magnet/lib/python3.9/multiprocessing/pool.py\", line 114, in worker\n",
      "    task = get()\n",
      "  File \"/home/Users/yl181/miniconda3/envs/magnet/lib/python3.9/multiprocessing/queues.py\", line 364, in get\n",
      "    with self._rlock:\n",
      "  File \"/home/Users/yl181/miniconda3/envs/magnet/lib/python3.9/multiprocessing/pool.py\", line 114, in worker\n",
      "    task = get()\n",
      "  File \"/home/Users/yl181/miniconda3/envs/magnet/lib/python3.9/multiprocessing/queues.py\", line 364, in get\n",
      "    with self._rlock:\n",
      "  File \"/home/Users/yl181/miniconda3/envs/magnet/lib/python3.9/multiprocessing/queues.py\", line 364, in get\n",
      "    with self._rlock:\n",
      "  File \"/home/Users/yl181/miniconda3/envs/magnet/lib/python3.9/multiprocessing/pool.py\", line 114, in worker\n",
      "    task = get()\n",
      "  File \"/home/Users/yl181/miniconda3/envs/magnet/lib/python3.9/multiprocessing/pool.py\", line 114, in worker\n",
      "    task = get()\n",
      "  File \"/home/Users/yl181/miniconda3/envs/magnet/lib/python3.9/multiprocessing/queues.py\", line 364, in get\n",
      "    with self._rlock:\n",
      "  File \"/home/Users/yl181/miniconda3/envs/magnet/lib/python3.9/multiprocessing/queues.py\", line 364, in get\n",
      "    with self._rlock:\n",
      "  File \"/home/Users/yl181/miniconda3/envs/magnet/lib/python3.9/multiprocessing/pool.py\", line 114, in worker\n",
      "    task = get()\n",
      "  File \"/home/Users/yl181/miniconda3/envs/magnet/lib/python3.9/multiprocessing/pool.py\", line 114, in worker\n",
      "    task = get()\n",
      "  File \"/home/Users/yl181/miniconda3/envs/magnet/lib/python3.9/multiprocessing/pool.py\", line 114, in worker\n",
      "    task = get()\n",
      "  File \"/home/Users/yl181/miniconda3/envs/magnet/lib/python3.9/multiprocessing/queues.py\", line 364, in get\n",
      "    with self._rlock:\n",
      "  File \"/home/Users/yl181/miniconda3/envs/magnet/lib/python3.9/multiprocessing/pool.py\", line 114, in worker\n",
      "    task = get()\n",
      "  File \"/home/Users/yl181/miniconda3/envs/magnet/lib/python3.9/multiprocessing/queues.py\", line 364, in get\n",
      "    with self._rlock:\n",
      "  File \"/home/Users/yl181/miniconda3/envs/magnet/lib/python3.9/multiprocessing/queues.py\", line 364, in get\n",
      "    with self._rlock:\n",
      "  File \"/home/Users/yl181/miniconda3/envs/magnet/lib/python3.9/multiprocessing/queues.py\", line 364, in get\n",
      "    with self._rlock:\n",
      "  File \"/home/Users/yl181/miniconda3/envs/magnet/lib/python3.9/multiprocessing/pool.py\", line 114, in worker\n",
      "    task = get()\n",
      "  File \"/home/Users/yl181/miniconda3/envs/magnet/lib/python3.9/multiprocessing/pool.py\", line 114, in worker\n",
      "    task = get()\n",
      "  File \"/home/Users/yl181/miniconda3/envs/magnet/lib/python3.9/multiprocessing/queues.py\", line 364, in get\n",
      "    with self._rlock:\n",
      "  File \"/home/Users/yl181/miniconda3/envs/magnet/lib/python3.9/multiprocessing/queues.py\", line 364, in get\n",
      "    with self._rlock:\n",
      "  File \"/home/Users/yl181/miniconda3/envs/magnet/lib/python3.9/multiprocessing/queues.py\", line 364, in get\n",
      "    with self._rlock:\n",
      "  File \"/home/Users/yl181/miniconda3/envs/magnet/lib/python3.9/multiprocessing/pool.py\", line 114, in worker\n",
      "    task = get()\n",
      "  File \"/home/Users/yl181/miniconda3/envs/magnet/lib/python3.9/multiprocessing/queues.py\", line 364, in get\n",
      "    with self._rlock:\n",
      "  File \"/home/Users/yl181/miniconda3/envs/magnet/lib/python3.9/multiprocessing/queues.py\", line 364, in get\n",
      "    with self._rlock:\n",
      "  File \"/home/Users/yl181/miniconda3/envs/magnet/lib/python3.9/multiprocessing/queues.py\", line 364, in get\n",
      "    with self._rlock:\n",
      "  File \"/home/Users/yl181/miniconda3/envs/magnet/lib/python3.9/multiprocessing/queues.py\", line 364, in get\n",
      "    with self._rlock:\n",
      "  File \"/home/Users/yl181/miniconda3/envs/magnet/lib/python3.9/multiprocessing/queues.py\", line 364, in get\n",
      "    with self._rlock:\n",
      "  File \"/home/Users/yl181/miniconda3/envs/magnet/lib/python3.9/multiprocessing/queues.py\", line 364, in get\n",
      "    with self._rlock:\n",
      "  File \"/home/Users/yl181/miniconda3/envs/magnet/lib/python3.9/multiprocessing/queues.py\", line 364, in get\n",
      "    with self._rlock:\n",
      "  File \"/home/Users/yl181/miniconda3/envs/magnet/lib/python3.9/multiprocessing/queues.py\", line 364, in get\n",
      "    with self._rlock:\n",
      "  File \"/home/Users/yl181/miniconda3/envs/magnet/lib/python3.9/multiprocessing/queues.py\", line 364, in get\n",
      "    with self._rlock:\n",
      "  File \"/home/Users/yl181/miniconda3/envs/magnet/lib/python3.9/multiprocessing/queues.py\", line 364, in get\n",
      "    with self._rlock:\n",
      "  File \"/home/Users/yl181/miniconda3/envs/magnet/lib/python3.9/multiprocessing/synchronize.py\", line 95, in __enter__\n",
      "    return self._semlock.__enter__()\n",
      "  File \"/home/Users/yl181/miniconda3/envs/magnet/lib/python3.9/multiprocessing/queues.py\", line 364, in get\n",
      "    with self._rlock:\n",
      "  File \"/home/Users/yl181/miniconda3/envs/magnet/lib/python3.9/multiprocessing/queues.py\", line 364, in get\n",
      "    with self._rlock:\n",
      "  File \"/home/Users/yl181/miniconda3/envs/magnet/lib/python3.9/multiprocessing/queues.py\", line 364, in get\n",
      "    with self._rlock:\n",
      "  File \"/home/Users/yl181/miniconda3/envs/magnet/lib/python3.9/multiprocessing/queues.py\", line 364, in get\n",
      "    with self._rlock:\n",
      "  File \"/home/Users/yl181/miniconda3/envs/magnet/lib/python3.9/multiprocessing/synchronize.py\", line 95, in __enter__\n",
      "    return self._semlock.__enter__()\n",
      "  File \"/home/Users/yl181/miniconda3/envs/magnet/lib/python3.9/multiprocessing/synchronize.py\", line 95, in __enter__\n",
      "    return self._semlock.__enter__()\n",
      "  File \"/home/Users/yl181/miniconda3/envs/magnet/lib/python3.9/multiprocessing/queues.py\", line 364, in get\n",
      "    with self._rlock:\n",
      "  File \"/home/Users/yl181/miniconda3/envs/magnet/lib/python3.9/multiprocessing/synchronize.py\", line 95, in __enter__\n",
      "    return self._semlock.__enter__()\n",
      "  File \"/home/Users/yl181/miniconda3/envs/magnet/lib/python3.9/multiprocessing/queues.py\", line 364, in get\n",
      "    with self._rlock:\n",
      "  File \"/home/Users/yl181/miniconda3/envs/magnet/lib/python3.9/multiprocessing/synchronize.py\", line 95, in __enter__\n",
      "    return self._semlock.__enter__()\n",
      "  File \"/home/Users/yl181/miniconda3/envs/magnet/lib/python3.9/multiprocessing/synchronize.py\", line 95, in __enter__\n",
      "    return self._semlock.__enter__()\n",
      "  File \"/home/Users/yl181/miniconda3/envs/magnet/lib/python3.9/multiprocessing/queues.py\", line 364, in get\n",
      "    with self._rlock:\n",
      "  File \"/home/Users/yl181/miniconda3/envs/magnet/lib/python3.9/multiprocessing/queues.py\", line 364, in get\n",
      "    with self._rlock:\n",
      "  File \"/home/Users/yl181/miniconda3/envs/magnet/lib/python3.9/multiprocessing/synchronize.py\", line 95, in __enter__\n",
      "    return self._semlock.__enter__()\n",
      "  File \"/home/Users/yl181/miniconda3/envs/magnet/lib/python3.9/multiprocessing/synchronize.py\", line 95, in __enter__\n",
      "    return self._semlock.__enter__()\n",
      "  File \"/home/Users/yl181/miniconda3/envs/magnet/lib/python3.9/multiprocessing/queues.py\", line 364, in get\n",
      "    with self._rlock:\n",
      "  File \"/home/Users/yl181/miniconda3/envs/magnet/lib/python3.9/multiprocessing/queues.py\", line 364, in get\n",
      "    with self._rlock:\n",
      "  File \"/home/Users/yl181/miniconda3/envs/magnet/lib/python3.9/multiprocessing/queues.py\", line 364, in get\n",
      "    with self._rlock:\n",
      "  File \"/home/Users/yl181/miniconda3/envs/magnet/lib/python3.9/multiprocessing/queues.py\", line 364, in get\n",
      "    with self._rlock:\n",
      "  File \"/home/Users/yl181/miniconda3/envs/magnet/lib/python3.9/multiprocessing/synchronize.py\", line 95, in __enter__\n",
      "    return self._semlock.__enter__()\n",
      "  File \"/home/Users/yl181/miniconda3/envs/magnet/lib/python3.9/multiprocessing/synchronize.py\", line 95, in __enter__\n",
      "    return self._semlock.__enter__()\n",
      "  File \"/home/Users/yl181/miniconda3/envs/magnet/lib/python3.9/multiprocessing/synchronize.py\", line 95, in __enter__\n",
      "    return self._semlock.__enter__()\n",
      "  File \"/home/Users/yl181/miniconda3/envs/magnet/lib/python3.9/multiprocessing/synchronize.py\", line 95, in __enter__\n",
      "    return self._semlock.__enter__()\n",
      "  File \"/home/Users/yl181/miniconda3/envs/magnet/lib/python3.9/multiprocessing/queues.py\", line 364, in get\n",
      "    with self._rlock:\n",
      "  File \"/home/Users/yl181/miniconda3/envs/magnet/lib/python3.9/multiprocessing/queues.py\", line 364, in get\n",
      "    with self._rlock:\n",
      "  File \"/home/Users/yl181/miniconda3/envs/magnet/lib/python3.9/multiprocessing/synchronize.py\", line 95, in __enter__\n",
      "    return self._semlock.__enter__()\n",
      "  File \"/home/Users/yl181/miniconda3/envs/magnet/lib/python3.9/multiprocessing/synchronize.py\", line 95, in __enter__\n",
      "    return self._semlock.__enter__()\n",
      "  File \"/home/Users/yl181/miniconda3/envs/magnet/lib/python3.9/multiprocessing/synchronize.py\", line 95, in __enter__\n",
      "    return self._semlock.__enter__()\n",
      "  File \"/home/Users/yl181/miniconda3/envs/magnet/lib/python3.9/multiprocessing/queues.py\", line 364, in get\n",
      "    with self._rlock:\n",
      "  File \"/home/Users/yl181/miniconda3/envs/magnet/lib/python3.9/multiprocessing/synchronize.py\", line 95, in __enter__\n",
      "    return self._semlock.__enter__()\n",
      "  File \"/home/Users/yl181/miniconda3/envs/magnet/lib/python3.9/multiprocessing/synchronize.py\", line 95, in __enter__\n",
      "    return self._semlock.__enter__()\n",
      "  File \"/home/Users/yl181/miniconda3/envs/magnet/lib/python3.9/multiprocessing/synchronize.py\", line 95, in __enter__\n",
      "    return self._semlock.__enter__()\n",
      "  File \"/home/Users/yl181/miniconda3/envs/magnet/lib/python3.9/multiprocessing/synchronize.py\", line 95, in __enter__\n",
      "    return self._semlock.__enter__()\n",
      "  File \"/home/Users/yl181/miniconda3/envs/magnet/lib/python3.9/multiprocessing/synchronize.py\", line 95, in __enter__\n",
      "    return self._semlock.__enter__()\n",
      "  File \"/home/Users/yl181/miniconda3/envs/magnet/lib/python3.9/multiprocessing/synchronize.py\", line 95, in __enter__\n",
      "    return self._semlock.__enter__()\n",
      "  File \"/home/Users/yl181/miniconda3/envs/magnet/lib/python3.9/multiprocessing/synchronize.py\", line 95, in __enter__\n",
      "    return self._semlock.__enter__()\n",
      "  File \"/home/Users/yl181/miniconda3/envs/magnet/lib/python3.9/multiprocessing/synchronize.py\", line 95, in __enter__\n",
      "    return self._semlock.__enter__()\n",
      "  File \"/home/Users/yl181/miniconda3/envs/magnet/lib/python3.9/multiprocessing/synchronize.py\", line 95, in __enter__\n",
      "    return self._semlock.__enter__()\n",
      "  File \"/home/Users/yl181/miniconda3/envs/magnet/lib/python3.9/multiprocessing/synchronize.py\", line 95, in __enter__\n",
      "    return self._semlock.__enter__()\n",
      "KeyboardInterrupt\n",
      "  File \"/home/Users/yl181/miniconda3/envs/magnet/lib/python3.9/multiprocessing/synchronize.py\", line 95, in __enter__\n",
      "    return self._semlock.__enter__()\n",
      "  File \"/home/Users/yl181/miniconda3/envs/magnet/lib/python3.9/multiprocessing/synchronize.py\", line 95, in __enter__\n",
      "    return self._semlock.__enter__()\n",
      "  File \"/home/Users/yl181/miniconda3/envs/magnet/lib/python3.9/multiprocessing/synchronize.py\", line 95, in __enter__\n",
      "    return self._semlock.__enter__()\n",
      "  File \"/home/Users/yl181/miniconda3/envs/magnet/lib/python3.9/multiprocessing/synchronize.py\", line 95, in __enter__\n",
      "    return self._semlock.__enter__()\n",
      "KeyboardInterrupt\n",
      "KeyboardInterrupt\n",
      "  File \"/home/Users/yl181/miniconda3/envs/magnet/lib/python3.9/multiprocessing/synchronize.py\", line 95, in __enter__\n",
      "    return self._semlock.__enter__()\n",
      "KeyboardInterrupt\n",
      "  File \"/home/Users/yl181/miniconda3/envs/magnet/lib/python3.9/multiprocessing/synchronize.py\", line 95, in __enter__\n",
      "    return self._semlock.__enter__()\n",
      "KeyboardInterrupt\n",
      "KeyboardInterrupt\n",
      "  File \"/home/Users/yl181/miniconda3/envs/magnet/lib/python3.9/multiprocessing/synchronize.py\", line 95, in __enter__\n",
      "    return self._semlock.__enter__()\n",
      "  File \"/home/Users/yl181/miniconda3/envs/magnet/lib/python3.9/multiprocessing/synchronize.py\", line 95, in __enter__\n",
      "    return self._semlock.__enter__()\n",
      "KeyboardInterrupt\n",
      "KeyboardInterrupt\n",
      "  File \"/home/Users/yl181/miniconda3/envs/magnet/lib/python3.9/multiprocessing/synchronize.py\", line 95, in __enter__\n",
      "    return self._semlock.__enter__()\n",
      "  File \"/home/Users/yl181/miniconda3/envs/magnet/lib/python3.9/multiprocessing/synchronize.py\", line 95, in __enter__\n",
      "    return self._semlock.__enter__()\n",
      "  File \"/home/Users/yl181/miniconda3/envs/magnet/lib/python3.9/multiprocessing/synchronize.py\", line 95, in __enter__\n",
      "    return self._semlock.__enter__()\n",
      "  File \"/home/Users/yl181/miniconda3/envs/magnet/lib/python3.9/multiprocessing/synchronize.py\", line 95, in __enter__\n",
      "    return self._semlock.__enter__()\n",
      "KeyboardInterrupt\n",
      "KeyboardInterrupt\n",
      "KeyboardInterrupt\n",
      "KeyboardInterrupt\n",
      "  File \"/home/Users/yl181/miniconda3/envs/magnet/lib/python3.9/multiprocessing/synchronize.py\", line 95, in __enter__\n",
      "    return self._semlock.__enter__()\n",
      "  File \"/home/Users/yl181/miniconda3/envs/magnet/lib/python3.9/multiprocessing/synchronize.py\", line 95, in __enter__\n",
      "    return self._semlock.__enter__()\n",
      "KeyboardInterrupt\n",
      "KeyboardInterrupt\n",
      "KeyboardInterrupt\n",
      "  File \"/home/Users/yl181/miniconda3/envs/magnet/lib/python3.9/multiprocessing/synchronize.py\", line 95, in __enter__\n",
      "    return self._semlock.__enter__()\n",
      "KeyboardInterrupt\n",
      "KeyboardInterrupt\n",
      "KeyboardInterrupt\n",
      "KeyboardInterrupt\n",
      "KeyboardInterrupt\n",
      "KeyboardInterrupt\n",
      "KeyboardInterrupt\n",
      "KeyboardInterrupt\n",
      "KeyboardInterrupt\n",
      "KeyboardInterrupt\n",
      "KeyboardInterrupt\n",
      "KeyboardInterrupt\n",
      "KeyboardInterrupt\n",
      "KeyboardInterrupt\n",
      "KeyboardInterrupt\n",
      "KeyboardInterrupt\n",
      "KeyboardInterrupt\n",
      "KeyboardInterrupt\n",
      "KeyboardInterrupt\n",
      "KeyboardInterrupt\n",
      "KeyboardInterrupt\n",
      "KeyboardInterrupt\n",
      "KeyboardInterrupt\n",
      "KeyboardInterrupt\n",
      "KeyboardInterrupt\n"
     ]
    }
   ],
   "source": [
    "for cmd_args in cmd_records:\n",
    "    run_magnet(cmd_args)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6a8b5e56-f23f-40ae-97e0-4140cb0c4191",
   "metadata": {},
   "source": [
    "## select zymo sample"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 101,
   "id": "f603751c-c87b-402e-b40a-388d07351b9d",
   "metadata": {},
   "outputs": [],
   "source": [
    "mob_output_path = '/home/Users/ns58/Mob-experiments/Mob-outputs/Zymo-EVEN-subsample'\n",
    "input_fastq_path = '/home/Users/ns58/Mob-experiments/Data-grouped/Zymo-EVEN-p0.01'\n",
    "magnet_output = '/home/Users/yl181/memu/dev_magnet_v222_zymo_even_strains'\n",
    "\n",
    "if not os.path.exists(magnet_output):\n",
    "    os.mkdir(magnet_output)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 102,
   "id": "3557ac31-bfd3-46b2-8640-bf245745205b",
   "metadata": {},
   "outputs": [],
   "source": [
    "records = []\n",
    "for f in os.listdir(mob_output_path):\n",
    "    if f.startswith('Zymo-EVEN-') and 'relative_abundance.tsv' in f:\n",
    "        sample_id = f[:13]\n",
    "        input_fastq = os.path.join(input_fastq_path, f'{sample_id}.fastq')\n",
    "        if not os.path.exists(input_fastq):\n",
    "            print(sample_id)\n",
    "        classification_output = os.path.join(mob_output_path, f)\n",
    "        if not os.path.exists(classification_output):\n",
    "            print(sample_id)\n",
    "            \n",
    "        strain_mob_output = os.path.join(mob_output_path, f'{sample_id}_P_rgs_df.tsv')\n",
    "        \n",
    "        if sample_id == 'Zymo-EVEN-004':\n",
    "            break"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "08252c0f-5aa8-4406-893e-8bfcd3d82fb3",
   "metadata": {},
   "source": [
    "## select zymo gut sample"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 103,
   "id": "bf49c92f-4640-4307-a9d9-e62648b5d38d",
   "metadata": {},
   "outputs": [],
   "source": [
    "mob_output_path = '/home/Users/ns58/Mob-experiments/Mob-outputs/Zymo-GUT'\n",
    "input_fastq_path = '/home/Users/ns58/Mob-experiments/Data-grouped/Zymo-GUT'\n",
    "magnet_output = '/home/Users/yl181/memu/dev_magnet_v222_zymo_gut'\n",
    "\n",
    "if not os.path.exists(magnet_output):\n",
    "    os.mkdir(magnet_output)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 104,
   "id": "531f54ac-2314-498d-bf4c-23d35bf6645b",
   "metadata": {},
   "outputs": [],
   "source": [
    "records = []\n",
    "for f in os.listdir(mob_output_path):\n",
    "    if f.startswith('nof-Zymo-gut-') and 'relative_abundance.tsv' in f:\n",
    "        sample_id = f[13]\n",
    "        input_fastq = os.path.join(input_fastq_path, f'Zymo-GUT-00{sample_id}.fastq')\n",
    "        if not os.path.exists(input_fastq):\n",
    "            print(sample_id)\n",
    "        classification_output = os.path.join(mob_output_path, f)\n",
    "        if not os.path.exists(classification_output):\n",
    "            print(sample_id)\n",
    "            \n",
    "        strain_mob_output = os.path.join(mob_output_path, f'{sample_id}_P_rgs_df.tsv')\n",
    "        \n",
    "        if sample_id == '6':\n",
    "            print(sample_id)\n",
    "            break"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "23b34b00-8b4f-4a06-9792-9eb2cb5776a9",
   "metadata": {},
   "source": [
    "## Soil"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 105,
   "id": "182f521c-7ae3-4a61-a8fe-0a6c25fbe547",
   "metadata": {},
   "outputs": [],
   "source": [
    "mob_output_path = '/home/Users/ns58/Mob-experiments/Mob-outputs/Soil'\n",
    "melon_output_path = '/home/Users/ns58/Mob-experiments/Melon-outputs/Soil'\n",
    "input_fastq_path = '/home/Users/ns58/Mob-experiments/Data-grouped/Soil'\n",
    "magnet_output = '/home/Users/yl181/memu/dev_magnet_v222_soil'\n",
    "\n",
    "if not os.path.exists(magnet_output):\n",
    "    os.mkdir(magnet_output)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b7671ccd-efa6-4858-b859-522a1c708efe",
   "metadata": {},
   "source": [
    "### Mob"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 106,
   "id": "834a3b7a-03f1-46ee-8a76-4f8c79d175eb",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Soil-5\n"
     ]
    }
   ],
   "source": [
    "for f in os.listdir(mob_output_path):\n",
    "    if 'relative_abundance.tsv' in f:\n",
    "        sample_id = \"-\".join(f.split('-')[0:2])\n",
    "        input_fastq = os.path.join(input_fastq_path, f'{sample_id}.fastq')\n",
    "        if not os.path.exists(input_fastq):\n",
    "            print(\"input_fastq missing:\", sample_id)\n",
    "        classification_output = os.path.join(mob_output_path, f)\n",
    "        if not os.path.exists(classification_output):\n",
    "            print('classification_output missing:', sample_id)\n",
    "            \n",
    "        strain_mob_output = os.path.join(mob_output_path, f'{sample_id}_P_rgs_df.tsv')\n",
    "        \n",
    "        if sample_id == 'Soil-5':\n",
    "            print(sample_id)\n",
    "            break"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 107,
   "id": "c447c24d-21cb-4c88-a0b6-89b062e3df5a",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "'-c /home/Users/ns58/Mob-experiments/Mob-outputs/Soil/Soil-5-relative_abundance.tsv -i /home/Users/ns58/Mob-experiments/Data-grouped/Soil/Soil-5.fastq -m ont -a 12 --min-abundance 0.0 -o /home/Users/yl181/memu/dev_magnet_v222_soil/Soil-5_mob --threads 40'"
      ]
     },
     "execution_count": 107,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "cmd_args = ['-c', classification_output,\n",
    "            '-i', input_fastq,\n",
    "            '-m', 'ont',\n",
    "            '-a', str(12),\n",
    "            '--min-abundance', str(0.000),\n",
    "            '-o', os.path.join(magnet_output, f\"{sample_id}_mob\"),\n",
    "            '--threads', '40']\n",
    "\" \".join(cmd_args)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "255bb533-b388-461a-b7af-e7915ed1c138",
   "metadata": {},
   "source": [
    "### Melon"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 69,
   "id": "b09bf86b-d179-4606-a970-34e901c5fb9f",
   "metadata": {},
   "outputs": [],
   "source": [
    "melon_output_f = os.path.join(melon_output_path, sample_id, f'{sample_id}.tsv')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 70,
   "id": "aa70c795-29d0-429d-afc5-f3a64453885d",
   "metadata": {},
   "outputs": [],
   "source": [
    "melon_taxids = []\n",
    "for item in pd.read_csv(melon_output_f, sep='\\t')['species'].values:\n",
    "    melon_taxids.append(item.split('|')[0])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 71,
   "id": "b344088a-b74b-4f56-a7e8-46bb5a6c16e7",
   "metadata": {},
   "outputs": [],
   "source": [
    "with open(os.path.join(magnet_output, f'{sample_id}_melon_taxids.txt'), 'w') as melon_taxid_f:\n",
    "    for item in melon_taxids:\n",
    "        melon_taxid_f.write(f\"{item}\\n\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 72,
   "id": "900db572-1d30-4bd3-a37a-31fd2a850720",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "'-c /home/Users/yl181/memu/dev_magnet_v222_soil/Soil-5_melon_taxids.txt -i /home/Users/ns58/Mob-experiments/Data-grouped/Soil/Soil-5.fastq -m ont -o /home/Users/yl181/memu/dev_magnet_v222_soil/Soil-5_melon --threads 40'"
      ]
     },
     "execution_count": 72,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "cmd_args = ['-c', os.path.join(magnet_output, f'{sample_id}_melon_taxids.txt'),\n",
    "            '-i', input_fastq,\n",
    "            '-m', 'ont',\n",
    "            '-o', os.path.join(magnet_output, f'{sample_id}_melon'),\n",
    "            '--threads', '40']\n",
    "\" \".join(cmd_args)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f0ce7936-46ce-484f-a3ca-beb61ccdce0a",
   "metadata": {},
   "source": [
    "## Strain magnet"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "cb14db0d-1970-4d5d-8cfb-264b1dfb161e",
   "metadata": {},
   "outputs": [],
   "source": [
    "# top_k = 5"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "38d2cad8-238c-47ee-a883-e2d88921e5fc",
   "metadata": {},
   "outputs": [],
   "source": [
    "# valid_species_df = pd.read_csv(classification_output, sep='\\t')\n",
    "# valid_species_df = valid_species_df[valid_species_df['F'] > 0.00]\n",
    "# valid_species_taxids = set(valid_species_df['Target_ID'].values)\n",
    "\n",
    "# strain_df = pd.read_csv(strain_mob_output, sep='\\t')\n",
    "# filtered_strain_df = strain_df[strain_df['Target_ID'].isin(valid_species_taxids)]\n",
    "\n",
    "# taxid2accession_record = []\n",
    "# for idx, row in filtered_strain_df.iterrows():\n",
    "#     genome_list = literal_eval(row['Genome'])\n",
    "#     for genome in genome_list:\n",
    "#         accession = genome.split('|')[0]\n",
    "#         taxid2accession_record.append({'taxid': row['Target_ID'],\n",
    "#                                        'accession': accession})\n",
    "# taxid2accession_df = pd.DataFrame(taxid2accession_record)\n",
    "# taxid2accession_df = taxid2accession_df.groupby(['taxid', 'accession']).size().reset_index(name='counts')\n",
    "\n",
    "# selected_genomes = []\n",
    "# for taxid in valid_species_taxids:\n",
    "#     selected_df = taxid2accession_df[taxid2accession_df['taxid'] == taxid].sort_values(['counts'], ascending=False)\n",
    "#     # selected_genomes += list(selected_df.head(top_k)['accession'].values)\n",
    "#     min_count = selected_df.iloc[0]['counts']\n",
    "#     if len(selected_df[selected_df['counts'] >= min_count]) >= top_k:\n",
    "#         candidate_selected_genomes = list(selected_df[selected_df['counts'] >= min_count]['accession'].values)\n",
    "#         # if len(candidate_selected_genomes) > top_k:\n",
    "#         #     print(taxid, candidate_selected_genomes, set(candidate_selected_genomes).intersection(ground_truth))\n",
    "#     else:\n",
    "#         candidate_selected_genomes = list(selected_df.head(top_k)['accession'].values)\n",
    "\n",
    "#     for genome in candidate_selected_genomes:\n",
    "#         selected_genomes.append({\n",
    "#             'assembly': genome,\n",
    "#             'taxid': taxid\n",
    "#         })\n",
    "        \n",
    "# len(selected_genomes)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "56761c3b-22b3-48ab-a01e-bdad39ca6d37",
   "metadata": {},
   "outputs": [],
   "source": [
    "# pd.DataFrame(selected_genomes).set_index('assembly').to_csv(f'{sample_id}.top{top_k}-strains.csv')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "6765c1f1-2b7b-4aa2-a355-5ada544c42a9",
   "metadata": {},
   "outputs": [],
   "source": [
    "# cmd_args = ['-c', f'{sample_id}.top{top_k}-strains.csv',\n",
    "#             '-i', input_fastq,\n",
    "#             '-m', 'ont',\n",
    "#             '-o', os.path.join(magnet_output, sample_id),\n",
    "#             '--threads', '40',\n",
    "#             '--accession']"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f6c82ff4-d086-4af6-9feb-04adb7d7a0b5",
   "metadata": {},
   "source": [
    "## construct parameters"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "id": "013c11a6-4476-486f-a40f-908c57cc7c77",
   "metadata": {},
   "outputs": [],
   "source": [
    "args = parser.parse_args(cmd_args)\n",
    "\n",
    "input_tsv = args.classification\n",
    "input_fastq = args.fastq\n",
    "input_fastq2 = args.fastq2\n",
    "mode = args.mode\n",
    "working_directory = args.output\n",
    "\n",
    "taxid_col_idx = args.taxid_idx\n",
    "abundance_col_idx = args.abundance_idx\n",
    "min_abundance = args.min_abundance\n",
    "min_mapq = args.min_mapq\n",
    "min_coverage_score = args.min_covscore\n",
    "threads = args.threads\n",
    "valid_kingdom_str = args.kingdom\n",
    "\n",
    "valid_kingdom = set()\n",
    "for i in valid_kingdom_str.split(','):\n",
    "    valid_kingdom.add(int(i))\n",
    "\n",
    "if args.include_mag:\n",
    "    mag_flag = 'all'\n",
    "else:\n",
    "    mag_flag = 'exclude'\n",
    "\n",
    "accession_flag = args.accession\n",
    "call_subspecies = args.subspecies\n",
    "\n",
    "sep = '\\t'\n",
    "if str(input_tsv)[-3:] == 'csv':\n",
    "    sep = ','\n",
    "\n",
    "ncbi_taxa_db = NCBITaxa()\n",
    "if not os.path.exists(working_directory):\n",
    "    os.mkdir(working_directory)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "354da8a9-af00-4105-875f-94f579fda714",
   "metadata": {},
   "source": [
    "# Refseq Downloader"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "id": "260a1b51-5faa-4704-9c18-396cdafd5a17",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "min_abundance: 0.0\n",
      "num_species: 82\n"
     ]
    }
   ],
   "source": [
    "if accession_flag:\n",
    "    abundance_col_idx = None\n",
    "    min_abundance = 0\n",
    "    input_df, min_abundance = parsing_input_f(input_tsv, sep, taxid_col_idx, abundance_col_idx, min_abundance)\n",
    "    valid_taxids = list(input_df['tax_id'].values)\n",
    "else:\n",
    "    input_df, min_abundance = parsing_input_f(input_tsv, sep, taxid_col_idx, abundance_col_idx, min_abundance)\n",
    "    # make valid_kingdom a variable?\n",
    "    valid_taxids = filter_input_df(input_df, min_abundance, ncbi_taxa_db, valid_kingdom=valid_kingdom, ret_subspecies=call_subspecies)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2926b442-e94d-42e9-909a-c6bdc112dbf6",
   "metadata": {},
   "source": [
    "## ! Missing genome need to be informed."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "id": "b8af64cd-7cfb-4c96-a12e-6d4331ec149b",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "2923522    \t GCF_031142765.1 \t Enterococcus sp. IsoGale005    \t IsoGale0.. \t Contig\n",
      "2555396    \t GCF_013350625.1 \t Salmonella sp. SG220           \t SG220      \t Contig\n",
      "1715212    \t GCF_001837115.1 \t Staphylococcus sp. HMSC063H1.. \t HMSC063H.. \t Scaffold\n",
      "562        \t GCF_000008865.2 \t Escherichia coli O157:H7 str.. \t Sakai su.. \t Complete Genome\n",
      "1739314    \t GCF_001811595.1 \t Enterococcus sp. HMSC078F03    \t HMSC078F.. \t Scaffold\n",
      "2315836    \t GCF_008082385.1 \t Enterococcus sp. N041.A-2      \t N041.A-2   \t Contig\n",
      "2502216    \t GCF_005503545.1 \t Pseudomonas sp. 3PA37B6        \t 3PA37B6    \t Contig\n",
      "1613       \t GCF_029961225.1 \t Limosilactobacillus fermentu.. \t EFEL6800   \t Complete Genome\n",
      "2665553    \t GCF_014492625.1 \t Salmonella sp. S048_01045      \t S048_010.. \t Contig\n",
      "5207       \t GCA_020975435.1 \t Cryptococcus neoformans        \t C36        \t Complete Genome\n",
      "2582616    \t GCF_009902665.1 \t Salmonella sp. zj-f50          \t zj-f50     \t Scaffold\n",
      "2665564    \t GCF_014492365.1 \t Salmonella sp. S060_01291      \t S060_012.. \t Contig\n",
      "2665570    \t GCF_014492255.1 \t Salmonella sp. S069_01685      \t S069_016.. \t Contig\n",
      "1639       \t GCF_000196035.1 \t Listeria monocytogenes EGD-e   \t EGD-e      \t Complete Genome\n",
      "621        \t GCF_001027225.1 \t Shigella boydii                \t ATCC 921.. \t Complete Genome\n",
      "622        \t GCF_022354085.1 \t Shigella dysenteriae           \t SWHEFF_4.. \t Complete Genome\n",
      "623        \t GCF_000006925.2 \t Shigella flexneri 2a str. 30.. \t 301        \t Complete Genome\n",
      "624        \t GCF_013374815.1 \t Shigella sonnei                \t SE6-1      \t Complete Genome\n",
      "2665594    \t GCF_014491785.1 \t Salmonella sp. S102_03650      \t S102_036.. \t Contig\n",
      "2837629    \t Genome Not Found.\n",
      "2767485    \t GCF_017152085.1 \t Staphylococcus sp. EG-SA-2     \t EG-SA-2    \t Contig\n",
      "2665598    \t GCF_014491695.1 \t Salmonella sp. S108-09120      \t S108-091.. \t Contig\n",
      "2767487    \t GCF_017152015.1 \t Staphylococcus sp. EG-SA-8     \t EG-SA-8    \t Scaffold\n",
      "2767493    \t GCF_017151845.1 \t Staphylococcus sp. EG-SA-17    \t EG-SA-17   \t Scaffold\n",
      "2767494    \t GCF_017151825.1 \t Staphylococcus sp. EG-SA-18    \t EG-SA-18   \t Scaffold\n",
      "2767496    \t GCF_017151875.1 \t Staphylococcus sp. EG-SA-21    \t EG-SA-21   \t Scaffold\n",
      "2767497    \t GCF_017151775.1 \t Staphylococcus sp. EG-SA-23    \t EG-SA-23   \t Scaffold\n",
      "2767499    \t GCF_017151735.1 \t Staphylococcus sp. EG-SA-27    \t EG-SA-27   \t Contig\n",
      "1739405    \t GCF_001813275.1 \t Enterococcus sp. HMSC073E09    \t HMSC073E.. \t Scaffold\n",
      "2767501    \t GCF_017151725.1 \t Staphylococcus sp. EG-SA-32    \t EG-SA-32   \t Scaffold\n",
      "1608872    \t GCF_001838925.1 \t Staphylococcus sp. HMSC36F05   \t HMSC36F0.. \t Scaffold\n",
      "2819762    \t GCF_017572025.1 \t Salmonella sp. 3202011420180.. \t 32020114.. \t Scaffold\n",
      "2819778    \t GCF_017571715.1 \t Salmonella sp. 32040203-2019.. \t 32040203.. \t Scaffold\n",
      "237777     \t Genome Not Found.\n",
      "28901      \t GCF_000006945.2 \t Salmonella enterica subsp. e.. \t LT2        \t Complete Genome\n",
      "1280       \t GCF_000013425.1 \t Staphylococcus aureus subsp... \t NCTC 832.. \t Complete Genome\n",
      "2678528    \t GCF_028867355.1 \t Listeria sp. LM90SB2           \t LM90SB2    \t Complete Genome\n",
      "2608391    \t GCF_008385845.1 \t Staphylococcus sp. 775         \t 775        \t Scaffold\n",
      "2608392    \t GCF_008385815.1 \t Staphylococcus sp. 5313        \t 5313       \t Contig\n",
      "2608394    \t GCF_008385875.1 \t Staphylococcus sp. 7716        \t 7716       \t Contig\n",
      "2846989    \t GCF_018920115.1 \t Enterococcus sp. S115_ASV_20   \t S115_ASV.. \t Contig\n",
      "2608399    \t GCF_008385975.1 \t Staphylococcus sp. 52717       \t 52717      \t Contig\n",
      "2608400    \t GCF_008385995.1 \t Staphylococcus sp. 53017       \t 53017      \t Scaffold\n",
      "2846991    \t GCF_018920065.1 \t Enterococcus sp. S123_ASV_20   \t S123_ASV.. \t Contig\n",
      "2846994    \t GCF_018919985.1 \t Enterococcus sp. S131_ASV_20   \t S131_ASV.. \t Contig\n",
      "2847001    \t GCF_018919865.1 \t Enterococcus sp. S159_ASV_20   \t S159_ASV.. \t Contig\n",
      "2847003    \t GCF_018919805.1 \t Enterococcus sp. S163_ASV_20   \t S163_ASV.. \t Contig\n",
      "2847005    \t GCF_018919785.1 \t Enterococcus sp. S171_ASV_20   \t S171_ASV.. \t Contig\n",
      "287        \t GCF_000006765.1 \t Pseudomonas aeruginosa PAO1    \t PAO1       \t Complete Genome\n",
      "2847009    \t GCF_018919685.1 \t Enterococcus sp. S183_ASV_20   \t S183_ASV.. \t Contig\n",
      "3061043    \t GCF_030499395.1 \t Enterococcus sp. A1(2023)      \t A1(2023)   \t Contig\n",
      "3061047    \t GCF_030499465.1 \t Enterococcus sp. A2(2023)      \t A2(2023)   \t Contig\n",
      "4932       \t GCF_000146045.2 \t Saccharomyces cerevisiae S28.. \t S288C      \t Complete Genome\n",
      "1351       \t GCF_001598635.1 \t Enterococcus faecalis          \t LD33       \t Complete Genome\n",
      "2807624    \t GCF_016917335.1 \t Staphylococcus sp. SM9054      \t SM9054     \t Complete Genome\n",
      "2908490    \t GCF_022556975.1 \t Salmonella sp. YP54-1          \t YP54-1     \t Contig\n",
      "2908492    \t GCF_022556915.1 \t Salmonella sp. YP67-1          \t YP67-1     \t Contig\n",
      "1579344    \t GCF_001056095.1 \t Enterococcus sp. 255_ESPC      \t 255_ESPC   \t Scaffold\n",
      "1579349    \t GCF_001053935.1 \t Enterococcus sp. 1140_ESPC     \t 1140_ESP.. \t Contig\n",
      "2908506    \t GCF_022556635.1 \t Salmonella sp. YP256-1         \t YP256-1    \t Contig\n",
      "2908508    \t GCF_022556555.1 \t Salmonella sp. YT2-4-2         \t YT2-4-2    \t Contig\n",
      "2908511    \t GCF_022556535.1 \t Salmonella sp. YT2-20-3        \t YT2-20-3   \t Contig\n",
      "1715097    \t GCF_001836155.1 \t Staphylococcus sp. HMSC035F1.. \t HMSC035F.. \t Scaffold\n",
      "2804146    \t GCF_017892245.1 \t Staphylococcus sp. GDY8P96A-.. \t GDY8P96A.. \t Contig\n",
      "1715150    \t GCF_001836465.1 \t Staphylococcus sp. HMSC058E0.. \t HMSC058E.. \t Scaffold\n",
      "2923488    \t GCF_031144525.1 \t Enterococcus sp. FR008         \t FR008      \t Contig\n",
      "83428      \t Genome Not Found.\n",
      "2923492    \t GCF_031144375.1 \t Enterococcus sp. FR042         \t FR042      \t Contig\n",
      "2920934    \t GCF_023539285.1 \t Enterococcus sp. LX10          \t LX10       \t Complete Genome\n",
      "2923493    \t GCF_031144365.1 \t Enterococcus sp. FR062         \t FR062      \t Contig\n",
      "2923494    \t GCF_031144325.1 \t Enterococcus sp. FR063         \t FR063      \t Contig\n",
      "2923495    \t GCF_031144215.1 \t Enterococcus sp. FR068         \t FR068      \t Contig\n",
      "2923496    \t GCF_031144155.1 \t Enterococcus sp. FR088         \t FR088      \t Contig\n",
      "2923497    \t GCF_031144105.1 \t Enterococcus sp. FR111         \t FR111      \t Contig\n",
      "2923498    \t GCF_031144145.1 \t Enterococcus sp. FR112         \t FR112      \t Contig\n",
      "2923500    \t GCF_031143985.1 \t Enterococcus sp. FR128         \t FR128      \t Contig\n",
      "2923502    \t GCF_031143855.1 \t Enterococcus sp. FR143         \t FR143      \t Contig\n",
      "2923503    \t GCF_031143745.1 \t Enterococcus sp. FR159         \t FR159      \t Contig\n",
      "96241      \t GCF_006094475.1 \t Bacillus spizizenii ATCC 663.. \t ATCC 663.. \t Complete Genome\n",
      "2923507    \t GCF_031143485.1 \t Enterococcus sp. FR180         \t FR180      \t Contig\n",
      "2923508    \t GCF_031143475.1 \t Enterococcus sp. FR191         \t FR191      \t Contig\n",
      "2923509    \t GCF_031143435.1 \t Enterococcus sp. FR192         \t FR192      \t Contig\n"
     ]
    }
   ],
   "source": [
    "reference_metadata = prepare_reference_genomes(valid_taxids, working_directory, ncbi_taxa_db, accession_flag=accession_flag, mag_flag=mag_flag)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "id": "2417ead8-e8cb-45e0-bab9-dda7f166fec3",
   "metadata": {},
   "outputs": [],
   "source": [
    "downloaded_assemblies = reference_metadata[reference_metadata['Downloaded']]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9c8c1ab4-d832-45db-85cb-a06a1c20a7dc",
   "metadata": {},
   "source": [
    "# FastANI and Clustering"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "410b3668-ddc6-4dcd-9918-bb2d8b093104",
   "metadata": {},
   "source": [
    "## Run FastANI, Clustering, and Choose Representative Genomes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "id": "d1461843-93c6-4c1f-968d-653c1073762e",
   "metadata": {},
   "outputs": [],
   "source": [
    "reference_genome_path = os.path.join(working_directory, 'reference_genomes')\n",
    "\n",
    "fastani_path = os.path.join(working_directory, 'fastANI')\n",
    "if not os.path.exists(fastani_path):\n",
    "    os.mkdir(fastani_path)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "046a24e2-82a2-4c19-8dcd-5299d338f9e6",
   "metadata": {},
   "source": [
    "## With-in genus ANI Clustering (deprecated)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "id": "848827cc-b7b6-4fd5-867a-22a2862d8aa8",
   "metadata": {},
   "outputs": [],
   "source": [
    "# genus_taxids = []\n",
    "# for idx, row in downloaded_assemblies.iterrows():\n",
    "#     species_taxid = row['Taxonomy ID']\n",
    "#     genus_taxid = get_genus_taxid(species_taxid, ncbi_taxa_db, valid_kingdom)\n",
    "#     genus_taxids.append(genus_taxid)\n",
    "# downloaded_assemblies['Genus Taxid'] = genus_taxids"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "id": "3f61ac71-67f7-4c14-9fa6-3bf542ffe607",
   "metadata": {},
   "outputs": [],
   "source": [
    "# def get_genus_taxid(taxid, ncbi_taxa_db, valid_kingdom):    \n",
    "#     try:\n",
    "#         lineage = ncbi_taxa_db.get_lineage(taxid)\n",
    "#     except ValueError:\n",
    "#         return 0 \n",
    "    \n",
    "#     if bool(set(lineage) & valid_kingdom):\n",
    "#         taxid2rank_dict = ncbi_taxa_db.get_rank(lineage)\n",
    "#         for lineage_taxid in lineage:\n",
    "#             if taxid2rank_dict[lineage_taxid] == 'genus':\n",
    "#                 return lineage_taxid\n",
    "            \n",
    "#     return 0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "id": "b5b67190-98d3-45d2-a88b-934ad4bf514e",
   "metadata": {},
   "outputs": [],
   "source": [
    "# genus_count_df = pd.DataFrame(downloaded_assemblies.groupby(['Genus Taxid'])['Taxonomy ID'].count())\n",
    "\n",
    "# multi_species_genera = genus_count_df[genus_count_df['Taxonomy ID'] > 1].index.to_list()\n",
    "# single_species_genera = genus_count_df[genus_count_df['Taxonomy ID'] == 1].index.to_list()+[0]\n",
    "\n",
    "# genus2represetitive_mapping = dict()\n",
    "\n",
    "# for genus_taxid in multi_species_genera:\n",
    "#     if genus_taxid != 0:\n",
    "#         fastani_assemblies = downloaded_assemblies[downloaded_assemblies['Genus Taxid'] == genus_taxid]['Assembly Accession ID'].values\n",
    "#         genus_assemblie_list_f = os.path.join(fastani_path, f\"genus_{genus_taxid}.txt\")\n",
    "#         with open(genus_assemblie_list_f, \"w\") as rl_f:\n",
    "#             for accession in fastani_assemblies:\n",
    "#                 reference_genome = os.path.join(reference_genome_path, f'{accession}.fasta')\n",
    "#                 rl_f.write(f\"{reference_genome}\\n\")\n",
    "\n",
    "#         subprocess.run(['fastANI',\n",
    "#                         '--rl', genus_assemblie_list_f,\n",
    "#                         '--ql', genus_assemblie_list_f,\n",
    "#                         \"--matrix\",\n",
    "#                         '-o', os.path.join(fastani_path, f'genus_{genus_taxid}')],\n",
    "#                        check=True,\n",
    "#                        stdout=open(os.path.join(fastani_path, \"fastani.log\"), \"a\"),\n",
    "#                        stderr=open(os.path.join(fastani_path, \"fastani.err\"), \"a\"))\n",
    "        \n",
    "#         representative_genomes = find_representative_within_genus(genus_taxid, fastani_path, downloaded_assemblies)\n",
    "#         genus2represetitive_mapping[genus_taxid] = representative_genomes\n",
    "\n",
    "# representative_labels = []\n",
    "# for idx, row in downloaded_assemblies.iterrows():\n",
    "#     genus_idx = row['Genus Taxid']\n",
    "#     if genus_idx in single_species_genera:\n",
    "#         representative_labels.append(True)\n",
    "#     elif row['Assembly Accession ID'] in genus2represetitive_mapping[genus_idx]:\n",
    "#         representative_labels.append(True)\n",
    "#     else:\n",
    "#         representative_labels.append(False)\n",
    "        \n",
    "# downloaded_assemblies['Cluster Representative'] = representative_labels"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9b612105-0c83-4373-a4c9-a98aebf97797",
   "metadata": {},
   "source": [
    "## Cross Genus ANI Clustering"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "id": "9615efde-31eb-4769-92c9-4426d7f34308",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "CompletedProcess(args=['fastANI', '--rl', '/home/Users/yl181/memu/magnet_v222_zymo_even_full/Zymo-EVEN-p0.75/Zymo-Even-3/fastANI/assemblie_list.txt', '--ql', '/home/Users/yl181/memu/magnet_v222_zymo_even_full/Zymo-EVEN-p0.75/Zymo-Even-3/fastANI/assemblie_list.txt', '--matrix', '--threads', '40', '-o', '/home/Users/yl181/memu/magnet_v222_zymo_even_full/Zymo-EVEN-p0.75/Zymo-Even-3/fastANI/pairwise_ani'], returncode=0)"
      ]
     },
     "execution_count": 40,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "fastani_assemblies = downloaded_assemblies['Assembly Accession ID'].values\n",
    "assemblie_list_f = os.path.join(fastani_path, f\"assemblie_list.txt\")\n",
    "with open(assemblie_list_f, \"w\") as rl_f:\n",
    "    for accession in fastani_assemblies:\n",
    "        reference_genome = os.path.join(reference_genome_path, f'{accession}.fasta')\n",
    "        rl_f.write(f\"{reference_genome}\\n\")\n",
    "        \n",
    "subprocess.run(['fastANI',\n",
    "                '--rl', assemblie_list_f,\n",
    "                '--ql', assemblie_list_f,\n",
    "                \"--matrix\",\n",
    "                '--threads', str(threads),\n",
    "                '-o', os.path.join(fastani_path, f'pairwise_ani')],\n",
    "               check=True,\n",
    "               stdout=open(os.path.join(fastani_path, \"fastani.log\"), \"a\"),\n",
    "               stderr=open(os.path.join(fastani_path, \"fastani.err\"), \"a\"))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "id": "f415383c-a93c-4157-b22b-b2e5507ff4d0",
   "metadata": {},
   "outputs": [],
   "source": [
    "representative_genomes, member2representative = find_representative_genome(fastani_path, fastani_assemblies, downloaded_assemblies)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "id": "693b6c50-e691-4cbb-9007-7b9d772a836e",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "defaultdict(list,\n",
       "            {'GCF_001598635.1': ['GCF_031142765.1',\n",
       "              'GCF_001811595.1',\n",
       "              'GCF_008082385.1',\n",
       "              'GCF_001813275.1',\n",
       "              'GCF_018920115.1',\n",
       "              'GCF_018920065.1',\n",
       "              'GCF_018919985.1',\n",
       "              'GCF_018919865.1',\n",
       "              'GCF_018919805.1',\n",
       "              'GCF_018919785.1',\n",
       "              'GCF_018919685.1',\n",
       "              'GCF_030499395.1',\n",
       "              'GCF_030499465.1',\n",
       "              'GCF_001598635.1',\n",
       "              'GCF_001056095.1',\n",
       "              'GCF_001053935.1',\n",
       "              'GCF_031144525.1',\n",
       "              'GCF_031144375.1',\n",
       "              'GCF_023539285.1',\n",
       "              'GCF_031144365.1',\n",
       "              'GCF_031144325.1',\n",
       "              'GCF_031144215.1',\n",
       "              'GCF_031144155.1',\n",
       "              'GCF_031144105.1',\n",
       "              'GCF_031144145.1',\n",
       "              'GCF_031143985.1',\n",
       "              'GCF_031143855.1',\n",
       "              'GCF_031143745.1',\n",
       "              'GCF_031143485.1',\n",
       "              'GCF_031143475.1',\n",
       "              'GCF_031143435.1'],\n",
       "             'GCF_000006945.2': ['GCF_013350625.1',\n",
       "              'GCF_014492625.1',\n",
       "              'GCF_009902665.1',\n",
       "              'GCF_014492365.1',\n",
       "              'GCF_014492255.1',\n",
       "              'GCF_014491785.1',\n",
       "              'GCF_014491695.1',\n",
       "              'GCF_017572025.1',\n",
       "              'GCF_017571715.1',\n",
       "              'GCF_000006945.2',\n",
       "              'GCF_022556975.1',\n",
       "              'GCF_022556915.1',\n",
       "              'GCF_022556635.1',\n",
       "              'GCF_022556555.1',\n",
       "              'GCF_022556535.1'],\n",
       "             'GCF_000013425.1': ['GCF_001837115.1',\n",
       "              'GCF_017152085.1',\n",
       "              'GCF_017152015.1',\n",
       "              'GCF_017151845.1',\n",
       "              'GCF_017151825.1',\n",
       "              'GCF_017151875.1',\n",
       "              'GCF_017151775.1',\n",
       "              'GCF_017151735.1',\n",
       "              'GCF_017151725.1',\n",
       "              'GCF_001838925.1',\n",
       "              'GCF_000013425.1',\n",
       "              'GCF_008385845.1',\n",
       "              'GCF_008385815.1',\n",
       "              'GCF_008385875.1',\n",
       "              'GCF_008385975.1',\n",
       "              'GCF_008385995.1',\n",
       "              'GCF_016917335.1',\n",
       "              'GCF_001836155.1',\n",
       "              'GCF_017892245.1',\n",
       "              'GCF_001836465.1'],\n",
       "             'GCF_013374815.1': ['GCF_000008865.2',\n",
       "              'GCF_001027225.1',\n",
       "              'GCF_022354085.1',\n",
       "              'GCF_000006925.2',\n",
       "              'GCF_013374815.1'],\n",
       "             'GCF_000006765.1': ['GCF_005503545.1', 'GCF_000006765.1'],\n",
       "             'GCF_029961225.1': ['GCF_029961225.1'],\n",
       "             'GCA_020975435.1': ['GCA_020975435.1'],\n",
       "             'GCF_000196035.1': ['GCF_000196035.1', 'GCF_028867355.1'],\n",
       "             'GCF_000146045.2': ['GCF_000146045.2'],\n",
       "             'GCF_006094475.1': ['GCF_006094475.1']})"
      ]
     },
     "execution_count": 43,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "representative_genomes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 44,
   "id": "102ab77e-522e-4fe7-a971-6f592c6ecfd5",
   "metadata": {},
   "outputs": [],
   "source": [
    "representative_labels = []\n",
    "cluster_members = []\n",
    "for idx, row in downloaded_assemblies.iterrows():\n",
    "    accession = row['Assembly Accession ID']\n",
    "    if accession in representative_genomes.keys():\n",
    "        representative_labels.append(True)\n",
    "        cluster_members.append(','.join(representative_genomes[accession]))\n",
    "    else:\n",
    "        representative_labels.append(False)\n",
    "        cluster_members.append(','.join(representative_genomes[member2representative[accession]]))\n",
    "downloaded_assemblies['Cluster Representative'] = representative_labels\n",
    "downloaded_assemblies['Cluster Members'] = cluster_members\n",
    "\n",
    "representative_df = downloaded_assemblies[downloaded_assemblies['Cluster Representative']]\n",
    "\n",
    "seq2assembly_dict = get_seq2assembly_dict(working_directory, representative_df)\n",
    "reference_fasta = merge_reference_fasta(list(representative_df['Assembly Accession ID']), working_directory)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 46,
   "id": "6819db75-17dc-4ff1-94b1-58b2e3bf26fb",
   "metadata": {},
   "outputs": [],
   "source": [
    "def sort_samfile(assembly_id, aligner_output, output_dir, min_mapq, threads=20):\n",
    "    '''converting and sorting alignment files'''\n",
    "    bam_files = os.path.join(output_dir, \"bam_files\")\n",
    "    # sam_files = os.path.join(output_dir, \"sam_files\")\n",
    "    \n",
    "    if not os.path.exists(bam_files):\n",
    "        os.mkdir(bam_files)\n",
    "\n",
    "    # covert sam file to binary bam file\n",
    "    samtools_view_res = subprocess.Popen([\n",
    "        \"/home/Users/yl181/tools/samtools-1.19.2/samtools\",\n",
    "        \"view\",\n",
    "        \"-@\", str(threads),\n",
    "        \"--min-MQ\", str(min_mapq),\n",
    "        \"-bS\"], #os.path.join(sam_files, f\"{assembly_id}.sam\")\n",
    "        stdin=aligner_output.stdout,\n",
    "        stdout=subprocess.PIPE)\n",
    "\n",
    "    # sort the bam file \n",
    "    subprocess.run([\n",
    "        \"/home/Users/yl181/tools/samtools-1.19.2/samtools\",\n",
    "        \"sort\",\n",
    "        \"-@\", str(threads),\n",
    "        \"-o\", os.path.join(bam_files, f\"{assembly_id}.sorted.bam\"),\n",
    "        \"-O\", \"BAM\"],\n",
    "        stdin=samtools_view_res.stdout,\n",
    "        stderr=subprocess.DEVNULL,\n",
    "        check=True)\n",
    "\n",
    "    # indexing the sorted bam file \n",
    "    subprocess.run([\n",
    "        \"/home/Users/yl181/tools/samtools-1.19.2/samtools\",\n",
    "        \"index\",\n",
    "        os.path.join(bam_files, f\"{assembly_id}.sorted.bam\")],\n",
    "                    check=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 47,
   "id": "c4d64548-8c89-4970-ab7d-c0b4cb7f0e97",
   "metadata": {},
   "outputs": [],
   "source": [
    "if mode == 'ont':\n",
    "    aligner_output = run_minimap2(input_fastq, reference_fasta, 'merged', working_directory, threads=threads)\n",
    "if mode == 'illumina':\n",
    "    aligner_output = run_bowtie2(input_fastq, input_fastq2, reference_fasta, 'merged', working_directory, threads=threads)\n",
    "sort_samfile('merged', aligner_output, working_directory, min_mapq=0, threads=threads)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 48,
   "id": "3147627c-e15e-4a10-a07b-8b4865ef86c8",
   "metadata": {},
   "outputs": [
    {
     "ename": "CalledProcessError",
     "evalue": "Command '['samtools', 'coverage', '/home/Users/yl181/memu/magnet_v222_zymo_even_full/Zymo-EVEN-p0.75/Zymo-Even-3/bam_files/merged.sorted.bam', '--ff', 'UNMAP,QCFAIL,DUP', '-q', '0']' died with <Signals.SIGSEGV: 11>.",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mRemoteTraceback\u001b[0m                           Traceback (most recent call last)",
      "\u001b[0;31mRemoteTraceback\u001b[0m: \n\"\"\"\nTraceback (most recent call last):\n  File \"/home/Users/yl181/miniconda3/envs/magnet/lib/python3.9/multiprocessing/pool.py\", line 125, in worker\n    result = (True, func(*args, **kwds))\n  File \"/home/Users/yl181/miniconda3/envs/magnet/lib/python3.9/multiprocessing/pool.py\", line 51, in starmapstar\n    return list(itertools.starmap(args[0], args[1]))\n  File \"/tmp/ipykernel_1529367/2249496809.py\", line 19, in samtools_calculate_coverage\n    subprocess.run(command,\n  File \"/home/Users/yl181/miniconda3/envs/magnet/lib/python3.9/subprocess.py\", line 528, in run\n    raise CalledProcessError(retcode, process.args,\nsubprocess.CalledProcessError: Command '['samtools', 'coverage', '/home/Users/yl181/memu/magnet_v222_zymo_even_full/Zymo-EVEN-p0.75/Zymo-Even-3/bam_files/merged.sorted.bam', '--ff', 'UNMAP,QCFAIL,DUP', '-q', '0']' died with <Signals.SIGSEGV: 11>.\n\"\"\"",
      "\nThe above exception was the direct cause of the following exception:\n",
      "\u001b[0;31mCalledProcessError\u001b[0m                        Traceback (most recent call last)",
      "Cell \u001b[0;32mIn [48], line 6\u001b[0m\n\u001b[1;32m      3\u001b[0m     os\u001b[38;5;241m.\u001b[39mmkdir(coverage_files)\n\u001b[1;32m      5\u001b[0m pool \u001b[38;5;241m=\u001b[39m Pool(processes\u001b[38;5;241m=\u001b[39mthreads)\n\u001b[0;32m----> 6\u001b[0m \u001b[43mpool\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mstarmap\u001b[49m\u001b[43m(\u001b[49m\u001b[43msamtools_calculate_coverage\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[38;5;28;43mzip\u001b[39;49m\u001b[43m(\u001b[49m\u001b[43mrepeat\u001b[49m\u001b[43m(\u001b[49m\u001b[43mworking_directory\u001b[49m\u001b[43m)\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43m[\u001b[49m\u001b[38;5;28;43;01mTrue\u001b[39;49;00m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[38;5;28;43;01mFalse\u001b[39;49;00m\u001b[43m]\u001b[49m\u001b[43m)\u001b[49m\u001b[43m)\u001b[49m\n\u001b[1;32m      7\u001b[0m pool\u001b[38;5;241m.\u001b[39mclose()\n\u001b[1;32m      8\u001b[0m pool\u001b[38;5;241m.\u001b[39mjoin()\n",
      "File \u001b[0;32m~/miniconda3/envs/magnet/lib/python3.9/multiprocessing/pool.py:372\u001b[0m, in \u001b[0;36mPool.starmap\u001b[0;34m(self, func, iterable, chunksize)\u001b[0m\n\u001b[1;32m    366\u001b[0m \u001b[38;5;28;01mdef\u001b[39;00m \u001b[38;5;21mstarmap\u001b[39m(\u001b[38;5;28mself\u001b[39m, func, iterable, chunksize\u001b[38;5;241m=\u001b[39m\u001b[38;5;28;01mNone\u001b[39;00m):\n\u001b[1;32m    367\u001b[0m     \u001b[38;5;124;03m'''\u001b[39;00m\n\u001b[1;32m    368\u001b[0m \u001b[38;5;124;03m    Like `map()` method but the elements of the `iterable` are expected to\u001b[39;00m\n\u001b[1;32m    369\u001b[0m \u001b[38;5;124;03m    be iterables as well and will be unpacked as arguments. Hence\u001b[39;00m\n\u001b[1;32m    370\u001b[0m \u001b[38;5;124;03m    `func` and (a, b) becomes func(a, b).\u001b[39;00m\n\u001b[1;32m    371\u001b[0m \u001b[38;5;124;03m    '''\u001b[39;00m\n\u001b[0;32m--> 372\u001b[0m     \u001b[38;5;28;01mreturn\u001b[39;00m \u001b[38;5;28;43mself\u001b[39;49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43m_map_async\u001b[49m\u001b[43m(\u001b[49m\u001b[43mfunc\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43miterable\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mstarmapstar\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mchunksize\u001b[49m\u001b[43m)\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mget\u001b[49m\u001b[43m(\u001b[49m\u001b[43m)\u001b[49m\n",
      "File \u001b[0;32m~/miniconda3/envs/magnet/lib/python3.9/multiprocessing/pool.py:771\u001b[0m, in \u001b[0;36mApplyResult.get\u001b[0;34m(self, timeout)\u001b[0m\n\u001b[1;32m    769\u001b[0m     \u001b[38;5;28;01mreturn\u001b[39;00m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39m_value\n\u001b[1;32m    770\u001b[0m \u001b[38;5;28;01melse\u001b[39;00m:\n\u001b[0;32m--> 771\u001b[0m     \u001b[38;5;28;01mraise\u001b[39;00m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39m_value\n",
      "\u001b[0;31mCalledProcessError\u001b[0m: Command '['samtools', 'coverage', '/home/Users/yl181/memu/magnet_v222_zymo_even_full/Zymo-EVEN-p0.75/Zymo-Even-3/bam_files/merged.sorted.bam', '--ff', 'UNMAP,QCFAIL,DUP', '-q', '0']' died with <Signals.SIGSEGV: 11>."
     ]
    }
   ],
   "source": [
    "coverage_files = os.path.join(working_directory, \"coverage_files\")\n",
    "if not os.path.exists(coverage_files):\n",
    "    os.mkdir(coverage_files)\n",
    "    \n",
    "pool = Pool(processes=threads)\n",
    "pool.starmap(samtools_calculate_coverage, zip(repeat(working_directory), [True, False]))\n",
    "pool.close()\n",
    "pool.join()\n",
    "    \n",
    "\n",
    "representative_df = alignment_summary(representative_df,\n",
    "                                      working_directory,\n",
    "                                      seq2assembly_dict,\n",
    "                                      include_supp=True)\n",
    "\n",
    "representative_df = alignment_summary(representative_df,\n",
    "                                      working_directory,\n",
    "                                      seq2assembly_dict,\n",
    "                                      include_supp=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 130,
   "id": "c0185772-db4f-4b48-a879-97a5f319dc95",
   "metadata": {},
   "outputs": [],
   "source": [
    "def samtools_merged_consensus(output_directory, threads):\n",
    "    merged_bam = os.path.join(output_directory, 'bam_files', 'merged.sorted.bam')\n",
    "    subprocess.run(['samtools', 'consensus', \n",
    "                    '--show-ins', 'no', \n",
    "                    '--show-del', 'yes',\n",
    "                    '--min-MQ', str(20),\n",
    "                    '-a',\n",
    "                    '--mode', \"simple\",\n",
    "                    '--threads', str(threads),\n",
    "                    merged_bam, \n",
    "                    '-o', os.path.join(output_directory, 'merged_consensus.fasta')],\n",
    "                  check=True)\n",
    "    \n",
    "    consensus_record_dict = SeqIO.to_dict(SeqIO.parse(os.path.join(output_directory, 'merged_consensus.fasta'), \"fasta\"))\n",
    "    return consensus_record_dict"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 131,
   "id": "2e1a82cf-cff7-4ca8-b7e2-bfba1975e2af",
   "metadata": {},
   "outputs": [],
   "source": [
    "consensus_record_dict = samtools_merged_consensus(working_directory, threads)\n",
    "representative_df = ani_summary(representative_df, consensus_record_dict, working_directory, threads)\n",
    "representative_df = call_present_absent(representative_df, min_coverage_score)\n",
    "representative_df.sort_values(['Primary Score'], ascending=False).to_csv(os.path.join(working_directory, \n",
    "                                                                                      'cluster_representative.csv'), \n",
    "                                                                         index=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 132,
   "id": "5572b780-8e8c-4382-95fb-4cf4f98532e3",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Index(['Taxonomy ID', 'Assembly Accession ID', 'Source Database',\n",
       "       'Is Representative', 'Assembly Level', 'Organism of Assembly', 'Strain',\n",
       "       'Total Length', 'Downloaded', 'Species', 'Cluster Representative',\n",
       "       'Cluster Members', 'Secondary Breadth', 'Secondary Expected',\n",
       "       'Secondary Score', 'Secondary Depth', 'Primary Breadth',\n",
       "       'Primary Expected', 'Primary Score', 'Primary Depth', 'Consensus ANI',\n",
       "       'Combined PS and ANI (Sqrt(ANI)xPSx100)', 'Presence/Absence'],\n",
       "      dtype='object')"
      ]
     },
     "execution_count": 132,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "representative_df.columns"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 133,
   "id": "08c38907-c2c7-481a-889f-c296558d0aa3",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>Taxonomy ID</th>\n",
       "      <th>Assembly Accession ID</th>\n",
       "      <th>Source Database</th>\n",
       "      <th>Is Representative</th>\n",
       "      <th>Assembly Level</th>\n",
       "      <th>Organism of Assembly</th>\n",
       "      <th>Strain</th>\n",
       "      <th>Total Length</th>\n",
       "      <th>Downloaded</th>\n",
       "      <th>Species</th>\n",
       "      <th>...</th>\n",
       "      <th>Secondary Expected</th>\n",
       "      <th>Secondary Score</th>\n",
       "      <th>Secondary Depth</th>\n",
       "      <th>Primary Breadth</th>\n",
       "      <th>Primary Expected</th>\n",
       "      <th>Primary Score</th>\n",
       "      <th>Primary Depth</th>\n",
       "      <th>Consensus ANI</th>\n",
       "      <th>Combined PS and ANI (Sqrt(ANI)xPSx100)</th>\n",
       "      <th>Presence/Absence</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>3056648</td>\n",
       "      <td>GCF_030413175.1</td>\n",
       "      <td>SOURCE_DATABASE_REFSEQ</td>\n",
       "      <td>True</td>\n",
       "      <td>Complete Genome</td>\n",
       "      <td>Bradyrhizobium roseus</td>\n",
       "      <td>S12-14-2</td>\n",
       "      <td>7319411.0</td>\n",
       "      <td>True</td>\n",
       "      <td>Bradyrhizobium sp. S12-14-2</td>\n",
       "      <td>...</td>\n",
       "      <td>0.143660</td>\n",
       "      <td>0.870010</td>\n",
       "      <td>1.240736</td>\n",
       "      <td>0.066733</td>\n",
       "      <td>0.073593</td>\n",
       "      <td>0.906780</td>\n",
       "      <td>1.145335</td>\n",
       "      <td>0.827142</td>\n",
       "      <td>82.47</td>\n",
       "      <td>Present</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>2954771</td>\n",
       "      <td>GCF_031202265.1</td>\n",
       "      <td>SOURCE_DATABASE_REFSEQ</td>\n",
       "      <td>False</td>\n",
       "      <td>Complete Genome</td>\n",
       "      <td>Bradyrhizobium sp. Ash2021</td>\n",
       "      <td>Ash2021</td>\n",
       "      <td>9704287.0</td>\n",
       "      <td>True</td>\n",
       "      <td>Bradyrhizobium sp. Ash2021</td>\n",
       "      <td>...</td>\n",
       "      <td>0.230135</td>\n",
       "      <td>0.841104</td>\n",
       "      <td>1.351045</td>\n",
       "      <td>0.107387</td>\n",
       "      <td>0.118962</td>\n",
       "      <td>0.902699</td>\n",
       "      <td>1.179293</td>\n",
       "      <td>0.803551</td>\n",
       "      <td>80.92</td>\n",
       "      <td>Present</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>3063329</td>\n",
       "      <td>GCF_030546645.1</td>\n",
       "      <td>SOURCE_DATABASE_REFSEQ</td>\n",
       "      <td>False</td>\n",
       "      <td>Contig</td>\n",
       "      <td>Flavitalea sp. BT771</td>\n",
       "      <td>BT771</td>\n",
       "      <td>8732077.0</td>\n",
       "      <td>True</td>\n",
       "      <td>Flavitalea sp. BT771</td>\n",
       "      <td>...</td>\n",
       "      <td>0.370252</td>\n",
       "      <td>0.795410</td>\n",
       "      <td>1.570058</td>\n",
       "      <td>0.248169</td>\n",
       "      <td>0.303316</td>\n",
       "      <td>0.818186</td>\n",
       "      <td>1.456170</td>\n",
       "      <td>0.799718</td>\n",
       "      <td>73.17</td>\n",
       "      <td>Present</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>5</th>\n",
       "      <td>2219043</td>\n",
       "      <td>GCF_003574315.2</td>\n",
       "      <td>SOURCE_DATABASE_REFSEQ</td>\n",
       "      <td>True</td>\n",
       "      <td>Complete Genome</td>\n",
       "      <td>Capsulimonas corticalis</td>\n",
       "      <td>AX-7</td>\n",
       "      <td>7645051.0</td>\n",
       "      <td>True</td>\n",
       "      <td>Capsulimonas corticalis</td>\n",
       "      <td>...</td>\n",
       "      <td>0.085282</td>\n",
       "      <td>0.816503</td>\n",
       "      <td>1.280021</td>\n",
       "      <td>0.050279</td>\n",
       "      <td>0.055528</td>\n",
       "      <td>0.905472</td>\n",
       "      <td>1.136095</td>\n",
       "      <td>0.790988</td>\n",
       "      <td>80.53</td>\n",
       "      <td>Present</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>7</th>\n",
       "      <td>722472</td>\n",
       "      <td>GCF_900141755.1</td>\n",
       "      <td>SOURCE_DATABASE_REFSEQ</td>\n",
       "      <td>False</td>\n",
       "      <td>Chromosome</td>\n",
       "      <td>Bradyrhizobium lablabi</td>\n",
       "      <td>GAS499</td>\n",
       "      <td>7910099.0</td>\n",
       "      <td>True</td>\n",
       "      <td>Bradyrhizobium lablabi</td>\n",
       "      <td>...</td>\n",
       "      <td>0.207450</td>\n",
       "      <td>0.871973</td>\n",
       "      <td>1.285189</td>\n",
       "      <td>0.099097</td>\n",
       "      <td>0.109425</td>\n",
       "      <td>0.905620</td>\n",
       "      <td>1.169292</td>\n",
       "      <td>0.796864</td>\n",
       "      <td>80.84</td>\n",
       "      <td>Present</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>...</th>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>139</th>\n",
       "      <td>2480626</td>\n",
       "      <td>GCF_007827045.1</td>\n",
       "      <td>SOURCE_DATABASE_REFSEQ</td>\n",
       "      <td>False</td>\n",
       "      <td>Scaffold</td>\n",
       "      <td>Trebonia kvetii</td>\n",
       "      <td>15TR583</td>\n",
       "      <td>9254646.0</td>\n",
       "      <td>True</td>\n",
       "      <td>Trebonia kvetii</td>\n",
       "      <td>...</td>\n",
       "      <td>0.628096</td>\n",
       "      <td>0.733255</td>\n",
       "      <td>2.147493</td>\n",
       "      <td>0.403171</td>\n",
       "      <td>0.541646</td>\n",
       "      <td>0.744344</td>\n",
       "      <td>1.934744</td>\n",
       "      <td>0.852379</td>\n",
       "      <td>68.72</td>\n",
       "      <td>Present</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>140</th>\n",
       "      <td>652787</td>\n",
       "      <td>GCF_900105165.1</td>\n",
       "      <td>SOURCE_DATABASE_REFSEQ</td>\n",
       "      <td>False</td>\n",
       "      <td>Chromosome</td>\n",
       "      <td>Mucilaginibacter mallensis</td>\n",
       "      <td>MP1X4</td>\n",
       "      <td>6014056.0</td>\n",
       "      <td>True</td>\n",
       "      <td>Mucilaginibacter mallensis</td>\n",
       "      <td>...</td>\n",
       "      <td>0.279387</td>\n",
       "      <td>0.931456</td>\n",
       "      <td>1.258843</td>\n",
       "      <td>0.232672</td>\n",
       "      <td>0.243372</td>\n",
       "      <td>0.956035</td>\n",
       "      <td>1.198360</td>\n",
       "      <td>0.820509</td>\n",
       "      <td>86.60</td>\n",
       "      <td>Present</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>141</th>\n",
       "      <td>1792502</td>\n",
       "      <td>GCF_014642875.1</td>\n",
       "      <td>SOURCE_DATABASE_REFSEQ</td>\n",
       "      <td>False</td>\n",
       "      <td>Scaffold</td>\n",
       "      <td>Puia dinghuensis</td>\n",
       "      <td>CGMCC 1.15448</td>\n",
       "      <td>7182087.0</td>\n",
       "      <td>True</td>\n",
       "      <td>Puia dinghuensis</td>\n",
       "      <td>...</td>\n",
       "      <td>0.424826</td>\n",
       "      <td>0.788081</td>\n",
       "      <td>1.651759</td>\n",
       "      <td>0.287182</td>\n",
       "      <td>0.355743</td>\n",
       "      <td>0.807273</td>\n",
       "      <td>1.530674</td>\n",
       "      <td>0.798597</td>\n",
       "      <td>72.14</td>\n",
       "      <td>Present</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>143</th>\n",
       "      <td>2080758</td>\n",
       "      <td>GCF_004661485.1</td>\n",
       "      <td>SOURCE_DATABASE_REFSEQ</td>\n",
       "      <td>False</td>\n",
       "      <td>Contig</td>\n",
       "      <td>Methylotenera oryzisoli</td>\n",
       "      <td>La3113</td>\n",
       "      <td>2664030.0</td>\n",
       "      <td>True</td>\n",
       "      <td>Methylotenera oryzisoli</td>\n",
       "      <td>...</td>\n",
       "      <td>0.111636</td>\n",
       "      <td>0.830866</td>\n",
       "      <td>1.275716</td>\n",
       "      <td>0.078677</td>\n",
       "      <td>0.090951</td>\n",
       "      <td>0.865047</td>\n",
       "      <td>1.211409</td>\n",
       "      <td>0.764349</td>\n",
       "      <td>75.63</td>\n",
       "      <td>Present</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>145</th>\n",
       "      <td>244734</td>\n",
       "      <td>GCF_008932115.1</td>\n",
       "      <td>SOURCE_DATABASE_REFSEQ</td>\n",
       "      <td>True</td>\n",
       "      <td>Complete Genome</td>\n",
       "      <td>Bradyrhizobium betae</td>\n",
       "      <td>PL7HG1</td>\n",
       "      <td>7419402.0</td>\n",
       "      <td>True</td>\n",
       "      <td>Bradyrhizobium betae</td>\n",
       "      <td>...</td>\n",
       "      <td>0.154391</td>\n",
       "      <td>0.900921</td>\n",
       "      <td>1.205537</td>\n",
       "      <td>0.063496</td>\n",
       "      <td>0.068461</td>\n",
       "      <td>0.927480</td>\n",
       "      <td>1.116716</td>\n",
       "      <td>0.825432</td>\n",
       "      <td>84.26</td>\n",
       "      <td>Present</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>99 rows × 23 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "     Taxonomy ID Assembly Accession ID         Source Database  \\\n",
       "1        3056648       GCF_030413175.1  SOURCE_DATABASE_REFSEQ   \n",
       "3        2954771       GCF_031202265.1  SOURCE_DATABASE_REFSEQ   \n",
       "4        3063329       GCF_030546645.1  SOURCE_DATABASE_REFSEQ   \n",
       "5        2219043       GCF_003574315.2  SOURCE_DATABASE_REFSEQ   \n",
       "7         722472       GCF_900141755.1  SOURCE_DATABASE_REFSEQ   \n",
       "..           ...                   ...                     ...   \n",
       "139      2480626       GCF_007827045.1  SOURCE_DATABASE_REFSEQ   \n",
       "140       652787       GCF_900105165.1  SOURCE_DATABASE_REFSEQ   \n",
       "141      1792502       GCF_014642875.1  SOURCE_DATABASE_REFSEQ   \n",
       "143      2080758       GCF_004661485.1  SOURCE_DATABASE_REFSEQ   \n",
       "145       244734       GCF_008932115.1  SOURCE_DATABASE_REFSEQ   \n",
       "\n",
       "    Is Representative   Assembly Level        Organism of Assembly  \\\n",
       "1                True  Complete Genome       Bradyrhizobium roseus   \n",
       "3               False  Complete Genome  Bradyrhizobium sp. Ash2021   \n",
       "4               False           Contig        Flavitalea sp. BT771   \n",
       "5                True  Complete Genome     Capsulimonas corticalis   \n",
       "7               False       Chromosome      Bradyrhizobium lablabi   \n",
       "..                ...              ...                         ...   \n",
       "139             False         Scaffold             Trebonia kvetii   \n",
       "140             False       Chromosome  Mucilaginibacter mallensis   \n",
       "141             False         Scaffold            Puia dinghuensis   \n",
       "143             False           Contig     Methylotenera oryzisoli   \n",
       "145              True  Complete Genome        Bradyrhizobium betae   \n",
       "\n",
       "            Strain  Total Length  Downloaded                      Species  \\\n",
       "1         S12-14-2     7319411.0        True  Bradyrhizobium sp. S12-14-2   \n",
       "3          Ash2021     9704287.0        True   Bradyrhizobium sp. Ash2021   \n",
       "4            BT771     8732077.0        True         Flavitalea sp. BT771   \n",
       "5             AX-7     7645051.0        True      Capsulimonas corticalis   \n",
       "7           GAS499     7910099.0        True       Bradyrhizobium lablabi   \n",
       "..             ...           ...         ...                          ...   \n",
       "139        15TR583     9254646.0        True              Trebonia kvetii   \n",
       "140          MP1X4     6014056.0        True   Mucilaginibacter mallensis   \n",
       "141  CGMCC 1.15448     7182087.0        True             Puia dinghuensis   \n",
       "143         La3113     2664030.0        True      Methylotenera oryzisoli   \n",
       "145         PL7HG1     7419402.0        True         Bradyrhizobium betae   \n",
       "\n",
       "     ...  Secondary Expected Secondary Score  Secondary Depth  \\\n",
       "1    ...            0.143660        0.870010         1.240736   \n",
       "3    ...            0.230135        0.841104         1.351045   \n",
       "4    ...            0.370252        0.795410         1.570058   \n",
       "5    ...            0.085282        0.816503         1.280021   \n",
       "7    ...            0.207450        0.871973         1.285189   \n",
       "..   ...                 ...             ...              ...   \n",
       "139  ...            0.628096        0.733255         2.147493   \n",
       "140  ...            0.279387        0.931456         1.258843   \n",
       "141  ...            0.424826        0.788081         1.651759   \n",
       "143  ...            0.111636        0.830866         1.275716   \n",
       "145  ...            0.154391        0.900921         1.205537   \n",
       "\n",
       "     Primary Breadth  Primary Expected  Primary Score  Primary Depth  \\\n",
       "1           0.066733          0.073593       0.906780       1.145335   \n",
       "3           0.107387          0.118962       0.902699       1.179293   \n",
       "4           0.248169          0.303316       0.818186       1.456170   \n",
       "5           0.050279          0.055528       0.905472       1.136095   \n",
       "7           0.099097          0.109425       0.905620       1.169292   \n",
       "..               ...               ...            ...            ...   \n",
       "139         0.403171          0.541646       0.744344       1.934744   \n",
       "140         0.232672          0.243372       0.956035       1.198360   \n",
       "141         0.287182          0.355743       0.807273       1.530674   \n",
       "143         0.078677          0.090951       0.865047       1.211409   \n",
       "145         0.063496          0.068461       0.927480       1.116716   \n",
       "\n",
       "     Consensus ANI  Combined PS and ANI (Sqrt(ANI)xPSx100)  Presence/Absence  \n",
       "1         0.827142                                   82.47           Present  \n",
       "3         0.803551                                   80.92           Present  \n",
       "4         0.799718                                   73.17           Present  \n",
       "5         0.790988                                   80.53           Present  \n",
       "7         0.796864                                   80.84           Present  \n",
       "..             ...                                     ...               ...  \n",
       "139       0.852379                                   68.72           Present  \n",
       "140       0.820509                                   86.60           Present  \n",
       "141       0.798597                                   72.14           Present  \n",
       "143       0.764349                                   75.63           Present  \n",
       "145       0.825432                                   84.26           Present  \n",
       "\n",
       "[99 rows x 23 columns]"
      ]
     },
     "execution_count": 133,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "representative_df[representative_df['Presence/Absence'] == 'Present']"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "128b40fc-2893-4d9c-b4d5-739eab53a260",
   "metadata": {},
   "source": [
    "# With-in Cluster Classification"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 428,
   "id": "cfa3d3bb-12ed-4a6e-bbf3-62321de702cb",
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_expected_coverage(genome_length, reads_mapped, genome_totol_count):\n",
    "    mean_mapping_length = genome_totol_count/reads_mapped\n",
    "    \n",
    "    N = genome_length/mean_mapping_length\n",
    "    x = reads_mapped\n",
    "    \n",
    "    expected_M = N*(1-((1-1/N)**x))\n",
    "    variance = N*((1-1/N)**x) + (N**2)*(1-1/N)*((1-2/N)**x)-(N**2)*((1-1/N)**(2*x))\n",
    "    \n",
    "    expected_coverage = expected_M/N\n",
    "    try:\n",
    "        std = math.sqrt(variance)\n",
    "    except ValueError:\n",
    "        std = 0\n",
    "    return expected_coverage, std"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 424,
   "id": "90df11b9-a987-447c-996b-382ce8595625",
   "metadata": {},
   "outputs": [],
   "source": [
    "def calculate_depth(assembly_id, taxa_records):\n",
    "    genome_length = taxa_records[assembly_id]['genome_length']\n",
    "    covbases = taxa_records[assembly_id]['covbases']\n",
    "    genome_totol_count = taxa_records[assembly_id]['genome_totol_count']\n",
    "    reads_mapped = taxa_records[assembly_id]['reads_mapped']\n",
    "    \n",
    "    if genome_totol_count == 0 or reads_mapped == 0:\n",
    "        breadth_coverage = 0\n",
    "        depth_coverage = 0\n",
    "        expected_breadth_coverage = 0\n",
    "    else:\n",
    "        breadth_coverage = covbases/genome_length\n",
    "        depth_coverage = genome_totol_count/covbases\n",
    "        expected_breadth_coverage, std = get_expected_coverage(genome_length, reads_mapped, genome_totol_count)\n",
    "    \n",
    "    return breadth_coverage, depth_coverage, expected_breadth_coverage"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 425,
   "id": "e5d23879-6a3d-4562-b437-a3b4b7cf350a",
   "metadata": {},
   "outputs": [],
   "source": [
    "def alignment_summary(downloaded_assemblies, output_directory, seq2assembly_dict, include_supp=True):\n",
    "    if include_supp:\n",
    "        coverage_file_name = 'secondary_coverage.tsv'\n",
    "        column_prefix = 'Secondary'\n",
    "    else:\n",
    "        coverage_file_name = 'primary_coverage.tsv'\n",
    "        column_prefix = 'Primary'\n",
    "    \n",
    "    columns = [f'{column_prefix} Breadth',\n",
    "               f'{column_prefix} Expected',\n",
    "               f'{column_prefix} Score',\n",
    "               f'{column_prefix} Depth']\n",
    "    \n",
    "    coverage_df = pd.read_csv(os.path.join(output_directory,\n",
    "                                           'coverage_files',\n",
    "                                           coverage_file_name),\n",
    "                              sep='\\t',\n",
    "                              header=None,\n",
    "                              names=['rname','startpos','endpos','numreads','covbases','coverage','meandepth','meanbaseq','meanmapq'])\n",
    "    \n",
    "    taxa_records = defaultdict(lambda: defaultdict(int))\n",
    "    for idx, row in coverage_df.iterrows():\n",
    "        taxa_reference = seq2assembly_dict[row['rname']]\n",
    "        taxa_records[taxa_reference]['genome_length'] += row['endpos']\n",
    "        taxa_records[taxa_reference]['reads_mapped'] += row['numreads']\n",
    "        taxa_records[taxa_reference]['genome_totol_count'] += int(row['meandepth'] * row['endpos'])\n",
    "        taxa_records[taxa_reference]['covbases'] += row['covbases']\n",
    "    \n",
    "    breadth_coverage_list = []\n",
    "    depth_coverage_list = []\n",
    "    expected_breadth_coverage_list = []\n",
    "    coverage_score = []\n",
    "    for assembly_id in downloaded_assemblies['Assembly Accession ID']:\n",
    "        breadth_coverage, depth_coverage, expected_breadth_coverage = calculate_depth(assembly_id, taxa_records)\n",
    "        breadth_coverage_list.append(breadth_coverage)\n",
    "        depth_coverage_list.append(depth_coverage)\n",
    "        expected_breadth_coverage_list.append(expected_breadth_coverage)\n",
    "        if expected_breadth_coverage != 0:\n",
    "            coverage_score.append(min(breadth_coverage/expected_breadth_coverage, 1))\n",
    "        else:\n",
    "            coverage_score.append(0)\n",
    "        \n",
    "    downloaded_assemblies[columns[0]] = breadth_coverage_list\n",
    "    downloaded_assemblies[columns[1]] = expected_breadth_coverage_list\n",
    "    downloaded_assemblies[columns[2]] = coverage_score\n",
    "    downloaded_assemblies[columns[3]] = depth_coverage_list\n",
    "    \n",
    "    downloaded_assemblies.to_csv(os.path.join(output_directory, 'alignment.csv'), index=False)\n",
    "    \n",
    "    return downloaded_assemblies"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 383,
   "id": "43a8d65c-9516-409e-b582-04e16d5b3b82",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "['GCF_000013425.1', 'GCF_001838925.1', 'GCF_016893225.1']\n",
      "['GCF_001598635.1', 'GCF_001053935.1', 'GCF_031143435.1', 'GCF_018919785.1']\n"
     ]
    }
   ],
   "source": [
    "counter = 0\n",
    "for idx, row in representative_df.iterrows():\n",
    "    cluster_members = row['Cluster Members'].split(\",\")\n",
    "    \n",
    "    if len(cluster_members) > 1:\n",
    "        print(cluster_members)\n",
    "        counter += 1\n",
    "    \n",
    "    if counter == 2:\n",
    "        break"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 384,
   "id": "49ab1e5a-0634-4837-bded-511fad3f44bb",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "['GCF_001598635.1', 'GCF_001053935.1', 'GCF_031143435.1', 'GCF_018919785.1']"
      ]
     },
     "execution_count": 384,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "cluster_members"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 421,
   "id": "9fe8b69d-9c70-49e6-9148-fcf554df4d2e",
   "metadata": {},
   "outputs": [],
   "source": [
    "members_df = downloaded_assemblies[downloaded_assemblies['Assembly Accession ID'].isin(cluster_members)]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 422,
   "id": "4f3cbaaf-6eb5-400c-8489-7530a588ddff",
   "metadata": {},
   "outputs": [],
   "source": [
    "seq2assembly_dict = get_seq2assembly_dict(working_directory, members_df)\n",
    "reference_fasta = merge_reference_fasta(list(members_df['Assembly Accession ID']), working_directory)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 423,
   "id": "9550eb97-4c28-4dff-8251-8ba18930d424",
   "metadata": {},
   "outputs": [],
   "source": [
    "if mode == 'ont':\n",
    "    aligner_output = run_minimap2(input_fastq, reference_fasta, 'merged', working_directory, threads=threads)\n",
    "if mode == 'illumina':\n",
    "    aligner_output = run_bowtie2(input_fastq, input_fastq2, reference_fasta, 'merged', working_directory, threads=threads)\n",
    "sort_samfile('merged', aligner_output, working_directory, min_mapq=0, threads=threads)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 431,
   "id": "f7d8d68f-0ab4-4194-89f0-3edded1475e0",
   "metadata": {},
   "outputs": [],
   "source": [
    "consensus_record_dict = samtools_merged_consensus(working_directory, threads)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 432,
   "id": "49c7305c-d423-4fe3-8ae2-be4fed1f63a6",
   "metadata": {},
   "outputs": [],
   "source": [
    "members_df = ani_summary(members_df, consensus_record_dict, working_directory, threads)\n",
    "members_df = call_present_absent(members_df, min_coverage_score)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 533,
   "id": "2e015a4f-9fac-42b1-986d-ffc47f2fe3ca",
   "metadata": {},
   "outputs": [],
   "source": [
    "def samtools_calculate_coverage(output_dir, include_supp=False):\n",
    "    coverage_files = os.path.join(output_dir, \"coverage_files\")\n",
    "    bam_files = os.path.join(output_dir, \"bam_files\")\n",
    "\n",
    "    command = [\"samtools\",\n",
    "               \"coverage\",\n",
    "               \"--no-header\",\n",
    "               os.path.join(bam_files, f\"merged.sorted.bam\")]\n",
    "    \n",
    "    if include_supp:\n",
    "        # samtools coverage --ff UNMAP,QCFAIL,DUP -q 0 merged.sorted.bam > secondary_coverage.tsv\n",
    "        coverage_file = os.path.join(coverage_files, f\"secondary_coverage.tsv\")\n",
    "        command += ['--ff', 'UNMAP,QCFAIL,DUP', '-q', str(0)]\n",
    "    else:\n",
    "        # samtools coverage merged.sorted.bam -q 20 > primary_coverage.tsv\n",
    "        coverage_file = os.path.join(coverage_files, f\"primary_coverage.tsv\")\n",
    "        command += ['-q', str(20)]\n",
    "        # command += ['--ff', 'UNMAP,QCFAIL,DUP', '-q', str(20)]\n",
    "\n",
    "    subprocess.run(command,\n",
    "                    check=True,\n",
    "                    stdout=open(coverage_file, \"w\"))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 534,
   "id": "a034d173-52ed-490e-b8fc-347e77761d8b",
   "metadata": {},
   "outputs": [],
   "source": [
    "coverage_files = os.path.join(working_directory, \"coverage_files\")\n",
    "if not os.path.exists(coverage_files):\n",
    "    os.mkdir(coverage_files)\n",
    "    \n",
    "pool = Pool(processes=threads)\n",
    "pool.starmap(samtools_calculate_coverage, zip(repeat(working_directory), [True, False]))\n",
    "pool.close()\n",
    "pool.join()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 535,
   "id": "9784e9b1-52e6-4691-8bf0-4ccd1150e833",
   "metadata": {},
   "outputs": [],
   "source": [
    "members_df = alignment_summary(members_df,\n",
    "                                      working_directory,\n",
    "                                      seq2assembly_dict,\n",
    "                                      include_supp=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 536,
   "id": "7b529ca1-207f-4e5a-9ae2-3c71957d4b04",
   "metadata": {},
   "outputs": [],
   "source": [
    "members_df = alignment_summary(members_df,\n",
    "                                      working_directory,\n",
    "                                      seq2assembly_dict,\n",
    "                                      include_supp=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 537,
   "id": "ae27af14-1374-4f7e-a61e-d2970151d9c5",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>Assembly Accession ID</th>\n",
       "      <th>Species</th>\n",
       "      <th>Assembly Level</th>\n",
       "      <th>Cluster Representative</th>\n",
       "      <th>Total Length</th>\n",
       "      <th>Secondary Breadth</th>\n",
       "      <th>Primary Breadth</th>\n",
       "      <th>Secondary Expected</th>\n",
       "      <th>Primary Expected</th>\n",
       "      <th>Secondary Score</th>\n",
       "      <th>Primary Score</th>\n",
       "      <th>Consensus ANI</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>GCF_001598635.1</td>\n",
       "      <td>Enterococcus faecalis</td>\n",
       "      <td>Complete Genome</td>\n",
       "      <td>True</td>\n",
       "      <td>2803429.0</td>\n",
       "      <td>0.876527</td>\n",
       "      <td>0.238773</td>\n",
       "      <td>0.996414</td>\n",
       "      <td>0.411209</td>\n",
       "      <td>0.879681</td>\n",
       "      <td>0.580661</td>\n",
       "      <td>0.938553</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>13</th>\n",
       "      <td>GCF_001053935.1</td>\n",
       "      <td>Enterococcus sp. 1140_ESPC</td>\n",
       "      <td>Contig</td>\n",
       "      <td>False</td>\n",
       "      <td>2713742.0</td>\n",
       "      <td>0.890697</td>\n",
       "      <td>0.170908</td>\n",
       "      <td>0.990717</td>\n",
       "      <td>0.346761</td>\n",
       "      <td>0.899043</td>\n",
       "      <td>0.492870</td>\n",
       "      <td>0.932959</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>14</th>\n",
       "      <td>GCF_031143435.1</td>\n",
       "      <td>Enterococcus sp. FR192</td>\n",
       "      <td>Contig</td>\n",
       "      <td>False</td>\n",
       "      <td>3090693.0</td>\n",
       "      <td>0.802494</td>\n",
       "      <td>0.268412</td>\n",
       "      <td>0.984906</td>\n",
       "      <td>0.606032</td>\n",
       "      <td>0.814792</td>\n",
       "      <td>0.442900</td>\n",
       "      <td>0.952038</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>16</th>\n",
       "      <td>GCF_018919785.1</td>\n",
       "      <td>Enterococcus sp. S171_ASV_20</td>\n",
       "      <td>Contig</td>\n",
       "      <td>False</td>\n",
       "      <td>3284381.0</td>\n",
       "      <td>0.756196</td>\n",
       "      <td>0.156462</td>\n",
       "      <td>0.978218</td>\n",
       "      <td>0.266662</td>\n",
       "      <td>0.773034</td>\n",
       "      <td>0.586741</td>\n",
       "      <td>0.929656</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "   Assembly Accession ID                       Species   Assembly Level  \\\n",
       "4        GCF_001598635.1         Enterococcus faecalis  Complete Genome   \n",
       "13       GCF_001053935.1    Enterococcus sp. 1140_ESPC           Contig   \n",
       "14       GCF_031143435.1        Enterococcus sp. FR192           Contig   \n",
       "16       GCF_018919785.1  Enterococcus sp. S171_ASV_20           Contig   \n",
       "\n",
       "    Cluster Representative  Total Length  Secondary Breadth  Primary Breadth  \\\n",
       "4                     True     2803429.0           0.876527         0.238773   \n",
       "13                   False     2713742.0           0.890697         0.170908   \n",
       "14                   False     3090693.0           0.802494         0.268412   \n",
       "16                   False     3284381.0           0.756196         0.156462   \n",
       "\n",
       "    Secondary Expected  Primary Expected  Secondary Score  Primary Score  \\\n",
       "4             0.996414          0.411209         0.879681       0.580661   \n",
       "13            0.990717          0.346761         0.899043       0.492870   \n",
       "14            0.984906          0.606032         0.814792       0.442900   \n",
       "16            0.978218          0.266662         0.773034       0.586741   \n",
       "\n",
       "    Consensus ANI  \n",
       "4        0.938553  \n",
       "13       0.932959  \n",
       "14       0.952038  \n",
       "16       0.929656  "
      ]
     },
     "execution_count": 537,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "members_df[['Assembly Accession ID', 'Species', 'Assembly Level', 'Cluster Representative', 'Total Length', \n",
    "            'Secondary Breadth', 'Primary Breadth', 'Secondary Expected', 'Primary Expected', 'Secondary Score', 'Primary Score', 'Consensus ANI']]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 521,
   "id": "263a5180-fd8d-4b78-a9b6-986d89eaa87c",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Index(['Taxonomy ID', 'Assembly Accession ID', 'Source Database',\n",
       "       'Is Representative', 'Assembly Level', 'Organism of Assembly', 'Strain',\n",
       "       'Total Length', 'Downloaded', 'Species', 'Genus Taxid',\n",
       "       'Cluster Representative', 'Cluster Members', 'Secondary Breadth',\n",
       "       'Secondary Expected', 'Secondary Score', 'Secondary Depth',\n",
       "       'Primary Breadth', 'Primary Expected', 'Primary Score', 'Primary Depth',\n",
       "       'Consensus ANI', 'Combined PS and ANI (Sqrt(ANI)xPSx100)',\n",
       "       'Presence/Absence'],\n",
       "      dtype='object')"
      ]
     },
     "execution_count": 521,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "members_df.columns"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "70fbb603-c5da-432d-9bab-86f1623fa48d",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.13"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}