notebooks/GlucoseGalactoseExp/GlucoseGalactoseExp_process_mito_fragments.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "fefbe9b6",
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "import numpy as np\n",
    "import matplotlib\n",
    "import matplotlib.pyplot as plt\n",
    "import os"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "42d27575",
   "metadata": {},
   "source": [
    "# Get mito genes location"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e1f2edc7",
   "metadata": {},
   "source": [
    "We subsetted the hg38 reference genome GTF file to obtain the location of chrM genes."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "id": "0f814b36",
   "metadata": {},
   "outputs": [],
   "source": [
    "# relative to repo\n",
    "mito_genes_path = '../../misc/mito_genes.csv'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "id": "17bb34d3",
   "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>seqname</th>\n",
       "      <th>source</th>\n",
       "      <th>feature</th>\n",
       "      <th>start</th>\n",
       "      <th>end</th>\n",
       "      <th>score</th>\n",
       "      <th>strand</th>\n",
       "      <th>frame</th>\n",
       "      <th>gene_id</th>\n",
       "      <th>gene_version</th>\n",
       "      <th>...</th>\n",
       "      <th>transcript_type</th>\n",
       "      <th>transcript_name</th>\n",
       "      <th>transcript_support_level</th>\n",
       "      <th>havana_transcript</th>\n",
       "      <th>exon_number</th>\n",
       "      <th>exon_id</th>\n",
       "      <th>exon_version</th>\n",
       "      <th>protein_id</th>\n",
       "      <th>ccdsid</th>\n",
       "      <th>ont</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>2765275</th>\n",
       "      <td>chrM</td>\n",
       "      <td>ENSEMBL</td>\n",
       "      <td>gene</td>\n",
       "      <td>3307</td>\n",
       "      <td>4262</td>\n",
       "      <td>NaN</td>\n",
       "      <td>+</td>\n",
       "      <td>0</td>\n",
       "      <td>ENSG00000198888</td>\n",
       "      <td>2</td>\n",
       "      <td>...</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2765279</th>\n",
       "      <td>chrM</td>\n",
       "      <td>ENSEMBL</td>\n",
       "      <td>gene</td>\n",
       "      <td>4470</td>\n",
       "      <td>5511</td>\n",
       "      <td>NaN</td>\n",
       "      <td>+</td>\n",
       "      <td>0</td>\n",
       "      <td>ENSG00000198763</td>\n",
       "      <td>3</td>\n",
       "      <td>...</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2765283</th>\n",
       "      <td>chrM</td>\n",
       "      <td>ENSEMBL</td>\n",
       "      <td>gene</td>\n",
       "      <td>5904</td>\n",
       "      <td>7445</td>\n",
       "      <td>NaN</td>\n",
       "      <td>+</td>\n",
       "      <td>0</td>\n",
       "      <td>ENSG00000198804</td>\n",
       "      <td>2</td>\n",
       "      <td>...</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2765288</th>\n",
       "      <td>chrM</td>\n",
       "      <td>ENSEMBL</td>\n",
       "      <td>gene</td>\n",
       "      <td>7586</td>\n",
       "      <td>8269</td>\n",
       "      <td>NaN</td>\n",
       "      <td>+</td>\n",
       "      <td>0</td>\n",
       "      <td>ENSG00000198712</td>\n",
       "      <td>1</td>\n",
       "      <td>...</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2765295</th>\n",
       "      <td>chrM</td>\n",
       "      <td>ENSEMBL</td>\n",
       "      <td>gene</td>\n",
       "      <td>8366</td>\n",
       "      <td>8572</td>\n",
       "      <td>NaN</td>\n",
       "      <td>+</td>\n",
       "      <td>0</td>\n",
       "      <td>ENSG00000228253</td>\n",
       "      <td>1</td>\n",
       "      <td>...</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2765302</th>\n",
       "      <td>chrM</td>\n",
       "      <td>ENSEMBL</td>\n",
       "      <td>gene</td>\n",
       "      <td>8527</td>\n",
       "      <td>9207</td>\n",
       "      <td>NaN</td>\n",
       "      <td>+</td>\n",
       "      <td>0</td>\n",
       "      <td>ENSG00000198899</td>\n",
       "      <td>2</td>\n",
       "      <td>...</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2765309</th>\n",
       "      <td>chrM</td>\n",
       "      <td>ENSEMBL</td>\n",
       "      <td>gene</td>\n",
       "      <td>9207</td>\n",
       "      <td>9990</td>\n",
       "      <td>NaN</td>\n",
       "      <td>+</td>\n",
       "      <td>0</td>\n",
       "      <td>ENSG00000198938</td>\n",
       "      <td>2</td>\n",
       "      <td>...</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2765314</th>\n",
       "      <td>chrM</td>\n",
       "      <td>ENSEMBL</td>\n",
       "      <td>gene</td>\n",
       "      <td>10059</td>\n",
       "      <td>10404</td>\n",
       "      <td>NaN</td>\n",
       "      <td>+</td>\n",
       "      <td>0</td>\n",
       "      <td>ENSG00000198840</td>\n",
       "      <td>2</td>\n",
       "      <td>...</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2765318</th>\n",
       "      <td>chrM</td>\n",
       "      <td>ENSEMBL</td>\n",
       "      <td>gene</td>\n",
       "      <td>10470</td>\n",
       "      <td>10766</td>\n",
       "      <td>NaN</td>\n",
       "      <td>+</td>\n",
       "      <td>0</td>\n",
       "      <td>ENSG00000212907</td>\n",
       "      <td>2</td>\n",
       "      <td>...</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2765325</th>\n",
       "      <td>chrM</td>\n",
       "      <td>ENSEMBL</td>\n",
       "      <td>gene</td>\n",
       "      <td>10760</td>\n",
       "      <td>12137</td>\n",
       "      <td>NaN</td>\n",
       "      <td>+</td>\n",
       "      <td>0</td>\n",
       "      <td>ENSG00000198886</td>\n",
       "      <td>2</td>\n",
       "      <td>...</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2765330</th>\n",
       "      <td>chrM</td>\n",
       "      <td>ENSEMBL</td>\n",
       "      <td>gene</td>\n",
       "      <td>12337</td>\n",
       "      <td>14148</td>\n",
       "      <td>NaN</td>\n",
       "      <td>+</td>\n",
       "      <td>0</td>\n",
       "      <td>ENSG00000198786</td>\n",
       "      <td>2</td>\n",
       "      <td>...</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2765336</th>\n",
       "      <td>chrM</td>\n",
       "      <td>ENSEMBL</td>\n",
       "      <td>gene</td>\n",
       "      <td>14149</td>\n",
       "      <td>14673</td>\n",
       "      <td>NaN</td>\n",
       "      <td>-</td>\n",
       "      <td>0</td>\n",
       "      <td>ENSG00000198695</td>\n",
       "      <td>2</td>\n",
       "      <td>...</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2765341</th>\n",
       "      <td>chrM</td>\n",
       "      <td>ENSEMBL</td>\n",
       "      <td>gene</td>\n",
       "      <td>14747</td>\n",
       "      <td>15887</td>\n",
       "      <td>NaN</td>\n",
       "      <td>+</td>\n",
       "      <td>0</td>\n",
       "      <td>ENSG00000198727</td>\n",
       "      <td>2</td>\n",
       "      <td>...</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>13 rows × 28 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "        seqname   source feature  start    end  score strand  frame  \\\n",
       "2765275    chrM  ENSEMBL    gene   3307   4262    NaN      +      0   \n",
       "2765279    chrM  ENSEMBL    gene   4470   5511    NaN      +      0   \n",
       "2765283    chrM  ENSEMBL    gene   5904   7445    NaN      +      0   \n",
       "2765288    chrM  ENSEMBL    gene   7586   8269    NaN      +      0   \n",
       "2765295    chrM  ENSEMBL    gene   8366   8572    NaN      +      0   \n",
       "2765302    chrM  ENSEMBL    gene   8527   9207    NaN      +      0   \n",
       "2765309    chrM  ENSEMBL    gene   9207   9990    NaN      +      0   \n",
       "2765314    chrM  ENSEMBL    gene  10059  10404    NaN      +      0   \n",
       "2765318    chrM  ENSEMBL    gene  10470  10766    NaN      +      0   \n",
       "2765325    chrM  ENSEMBL    gene  10760  12137    NaN      +      0   \n",
       "2765330    chrM  ENSEMBL    gene  12337  14148    NaN      +      0   \n",
       "2765336    chrM  ENSEMBL    gene  14149  14673    NaN      -      0   \n",
       "2765341    chrM  ENSEMBL    gene  14747  15887    NaN      +      0   \n",
       "\n",
       "                 gene_id  gene_version  ... transcript_type transcript_name  \\\n",
       "2765275  ENSG00000198888             2  ...             NaN             NaN   \n",
       "2765279  ENSG00000198763             3  ...             NaN             NaN   \n",
       "2765283  ENSG00000198804             2  ...             NaN             NaN   \n",
       "2765288  ENSG00000198712             1  ...             NaN             NaN   \n",
       "2765295  ENSG00000228253             1  ...             NaN             NaN   \n",
       "2765302  ENSG00000198899             2  ...             NaN             NaN   \n",
       "2765309  ENSG00000198938             2  ...             NaN             NaN   \n",
       "2765314  ENSG00000198840             2  ...             NaN             NaN   \n",
       "2765318  ENSG00000212907             2  ...             NaN             NaN   \n",
       "2765325  ENSG00000198886             2  ...             NaN             NaN   \n",
       "2765330  ENSG00000198786             2  ...             NaN             NaN   \n",
       "2765336  ENSG00000198695             2  ...             NaN             NaN   \n",
       "2765341  ENSG00000198727             2  ...             NaN             NaN   \n",
       "\n",
       "         transcript_support_level havana_transcript  exon_number  exon_id  \\\n",
       "2765275                       NaN               NaN          NaN      NaN   \n",
       "2765279                       NaN               NaN          NaN      NaN   \n",
       "2765283                       NaN               NaN          NaN      NaN   \n",
       "2765288                       NaN               NaN          NaN      NaN   \n",
       "2765295                       NaN               NaN          NaN      NaN   \n",
       "2765302                       NaN               NaN          NaN      NaN   \n",
       "2765309                       NaN               NaN          NaN      NaN   \n",
       "2765314                       NaN               NaN          NaN      NaN   \n",
       "2765318                       NaN               NaN          NaN      NaN   \n",
       "2765325                       NaN               NaN          NaN      NaN   \n",
       "2765330                       NaN               NaN          NaN      NaN   \n",
       "2765336                       NaN               NaN          NaN      NaN   \n",
       "2765341                       NaN               NaN          NaN      NaN   \n",
       "\n",
       "         exon_version  protein_id  ccdsid  ont  \n",
       "2765275           NaN         NaN     NaN  NaN  \n",
       "2765279           NaN         NaN     NaN  NaN  \n",
       "2765283           NaN         NaN     NaN  NaN  \n",
       "2765288           NaN         NaN     NaN  NaN  \n",
       "2765295           NaN         NaN     NaN  NaN  \n",
       "2765302           NaN         NaN     NaN  NaN  \n",
       "2765309           NaN         NaN     NaN  NaN  \n",
       "2765314           NaN         NaN     NaN  NaN  \n",
       "2765318           NaN         NaN     NaN  NaN  \n",
       "2765325           NaN         NaN     NaN  NaN  \n",
       "2765330           NaN         NaN     NaN  NaN  \n",
       "2765336           NaN         NaN     NaN  NaN  \n",
       "2765341           NaN         NaN     NaN  NaN  \n",
       "\n",
       "[13 rows x 28 columns]"
      ]
     },
     "execution_count": 24,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Load mito genes\n",
    "mito_genes = pd.read_csv(f'{mito_genes_path}', index_col=0)\n",
    "mito_genes"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "916a26a6",
   "metadata": {},
   "source": [
    "# 00 Getting ATAC Fragments"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "adfe6bf8",
   "metadata": {},
   "source": [
    "**NOTE:** We obtained our mito fragments from the fastq files by running Cellranger ARC for multiomic data. \n",
    "\n",
    "**However, we made the following changes to run Cellranger ARC:**\n",
    "- Used a custom CellRanger ARC reference genome. We blacklisted nuclear regions that were homologous to the mitochondrial genome in our reference genome. We followed the steps here (https://github.com/caleblareau/mitoblacklist). This prevents mito reads from being multimapped and thrown out and increased the coverage of our mitochondrial reads.\n",
    "- Move `chrM` from `non_nuclear_contigs` to `primary_contigs` in the CellRanger ARC configuration file. This ensures the mitochondrial reads are not discarded during processing.\n",
    "\n",
    "After processing, CellRanger ARC provides us with `atac_fragments.tsv.gz` as output. This file contains information on all the ATAC fragments sequenced. This file is used as input for the next steps of processing:"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "09f80db1",
   "metadata": {},
   "source": [
    "# 01 Subset Multiome ATAC fragments to just Mito"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9c9ac4d7",
   "metadata": {},
   "source": [
    "The steps in this file are essentially the same processing steps as we did for the 8HetClonesHashtagMix samples. We just use a for loop to process our 3 samples (Glucose + Pyruvate, Glucose - Pyruvate, Galactose)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "id": "4f5f43b3",
   "metadata": {},
   "outputs": [],
   "source": [
    "samples = [\n",
    "    \"YF-2721_48mix_25dGlu_3d_mito_hLLLdogma_multiome\", # Glucose + Pyruvate\n",
    "    \"YF-2721_48mix_10dGal_3d_mito_hLLLdogma_multiome\", # Galactose\n",
    "    \"YF-2721_48mix_25dGlu_NoPyr_3d_mito_hLLLdogma_multiome\", # Glucose - Pyruvate\n",
    "]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9ff85db0",
   "metadata": {},
   "source": [
    "Our file structure may be different from yours, just point to the correct files. What you will need is:\n",
    "- input path to `atac_fragments.tsv.gz` which is the output of Cellranger ARC.\n",
    "- output path to store results"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8161d4a6",
   "metadata": {},
   "outputs": [],
   "source": [
    "# edit to your file structure\n",
    "data_dir = '/data/peer/landm/Projects/sail-dogma/yi-clone-samples/data/'\n",
    "cache_dir = '/data/peer/landm/Projects/sail-dogma/yi-clone-samples/cache/'\n",
    "\n",
    "for s in samples:\n",
    "    sample_dir = f\"{data_dir}/{s}/cr-arc-results/\"\n",
    "    \n",
    "    # READ INPUT\n",
    "    fragments = pd.read_csv(f\"{sample_dir}atac_fragments.tsv.gz\", compression='gzip', header = None, sep='\\t', quotechar='\"', comment = '#')\n",
    "    mito_fragments = fragments[fragments[0] == 'chrM']\n",
    "    print(f'{mito_fragments.shape[0]} mitochondrial fragments')\n",
    "\n",
    "    # Annotate fragments by where they fall on the genome\n",
    "    mito_fragments['gene'] = 'chrM'\n",
    "    mito_fragments.columns = ['chromosome', 'start', 'end', 'barcode', 'umis', 'gene']\n",
    "    for row in mito_genes.iterrows():\n",
    "        row = row[1]\n",
    "        start = row['start']\n",
    "        end = row['end']\n",
    "\n",
    "        mito_fragments.loc[(mito_fragments['start']>=start) & (mito_fragments['end'] < end), 'gene'] = row['gene_name']\n",
    "    # mito_fragments.head()\n",
    "\n",
    "    mito_dir = f\"{cache_dir}{s}/mito/\"\n",
    "    if not os.path.exists(mito_dir):\n",
    "        os.mkdir(mito_dir)\n",
    "        \n",
    "    # WRITE OUTPUT \n",
    "    mito_fragments.to_csv(f'{mito_dir}mito_fragments.csv')\n",
    "\n",
    "    # Get counts for each gene per cell barcode\n",
    "    mito_counts = mito_fragments[['barcode', 'gene']].groupby('barcode').value_counts().unstack(fill_value = 0)\n",
    "    mito_counts.head()\n",
    "    mito_counts.to_csv(f'{mito_dir}mito_fragments_counts.csv')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cde3001e",
   "metadata": {},
   "source": [
    "## Split frags greater than 3000 bp\n",
    "\n",
    "From other analysis, we know that frags greater than 3000 bp mostly span the deleted region and are therefore representative of endjoined mito. Thus, we split the fragments into two with a padding of 50. \n",
    "\n",
    "For example,\n",
    "\n",
    "Fragment `7692 - 11869` --> `7692 - 7742` and `11819 - 11869`\n",
    "\n",
    "If you do not need to remove any mitochondrial fragments, can SKIP this step."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "id": "f9dd76cb",
   "metadata": {},
   "outputs": [],
   "source": [
    "def split_long_frag_df(df, frag_threshold, padding):\n",
    "    assert padding >= 0, \"padding must be positive\"\n",
    "    \n",
    "    long_frags_df = df[(df['end'] - df['start']) >= frag_threshold]\n",
    "    start_frags_df =  long_frags_df.copy()\n",
    "    start_frags_df['end'] = start_frags_df['start'] + padding\n",
    "    end_frags_df =  long_frags_df.copy()\n",
    "    end_frags_df['start'] = end_frags_df['end'] - padding\n",
    "    \n",
    "    \n",
    "    return pd.concat([start_frags_df, end_frags_df]).sort_values('start')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "41508017",
   "metadata": {},
   "outputs": [],
   "source": [
    "frag_threshold = 3000\n",
    "padding = 50\n",
    "\n",
    "for s in samples:\n",
    "    mito_dir = f\"{cache_dir}{s}/mito/\"\n",
    "    \n",
    "    # READ IN INPUT (output of above code)\n",
    "    unfilt_path = f'{mito_dir}mito_fragments.csv'\n",
    "    \n",
    "    # OUTPUT path\n",
    "    filt_path = f\"{mito_dir}mito_fragments_split_{frag_threshold}_bp_frags.csv\"\n",
    "    \n",
    "    df = pd.read_csv(unfilt_path, index_col = 0)\n",
    "    no_long_frags_df = df[(df['end'] - df['start']) < frag_threshold]\n",
    "    split_frags_df = split_long_frag_df(df, frag_threshold, padding)\n",
    "    new_frags_df = pd.concat([no_long_frags_df, split_frags_df]).sort_values('start')\n",
    "    \n",
    "    new_frags_df.to_csv(filt_path)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7517c840",
   "metadata": {},
   "source": [
    "# 02 Generate Bulk Mitochondrial Coverage Plots\n",
    "\n",
    "See `8HetClonesHashtagMix_process_mito_fragments.ipynb`, same code just adapt inputs and outputs. Since not essential and only for visualization purposes, we do not repeat in this notebook. "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6d7c874b",
   "metadata": {},
   "source": [
    "# 03 Compute Single-cell coverage (takes longer to run)\n",
    "\n",
    "For each cell, we report the number of unique fragments (different START and END) that overlap each base pair."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c6629916",
   "metadata": {},
   "outputs": [],
   "source": [
    "# relative to repo\n",
    "cov_script = \"../../scripts/single-cell-coverage.py\"\n",
    "mito_genes_path = '../../misc/mito_genes.csv'\n",
    "\n",
    "samples = [\n",
    "    \"YF-2721_48mix_25dGlu_3d_mito_hLLLdogma_multiome\",\n",
    "    \"YF-2721_48mix_10dGal_3d_mito_hLLLdogma_multiome\",\n",
    "    \"YF-2721_48mix_25dGlu_NoPyr_3d_mito_hLLLdogma_multiome\",\n",
    "]\n",
    "\n",
    "for s in samples:\n",
    "    mito_frags_path = f'{cache_dir}{s}/mito/mito_fragments_split_3000_bp_frags.csv'\n",
    "    barcodes_path = f\"{cache_dir}{s}/mito/barcodes.tsv\"\n",
    "    out_dir = f'{cache_dir}{s}/mito/single-cell-coverage_split_3000_bp_frags/'\n",
    "    if not os.path.exists(out_dir):\n",
    "        os.mkdir(out_dir)\n",
    "    \n",
    "    memory = 10\n",
    "    conda_env = 'scanpy'\n",
    "    options = dict()\n",
    "    options['J'] = f'{s}_coverage'\n",
    "    options['n'] = f'5'\n",
    "    options['R'] = f'rusage[mem={memory}]'\n",
    "    options['W'] = f'120:00'\n",
    "    options['o'] = f'{s}.stdout'\n",
    "    options['eo'] = f'{s}.stdout'\n",
    "\n",
    "\n",
    "    if os.path.exists(options['o']):\n",
    "        os.remove(options['o'])\n",
    "\n",
    "    option_string = ' '.join([f'-{flag} {val}' for flag, val in options.items()])\n",
    "    job_string = f\"source ~/.bashrc; conda activate {conda_env};\"\n",
    "    \n",
    "    # THIS IS THE ACTUAL COMMAND LINE TO RUN THE SCRIPT (everything else is to allow it to run on MSKCC servers)\n",
    "    run_script_string = f\"python {cov_script} --mito_genes_path {mito_genes_path} --mito_frags_path {mito_frags_path} --barcodes_path {barcodes_path} --out_dir {out_dir}\"\n",
    "\n",
    "    command = f'bsub {option_string} \\\"{job_string} {run_script_string}\"'\n",
    "    os.system(command)\n",
    "    print(command + '\\n')"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "scanpy",
   "language": "python",
   "name": "scanpy"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.17"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}