docs/source/notebooks/preprocess.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Pipeline through notebooks\n",
    "### Preprocess output directory "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "3.7.10 (default, Feb 26 2021, 18:47:35) \n",
      "[GCC 7.3.0]\n"
     ]
    }
   ],
   "source": [
    "import sys\n",
    "print(sys.version)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Import the needed packages"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "from xml.etree import ElementTree as et\n",
    "from dataclasses import dataclass \n",
    "from pathlib import Path\n",
    "from openpyxl import load_workbook\n",
    "import pandas as pd\n",
    "import sys\n",
    "import os\n",
    "import csv\n",
    "from typing import Dict, List, Tuple\n",
    "\n",
    "sys.path.append('/home/hsarkar/code/slideseq-tools/src/')\n",
    "\n",
    "\n",
    "from slideseq.metadata import Manifest\n",
    "from slideseq.metadata import Manifest, split_sample_lanes, validate_run_df\n",
    "from slideseq.config import get_config\n",
    "import slideseq.util.constants as constants\n",
    "from slideseq.util import run_command\n",
    "\n",
    "import logging\n",
    "log = logging.getLogger(__name__)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Classes\n",
    "We followed the same structure created for the original slide-seqtools repository. The classes and the containing \n",
    "datastructure are repeatedly used for parsing the excel file containing information about the experiment. The main goal\n",
    "of the associated dataframe is to iterate over all the directories within different libraries to create inermediate files\n",
    "within the output directory."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_flowcell(run_info: et.ElementTree) -> str:\n",
    "    \"\"\"\n",
    "    Get the flowcell name from RunInfo.xml. This should be more reliable than\n",
    "    trying to parse the run directory, which might have been renamed.\n",
    "\n",
    "    :param run_info: ElementTree representing RunInfo.xml\n",
    "    :return: The flowcell name for this run\n",
    "    \"\"\"\n",
    "    flowcell = run_info.find(\"./Run/Flowcell\").text\n",
    "    return flowcell\n",
    "\n",
    "\n",
    "def get_read_structure(run_info: et.ElementTree) -> str:\n",
    "    \"\"\"\n",
    "    Get read structure from RunInfo.xml. Assumes one index sequence only,\n",
    "    will warn if two are present and ignore the second.\n",
    "\n",
    "    :param run_info: ElementTree representing RunInfo.xml\n",
    "    :return: Formatting string representing the read structure\n",
    "    \"\"\"\n",
    "    read_elems = run_info.findall(\"./Run/Reads/Read[@NumCycles][@Number]\")\n",
    "    read_elems.sort(key=lambda el: int(el.get(\"Number\")))\n",
    "\n",
    "    if len(read_elems) == 4:\n",
    "        # two index reads. We will just ignore the second index\n",
    "        log.warning(\n",
    "            \"This sequencing run has two index reads, we are ignoring the second one\"\n",
    "        )\n",
    "        return \"{}T{}B{}S{}T\".format(*(el.get(\"NumCycles\") for el in read_elems))\n",
    "    elif len(read_elems) != 3:\n",
    "        raise ValueError(f\"Expected three reads, got {len(read_elems)}\")\n",
    "\n",
    "    return \"{}T{}B{}T\".format(*(el.get(\"NumCycles\") for el in read_elems))\n",
    "\n",
    "def get_lanes(run_info: et.ElementTree) -> range:\n",
    "    \"\"\"\n",
    "    Return the lanes for this run, as a range\n",
    "    :param run_info: ElementTree representing RunInfo.xml\n",
    "    :return: range object for the lanes of the run\n",
    "    \"\"\"\n",
    "    lane_count = int(run_info.find(\"./Run/FlowcellLayout[@LaneCount]\").get(\"LaneCount\"))\n",
    "    return range(1, lane_count + 1)\n",
    "\n",
    "\n",
    "@dataclass\n",
    "class RunInfo:\n",
    "    \"\"\"A dataclass to represent RunInfo.xml for a sequencing run (a.k.a. flowcell)\"\"\"\n",
    "\n",
    "    run_dir: Path\n",
    "    flowcell: str\n",
    "    lanes: range\n",
    "    read_structure: str\n",
    "\n",
    "    @property\n",
    "    def demux_log(self):\n",
    "        return f\"demultiplex.{self.flowcell}.L00$TASK_ID.log\"\n",
    "\n",
    "    @property\n",
    "    def basecall_dir(self):\n",
    "        return self.run_dir / \"Data\" / \"Intensities\" / \"BaseCalls\"\n",
    "\n",
    "    def alignment_log(self, lane: int):\n",
    "        return f\"alignment.{self.flowcell}.L{lane:03d}.$TASK_ID.log\"\n",
    "\n",
    "\n",
    "def get_run_info(run_dir: Path) -> RunInfo:\n",
    "    # open the RunInfo.xml file and parse it with element tree\n",
    "    with (run_dir / \"RunInfo.xml\").open() as f:\n",
    "        run_info = et.parse(f)\n",
    "\n",
    "    return RunInfo(\n",
    "        run_dir,\n",
    "        get_flowcell(run_info),\n",
    "        get_lanes(run_info),\n",
    "        get_read_structure(run_info),\n",
    "    )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "run_info = get_run_info(\n",
    "    Path('/home/hsarkar/code/slide-seq-data/P25255/220329_A00689_0513_BHYK23DRXY')\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Run infomation shows that there is one flow-cell and the files are present in the `run_dir`. A library contains two lanes."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "RunInfo(run_dir=PosixPath('/home/hsarkar/code/slide-seq-data/P25255/220329_A00689_0513_BHYK23DRXY'), flowcell='HYK23DRXY', lanes=range(1, 3), read_structure='42T8B41T')"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "run_info"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We now construct the required sheet for running further steps in slide-seq pipeline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "ws = load_workbook('/home/hsarkar/notebooks/slide_seq_analysis/2023/example_metadata.xlsx').active\n",
    "from itertools import islice\n",
    "data = ws.values\n",
    "cols = next(data)\n",
    "data = list(data)\n",
    "#idx = [r[0] for r in data]\n",
    "#data = (islice(r, 1, None) for r in data)\n",
    "df = pd.DataFrame(data, columns=cols)\n",
    "run_df = df[constants.METADATA_COLS]\n",
    "run_df.columns = [c.lower() for c in constants.METADATA_COLS]\n",
    "import numpy as np\n",
    "run_df = run_df.fillna(value=0)\n",
    "run_df = run_df.astype(constants.METADATA_TYPES)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "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>library</th>\n",
       "      <th>date</th>\n",
       "      <th>flowcell</th>\n",
       "      <th>run_name</th>\n",
       "      <th>bclpath</th>\n",
       "      <th>lane</th>\n",
       "      <th>sample_barcode</th>\n",
       "      <th>bead_structure</th>\n",
       "      <th>reference</th>\n",
       "      <th>run_barcodematching</th>\n",
       "      <th>locus_function_list</th>\n",
       "      <th>start_sequence</th>\n",
       "      <th>base_quality</th>\n",
       "      <th>min_transcripts_per_cell</th>\n",
       "      <th>email</th>\n",
       "      <th>puckcaller_path</th>\n",
       "      <th>bead_type</th>\n",
       "      <th>gen_read1_plot</th>\n",
       "      <th>gen_downsampling</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>P25255_1001</td>\n",
       "      <td>2022-28-03</td>\n",
       "      <td>Kharchenko</td>\n",
       "      <td>1</td>\n",
       "      <td>/home/hsarkar/code/slide-seq-data/P25255/22032...</td>\n",
       "      <td>1</td>\n",
       "      <td>TAAGGCGA</td>\n",
       "      <td>8C18X6C1X8M1X</td>\n",
       "      <td>/home/hsarkar/code/slide-seq-data/slide-seq-re...</td>\n",
       "      <td>False</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>hiraksarkar.cs@gmail.com</td>\n",
       "      <td>/home/hsarkar/code/slide-seq-data/P25255/Barco...</td>\n",
       "      <td>0</td>\n",
       "      <td>False</td>\n",
       "      <td>False</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>P25255_1002</td>\n",
       "      <td>2022-28-03</td>\n",
       "      <td>Kharchenko</td>\n",
       "      <td>1</td>\n",
       "      <td>/home/hsarkar/code/slide-seq-data/P25255/22032...</td>\n",
       "      <td>1</td>\n",
       "      <td>CGTACTAG</td>\n",
       "      <td>8C18X6C1X8M1X</td>\n",
       "      <td>/home/hsarkar/code/slide-seq-data/slide-seq-re...</td>\n",
       "      <td>False</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>hiraksarkar.cs@gmail.com</td>\n",
       "      <td>/home/hsarkar/code/slide-seq-data/P25255/Barco...</td>\n",
       "      <td>0</td>\n",
       "      <td>False</td>\n",
       "      <td>False</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>P25255_1003</td>\n",
       "      <td>2022-28-03</td>\n",
       "      <td>Kharchenko</td>\n",
       "      <td>1</td>\n",
       "      <td>/home/hsarkar/code/slide-seq-data/P25255/22032...</td>\n",
       "      <td>1</td>\n",
       "      <td>AGGCAGAA</td>\n",
       "      <td>8C18X6C1X8M1X</td>\n",
       "      <td>/home/hsarkar/code/slide-seq-data/slide-seq-re...</td>\n",
       "      <td>False</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>hiraksarkar.cs@gmail.com</td>\n",
       "      <td>/home/hsarkar/code/slide-seq-data/P25255/Barco...</td>\n",
       "      <td>0</td>\n",
       "      <td>False</td>\n",
       "      <td>False</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>P25255_1004</td>\n",
       "      <td>2022-28-03</td>\n",
       "      <td>Kharchenko</td>\n",
       "      <td>1</td>\n",
       "      <td>/home/hsarkar/code/slide-seq-data/P25255/22032...</td>\n",
       "      <td>1</td>\n",
       "      <td>TCCTGAGC</td>\n",
       "      <td>8C18X6C1X8M1X</td>\n",
       "      <td>/home/hsarkar/code/slide-seq-data/slide-seq-re...</td>\n",
       "      <td>False</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>hiraksarkar.cs@gmail.com</td>\n",
       "      <td>/home/hsarkar/code/slide-seq-data/P25255/Barco...</td>\n",
       "      <td>0</td>\n",
       "      <td>False</td>\n",
       "      <td>False</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>P25255_1001</td>\n",
       "      <td>2022-28-03</td>\n",
       "      <td>Kharchenko</td>\n",
       "      <td>1</td>\n",
       "      <td>/home/hsarkar/code/slide-seq-data/P25255/22032...</td>\n",
       "      <td>2</td>\n",
       "      <td>TAAGGCGA</td>\n",
       "      <td>8C18X6C1X8M1X</td>\n",
       "      <td>/home/hsarkar/code/slide-seq-data/slide-seq-re...</td>\n",
       "      <td>False</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>hiraksarkar.cs@gmail.com</td>\n",
       "      <td>/home/hsarkar/code/slide-seq-data/P25255/Barco...</td>\n",
       "      <td>0</td>\n",
       "      <td>False</td>\n",
       "      <td>False</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>5</th>\n",
       "      <td>P25255_1002</td>\n",
       "      <td>2022-28-03</td>\n",
       "      <td>Kharchenko</td>\n",
       "      <td>1</td>\n",
       "      <td>/home/hsarkar/code/slide-seq-data/P25255/22032...</td>\n",
       "      <td>2</td>\n",
       "      <td>CGTACTAG</td>\n",
       "      <td>8C18X6C1X8M1X</td>\n",
       "      <td>/home/hsarkar/code/slide-seq-data/slide-seq-re...</td>\n",
       "      <td>False</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>hiraksarkar.cs@gmail.com</td>\n",
       "      <td>/home/hsarkar/code/slide-seq-data/P25255/Barco...</td>\n",
       "      <td>0</td>\n",
       "      <td>False</td>\n",
       "      <td>False</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>6</th>\n",
       "      <td>P25255_1003</td>\n",
       "      <td>2022-28-03</td>\n",
       "      <td>Kharchenko</td>\n",
       "      <td>1</td>\n",
       "      <td>/home/hsarkar/code/slide-seq-data/P25255/22032...</td>\n",
       "      <td>2</td>\n",
       "      <td>AGGCAGAA</td>\n",
       "      <td>8C18X6C1X8M1X</td>\n",
       "      <td>/home/hsarkar/code/slide-seq-data/slide-seq-re...</td>\n",
       "      <td>False</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>hiraksarkar.cs@gmail.com</td>\n",
       "      <td>/home/hsarkar/code/slide-seq-data/P25255/Barco...</td>\n",
       "      <td>0</td>\n",
       "      <td>False</td>\n",
       "      <td>False</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>7</th>\n",
       "      <td>P25255_1004</td>\n",
       "      <td>2022-28-03</td>\n",
       "      <td>Kharchenko</td>\n",
       "      <td>1</td>\n",
       "      <td>/home/hsarkar/code/slide-seq-data/P25255/22032...</td>\n",
       "      <td>2</td>\n",
       "      <td>TCCTGAGC</td>\n",
       "      <td>8C18X6C1X8M1X</td>\n",
       "      <td>/home/hsarkar/code/slide-seq-data/slide-seq-re...</td>\n",
       "      <td>False</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>hiraksarkar.cs@gmail.com</td>\n",
       "      <td>/home/hsarkar/code/slide-seq-data/P25255/Barco...</td>\n",
       "      <td>0</td>\n",
       "      <td>False</td>\n",
       "      <td>False</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "       library        date    flowcell  run_name  \\\n",
       "0  P25255_1001  2022-28-03  Kharchenko         1   \n",
       "1  P25255_1002  2022-28-03  Kharchenko         1   \n",
       "2  P25255_1003  2022-28-03  Kharchenko         1   \n",
       "3  P25255_1004  2022-28-03  Kharchenko         1   \n",
       "4  P25255_1001  2022-28-03  Kharchenko         1   \n",
       "5  P25255_1002  2022-28-03  Kharchenko         1   \n",
       "6  P25255_1003  2022-28-03  Kharchenko         1   \n",
       "7  P25255_1004  2022-28-03  Kharchenko         1   \n",
       "\n",
       "                                             bclpath  lane sample_barcode  \\\n",
       "0  /home/hsarkar/code/slide-seq-data/P25255/22032...     1       TAAGGCGA   \n",
       "1  /home/hsarkar/code/slide-seq-data/P25255/22032...     1       CGTACTAG   \n",
       "2  /home/hsarkar/code/slide-seq-data/P25255/22032...     1       AGGCAGAA   \n",
       "3  /home/hsarkar/code/slide-seq-data/P25255/22032...     1       TCCTGAGC   \n",
       "4  /home/hsarkar/code/slide-seq-data/P25255/22032...     2       TAAGGCGA   \n",
       "5  /home/hsarkar/code/slide-seq-data/P25255/22032...     2       CGTACTAG   \n",
       "6  /home/hsarkar/code/slide-seq-data/P25255/22032...     2       AGGCAGAA   \n",
       "7  /home/hsarkar/code/slide-seq-data/P25255/22032...     2       TCCTGAGC   \n",
       "\n",
       "  bead_structure                                          reference  \\\n",
       "0  8C18X6C1X8M1X  /home/hsarkar/code/slide-seq-data/slide-seq-re...   \n",
       "1  8C18X6C1X8M1X  /home/hsarkar/code/slide-seq-data/slide-seq-re...   \n",
       "2  8C18X6C1X8M1X  /home/hsarkar/code/slide-seq-data/slide-seq-re...   \n",
       "3  8C18X6C1X8M1X  /home/hsarkar/code/slide-seq-data/slide-seq-re...   \n",
       "4  8C18X6C1X8M1X  /home/hsarkar/code/slide-seq-data/slide-seq-re...   \n",
       "5  8C18X6C1X8M1X  /home/hsarkar/code/slide-seq-data/slide-seq-re...   \n",
       "6  8C18X6C1X8M1X  /home/hsarkar/code/slide-seq-data/slide-seq-re...   \n",
       "7  8C18X6C1X8M1X  /home/hsarkar/code/slide-seq-data/slide-seq-re...   \n",
       "\n",
       "   run_barcodematching  locus_function_list  start_sequence  base_quality  \\\n",
       "0                False                    0               0             0   \n",
       "1                False                    0               0             0   \n",
       "2                False                    0               0             0   \n",
       "3                False                    0               0             0   \n",
       "4                False                    0               0             0   \n",
       "5                False                    0               0             0   \n",
       "6                False                    0               0             0   \n",
       "7                False                    0               0             0   \n",
       "\n",
       "   min_transcripts_per_cell                     email  \\\n",
       "0                         0  hiraksarkar.cs@gmail.com   \n",
       "1                         0  hiraksarkar.cs@gmail.com   \n",
       "2                         0  hiraksarkar.cs@gmail.com   \n",
       "3                         0  hiraksarkar.cs@gmail.com   \n",
       "4                         0  hiraksarkar.cs@gmail.com   \n",
       "5                         0  hiraksarkar.cs@gmail.com   \n",
       "6                         0  hiraksarkar.cs@gmail.com   \n",
       "7                         0  hiraksarkar.cs@gmail.com   \n",
       "\n",
       "                                     puckcaller_path  bead_type  \\\n",
       "0  /home/hsarkar/code/slide-seq-data/P25255/Barco...          0   \n",
       "1  /home/hsarkar/code/slide-seq-data/P25255/Barco...          0   \n",
       "2  /home/hsarkar/code/slide-seq-data/P25255/Barco...          0   \n",
       "3  /home/hsarkar/code/slide-seq-data/P25255/Barco...          0   \n",
       "4  /home/hsarkar/code/slide-seq-data/P25255/Barco...          0   \n",
       "5  /home/hsarkar/code/slide-seq-data/P25255/Barco...          0   \n",
       "6  /home/hsarkar/code/slide-seq-data/P25255/Barco...          0   \n",
       "7  /home/hsarkar/code/slide-seq-data/P25255/Barco...          0   \n",
       "\n",
       "   gen_read1_plot  gen_downsampling  \n",
       "0           False             False  \n",
       "1           False             False  \n",
       "2           False             False  \n",
       "3           False             False  \n",
       "4           False             False  \n",
       "5           False             False  \n",
       "6           False             False  \n",
       "7           False             False  "
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "run_df"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The above data-frame are manually curated. One needs to check with the actual file-structure\n",
    "before creating such a worksheet. Traditionally, the researchers at Broad work with google sheets and fetch data with \n",
    "the sheet id. \n",
    "\n",
    "Now we would describe the generation process step by step."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Step 1 (Preparing for demultiplexing)\n",
    "\n",
    "Create folders and files for carrying out demultilexing"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "from slideseq.metadata import Manifest\n",
    "from typing import Dict, List, Tuple\n",
    "import logging\n",
    "log = logging.getLogger(__name__)\n",
    "import csv\n",
    "\n",
    "def gen_barcode_file(manifest: Manifest, flowcell: str, lane: int, output_file: Path):\n",
    "    with output_file.open(\"w\") as out:\n",
    "        wtr = csv.writer(out, delimiter=\"\\t\")\n",
    "        wtr.writerow((\"barcode_sequence_1\", \"library_name\", \"barcode_name\"))\n",
    "\n",
    "        for library in manifest.libraries:\n",
    "            if (flowcell, lane) in library.samples:\n",
    "                for barcode in library.samples[flowcell, lane]:\n",
    "                    # we don't write out barcode_name but the column is required\n",
    "                    wtr.writerow((barcode, library.name, \"\"))\n",
    "\n",
    "\n",
    "def gen_library_params(manifest: Manifest, flowcell: str, lane: int, output_file: Path):\n",
    "    with output_file.open(\"w\") as out:\n",
    "        wtr = csv.writer(out, delimiter=\"\\t\")\n",
    "        wtr.writerow((\"OUTPUT\", \"SAMPLE_ALIAS\", \"LIBRARY_NAME\", \"BARCODE_1\"))\n",
    "\n",
    "        for sample in manifest.samples:\n",
    "            if sample.flowcell == flowcell and sample.lane == lane:\n",
    "                # output the uBAM directly to library directory\n",
    "                sample.lane_dir.mkdir(exist_ok=True, parents=True)\n",
    "\n",
    "                for barcode, output_ubam in zip(sample.barcodes, sample.barcode_ubams):\n",
    "                    wtr.writerow((output_ubam, sample.name, sample.name, barcode))\n",
    "\n",
    "\n",
    "def prepare_demux(run_info_list: List[RunInfo], manifest: Manifest):\n",
    "    \"\"\"create a bunch of directories, and write some input files for picard\"\"\"\n",
    "    # Create directories\n",
    "    log.info(\n",
    "        f\"Creating directories in {manifest.workflow_dir} and {manifest.library_dir}\"\n",
    "    )\n",
    "\n",
    "    for run_info in run_info_list:\n",
    "        for lane in run_info.lanes:\n",
    "            output_lane_dir = manifest.workflow_dir / run_info.flowcell / f\"L{lane:03d}\"\n",
    "\n",
    "            output_lane_dir.mkdir(exist_ok=True, parents=True)\n",
    "            (output_lane_dir / \"barcodes\").mkdir(exist_ok=True)\n",
    "\n",
    "            # Generate barcode_params.txt that is needed by ExtractIlluminaBarcodes\n",
    "            gen_barcode_file(\n",
    "                manifest,\n",
    "                run_info.flowcell,\n",
    "                lane,\n",
    "                output_lane_dir / \"barcode_params.txt\",\n",
    "            )\n",
    "\n",
    "            # Generate library_params that is needed by IlluminaBasecallsToSam\n",
    "            gen_library_params(\n",
    "                manifest,\n",
    "                run_info.flowcell,\n",
    "                lane,\n",
    "                output_lane_dir / \"library_params.txt\",\n",
    "            )\n",
    "\n",
    "\n",
    "def validate_demux(manifest: Manifest):\n",
    "    \"\"\"verify that `prepare_demux` was run previously\"\"\"\n",
    "    if not manifest.workflow_dir.exists():\n",
    "        log.error(f\"{manifest.workflow_dir} does not exist\")\n",
    "        return False\n",
    "\n",
    "    for flowcell_dir in manifest.flowcell_dirs:\n",
    "        run_info = get_run_info(flowcell_dir)\n",
    "\n",
    "        # Create directories\n",
    "        log.info(f\"Checking directories in {manifest.workflow_dir / run_info.flowcell}\")\n",
    "        for lane in run_info.lanes:\n",
    "            output_lane_dir = manifest.workflow_dir / run_info.flowcell / f\"L{lane:03d}\"\n",
    "\n",
    "            for p in (\n",
    "                output_lane_dir,\n",
    "                output_lane_dir / \"barcodes\",\n",
    "                output_lane_dir / \"barcode_params.txt\",\n",
    "                output_lane_dir / \"library_params.txt\",\n",
    "            ):\n",
    "                if not p.exists():\n",
    "                    log.error(f\"{p} does not exist, demux looks incomplete\")\n",
    "                    return False\n",
    "\n",
    "    return True\n",
    "\n",
    "\n",
    "def validate_alignment(manifest: Manifest, n_libraries: int):\n",
    "    \"\"\"verify that alignment was run and output is present\"\"\"\n",
    "\n",
    "    for i in range(n_libraries):\n",
    "        library = manifest.get_library(i)\n",
    "\n",
    "        for p_list in (\n",
    "            library.polya_filtering_summaries,\n",
    "            library.star_logs,\n",
    "            library.alignment_pickles,\n",
    "            library.processed_bams,\n",
    "        ):\n",
    "            for p in p_list:\n",
    "                if not p.exists():\n",
    "                    log.error(f\"{p} does not exist, alignment looks incomplete\")\n",
    "                    return False\n",
    "    else:\n",
    "        return True\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "run_name = '1'\n",
    "output_dir = Path('/home/hsarkar/code/slide-seq-data/P25255/processed_data')/\\\n",
    "    Path(run_name)\n",
    "flowcell_dirs = sorted(Path(fd) for fd in set(run_df.bclpath))\n",
    "manifest_file = output_dir / \"manifest.yaml\"\n",
    "metadata_file = output_dir / \"metadata.csv\"\n",
    "\n",
    "run_info_list = [ run_info ]\n",
    "config = get_config()\n",
    "\n",
    "manifest = Manifest(\n",
    "            run_name=run_name,\n",
    "            flowcell_dirs=flowcell_dirs,\n",
    "            workflow_dir=output_dir,\n",
    "            library_dir=Path('/home/hsarkar/code/slide-seq-data/P25255/processed_data/libraries'),\n",
    "            metadata_file=metadata_file,\n",
    "            metadata=split_sample_lanes(run_df, run_info_list),\n",
    "            email_addresses=sorted(\n",
    "                set(e.strip() for v in run_df.email for e in v.split(\",\"))\n",
    "            ),\n",
    "        )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "At this point we have read the config file stored within the `slideseq-tools` directory, curtailed to our need"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Config(picard=PosixPath('/home/hsarkar/code/slideseq-tools/soft/picard/build/libs/picard.jar'), dropseq_dir=PosixPath('/home/hsarkar/code/slideseq-tools/soft/Drop-seq_tools-2.5.1'), reference_dir=PosixPath('/broad/macosko/reference'), workflow_dir=PosixPath('/broad/macosko/data/workflows/flowcell'), library_dir=PosixPath('/broad/macosko/data/libraries'), gsecret_name='projects/velina-208320/secrets/sequencing-credentials/versions/latest', gsheet_id='1kwnKrkbl80LyE9lND0UZZJXipL4yfBbGjkTe6hcwJic', worksheet='Experiment Log', gs_path=PosixPath('macosko_data/libraries'))"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "config"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Manifest(run_name='1', flowcell_dirs=[PosixPath('/home/hsarkar/code/slide-seq-data/P25255/220329_A00689_0513_BHYK23DRXY')], workflow_dir=PosixPath('/home/hsarkar/code/slide-seq-data/P25255/processed_data/1'), library_dir=PosixPath('/home/hsarkar/code/slide-seq-data/P25255/processed_data/libraries'), metadata_file=PosixPath('/home/hsarkar/code/slide-seq-data/P25255/processed_data/1/metadata.csv'), metadata=       library        date   flowcell  run_name  \\\n",
       "0  P25255_1001  2022-28-03  HYK23DRXY         1   \n",
       "1  P25255_1002  2022-28-03  HYK23DRXY         1   \n",
       "2  P25255_1003  2022-28-03  HYK23DRXY         1   \n",
       "3  P25255_1004  2022-28-03  HYK23DRXY         1   \n",
       "4  P25255_1001  2022-28-03  HYK23DRXY         1   \n",
       "5  P25255_1002  2022-28-03  HYK23DRXY         1   \n",
       "6  P25255_1003  2022-28-03  HYK23DRXY         1   \n",
       "7  P25255_1004  2022-28-03  HYK23DRXY         1   \n",
       "\n",
       "                                             bclpath  lane sample_barcode  \\\n",
       "0  /home/hsarkar/code/slide-seq-data/P25255/22032...     1       TAAGGCGA   \n",
       "1  /home/hsarkar/code/slide-seq-data/P25255/22032...     1       CGTACTAG   \n",
       "2  /home/hsarkar/code/slide-seq-data/P25255/22032...     1       AGGCAGAA   \n",
       "3  /home/hsarkar/code/slide-seq-data/P25255/22032...     1       TCCTGAGC   \n",
       "4  /home/hsarkar/code/slide-seq-data/P25255/22032...     2       TAAGGCGA   \n",
       "5  /home/hsarkar/code/slide-seq-data/P25255/22032...     2       CGTACTAG   \n",
       "6  /home/hsarkar/code/slide-seq-data/P25255/22032...     2       AGGCAGAA   \n",
       "7  /home/hsarkar/code/slide-seq-data/P25255/22032...     2       TCCTGAGC   \n",
       "\n",
       "  bead_structure                                          reference  \\\n",
       "0  8C18X6C1X8M1X  /home/hsarkar/code/slide-seq-data/slide-seq-re...   \n",
       "1  8C18X6C1X8M1X  /home/hsarkar/code/slide-seq-data/slide-seq-re...   \n",
       "2  8C18X6C1X8M1X  /home/hsarkar/code/slide-seq-data/slide-seq-re...   \n",
       "3  8C18X6C1X8M1X  /home/hsarkar/code/slide-seq-data/slide-seq-re...   \n",
       "4  8C18X6C1X8M1X  /home/hsarkar/code/slide-seq-data/slide-seq-re...   \n",
       "5  8C18X6C1X8M1X  /home/hsarkar/code/slide-seq-data/slide-seq-re...   \n",
       "6  8C18X6C1X8M1X  /home/hsarkar/code/slide-seq-data/slide-seq-re...   \n",
       "7  8C18X6C1X8M1X  /home/hsarkar/code/slide-seq-data/slide-seq-re...   \n",
       "\n",
       "   run_barcodematching  locus_function_list  start_sequence  base_quality  \\\n",
       "0                False                    0               0             0   \n",
       "1                False                    0               0             0   \n",
       "2                False                    0               0             0   \n",
       "3                False                    0               0             0   \n",
       "4                False                    0               0             0   \n",
       "5                False                    0               0             0   \n",
       "6                False                    0               0             0   \n",
       "7                False                    0               0             0   \n",
       "\n",
       "   min_transcripts_per_cell                     email  \\\n",
       "0                         0  hiraksarkar.cs@gmail.com   \n",
       "1                         0  hiraksarkar.cs@gmail.com   \n",
       "2                         0  hiraksarkar.cs@gmail.com   \n",
       "3                         0  hiraksarkar.cs@gmail.com   \n",
       "4                         0  hiraksarkar.cs@gmail.com   \n",
       "5                         0  hiraksarkar.cs@gmail.com   \n",
       "6                         0  hiraksarkar.cs@gmail.com   \n",
       "7                         0  hiraksarkar.cs@gmail.com   \n",
       "\n",
       "                                     puckcaller_path  bead_type  \\\n",
       "0  /home/hsarkar/code/slide-seq-data/P25255/Barco...          0   \n",
       "1  /home/hsarkar/code/slide-seq-data/P25255/Barco...          0   \n",
       "2  /home/hsarkar/code/slide-seq-data/P25255/Barco...          0   \n",
       "3  /home/hsarkar/code/slide-seq-data/P25255/Barco...          0   \n",
       "4  /home/hsarkar/code/slide-seq-data/P25255/Barco...          0   \n",
       "5  /home/hsarkar/code/slide-seq-data/P25255/Barco...          0   \n",
       "6  /home/hsarkar/code/slide-seq-data/P25255/Barco...          0   \n",
       "7  /home/hsarkar/code/slide-seq-data/P25255/Barco...          0   \n",
       "\n",
       "   gen_read1_plot  gen_downsampling  \n",
       "0           False             False  \n",
       "1           False             False  \n",
       "2           False             False  \n",
       "3           False             False  \n",
       "4           False             False  \n",
       "5           False             False  \n",
       "6           False             False  \n",
       "7           False             False  , email_addresses=['hiraksarkar.cs@gmail.com'])"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "manifest"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Create folders using the manifest files"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "prepare_demux(run_info_list, manifest)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Step 2 Demultiplex (Has to be run in command line) \n",
    "`demultiplex.sh` should contain the locations that can be determined from the worksheet and the manifest file\n",
    "```\n",
    "TMP_DIR=\"/home/hsarkar/code/slide-seq-data/P25255/processed_data/tmp\"\n",
    "PICARD_JAR=\"/home/hsarkar/code/slideseq-tools/soft/picard/build/libs/picard.jar\"\n",
    "SGE_TASK_ID=\"2\"\n",
    "BASECALLS_DIR=\"/home/hsarkar/code/slide-seq-data/P25255/220329_A00689_0513_BHYK23DRXY/Data/Intensities/BaseCalls\"\n",
    "READ_STRUCTURE=\"42T8B41T\"\n",
    "OUTPUT_DIR=\"/home/hsarkar/code/slide-seq-data/P25255/processed_data/1\"\n",
    "FLOWCELL=\"HYK23DRXY\"\n",
    "```\n",
    "Run the following script\n",
    "```\n",
    "./scripts/demultiplex.sh\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can see the created files here"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[0,HYK23DRXY,1]\n",
      "/home/hsarkar/code/slide-seq-data/P25255/processed_data/libraries/2022-28-03_P25255_1001/HYK23DRXY/L001/P25255_1001.TAAGGCGA.unmapped.bam\n",
      "[0,HYK23DRXY,2]\n",
      "/home/hsarkar/code/slide-seq-data/P25255/processed_data/libraries/2022-28-03_P25255_1001/HYK23DRXY/L002/P25255_1001.TAAGGCGA.unmapped.bam\n",
      "[1,HYK23DRXY,1]\n",
      "/home/hsarkar/code/slide-seq-data/P25255/processed_data/libraries/2022-28-03_P25255_1002/HYK23DRXY/L001/P25255_1002.CGTACTAG.unmapped.bam\n",
      "[1,HYK23DRXY,2]\n",
      "/home/hsarkar/code/slide-seq-data/P25255/processed_data/libraries/2022-28-03_P25255_1002/HYK23DRXY/L002/P25255_1002.CGTACTAG.unmapped.bam\n",
      "[2,HYK23DRXY,1]\n",
      "/home/hsarkar/code/slide-seq-data/P25255/processed_data/libraries/2022-28-03_P25255_1003/HYK23DRXY/L001/P25255_1003.AGGCAGAA.unmapped.bam\n",
      "[2,HYK23DRXY,2]\n",
      "/home/hsarkar/code/slide-seq-data/P25255/processed_data/libraries/2022-28-03_P25255_1003/HYK23DRXY/L002/P25255_1003.AGGCAGAA.unmapped.bam\n",
      "[3,HYK23DRXY,1]\n",
      "/home/hsarkar/code/slide-seq-data/P25255/processed_data/libraries/2022-28-03_P25255_1004/HYK23DRXY/L001/P25255_1004.TCCTGAGC.unmapped.bam\n",
      "[3,HYK23DRXY,2]\n",
      "/home/hsarkar/code/slide-seq-data/P25255/processed_data/libraries/2022-28-03_P25255_1004/HYK23DRXY/L002/P25255_1004.TCCTGAGC.unmapped.bam\n"
     ]
    }
   ],
   "source": [
    "n_libraries = len(list(manifest.libraries))\n",
    "for library_index in range(n_libraries):\n",
    "    for run_info in run_info_list:\n",
    "        for lane in run_info.lanes:\n",
    "\n",
    "            print('[{},{},{}]'.format(\n",
    "                library_index,run_info.flowcell, lane))\n",
    "            sample = manifest.get_sample(\n",
    "                library_index, run_info.flowcell, lane)\n",
    "\n",
    "            \n",
    "            for barcode_ubam in sample.barcode_ubams:\n",
    "                print(barcode_ubam)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "After demultiplexing it multiple unmapped bamfiles will be created for each barcode within lane. If there is one `sample_barcode` for each lane then it will only create one file.\n",
    "\n",
    "## Step 3 Merge the barcode.unmapped ubams and "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "n_libraries = len(list(manifest.libraries))\n",
    "for library_index in range(n_libraries):\n",
    "    for run_info in run_info_list:\n",
    "        for lane in run_info.lanes:\n",
    "\n",
    "            print('[{},{},{}]'.format(\n",
    "                library_index,run_info.flowcell, lane))\n",
    "            sample = manifest.get_sample(\n",
    "                library_index, run_info.flowcell, lane)\n",
    "\n",
    "            \n",
    "            for barcode_ubam in sample.barcode_ubams:\n",
    "                print(barcode_ubam)\n",
    "            bead_structure = sample.get_bead_structure()\n",
    "            xc_range = \":\".join(f\"{i}-{j}\" for c, i, j in bead_structure if c == \"C\")\n",
    "            xm_range = \":\".join(f\"{i}-{j}\" for c, i, j in bead_structure if c == \"M\")\n",
    "            \n",
    "            if not sample.raw_ubam.exists():\n",
    "                if len(sample.barcode_ubams) > 1:\n",
    "                    cmd = config.picard_cmd(\"MergeSamFiles\", manifest.tmp_dir)\n",
    "                    cmd.extend(\n",
    "                        [\n",
    "                            \"--OUTPUT\",\n",
    "                            sample.raw_ubam,\n",
    "                            \"--SORT_ORDER\",\n",
    "                            \"unsorted\",\n",
    "                            \"--ASSUME_SORTED\",\n",
    "                            \"true\",\n",
    "                        ]\n",
    "                    )\n",
    "\n",
    "                    for ubam_file in sample.barcode_ubams:\n",
    "                        cmd.extend([\"--INPUT\", ubam_file])\n",
    "                    run_command(cmd, \"MergeBamFiles\", sample)\n",
    "                else:\n",
    "                    if sample.raw_ubam.exists():\n",
    "                        print(sample.raw_ubam)\n",
    "                    else:\n",
    "                        os.rename(sample.barcode_ubams[0], sample.raw_ubam)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Step 4 recover whitelisted barcodes (Takes long run from command prompt)\n",
    "\n",
    "### Step 4a: Extract the barcodes (42bp) from the bam files"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[0,HYK23DRXY,1]\n",
      "[0,HYK23DRXY,2]\n",
      "[1,HYK23DRXY,1]\n",
      "[1,HYK23DRXY,2]\n",
      "[2,HYK23DRXY,1]\n",
      "[2,HYK23DRXY,2]\n",
      "[3,HYK23DRXY,1]\n",
      "[3,HYK23DRXY,2]\n"
     ]
    }
   ],
   "source": [
    "fp = open('../../../scripts/write_shell.txt','w')\n",
    "for library_index in range(n_libraries):\n",
    "    for run_info in run_info_list:\n",
    "        for lane in run_info.lanes:\n",
    "            print('[{},{},{}]'.format(\n",
    "                library_index,run_info.flowcell, lane))\n",
    "            sample = manifest.get_sample(\n",
    "                library_index, run_info.flowcell, lane)\n",
    "            if sample is None:\n",
    "                print('sample not present')\n",
    "            \n",
    "            awk_command = \"awk 'NR%2 {{print > \\\"{}\\\"}}'\".format(str(sample.raw_barcode))\n",
    "            cmd = [\n",
    "                \"samtools\",\n",
    "                \"view\",\n",
    "                str(sample.raw_ubam),\n",
    "                \"|\",\n",
    "                \"cut -f10 |\",\n",
    "                awk_command\n",
    "            ]\n",
    "            fp.write(\"{}\\n\".format(' '.join(cmd)))\n",
    "fp.close()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Run the commands\n",
    "```\n",
    "parallel -j 8 < /home/hsarkar/code/slideseq-tools/write_shell.txt\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Step 4b: Create whitelist (Ideally one should borrow information from image analysis)\n",
    "\n",
    "**All the follwoing commands should be run from command line and not notebook**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The barcodes will be written using the following script\n",
    "\n",
    "```\n",
    "python scripts/createwhitelist.py\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Step 5: Create the refined bam files\n",
    "\n",
    "```\n",
    "python scripts/recoverbam.py\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Step 6: Run the alignment step\n",
    "\n",
    "```\n",
    "python scripts/alignment.py\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Step 7: Run processing pipeline\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "```\n",
    "python scripts/processing.py\n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "/home/hsarkar/code/slide-seq-pipeline/docs/source/notebooks\n"
     ]
    }
   ],
   "source": [
    "!pwd"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": []
  }
 ],
 "metadata": {
  "interpreter": {
   "hash": "78b42e468616091eed2c9f1a76c2fa153b13dae9da70b8ad40e386413158e954"
  },
  "kernelspec": {
   "display_name": "Python 3.7.10 ('jupyter-lab')",
   "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.7.10"
  },
  "orig_nbformat": 4
 },
 "nbformat": 4,
 "nbformat_minor": 2
}