2022_sampling_pairs.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Sampling to equal read depth from pairs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Import the packages we will use\n",
    "%matplotlib inline\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import h5py\n",
    "import bioframe\n",
    "from matplotlib.gridspec import GridSpec\n",
    "from matplotlib.gridspec import GridSpecFromSubplotSpec\n",
    "import matplotlib.colors as colors\n",
    "import cooler\n",
    "import seaborn as sns\n",
    "from matplotlib.colors import ListedColormap\n",
    "import os\n",
    "import glob\n",
    "import pairtools\n",
    "import random"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {},
   "outputs": [],
   "source": [
    "#For many libraries\n",
    "inDir = \"/nl/umw_job_dekker/users/as38w/JZ2022_git/pairs/\"\n",
    "outDir = \"/nl/umw_job_dekker/users/as38w/JZ2022_git/coolers/\"\n",
    "\n",
    "conditions = [\n",
    "    'DMSO-ND',\n",
    "    'DMSO-4OHT'\n",
    "        \n",
    "]\n",
    "\n",
    "long_names = {\n",
    "    'DMSO-ND':'JZ-MEF-AsiS1-ND-dmso6h-HiCD2R1',\n",
    "    'DMSO-4OHT':'JZ-MEF-AsiS1-4OHT-dmso6h-HiCD2R1'\n",
    "}\n",
    "\n",
    "pairsPath = {}\n",
    "for cond in conditions[0:2]:\n",
    "    pairsPath[cond] = '{}/{}__mm10.mm10.nodups.pairs.gz'.format(inDir, long_names[cond])\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{'DMSO-ND': '/nl/umw_job_dekker/users/as38w/JZ2022_git/pairs//JZ-MEF-AsiS1-ND-dmso6h-HiCD2R1__mm10.mm10.nodups.pairs.gz',\n",
       " 'DMSO-4OHT': '/nl/umw_job_dekker/users/as38w/JZ2022_git/pairs//JZ-MEF-AsiS1-4OHT-dmso6h-HiCD2R1__mm10.mm10.nodups.pairs.gz'}"
      ]
     },
     "execution_count": 37,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "pairsPath"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "metadata": {},
   "outputs": [],
   "source": [
    "noheaderPath = {}\n",
    "for cond in conditions:\n",
    "    noheaderPath[cond] = '{}/{}.mm10.nodups.noheader.pairs'.format(outDir, long_names[cond])\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "INFO: Total memory requested is 48000 MB (4 cores x 12000 MB)\n",
      "Job <1543594> is submitted to queue <short>.\n",
      "INFO: Total memory requested is 48000 MB (4 cores x 12000 MB)\n",
      "Job <1543595> is submitted to queue <short>.\n"
     ]
    }
   ],
   "source": [
    "#Remove header and unzip noheader file\n",
    "for cond in conditions:\n",
    "    infile = pairsPath[cond]\n",
    "    noheaderfile = noheaderPath[cond]\n",
    "    !bsub -q short -W 04:00 -n 4 -R \"span[hosts=1]\" -R rusage[mem=12000] -R \"select[rh=8]\"\\\n",
    "        \"zgrep -o '^[^#]*' $infile > $noheaderfile\"\n",
    "        "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "metadata": {},
   "outputs": [],
   "source": [
    "sampledPaths = {}\n",
    "\n",
    "for cond in conditions:\n",
    "    sampledPaths[cond] = '{}/{}.mm10.nodups.noheader_sampled_64509240M.pairs'.format(outDir, long_names[cond])\n",
    " "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "INFO: Total memory requested is 96000 MB (4 cores x 24000 MB)\n",
      "Job <1543596> is submitted to queue <short>.\n",
      "INFO: Total memory requested is 96000 MB (4 cores x 24000 MB)\n",
      "Job <1543597> is submitted to queue <short>.\n"
     ]
    }
   ],
   "source": [
    "#shuffle\n",
    "for cond in conditions:\n",
    "    noheaderfile = noheaderPath[cond]\n",
    "    shuffile = sampledPaths[cond]\n",
    "    !bsub -q short -W 04:00 -n 4 -R \"span[hosts=1]\" -R rusage[mem=24000] -R \"select[rh=8]\"\\\n",
    "        \"shuf -n 64509240 -o $shuffile $noheaderfile\"\n",
    "        "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 47,
   "metadata": {},
   "outputs": [],
   "source": [
    "#add header before sorting\n",
    "headerPaths = {}\n",
    "\n",
    "for cond in conditions:\n",
    "    headerPaths[cond] = '{}/{}.mm10.nodups.header_sampled_65M.pairs'.format(outDir, long_names[cond])\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 48,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "INFO: Total memory requested is 48000 MB (4 cores x 12000 MB)\n",
      "Job <1543602> is submitted to queue <short>.\n",
      "INFO: Total memory requested is 48000 MB (4 cores x 12000 MB)\n",
      "Job <1543603> is submitted to queue <short>.\n"
     ]
    }
   ],
   "source": [
    "#add header to sampled files ... #IMPORTANT: need the chromsizes (1-22, X, Y, M) and columns line.  I made this into one file:\n",
    "headerfile = '/nl/umw_job_dekker/users/as38w/JZ-MEF-2020/sampled_2022/pairsheader_mm10'\n",
    "for cond in conditions:\n",
    "    in_fname = sampledPaths[cond]\n",
    "    tmp_fname = headerPaths[cond]\n",
    "    !bsub -q short -W 01:00 -n 4 -R \"span[hosts=1]\" -R rusage[mem=12000] -R \"select[rh=8]\"\\\n",
    "        \"cat $headerfile $in_fname > $tmp_fname\"\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 49,
   "metadata": {},
   "outputs": [],
   "source": [
    "sortedPaths = {}\n",
    "\n",
    "for cond in conditions:\n",
    "    sortedPaths[cond] = '{}/{}.mm10.nodups.sorted.sampled_65M.pairs'.format(outDir, long_names[cond])\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 50,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Job <1543604> is submitted to queue <short>.\n",
      "Job <1543605> is submitted to queue <short>.\n"
     ]
    }
   ],
   "source": [
    "#sort sampled files plus header\n",
    "for cond in conditions:\n",
    "    in_fname = headerPaths[cond]\n",
    "    out_fname = sortedPaths[cond]\n",
    "    !bsub -q short -W 04:00 -n 2 -R span[hosts=1] -R select[ib] -R \"select[rh=8]\" -R rusage[mem=12000] \"source activate coolerenv_manual_master; pairtools sort -o $out_fname $in_fname\"\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 51,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "INFO: Total memory requested is 96000 MB (4 cores x 24000 MB)\n",
      "Job <1543606> is submitted to queue <short>.\n",
      "INFO: Total memory requested is 96000 MB (4 cores x 24000 MB)\n",
      "Job <1543607> is submitted to queue <short>.\n"
     ]
    }
   ],
   "source": [
    "#rezip using bgzip on sampled files\n",
    "for cond in conditions:\n",
    "    in_fname = sortedPaths[cond]\n",
    "    !bsub -q short -W 04:00 -n 4 -R \"span[hosts=1]\" -R rusage[mem=24000] -R \"select[rh=8]\"\\\n",
    "        \"bgzip $in_fname\"\n",
    "        "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 53,
   "metadata": {},
   "outputs": [],
   "source": [
    "#For many libraries\n",
    "inDir = \"/nl/umw_job_dekker/users/as38w/JZ2022_git/pairs/\"\n",
    "outDir = \"/nl/umw_job_dekker/users/as38w/JZ2022_git/coolers/\"\n",
    "\n",
    "newconditions = [\n",
    "    'DMSO-ND',\n",
    "    'DMSO-4OHT'\n",
    "        \n",
    "]\n",
    "\n",
    "newlong_names = {\n",
    "    'DMSO-ND':'JZ-MEF-AsiS1-ND-dmso6h-HiCD2R1',\n",
    "    'DMSO-4OHT':'JZ-MEF-AsiS1-4OHT-dmso6h-HiCD2R1'\n",
    "}\n",
    "\n",
    "newPairsPath = {}\n",
    "for cond in newconditions:\n",
    "    newPairsPath[cond] = '{}/{}.mm10.nodups.sorted.sampled_65M.pairs.gz'.format(outDir, newlong_names[cond])\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 54,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "INFO: Total memory requested is 96000 MB (4 cores x 24000 MB)\n",
      "Job <1543608> is submitted to queue <short>.\n",
      "INFO: Total memory requested is 96000 MB (4 cores x 24000 MB)\n",
      "Job <1543609> is submitted to queue <short>.\n"
     ]
    }
   ],
   "source": [
    "#index the pairs files after sorting\n",
    "for cond in newconditions:\n",
    "    in_fname = newPairsPath[cond]\n",
    "    !bsub -q short -W 04:00 -n 4 -R \"span[hosts=1]\" -R rusage[mem=24000] -R \"select[rh=8]\"\\\n",
    "        \"pairix -p pairs $in_fname\"\n",
    "        "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 55,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "INFO: Total memory requested is 48000 MB (4 cores x 12000 MB)\n",
      "Job <1543610> is submitted to queue <short>.\n",
      "INFO: Total memory requested is 48000 MB (4 cores x 12000 MB)\n",
      "Job <1543611> is submitted to queue <short>.\n"
     ]
    }
   ],
   "source": [
    "#stats\n",
    "for cond in newconditions:\n",
    "    in_fname = newPairsPath[cond]\n",
    "    out_fname = '{}/{}.mm10.nodups.sorted.sampled_65M.pairs.stats.tsv'.format(outDir, newlong_names[cond])\n",
    "    !bsub -q short -W 03:00 -n 4 -R \"span[hosts=1]\" -R rusage[mem=12000] -R \"select[rh=8]\"\\\n",
    "        \"pairtools stats -o $out_fname $in_fname\"\n",
    "       "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 56,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "INFO: Total memory requested is 48000 MB (4 cores x 12000 MB)\n",
      "Job <1543612> is submitted to queue <short>.\n"
     ]
    }
   ],
   "source": [
    "chromsizes=\"/nl/umw_job_dekker/cshare/reference/sorted_chromsizes/mm10.reduced.chrom.sizes\"\n",
    "binfile=\"/nl/umw_job_dekker/users/as38w/JZ2022_git/coolers/mm10_1000_bins\"\n",
    "\n",
    "!bsub -q short -W 03:00 -n 4 -R \"span[hosts=1]\" -R \"select[rh=8]\" -R rusage[mem=12000] \"source activate coolerenv_manual_master; cooler makebins -o $binfile $chromsizes 1000\"\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 57,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "INFO: Total memory requested is 48000 MB (4 cores x 12000 MB)\n",
      "Job <1543614> is submitted to queue <short>.\n",
      "INFO: Total memory requested is 48000 MB (4 cores x 12000 MB)\n",
      "Job <1543615> is submitted to queue <short>.\n"
     ]
    }
   ],
   "source": [
    "#make 1kb sampled coolers\n",
    "binfile = '/nl/umw_job_dekker/users/as38w/JZ2022_git/coolers/mm10_1000_bins'\n",
    "\n",
    "for cond in newconditions:\n",
    "    in_fname = newPairsPath[cond]\n",
    "    coolfile = '{}/{}.mm10.nodups.65M.1000.cool'.format(outDir, newlong_names[cond])\n",
    "    !bsub -q short -W 03:00 -n 4 -R \"span[hosts=1]\" -R \"select[rh=8]\" -R rusage[mem=12000] \"source activate coolerenv_manual_master; cooler cload pairs --assembly mm10 -c1 2 -p1 3 -c2 4 -p2 5 $binfile $in_fname $coolfile\"\n",
    "        "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 59,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "INFO: Total memory requested is 48000 MB (4 cores x 12000 MB)\n",
      "Job <1543641> is submitted to queue <short>.\n",
      "INFO: Total memory requested is 48000 MB (4 cores x 12000 MB)\n",
      "Job <1543642> is submitted to queue <short>.\n"
     ]
    }
   ],
   "source": [
    "#balance single resolution coolers\n",
    "for cond in newconditions:\n",
    "    coolfile = '{}/{}.mm10.nodups.65M.1000.cool'.format(outDir, newlong_names[cond])\n",
    "    !bsub -q short -W 03:00 -n 4 -R \"span[hosts=1]\" -R \"select[rh=8]\" -R rusage[mem=12000] \"source activate coolerenv_manual_master; cooler balance -p 8 $coolfile\"\n",
    "                            "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 60,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "INFO: Total memory requested is 48000 MB (4 cores x 12000 MB)\n",
      "Job <1543643> is submitted to queue <short>.\n",
      "INFO: Total memory requested is 48000 MB (4 cores x 12000 MB)\n",
      "Job <1543644> is submitted to queue <short>.\n"
     ]
    }
   ],
   "source": [
    "#zoomify and balance 1kb coolers into mcools\n",
    "res_str = '10000000,5000000,2500000,1000000,500000,250000,100000,50000,25000,10000,5000,2000,1000'\n",
    "\n",
    "#sampled files\n",
    "for cond in newconditions:\n",
    "    coolfile = '{}/{}.mm10.nodups.65M.1000.cool'.format(outDir, newlong_names[cond])\n",
    "    mcoolfile = '{}/{}.mm10.nodups.65M.1000.mcool'.format(outDir, newlong_names[cond])\n",
    "    !bsub -q short -W 03:00 -n 4 -R \"span[hosts=1]\" -R \"select[rh=8]\" -R rusage[mem=12000] \"source activate coolerenv_updated_master; cooler zoomify -p 8 --balance --resolutions $res_str $coolfile -o $mcoolfile\"\n"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "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"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}