usecases/prepare_small_refgenome.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Prepare a small reference genome\n",
    "\n",
    "We cut out a small ref genome from the CHM13 ref genome and normalize all N's by replacing them by random DNA letters.\n",
    "\n",
    "Since NanoSim itself does not scale well due to locking, we sometimes run it in parallel. \n",
    "NanoSim automatically replaces N's by random letters. To ensure the same reference across all parallel instances, we replace the N letters.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "/Users/maximilianmordig/ont_project_all/ont_project\n"
     ]
    }
   ],
   "source": [
    "%load_ext autoreload\n",
    "%autoreload 2\n",
    "\n",
    "%cd ~/ont_project_all/ont_project\n",
    "%pwd\n",
    "from pathlib import Path\n",
    "\n",
    "import gzip\n",
    "from Bio import SeqIO\n",
    "from Bio.Seq import Seq\n",
    "from Bio import bgzf\n",
    "import pysam\n",
    "import itertools\n",
    "\n",
    "from simreaduntil.shared_utils.nanosim_parsing import case_convert_dna\n",
    "from simreaduntil.shared_utils.utils import tqdm_with_name"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "/Users/maximilianmordig/ont_project_all/tools/bin/minimap2\n"
     ]
    }
   ],
   "source": [
    "import os\n",
    "os.environ[\"PATH\"] = str(Path(\"~/ont_project_all/tools/bin\").expanduser()) + \":\" + os.environ[\"PATH\"]\n",
    "!which minimap2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "run_dir = Path(\"runs/enrich_usecase\")\n",
    "orig_ref_genome_path = run_dir / \"data/chm13v2.0.fa.gz\"\n",
    "ref_genome_path = run_dir / \"data/chm13v2.0_normalized1000000firsttwo.fa.gz\"\n",
    "\n",
    "run_dir.mkdir(exist_ok=True, parents=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# use either the aws cli or curl to download the reference genome\n",
    "!curl https://s3-us-west-2.amazonaws.com/human-pangenomics/T2T/CHM13/assemblies/analysis_set/chm13v2.0.fa.gz --output {orig_ref_genome_path}\n",
    "# or\n",
    "# !aws s3 --no-sign-request cp s3://human-pangenomics/T2T/CHM13/assemblies/analysis_set/chm13v2.0.fa.gz {orig_ref_genome_path}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def normalize_fasta_modified(fasta_gz):\n",
    "    # use SeqIO rather than pysam.FastaFile because it discards the additional chromosome name\n",
    "    with gzip.open(fasta_gz, mode=\"rt\") as f:\n",
    "        for (i, record) in enumerate(SeqIO.parse(f, \"fasta\")):\n",
    "            # yield (record.id, SeqIO.SeqRecord(seq=Seq(case_convert_dna(record.seq)), id=record.id, description=f\"{record.description} normalized\"))\n",
    "            yield (record.id, SeqIO.SeqRecord(seq=Seq(case_convert_dna(record.seq[1_000_000:2_000_000])), id=record.id, description=f\"{record.description} normalized\"))\n",
    "            if i>=1:\n",
    "                break\n",
    "            # print(record.format(\"fasta\"))\n",
    "        \n",
    "# bgzip compression needed for pysam\n",
    "with bgzf.open(ref_genome_path, \"wt\") as f:\n",
    "    SeqIO.write(tqdm_with_name(normalize_fasta_modified(orig_ref_genome_path)), f, \"fasta\")\n",
    "print(f\"Written to '{ref_genome_path}'\")\n",
    "\n",
    "print(f\"Chromosome names: {pysam.FastaFile(ref_genome_path).references}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# # take first chromosome from fasta file and write to new file\n",
    "# with bgzf.open(enrich_ref_genome_path, \"wt\") as f:\n",
    "#     with gzip.open(ref_genome_path, mode=\"rt\") as f_in:\n",
    "#         SeqIO.write(itertools.islice(SeqIO.parse(f_in, \"fasta\"), 1), f, \"fasta\")\n",
    "# print(f\"Chromosome names: {pysam.FastaFile(enrich_ref_genome_path).references}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%bash -e\n",
    "\n",
    "# cd '{run_dir}' # not working\n",
    "cd runs/enrich_usecase\n",
    "export PATH=~/ont_project_all/tools/bin:$PATH && which minimap2 # make minimap2 available\n",
    "# minimap2 -d data/chm13v2.0_normalized1000000first.mmi data/chm13v2.0_normalized1000000first.fa.gz\n",
    "minimap2 -d data/chm13v2.0_normalized1000000firsttwo.mmi data/chm13v2.0_normalized1000000firsttwo.fa.gz"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Enrich full genome"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "run_dir = Path(\"runs/enrich_usecase\")\n",
    "orig_ref_genome_path = run_dir / \"data/chm13v2.0.fa.gz\"\n",
    "ref_genome_path = run_dir / \"data/chm13v2.0_normalized.fa.gz\"\n",
    "\n",
    "run_dir.mkdir(exist_ok=True, parents=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "usage: normalize_fasta [-h] in_fasta_path out_fasta_path\n",
      "\n",
      "Normalize a FASTA file to upper case and replacing N letters\n",
      "\n",
      "positional arguments:\n",
      "  in_fasta_path   input FASTA file\n",
      "  out_fasta_path  output FASTA file\n",
      "\n",
      "options:\n",
      "  -h, --help      show this help message and exit\n",
      "2023-08-07 09:58:39,755 - in_fasta_path=runs/enrich_usecase/data/chm13v2.0.fa.gz --- utils.py:145 (print_args) INFO ##\n",
      "2023-08-07 09:58:39,755 - out_fasta_path=runs/enrich_usecase/data/chm13v2.0_normalized1.fa.gz --- utils.py:145 (print_args) INFO ##\n",
      "2023-08-07 09:58:39,755 - Detected .gz extension for file 'runs/enrich_usecase/data/chm13v2.0_normalized1.fa.gz' --- utils.py:107 (get_fasta_open_method) INFO ##\n",
      "0it [00:00, ?it/s]2023-08-07 09:58:39,788 - Detected .gz extension for file 'runs/enrich_usecase/data/chm13v2.0.fa.gz' --- utils.py:107 (get_fasta_open_method) INFO ##\n",
      "^C\n"
     ]
    }
   ],
   "source": [
    "# from simreaduntil.usecase_helpers.utils import normalize_fasta\n",
    "# normalize_fasta(orig_ref_genome_path, ref_genome_path)\n",
    "\n",
    "!normalize_fasta --help\n",
    "!normalize_fasta {orig_ref_genome_path} {ref_genome_path}\n",
    "\n",
    "# # todo: can be parallelized\n",
    "# with bgzf.open(ref_genome_path, \"wt\") as f:\n",
    "#     SeqIO.write(tqdm_with_name(normalize_fasta_gen(orig_ref_genome_path)), f, \"fasta\")\n",
    "# print(f\"Written to '{ref_genome_path}'\")\n",
    "\n",
    "# print(f\"Chromosome names: {pysam.FastaFile(ref_genome_path).references}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "ont_project_venv",
   "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.10.12"
  },
  "orig_nbformat": 4
 },
 "nbformat": 4,
 "nbformat_minor": 2
}