02_TADBit_HOMER_Overlap_SNPs.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "toc": true
   },
   "source": [
    "<h1>Table of Contents<span class=\"tocSkip\"></span></h1>\n",
    "<div class=\"toc\"><ul class=\"toc-item\"><li><span><a href=\"#TADBit-HiC-processing-parameters-and-commands\" data-toc-modified-id=\"TADBit-HiC-processing-parameters-and-commands-1\"><span class=\"toc-item-num\">1&nbsp;&nbsp;</span>TADBit HiC processing parameters and commands</a></span></li><li><span><a href=\"#TADBit-merging-7-replicas\" data-toc-modified-id=\"TADBit-merging-7-replicas-2\"><span class=\"toc-item-num\">2&nbsp;&nbsp;</span>TADBit merging 7 replicas</a></span><ul class=\"toc-item\"><li><span><a href=\"#Change-to-HOMER-format\" data-toc-modified-id=\"Change-to-HOMER-format-2.1\"><span class=\"toc-item-num\">2.1&nbsp;&nbsp;</span>Change to HOMER format</a></span></li></ul></li><li><span><a href=\"#HOMER-processing\" data-toc-modified-id=\"HOMER-processing-3\"><span class=\"toc-item-num\">3&nbsp;&nbsp;</span>HOMER processing</a></span></li><li><span><a href=\"#Overlapping-the-significant-interactions-with-SNPs\" data-toc-modified-id=\"Overlapping-the-significant-interactions-with-SNPs-4\"><span class=\"toc-item-num\">4&nbsp;&nbsp;</span>Overlapping the significant interactions with SNPs</a></span></li></ul></div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**With this script you can process the HiC data (from valid pairs) and overlap it with our 624 LMI selected SNPs to get their interacting counter parts.**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### TADBit HiC processing parameters and commands"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We processed our 7 replicas from pancreas HiC experiment using TADbit (Serra *et al.,* PLoS Comp Bio, 2017) to obtain one file with the valid pair read interactions for each replica."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Start here below processing the valid pairs from your samples:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### TADBit merging 7 replicas"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Count the valid pairs at each of the replicas\n",
    "parallel wc -l /home/jrodriguez/scratch/pancreas_HiC/pancreas{}/03_filtered_reads/valid* ::: {1..7}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "2269506 /home/jrodriguez/scratch/pancreas_HiC/pancreas5/03_filtered_reads/valid_r1-r2_intersection_e7dd1ea461.tsv\n",
    " \n",
    "2578229 /home/jrodriguez/scratch/pancreas_HiC/pancreas4/03_filtered_reads/valid_r1-r2_intersection_e7dd1ea461.tsv\n",
    "\n",
    "3759181 /home/jrodriguez/scratch/pancreas_HiC/pancreas6/03_filtered_reads/valid_r1-r2_intersection_e7dd1ea461.tsv \n",
    "\n",
    "5127448 /home/jrodriguez/scratch/pancreas_HiC/pancreas3/03_filtered_reads/valid_r1-r2_intersection_e7dd1ea461.tsv\n",
    "\n",
    "19815943 /home/jrodriguez/scratch/pancreas_HiC/pancreas7/03_filtered_reads/valid_r1-r2_intersection_e7dd1ea461.tsv\n",
    "\n",
    "24297763 /home/jrodriguez/scratch/pancreas_HiC/pancreas1/03_filtered_reads/valid_r1-r2_intersection_e7dd1ea461.tsv\n",
    "\n",
    "41226187 /home/jrodriguez/scratch/pancreas_HiC/pancreas2/03_filtered_reads/valid_r1-r2_intersection_e7dd1ea461.tsv"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**All summed up result in 99,074,257**\n",
    "Each valid pairs files has a header of 25 lines. 25x7=175\n",
    "\n",
    "99,074,257 - 175 = **99,074,082** valid pairs, which correspond to the number seen after the merging of the 7 replicas"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Concatenate the seven replicas\n",
    "cat {YOUR_PATH}/pancreas{1..7}/03*/valid* | grep -v '^#' > valid_r1-r2_7replicas.tsv"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Change to HOMER format"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Change TADbit valid pairs file to HOMER format\n",
    "awk -v OFS=$'\\t' '{print $1,$2,$3,$4,$8,$9,$10}' valid_r1-r2_7replicas.tsv \\| \n",
    "awk -F\"\\t\" '{ $4 = ($4 == \"0\" ? \"-\" : \"+\") } 1' OFS=\"\\t\" \\| \n",
    "awk -F\"\\t\" '{ $7 = ($7 == \"0\" ? \"-\" : \"+\") } 1' OFS=\"\\t\" > valid_r1-r2_7replicas_HOMER.tsv"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### HOMER processing"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "create a tag directory first"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "makeTagDirectory {YOUR_PATH}/hic_data_healthy_homer # name of your to-be-created tag directory\n",
    "    -format HiCsummary # format of the interactions files \n",
    "    valid_r1-r2_7replicas_HOMER.tsv # your file"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "create a background model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "analyzeHiC {YOUR_PATH}hic_data_healthy_homer \n",
    "    -res 40000 # resolution\n",
    "    -bgonly # generate a background model\n",
    "    -cpu 24"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "and finally get the significant interactions "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "analyzeHiC {YOUR_PATH}hic_data_healthy_homer \n",
    "    -res 40000 \n",
    "    -interactions significantInteractions_40000_pancreas_17_new # significant interactions output\n",
    "    -center # centering the interactions based on read average mapping\n",
    "    -maxDist 2000000 # maximum distance to find interactions\n",
    "    -nomatrix \n",
    "    -cpu 24"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We get a total of `41,832` significant interactions. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "By default, `HOMER` reports interactions with pv <= 0.001. Each interaction has its pv based on a test and is reported on the file that we get after the command above. We decided that we could further filter that output, setting the threshold based on Bonferroni correction, estimated from the number of tests performed. *(See main text for a more detailed explanation)* "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "log(1e-05) # = -11.51293 (This is an R command)\n",
    "# Filter the Z-score LogP field by those that have lower than -11.51293\n",
    "logp=-11.51293\n",
    "# Filtering the interactions obtained\n",
    "awk -F'\\t' -v b=${logp} ' $18 <= b ' significantInteractions_40000_pancreas_17_new > significantInteractions_40000_pancreas_17_new_binp1e5"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Overlapping the significant interactions with SNPs "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we select the both sides of the interactions..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "awk -F'\\t' -v OFS='\\t' '{print $3,$4,$5,$1}' significantInteractions_40000_pancreas_17_new_binp1e5 > interactions_40000_pancreas_17_site1_binp1e5.bed\n",
    "awk -F'\\t' -v OFS='\\t' '{print $9,$10,$11,$1}' significantInteractions_40000_pancreas_17_new_binp1e5 > interactions_40000_pancreas_17_site2_binp1e5.bed"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Run these commands to get the overlap of the `624 GWAS-LMI SNPs` with the interactions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "input=\"LMI_pv.snps.bed\" # the file with the GWAS-LMI SNPs\n",
    "binp=\"1e5\"\n",
    "mod=\"LMI_pv_snp_new\"\n",
    "cd {YOUR_PATH}/go_analysis_with_new_interactions/GO_snp ; "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Interactions\n",
    "intersectBed -a ../interactions_40000_pancreas_17_site2_binp1e5.bed -b ${input} -wa -wb | cut -f4 > site2_${mod}_binp${binp}.interactions\n",
    "intersectBed -a ../interactions_40000_pancreas_17_site1_binp1e5.bed -b ${input} -wa  -wb | cut -f4 > site1_${mod}_binp${binp}.interactions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "direct_bins=($(grep -Fwf site1_${mod}_binp${binp}.interactions ../significantInteractions_40000_pancreas_17_new_binp1e5 | awk -F'\\t' -v OFS=$'\\t' '{print $1,$2,$3,$4,$5,$6}';\n",
    "               grep -Fwf site2_${mod}_binp${binp}.interactions ../significantInteractions_40000_pancreas_17_new_binp1e5 | awk -F'\\t' -v OFS=$'\\t' '{print $1,$8,$9,$10,$11,$12}'))\n",
    "printf \"%s\\t%s\\t%s\\t%s\\t%s\\t%s\\n\" \"${direct_bins[@]}\" | sort -k1 > GWAS_bait_${mod}_bins_binp${binp}\n",
    "\n",
    "interacting_bins=($(grep -Fwf site1_${mod}_binp${binp}.interactions ../significantInteractions_40000_pancreas_17_new_binp1e5 | awk -F'\\t' -v OFS=$'\\t' '{print $1,$8,$9,$10,$11,$12}' ; \n",
    "                    grep -Fwf site2_${mod}_binp${binp}.interactions ../significantInteractions_40000_pancreas_17_new_binp1e5 | awk -F'\\t' -v OFS=$'\\t' '{print $1,$2,$3,$4,$5,$6}'))\n",
    "printf \"%s\\t%s\\t%s\\t%s\\t%s\\t%s\\n\" \"${interacting_bins[@]}\" | sort -k1 > GWAS_target_${mod}_bins_binp${binp}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "## Identify and remove from the target list those interactions in which bait and target are overlapping a same bin.\n",
    "## This will lead to targets (3D genomics) not identifying anything different \n",
    "## from the standard approach.\n",
    "## First, identifies them and then excludes them from the target list\n",
    "grep -vf <(comm -12 \\\n",
    "         <(cut -f2 GWAS_bait_LMI_pv_snp_new_bins_binp1e5 | sort | uniq ) \\ \n",
    "         <(cut -f2 GWAS_target_LMI_pv_snp_new_bins_binp1e5 | sort | uniq)) \\\n",
    "         GWAS_target_LMI_pv_snp_new_bins_binp1e5 > GWAS_target_LMI_pv_snp_new_bins_binp1e5.NoOverlapBaitTarget"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# change the format of the baits file to a bed to intersect with the SNPs\n",
    "awk '{print}' GWAS_bait_${mod}_bins_binp${binp} > bait_overlap_LMI_pv_snps.bed\n",
    "\n",
    "#### TO FINISH THE AWK COMMAND TO SELECT FIELDS "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# identifies the interactions that overlap SNPs on both sides\n",
    "comm -13 <(cut -f1 GWAS_target_LMI_pv_snp_new_bins_binp1e5.NoOverlapBaitTarget | sort | uniq) <(cut -f1 GWAS_target_LMI_pv_snp_new_bins_binp1e5 | sort | uniq) > both_bait_and_target_overlap.ints \n",
    "# removing the interactions we do not want\n",
    "grep -vf both_bait_and_target_overlap.ints <(intersectBed -a bait_overlap_LMI_pv_snps.bed -b LMI_pv.snps.bed -wa -wb | sort -k4,4) | cut -f5- | sort | uniq > 76_bait_regions.SNPs\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Create a combined table with baits, targets, snps and genes.**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The baits: \n",
    "\n",
    "*(is the same that we just got `76_bait_regions.SNPs`, but with more fields.)*"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "collapsed": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "14\t53952691\t53992691\tinteraction1007\t14\t53980388\t53980389\trs1959834\n",
      "18\t76361162\t76401162\tinteraction1267\t18\t76392789\t76392790\trs2007105\n",
      "11\t29314259\t29354259\tinteraction1371\t11\t29353459\t29353460\trs12226886\n",
      "14\t53916156\t53956156\tinteraction1422\t14\t53920178\t53920179\trs11849931\n",
      "14\t53916156\t53956156\tinteraction1422\t14\t53924495\t53924496\trs11850783\n",
      "14\t53916156\t53956156\tinteraction1422\t14\t53942112\t53942113\trs7153908\n",
      "20\t32436028\t32476028\tinteraction1448\t20\t32456862\t32456863\trs62209634\n",
      "22\t28477385\t28517385\tinteraction1478\t22\t28502867\t28502868\tchr22_28502867_A_C\n",
      "22\t28477385\t28517385\tinteraction1478\t22\t28504183\t28504184\tchr22_28504183_A_G\n",
      "22\t28477385\t28517385\tinteraction1478\t22\t28504258\t28504259\tchr22_28504258_A_G\n",
      "grep: write error: Broken pipe\n"
     ]
    }
   ],
   "source": [
    "grep -vf both_bait_and_target_overlap.ints \\ \n",
    "<(intersectBed -a bait_overlap_LMI_pv_snps.bed -b LMI_pv.snps.bed -wa -wb | sort -k4,4) > bait_overlaps.snps"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we need to cross this file above, which contains all the SNPs that overlap a bait counterpart, repeated as many times as SNPs are overlaping each interaction, with the file of individual targets, see below:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "collapsed": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "interaction1007\t14-53520000\t14\t53513149\t53553149\t+\n",
      "interaction1267\t18-76680000\t18\t76676433\t76716433\t+\n",
      "interaction1371\t11-27680000\t11\t27683551\t27723551\t+\n",
      "interaction1422\t14-53520000\t14\t53513872\t53553872\t+\n",
      "interaction1448\t20-32400000\t20\t32407377\t32447377\t+\n",
      "interaction1478\t22-28240000\t22\t28237254\t28277254\t+\n",
      "interaction1581\t12-132920000\t12\t132924371\t132964371\t+\n",
      "interaction1778\t18-30600000\t18\t30604852\t30644852\t+\n",
      "interaction1982\t9-108520000\t9\t108521161\t108561161\t+\n",
      "interaction2480\t20-32280000\t20\t32282049\t32322049\t+\n"
     ]
    }
   ],
   "source": [
    "# File with target for all these. Just checking :)\n",
    "head -n10 GWAS_target_LMI_pv_snp_new_bins_binp1e5.NoOverlapBaitTarget"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Join the files\n",
    "awk -v OFS=$'\\t' 'NR==FNR {a[$1]=$0; next } $4 in a { print $0, a[$4]}' \\ \n",
    "GWAS_target_LMI_pv_snp_new_bins_binp1e5.NoOverlapBaitTarget bait_overlaps.snps > all_baits_and_targets.txt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**## END ## **"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Bash",
   "language": "bash",
   "name": "bash"
  },
  "language_info": {
   "codemirror_mode": "shell",
   "file_extension": ".sh",
   "mimetype": "text/x-sh",
   "name": "bash"
  },
  "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": true,
   "toc_position": {},
   "toc_section_display": true,
   "toc_window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}