rewrite_info_cLAD.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 102,
   "metadata": {},
   "outputs": [],
   "source": [
    "# rewrite_info_cLAD.ipy\n",
    "\n",
    "import os\n",
    "import sys\n",
    "import math\n",
    "import scipy\n",
    "import random\n",
    "import warnings\n",
    "import scipy.stats\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import seaborn as sns\n",
    "import matplotlib as mpl\n",
    "from scipy import optimize\n",
    "from scipy.stats import norm\n",
    "from sklearn import datasets\n",
    "from scipy import stats as st\n",
    "import matplotlib.pyplot as plt\n",
    "from sklearn.metrics import pairwise_distances\n",
    "from sklearn.preprocessing import StandardScaler\n",
    "from collections import defaultdict\n",
    "\n",
    "\n",
    "warnings.filterwarnings('ignore')\n",
    "# np.set_printoptions(threshold='nan')\n",
    "np.set_printoptions(threshold=sys.maxsize)\n",
    "\n",
    "%matplotlib inline\n",
    "# Load data and select first column\n",
    "\n",
    "sns.set_style(\"white\")\n",
    "sns.set_context(\"talk\")\n",
    "resolution = 10000\n",
    "chr_name_list = [\"chr\"+str(i) for i in range(1,20)]\n",
    "given_resolution = 1e+4\n",
    "# cLAD_dropout = 0.01\n",
    "# cLAD_dropout = 0.05\n",
    "cLAD_dropout = 0.10\n",
    "# cLAD_dropout = 0.15\n",
    "# cLAD_dropout = 0.20\n",
    "\n",
    "cell_name_all=[17,18,19,20,25,26,27,28,29,38,39]\n",
    "cell_name_index = 29\n",
    "seed_sample = cell_name_index #cell index\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 103,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "dropout rate 0.1, num of cLAD, before dropout 3903, after dropout 3512\n"
     ]
    },
    {
     "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>chr</th>\n",
       "      <th>loci_head</th>\n",
       "      <th>loci_tail</th>\n",
       "      <th>3</th>\n",
       "      <th>4</th>\n",
       "      <th>5</th>\n",
       "      <th>head_resolution</th>\n",
       "      <th>tail_resolution</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>chr1</td>\n",
       "      <td>3010218</td>\n",
       "      <td>3087058</td>\n",
       "      <td>.</td>\n",
       "      <td>1</td>\n",
       "      <td>.</td>\n",
       "      <td>3010000</td>\n",
       "      <td>3080000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>chr1</td>\n",
       "      <td>3126730</td>\n",
       "      <td>3196670</td>\n",
       "      <td>.</td>\n",
       "      <td>1</td>\n",
       "      <td>.</td>\n",
       "      <td>3120000</td>\n",
       "      <td>3190000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>chr1</td>\n",
       "      <td>3220610</td>\n",
       "      <td>3282860</td>\n",
       "      <td>.</td>\n",
       "      <td>1</td>\n",
       "      <td>.</td>\n",
       "      <td>3220000</td>\n",
       "      <td>3280000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>chr1</td>\n",
       "      <td>3339298</td>\n",
       "      <td>3357024</td>\n",
       "      <td>.</td>\n",
       "      <td>1</td>\n",
       "      <td>.</td>\n",
       "      <td>3330000</td>\n",
       "      <td>3350000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>chr1</td>\n",
       "      <td>3371740</td>\n",
       "      <td>3389850</td>\n",
       "      <td>.</td>\n",
       "      <td>1</td>\n",
       "      <td>.</td>\n",
       "      <td>3370000</td>\n",
       "      <td>3380000</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "    chr  loci_head  loci_tail  3  4  5  head_resolution  tail_resolution\n",
       "0  chr1    3010218    3087058  .  1  .          3010000          3080000\n",
       "1  chr1    3126730    3196670  .  1  .          3120000          3190000\n",
       "2  chr1    3220610    3282860  .  1  .          3220000          3280000\n",
       "3  chr1    3339298    3357024  .  1  .          3330000          3350000\n",
       "4  chr1    3371740    3389850  .  1  .          3370000          3380000"
      ]
     },
     "execution_count": 103,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "LAD_dir = \"/home/ygli/gam_paper_data/lad_data/GSE17051_cLAD_regions_mm10.bed\"\n",
    "LAD_data = pd.read_csv(LAD_dir,sep='\\t',header=None)\n",
    "LAD_data.columns=['chr','loci_head','loci_tail','3','4','5']\n",
    "LAD_data['head_resolution'] = [ int(i.loci_head//given_resolution*given_resolution) for _,i in LAD_data.iterrows()]\n",
    "LAD_data['tail_resolution'] = [ int(i.loci_tail//given_resolution*given_resolution) for _,i in LAD_data.iterrows()]\n",
    "\n",
    "#dropout\n",
    "len_LAD = len(LAD_data)\n",
    "LAD_kept_num = int(len_LAD*(1-cLAD_dropout))\n",
    "np.random.seed(seed_sample)\n",
    "LAD_kept_index = random.sample(LAD_data.index, LAD_kept_num)\n",
    "LAD_data = LAD_data.loc[LAD_data.index.isin(LAD_kept_index)]\n",
    "LAD_data.reset_index(drop=True, inplace=True)\n",
    "print(\"dropout rate %s, num of cLAD, before dropout %s, after dropout %s\"%(cLAD_dropout, len_LAD,len(LAD_data)))\n",
    "LAD_data.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 104,
   "metadata": {},
   "outputs": [],
   "source": [
    "windows_sep_dir = \"/home/ygli/gam_paper_data/gam_seq_mapped4/bowtie_result/chrom_windows/\"\n",
    "# chr_name = \"chr1\"\n",
    "# windows_sep_file = \"mm10_\"+chr_name+\"_windows_10000.bed\"\n",
    "# windows_sep = pd.read_csv(windows_sep_dir+windows_sep_file, sep='\\t', header=None)\n",
    "# windows_sep.columns = ['chr','loci_head', 'loci_tail']\n",
    "# windows_sep.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 105,
   "metadata": {},
   "outputs": [],
   "source": [
    "# return the shortest dis from each element in listB to listA\n",
    "# return a dictionary with the same length as listB, with key as indexB_sort, value as the shorteds dis\n",
    "def get_intersection_shortest_distance(listA, listB, indexB):\n",
    "    listA_sort = sorted(listA)\n",
    "    listB_sort = sorted(listB)\n",
    "    indexB_sort = [x for (y,x) in sorted(zip([i[0] for i in listB], indexB))]\n",
    "    inter_index = []\n",
    "    inter_dis = defaultdict(list)\n",
    "    i,j = 0,0\n",
    "    while i < len(listA_sort) and j < len(listB_sort):\n",
    "#     while j < len(listB_sort):\n",
    "        l = max(listA_sort[i][0], listB_sort[j][0])\n",
    "        r = min(listA_sort[i][1], listB_sort[j][1])\n",
    "        if l <= r:\n",
    "            inter_index.append(indexB_sort[j])\n",
    "            inter_dis[indexB_sort[j]].append(0) #intersection\n",
    "        else:\n",
    "            inter_dis[indexB_sort[j]].append(l-r) #this two eles not intersection\n",
    "        if listA_sort[i][1] < listB_sort[j][1]:\n",
    "            i += 1\n",
    "        else:\n",
    "            j += 1\n",
    "    while j < len(listB_sort):\n",
    "        l = max(listA_sort[-1][0], listB_sort[j][0])\n",
    "        r = min(listA_sort[-1][1], listB_sort[j][1])\n",
    "        inter_dis[indexB_sort[j]].append(l-r)\n",
    "        j+=1\n",
    "    \n",
    "    inter_dis = {i:np.min(inter_dis[i]) for i in inter_dis}\n",
    "    return inter_dis\n",
    "\n",
    "# listA = [[4,7], [1,3], [5,9]]\n",
    "# # listB = [[20,50], [2,4], [5,8], [0,3]]\n",
    "# listB = [[1,3], [3,6], [6,9], [12,15]]\n",
    "# indexB = [1,2,4,0]\n",
    "# get_intersection_index(listA, listB, indexB)\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 106,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "/home/ygli/gam_paper_data/lad_data/GSE17051_cLAD_regions_mm10_10000/chr1_shortest_distance_dpot0.1_seed29.csv\n",
      "/home/ygli/gam_paper_data/lad_data/GSE17051_cLAD_regions_mm10_10000/chr2_shortest_distance_dpot0.1_seed29.csv\n",
      "/home/ygli/gam_paper_data/lad_data/GSE17051_cLAD_regions_mm10_10000/chr3_shortest_distance_dpot0.1_seed29.csv\n",
      "/home/ygli/gam_paper_data/lad_data/GSE17051_cLAD_regions_mm10_10000/chr4_shortest_distance_dpot0.1_seed29.csv\n",
      "/home/ygli/gam_paper_data/lad_data/GSE17051_cLAD_regions_mm10_10000/chr5_shortest_distance_dpot0.1_seed29.csv\n",
      "/home/ygli/gam_paper_data/lad_data/GSE17051_cLAD_regions_mm10_10000/chr6_shortest_distance_dpot0.1_seed29.csv\n",
      "/home/ygli/gam_paper_data/lad_data/GSE17051_cLAD_regions_mm10_10000/chr7_shortest_distance_dpot0.1_seed29.csv\n",
      "/home/ygli/gam_paper_data/lad_data/GSE17051_cLAD_regions_mm10_10000/chr8_shortest_distance_dpot0.1_seed29.csv\n",
      "/home/ygli/gam_paper_data/lad_data/GSE17051_cLAD_regions_mm10_10000/chr9_shortest_distance_dpot0.1_seed29.csv\n",
      "/home/ygli/gam_paper_data/lad_data/GSE17051_cLAD_regions_mm10_10000/chr10_shortest_distance_dpot0.1_seed29.csv\n",
      "/home/ygli/gam_paper_data/lad_data/GSE17051_cLAD_regions_mm10_10000/chr11_shortest_distance_dpot0.1_seed29.csv\n",
      "/home/ygli/gam_paper_data/lad_data/GSE17051_cLAD_regions_mm10_10000/chr12_shortest_distance_dpot0.1_seed29.csv\n",
      "/home/ygli/gam_paper_data/lad_data/GSE17051_cLAD_regions_mm10_10000/chr13_shortest_distance_dpot0.1_seed29.csv\n",
      "/home/ygli/gam_paper_data/lad_data/GSE17051_cLAD_regions_mm10_10000/chr14_shortest_distance_dpot0.1_seed29.csv\n",
      "/home/ygli/gam_paper_data/lad_data/GSE17051_cLAD_regions_mm10_10000/chr15_shortest_distance_dpot0.1_seed29.csv\n",
      "/home/ygli/gam_paper_data/lad_data/GSE17051_cLAD_regions_mm10_10000/chr16_shortest_distance_dpot0.1_seed29.csv\n",
      "/home/ygli/gam_paper_data/lad_data/GSE17051_cLAD_regions_mm10_10000/chr17_shortest_distance_dpot0.1_seed29.csv\n",
      "/home/ygli/gam_paper_data/lad_data/GSE17051_cLAD_regions_mm10_10000/chr18_shortest_distance_dpot0.1_seed29.csv\n",
      "/home/ygli/gam_paper_data/lad_data/GSE17051_cLAD_regions_mm10_10000/chr19_shortest_distance_dpot0.1_seed29.csv\n"
     ]
    }
   ],
   "source": [
    "save_dir = \"/home/ygli/gam_paper_data/lad_data/GSE17051_cLAD_regions_mm10_10000/\"\n",
    "# chr_name = \"chr1\"\n",
    "\n",
    "for chr_name in chr_name_list:\n",
    "# for chr_name in ['chr1']:\n",
    "    if cLAD_dropout==0:\n",
    "        save_file_name = chr_name+\"_shortest_distance.csv\"\n",
    "    else:\n",
    "        save_file_name = chr_name+\"_shortest_distance_dpot\"+str(cLAD_dropout)+\"_seed\"+str(seed_sample)+\".csv\"\n",
    "#     save_file_name = chr_name+\"_shortest_distance.csv\"\n",
    "    windows_sep_file = \"mm10_\"+chr_name+\"_windows_10000.bed\"\n",
    "    windows_sep = pd.read_csv(windows_sep_dir+windows_sep_file, sep='\\t', header=None)\n",
    "    windows_sep.columns = ['chr','loci_head', 'loci_tail']\n",
    "\n",
    "    listA = [[i.loci_head, i.loci_tail] for _,i in LAD_data[LAD_data.chr==chr_name].iterrows()]\n",
    "    listB = [[i.loci_head, i.loci_tail] for _,i in windows_sep.iterrows()]\n",
    "    indexB = [i for i,_ in windows_sep.iterrows()]\n",
    "    intersection_shortest_distance_chr = get_intersection_shortest_distance(listA, listB, indexB)\n",
    "    intersection_shortest_distance_chr_list = [intersection_shortest_distance_chr[i] for i in intersection_shortest_distance_chr]\n",
    "    head_lociB = [listB[i][0] for i in intersection_shortest_distance_chr]\n",
    "    tail_lociB = [listB[i][1] for i in intersection_shortest_distance_chr]\n",
    "    intersection_shortest_distance_chr_df = pd.DataFrame({\n",
    "        \"head_loci\": head_lociB,\n",
    "        \"tail_loci\": tail_lociB,\n",
    "        \"shortest_distance\": intersection_shortest_distance_chr_list,\n",
    "    })\n",
    "#     print(intersection_shortest_distance_chr_df.head())\n",
    "    intersection_shortest_distance_chr_df.to_csv(save_dir+save_file_name, columns=['head_loci', 'tail_loci', 'shortest_distance'], index=0)\n",
    "    print(save_dir+save_file_name)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 101,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(6144, 3)"
      ]
     },
     "execution_count": 101,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "intersection_shortest_distance_chr_df.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 66,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(6144, 3)"
      ]
     },
     "execution_count": 66,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "windows_sep.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.15"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}