{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "afbe84db",
   "metadata": {},
   "source": [
    "# Preprocessing"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f0892ec2",
   "metadata": {},
   "source": [
    "## Import pakcages"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "b7edb9aa",
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "import sys\n",
    "import copy\n",
    "import numpy as np"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "496bc92d",
   "metadata": {},
   "outputs": [],
   "source": [
    "import ConsTADs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "d203a565",
   "metadata": {},
   "outputs": [],
   "source": [
    "from ConsTADs import source\n",
    "from ConsTADs import Preprocess as pre\n",
    "from ConsTADs import TadSeparationLandscape as TSL\n",
    "from ConsTADs import IdentifyBoundaryRegions as IBR\n",
    "from ConsTADs import GetConsTADs as GCT\n",
    "from ConsTADs import ConsTADs_main as main\n",
    "from ConsTADs import PlotFunction as PF"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3600a815",
   "metadata": {},
   "source": [
    "# Build TAD separation landscape and get ConsTADs"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "dcb5f466",
   "metadata": {},
   "source": [
    "## paramaters setting"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "e0702214",
   "metadata": {},
   "outputs": [],
   "source": [
    "sample = 'GM12878_MboI'\n",
    "#sample = 'K562_MboI'\n",
    "# Note: the TAD_caller_result_add and mat_file should use your own results\n",
    "# Step 1,2 paramaters\n",
    "TAD_caller_result_add = '/home/dcdang/ConsTADs/TAD_caller_result/' + sample + '/chr2_method_result'\n",
    "target_chr = 'chr2'\n",
    "resolution = 50000\n",
    "chr_size = 243199373\n",
    "method_list = ['3DNetMod', 'CaTCH', 'CHDF', 'ClusterTAD', \n",
    "               'deDoc', 'DI', 'GMAP', 'HiCDB', 'HiCseg', \n",
    "               'HiTAD', 'ICFinder', 'IS', 'MSTD', 'OnTAD',\n",
    "               'Spectral','TopDom']\n",
    "\n",
    "# Step 3 paramaters\n",
    "mat_file = '/mnt/disk2/ddc/project/ConsTADs/hic_mat/' + 'Rao2014_' + sample + '_50kb_chr2_iced_sparse.txt'\n",
    "\n",
    "mat_type = 'sparse'\n",
    "mat_norm_check = False\n",
    "window_list = [100000, 150000, 200000, 300000, 400000, 500000, 600000, 800000, 1000000, 2000000] \n",
    "p_cut = 0.05\n",
    "high_score_cut = 5\n",
    "combine_dist = 2\n",
    "\n",
    "# Step 4 paramaters\n",
    "K = 3\n",
    "color_bd = ['#D65F4D', '#459457', '#4392C3']\n",
    "\n",
    "# Step 5 paramaters\n",
    "weight = 0.5"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5165b070",
   "metadata": {},
   "source": [
    "## Main process"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "380e0b48",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Help on function ConsTADs in module ConsTADs.ConsTADs_main:\n",
      "\n",
      "ConsTADs(TAD_caller_result_add, target_chr, resolution, chr_size, method_list, mat_file, mat_type, window_list, mat_norm_check=False, p_cut=0.05, high_score_cut=5, combine_dist=2, K=3, color_bd=['#D65F4D', '#459457', '#4392C3'], weight=0.5)\n",
      "    Parameters\n",
      "    ----------\n",
      "    TAD_caller_result_add : str\n",
      "        Path of fold containing results of multiple TAD callers.\n",
      "    target_chr : str\n",
      "        Symbol of target chromosome, eg: chr2.\n",
      "    resolution : int\n",
      "        Resolution of Hi-C contact map.\n",
      "    chr_size : int\n",
      "        Length of target chromosome in base pair.\n",
      "    method_list : list\n",
      "        List of TAD caller names.\n",
      "    mat_file : str\n",
      "        Path of file for Hi-C contact map of target chr.\n",
      "    mat_type : str\n",
      "        Type of Hi-C contact map file, dense or sparse.\n",
      "    window_list : list\n",
      "        List of window size for calculating multi-scale contrast p-value. These sizes are lengths in base pair. \n",
      "    mat_norm_check : bool, optional\n",
      "        Dose the Hi-C contact map are distance-dependent normalized. The default is False. \n",
      "    p_cut : float, optional\n",
      "        Cut off of contrast p-value. The default is 0.05.\n",
      "    high_score_cut : int, optional\n",
      "        Cut off of confident boundary score. The default is 5.\n",
      "    combine_dist : int, optional\n",
      "        The cut off of distance for boundary regions' combination. The default is 2.\n",
      "    K : int, optional\n",
      "        Number of boundary region types. The default is 3.\n",
      "    color_bd : list, optional\n",
      "        List of color for different boundary region types. The default is ['#D65F4D', '#459457', '#4392C3'].\n",
      "    weight : float, optional\n",
      "        weight to calculate best bin pair. \n",
      "        Balance affect of boundary score and domain value. The default is 0.5 for boundary score.\n",
      "    \n",
      "    Returns\n",
      "    -------\n",
      "    TAD_caller_result_all : dict\n",
      "        Dictionary containing the results of multiple TAD callers for the Hi-C contact map of target chromosme.\n",
      "    mat_dense : numpy.array\n",
      "        Hi-C contact map in dense format.\n",
      "    mat_norm : numpy.array\n",
      "        Distance-dependent normalized Hi-C contact map in dense format.\n",
      "    result_record : dict\n",
      "        Result of each oreration when building the TAD separation landscape.\n",
      "    w_best : int\n",
      "        Best window size for contrast p-value.\n",
      "    df_bd_insul_pvalue : pandas.DataFrame\n",
      "        Multi-scale contrast p-value for each bin along the chromosome.\n",
      "    df_pvalue_score_cor : pandas.DataFrame\n",
      "        Pearson correlation of multi-scale contrast p-value and boundary score profile.\n",
      "    df_boundary_region_with_types : pandas.DataFrame\n",
      "        DataFrame of Boundary regions and their types.\n",
      "    df_tad_cons : pandas.DataFrame\n",
      "        DataFrame of ConsTADs.\n",
      "    df_boundary_cons : pandas.DataFrame\n",
      "        DataFrame of ConsTADs boundaries.\n",
      "\n"
     ]
    }
   ],
   "source": [
    "## Information of the input and output of the main process of ConsTADs\n",
    "help(main.ConsTADs)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "4c4818ae",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\u001b[1mStep 1: Preprocess of TAD caller results...\u001b[0m\n",
      "Done!\n",
      "\u001b[1mStep 2: Boundary voting...\u001b[0m\n",
      "Done!\n",
      "\u001b[1mStep 3: Build TAD separation landscape...\u001b[0m\n",
      "Load Hi-C matrix...\n",
      "Normalize the Hi-C matrix...\n",
      "Calculate multi-scale contrast pvalue...\n",
      "Best window size:200kb-window\n",
      "Building TAD seperation landscape...\n",
      "Operation 1: Add\n",
      "Add score for 69 bins\n",
      "Operation 2: Filter\n",
      "Operation 3: Combine\n",
      "Done!\n",
      "\u001b[1mStep 4: Identify three types of boundary regions...\u001b[0m\n",
      "Length cut-off:5.0\n",
      "Score cut-off:5.0\n",
      "Boundary region types: Narrow-strong, Number: 239\n",
      "Boundary region types: Narrow-weak, Number: 269\n",
      "Boundary region types: Wide, Number: 93\n"
     ]
    },
    {
     "data": {
      "image/png": "[encoded data, 22944 characters]\n",
      "text/plain": [
       "<Figure size 360x302.4 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Done!\n",
      "\u001b[1mStep 5: Get ConsTADs based on TAD separation landscape...\u001b[0m\n",
      "Get 600 ConsTADs\n",
      "Done!\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "<Figure size 432x288 with 0 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "TAD_caller_result_all, mat_dense, mat_norm, result_record, w_best, df_bd_insul_pvalue, df_pvalue_score_cor, df_boundary_region_with_types, df_tad_cons, df_boundary_cons = main.ConsTADs(TAD_caller_result_add, \n",
    "             target_chr, resolution, chr_size, method_list, mat_file, mat_type, window_list, mat_norm_check = mat_norm_check, \n",
    "             p_cut = p_cut, high_score_cut = high_score_cut, combine_dist = combine_dist, K = K, color_bd = color_bd, weight = weight)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ecd7f5f9",
   "metadata": {},
   "source": [
    "## Results saving"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "ee3d9a4d",
   "metadata": {},
   "outputs": [],
   "source": [
    "save_add = '/mnt/disk2/ddc/project/ConsTADs/result_save_jupyter'\n",
    "save_file = save_add + '/' + sample + '_' + target_chr + '_' + str(int(resolution / 1000)) + 'kb' + '_ConsTADs_results_all.pkl'\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "196d59e4",
   "metadata": {},
   "outputs": [],
   "source": [
    "sample_res_all = {}\n",
    "sample_res_all['TAD_caller_results'] = TAD_caller_result_all\n",
    "sample_res_all['Hi-C_map'] = mat_dense\n",
    "sample_res_all['TSL'] = result_record\n",
    "sample_res_all['Contrast_pvalue'] = df_bd_insul_pvalue\n",
    "sample_res_all['boundary_region'] = df_boundary_region_with_types\n",
    "sample_res_all['ConsTADs'] = df_tad_cons\n",
    "sample_res_all['ConsTADs_bd'] = df_boundary_cons"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "id": "59cbd00b",
   "metadata": {},
   "outputs": [],
   "source": [
    "source.SaveData(save_file, sample_res_all)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "74c17dcd",
   "metadata": {},
   "source": [
    "# Results display"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "517a1875",
   "metadata": {},
   "source": [
    "## TAD separation landscape and multi-scale contrast P-value"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "80df313f",
   "metadata": {},
   "outputs": [],
   "source": [
    "Chr = target_chr\n",
    "bd_cell_score_original = result_record['Original']['TAD_score']\n",
    "bd_cell_score_refined = result_record['Combine']['TAD_score']\n",
    "df_bd_pvalue_result = df_bd_insul_pvalue\n",
    "w_list = ['200kb-window', '400kb-window', '800kb-window', '1000kb-window']\n",
    "p_cut = 0.05\n",
    "save_name = ''\n",
    "st = 750\n",
    "ed = 820\n",
    "fgsize = (5, 10)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "c1aa26fe",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 73664 characters]\n",
      "text/plain": [
       "<Figure size 360x720 with 9 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Note: The best window size for contrast P-value is colored in red. \n",
    "# Note: The TAD separation landscape (TSL) is colored in red. \n",
    "PF.DrawBdScoreMultiScalePvalue(mat_dense, Chr, st, ed, bd_cell_score_original, bd_cell_score_refined, df_bd_pvalue_result, \n",
    "                               w_list, resolution, save_name, p_cut, target_site = [], fgsize = fgsize, bin_size = 10)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "02639339",
   "metadata": {},
   "source": [
    "## ConsTADs with annotated boundary region types"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "b005b7b8",
   "metadata": {},
   "outputs": [],
   "source": [
    "Chr = target_chr\n",
    "bd_score_cell = result_record['Combine']['TAD_score']\n",
    "symbol_dic = {'No-bd':0, 'Narrow-weak':1, 'Narrow-strong':2, 'Wide':3}\n",
    "bd_symbol = PF.get_bd_type_symbol(df_boundary_region_with_types, bd_score_cell, symbol_dic)\n",
    "fgsize = (6, 7)\n",
    "st = 750\n",
    "ed = 820\n",
    "TAD_list = []\n",
    "for i in range(len(df_tad_cons)):\n",
    "    tad_st = df_tad_cons['bd_st'][i]\n",
    "    tad_ed = df_tad_cons['bd_ed'][i]\n",
    "    if tad_st >= st and tad_ed <= ed:\n",
    "        TAD_list.append((tad_st, tad_ed))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "ccb75a37",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 39600 characters]\n",
      "text/plain": [
       "<Figure size 432x504 with 5 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Note: The black frames represent the ConsTADs and the boundary regions types are shown below. \n",
    "#       Wide: blue, Narrow-strong: red, Narrow-weak: green\n",
    "PF.DrawConsTADsAndBdRegion(mat_dense, Chr, st, ed, TAD_list, bd_score_cell, bd_symbol, resolution, \n",
    "                           fgsize = fgsize, save_name = '', bin_size = 10)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fa6970b5",
   "metadata": {},
   "source": [
    "## Comparison of TAD separation landscape between cell types"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "e7daa276",
   "metadata": {},
   "outputs": [],
   "source": [
    "save_add = '/mnt/disk2/ddc/project/ConsTADs/result_save_jupyter'\n",
    "Chr = target_chr = 'chr2'\n",
    "resolution = 50000\n",
    "sample1 = 'GM12878_MboI'\n",
    "sample_res_all_1 = source.ReadSaveData(save_add + '/' + sample1 + '_' + target_chr + '_' + str(int(resolution / 1000)) + 'kb' + '_ConsTADs_results_all.pkl')\n",
    "sample2 = 'K562_MboI'\n",
    "sample_res_all_2 = source.ReadSaveData(save_add + '/' + sample2 + '_' + target_chr + '_' + str(int(resolution / 1000)) + 'kb' + '_ConsTADs_results_all.pkl')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "6c1776b1",
   "metadata": {},
   "outputs": [],
   "source": [
    "cell_1 = 'GM12878'\n",
    "bd_score_cell_1 = sample_res_all_1['TSL']['Combine']['TAD_score']\n",
    "df_tad_cons_1 = sample_res_all_1['ConsTADs']\n",
    "mat_dense1 = sample_res_all_1['Hi-C_map']\n",
    "\n",
    "cell_2 = 'K562'\n",
    "bd_score_cell_2 = sample_res_all_2['TSL']['Combine']['TAD_score']\n",
    "df_tad_cons_2 = sample_res_all_2['ConsTADs']\n",
    "mat_dense2 = sample_res_all_2['Hi-C_map']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "a5264f89",
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 44368 characters]\n",
      "text/plain": [
       "<Figure size 432x576 with 5 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "## Region around a GM12878-gained boundary region (K562-lost)\n",
    "st = 3646 - 30\n",
    "ed = 3646 + 30\n",
    "figsize = (6, 8)\n",
    "TAD_list_1 = PF.get_tad_list_in_target_ranges(st, ed, df_tad_cons_1, pos_type = 'cord', resolution = resolution)\n",
    "TAD_list_2 = PF.get_tad_list_in_target_ranges(st, ed, df_tad_cons_2, pos_type = 'cord', resolution = resolution)\n",
    "PF.DrawPairWiseMapCompare(cell_1, cell_2, mat_dense1, mat_dense2, Chr, st, ed, bd_score_cell_1, bd_score_cell_2, \n",
    "                          TAD_list_1, TAD_list_2, resolution, figsize = figsize, TAD_color_1 = 'black', TAD_color_2 = 'white', \n",
    "                          target_site = [], save_name = '', bin_size = 10) \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "5ce8cebe",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 47760 characters]\n",
      "text/plain": [
       "<Figure size 432x576 with 5 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "## Region around a GM12878-lost boundary region (K562-gained)\n",
    "st = 565 - 30\n",
    "ed = 565 + 30 \n",
    "\n",
    "figsize = (6, 8)\n",
    "TAD_list_1 = PF.get_tad_list_in_target_ranges(st, ed, df_tad_cons_1, pos_type = 'cord', resolution = resolution)\n",
    "TAD_list_2 = PF.get_tad_list_in_target_ranges(st, ed, df_tad_cons_2, pos_type = 'cord', resolution = resolution)\n",
    "PF.DrawPairWiseMapCompare(cell_1, cell_2, mat_dense1, mat_dense2, Chr, st, ed, bd_score_cell_1, bd_score_cell_2, \n",
    "                          TAD_list_1, TAD_list_2, resolution, figsize = figsize, TAD_color_1 = 'black', TAD_color_2 = 'white', \n",
    "                          target_site = [], save_name = '', bin_size = 10) \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2fd103da",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "567ab10f",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "ConsTADs",
   "language": "python",
   "name": "constads"
  },
  "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.7"
  },
  "toc": {
   "base_numbering": 1,
   "nav_menu": {},
   "number_sections": true,
   "sideBar": true,
   "skip_h1_title": false,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": false,
   "toc_position": {
    "height": "calc(100% - 180px)",
    "left": "10px",
    "top": "150px",
    "width": "303.837px"
   },
   "toc_section_display": true,
   "toc_window_display": true
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
