1.1.HDMI_Processing.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "1cc547b2-f0b2-4f18-aae5-fc7c23464632",
   "metadata": {},
   "source": [
    "# HDMI Processing\n",
    "\n",
    "This notebook is for the initial processing of the HDMI sequencing round.\n",
    "\n",
    "We will confirm the quality of the barcodes, extract the coordinates per barcode and save these to disk, as well as prepare some files for later analyses"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d9149578-469f-4682-a250-92ede3b9982a",
   "metadata": {},
   "outputs": [],
   "source": [
    "import pysam\n",
    "import collections as c\n",
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "import pickle\n",
    "import os\n",
    "\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0e39d4f9-5a4f-47b8-865a-5ddbd85509ba",
   "metadata": {},
   "source": [
    "Check overall nucleotide distribution of HDMIs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "34a4bd6f-4892-4c51-bef5-f20cbb00e0fd",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# Generated from demultiplex_per_tile.sh\n",
    "tiles_fastq_dir = '' # /path/to/Demultiplexed_per_tile/output"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0226de14-420a-432b-a5d3-e33cbd8bc863",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "tile_fastqs = [x for x in os.listdir(tiles_fastq_dir) if x.endswith('.fastq.gz')]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "05800056-d2d8-4ee3-81e7-f62f42a7aacc",
   "metadata": {},
   "source": [
    "# Barcode structure check\n",
    "\n",
    "We look at the reads from the first tile of the flowcell and ensure that all of the barcode reads together appear to have the expected sequence (NNVNBVNNVNNVNNVNNVNNVNNVNNVNNNNN)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e7858aab-6553-4212-bac4-2407eb254d72",
   "metadata": {},
   "outputs": [],
   "source": [
    "base_pairs = c.defaultdict(lambda: c.defaultdict(lambda: 0))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "dd33506f-fe6d-46fb-96cd-93bef4309ff5",
   "metadata": {
    "scrolled": true,
    "tags": []
   },
   "outputs": [],
   "source": [
    "# Grab reads from the first tile to check barcode structure\n",
    "with pysam.FastqFile(f'{tiles_fastq_dir}/{tile_fastqs[0]}') as fh:\n",
    "    for en, entry in enumerate(fh):\n",
    "        for n, base in enumerate(entry.sequence):\n",
    "            base_pairs[n][base] += 1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a2de2008-92c3-49ca-8d8b-60f18f7ce067",
   "metadata": {},
   "outputs": [],
   "source": [
    "bp_df = pd.DataFrame(base_pairs)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cf061c33-cb0e-44fa-bb4d-3b5d06e00014",
   "metadata": {},
   "outputs": [],
   "source": [
    "bp_df = bp_df.drop('N')\n",
    "bp_df.columns = [int(x) + 1 for x in bp_df.columns]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "66f24834-4d2a-42a2-99c0-7a6ec0a9de78",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "base_proportions = bp_df/bp_df.sum()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ccdbe954-fa09-4dba-93b0-f314125deb96",
   "metadata": {},
   "source": [
    "# NNVNBVNNVNNVNNVNNVNNVNNVNNVNNNNN"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "19f440b1-1fa6-456b-a902-0b0ecb245a2b",
   "metadata": {},
   "outputs": [],
   "source": [
    "base_proportions.T.plot(kind='line', figsize=(18, 12))\n",
    "plt.grid(True, which='major')\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "704dda45-f44d-47ae-9b7d-53dccbe692e2",
   "metadata": {},
   "source": [
    "# Coordinate Extraction\n",
    "\n",
    "Next, we load each read and check it against the expected patterns, discarding any which do not match.\n",
    "\n",
    "For those reads that pass filtering, we store the barcode and the x and y coordinates (local to the tile) in dictionaries and then to disk using pickle files.\n",
    "\n",
    "We also save a small number of the reads (10000 per tile) to a seperate file allowing fast identifications of tiles -> Nova-ST chips later."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c22e5da0-d2bf-41d8-b552-a0b563d27389",
   "metadata": {},
   "outputs": [],
   "source": [
    "tiles_dir = 'HDMI_Tiles_Data' # Path to the folder to output results"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cffe0572-467e-427c-b093-34845cee4529",
   "metadata": {},
   "outputs": [],
   "source": [
    "if not os.path.isdir(tiles_dir):\n",
    "    os.mkdir(tiles_dir)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4f871ad5-49b7-49bf-a7ba-84567fc3da3e",
   "metadata": {},
   "outputs": [],
   "source": [
    "bc_pattern = 'NNVNBVNNVNNVNNVNNVNNVNNVNNVNNNNN'\n",
    "filter_list = []\n",
    "for n, b in enumerate(bc_pattern):\n",
    "    if b == 'V':\n",
    "        filter_list.append((n, 'T'))\n",
    "    elif b == 'B':\n",
    "        filter_list.append((n, 'A'))\n",
    "    else:\n",
    "        continue\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "08bedb8a-66cf-4b4c-974f-f91d45cb9d91",
   "metadata": {},
   "outputs": [],
   "source": [
    "def is_read_bad(seq, filter_list, allowed_mismatches=0):\n",
    "    bad_count = 0\n",
    "    for n, not_base in filter_list:\n",
    "        if seq[n] == not_base:\n",
    "            bad_count += 1\n",
    "        if bad_count > allowed_mismatches:\n",
    "            return True\n",
    "    return False\n",
    "        "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b0fcdc2d-0a8a-4300-abe0-85befa36b4f3",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "def extract_coords(tileinfo, subset_size = 10000):\n",
    "    fastq_dir, tiles_dir, tile_id = tileinfo\n",
    "    with pysam.FastxFile(f\"{fastq_dir}/{tile_id}.fastq.gz\", persist=False) as fh:\n",
    "        total_bad = 0\n",
    "        tile_data = {}\n",
    "        subset = set()\n",
    "        for en, entry in enumerate(fh):\n",
    "            seq = entry.sequence\n",
    "            # Don't save barcodes which don't match the expected pattern\n",
    "            if is_read_bad(seq, filter_list, allowed_mismatches=0):\n",
    "                total_bad += 1\n",
    "                continue\n",
    "            _, x, y = entry.name.split(\":\")[-3:]\n",
    "            tile_data[seq] = (x, y)\n",
    "            if en < subset_size:\n",
    "                subset.add(seq)\n",
    "        with open(f\"{tiles_dir}/{tile_id}_barcodes.pickle\", \"wb\") as tile_fh:\n",
    "            pickle.dump(tile_data, tile_fh)\n",
    "        return {\n",
    "            \"tile\": tile_id,\n",
    "            \"data\": f\"{tiles_dir}/{tile_id}_barcodes.pickle\",\n",
    "            \"subset\": subset,\n",
    "            \"good_barcodes\": en - total_bad,\n",
    "            \"bad_barcodes\": total_bad,\n",
    "        }"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "75a4a537-4aae-45af-a7bc-da895c9983de",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "all_tiles_data = []\n",
    "for tile in tile_fastqs:\n",
    "    tile = tile.split(\".\")[0]\n",
    "    all_tiles_data.append((tiles_fastq_dir, tiles_dir, tile))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "925ca52e-de93-420b-9c4b-417587ef9d35",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "from multiprocessing import Pool"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e6680a5f-e221-4286-a297-5f4a27bcbdf7",
   "metadata": {
    "scrolled": true,
    "tags": []
   },
   "outputs": [],
   "source": [
    "with Pool(processes=70) as pool:\n",
    "    results = pool.map(extract_coords, all_tiles_data)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "14e16489",
   "metadata": {},
   "outputs": [],
   "source": [
    "total_bad_reads = 0\n",
    "total_good_reads = 0\n",
    "for tile in results:\n",
    "    total_bad_reads += tile['bad_barcodes']\n",
    "    total_good_reads += tile['good_barcodes']\n",
    "\n",
    "print('Final counts:')\n",
    "print(f\"Total good barcodes: {total_good_reads}\")\n",
    "print(f\"Total bad barcodes: {total_bad_reads}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "881df5c2",
   "metadata": {},
   "outputs": [],
   "source": [
    "print(f'{total_bad_reads / (total_good_reads + total_bad_reads) * 100:2f}% of barcodes were bad')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "40c5e700",
   "metadata": {},
   "outputs": [],
   "source": [
    "subset_datas = c.defaultdict(lambda: c.defaultdict(set))\n",
    "for tile in results:\n",
    "    tile_info = tile['tile']\n",
    "    lane = tile_info.split('_')[0]\n",
    "    surface = tile_info.split('_')[1][0]\n",
    "    swath = tile_info.split('_')[1][1]\n",
    "    tile_no = tile_info.split('_')[1][2:]\n",
    "    subset_datas[f'{lane}_{surface}{swath}'][tile_no].update(tile['subset'])\n",
    "\n",
    "\n",
    "for subset in subset_datas:\n",
    "    with open(f'{tiles_dir}/{subset}_subset_barcodes.pickle', 'wb') as fh:\n",
    "        pickle.dump(subset_datas[subset], fh)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python3.9.5 - BiomedBundle",
   "language": "Python",
   "name": "biomedbundle"
  },
  "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.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}