generate_union_Etv2_peakset.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Generate a union set of all Etv2 ChIP-seq peaks"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Loading required package: stats4\n",
      "\n",
      "Loading required package: BiocGenerics\n",
      "\n",
      "Loading required package: parallel\n",
      "\n",
      "\n",
      "Attaching package: ‘BiocGenerics’\n",
      "\n",
      "\n",
      "The following objects are masked from ‘package:parallel’:\n",
      "\n",
      "    clusterApply, clusterApplyLB, clusterCall, clusterEvalQ,\n",
      "    clusterExport, clusterMap, parApply, parCapply, parLapply,\n",
      "    parLapplyLB, parRapply, parSapply, parSapplyLB\n",
      "\n",
      "\n",
      "The following objects are masked from ‘package:stats’:\n",
      "\n",
      "    IQR, mad, sd, var, xtabs\n",
      "\n",
      "\n",
      "The following objects are masked from ‘package:base’:\n",
      "\n",
      "    anyDuplicated, append, as.data.frame, basename, cbind, colnames,\n",
      "    dirname, do.call, duplicated, eval, evalq, Filter, Find, get, grep,\n",
      "    grepl, intersect, is.unsorted, lapply, Map, mapply, match, mget,\n",
      "    order, paste, pmax, pmax.int, pmin, pmin.int, Position, rank,\n",
      "    rbind, Reduce, rownames, sapply, setdiff, sort, table, tapply,\n",
      "    union, unique, unsplit, which, which.max, which.min\n",
      "\n",
      "\n",
      "Loading required package: S4Vectors\n",
      "\n",
      "\n",
      "Attaching package: ‘S4Vectors’\n",
      "\n",
      "\n",
      "The following object is masked from ‘package:base’:\n",
      "\n",
      "    expand.grid\n",
      "\n",
      "\n",
      "Loading required package: IRanges\n",
      "\n",
      "Loading required package: GenomeInfoDb\n",
      "\n"
     ]
    }
   ],
   "source": [
    "library(GenomicRanges)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "bed_files <- c(\n",
    "  'MEF_Dox_d1_Etv2' = 'https://s3.msi.umn.edu/gongx030/datasets/dataset=Etv2PioneerChIPseq_version=20191203a/MEF_Dox_d1_Etv2_summits.bed',\n",
    "  'MEF_Dox_d2_Etv2' = 'https://s3.msi.umn.edu/gongx030/datasets/dataset=Etv2PioneerChIPseq_version=20191203a/MEF_Dox_d2_Etv2_summits.bed',\n",
    "  'MEF_Dox_d7_Etv2' = 'https://s3.msi.umn.edu/gongx030/datasets/dataset=Etv2PioneerChIPseq_version=20191203a/MEF_Dox_d7_Etv2_summits.bed',\n",
    "  'EB_Dox_3h_Etv2' = 'https://s3.msi.umn.edu/gongx030/datasets/dataset=Etv2PioneerChIPseq_version=20191203a/EB_Dox_3h_Etv2_summits.bed',\n",
    "  'EB_Dox_12h_Etv2' = 'https://s3.msi.umn.edu/gongx030/datasets/dataset=Etv2PioneerChIPseq_version=20191203a/EB_Dox_12h_Etv2_summits.bed'\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "gr_list <- lapply(1:length(bed_files), function(i){\n",
    "  x <- read.table(bed_files[i], header = FALSE, sep = '\\t')\n",
    "  gr <- GRanges(seqnames = x[, 1], range = IRanges(x[, 2], x[, 3]), peak_id = x[, 4], score = x[, 5])\n",
    "  gr$source <- names(bed_files)[i]\n",
    "  gr\n",
    "})"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "GRanges object with 162010 ranges and 4 metadata columns:\n",
       "           seqnames            ranges strand |\n",
       "              <Rle>         <IRanges>  <Rle> |\n",
       "       [1]     chr1   3037784-3037983      * |\n",
       "       [2]     chr1   3263790-3263989      * |\n",
       "       [3]     chr1   3344682-3344881      * |\n",
       "       [4]     chr1   3460640-3460839      * |\n",
       "       [5]     chr1   3490314-3490513      * |\n",
       "       ...      ...               ...    ... .\n",
       "  [162006]     chrY 90828936-90829135      * |\n",
       "  [162007]     chrY 90836259-90836458      * |\n",
       "  [162008]     chrY 90841624-90841823      * |\n",
       "  [162009]     chrY 90842563-90842762      * |\n",
       "  [162010]     chrY 90843810-90844009      * |\n",
       "                                                                                                                 peak_id\n",
       "                                                                                                                <factor>\n",
       "       [1]      /panfs/roc/scratch/gongx030/datasets/dataset=Etv2PioneerChIPseq_version=20191203a/MEF_Dox_d2_Etv2_peak_1\n",
       "       [2]       /panfs/roc/scratch/gongx030/datasets/dataset=Etv2PioneerChIPseq_version=20191203a/EB_Dox_3h_Etv2_peak_1\n",
       "       [3]      /panfs/roc/scratch/gongx030/datasets/dataset=Etv2PioneerChIPseq_version=20191203a/MEF_Dox_d2_Etv2_peak_2\n",
       "       [4]      /panfs/roc/scratch/gongx030/datasets/dataset=Etv2PioneerChIPseq_version=20191203a/MEF_Dox_d2_Etv2_peak_3\n",
       "       [5]      /panfs/roc/scratch/gongx030/datasets/dataset=Etv2PioneerChIPseq_version=20191203a/MEF_Dox_d2_Etv2_peak_4\n",
       "       ...                                                                                                           ...\n",
       "  [162006] /panfs/roc/scratch/gongx030/datasets/dataset=Etv2PioneerChIPseq_version=20191203a/MEF_Dox_d2_Etv2_peak_61360b\n",
       "  [162007]  /panfs/roc/scratch/gongx030/datasets/dataset=Etv2PioneerChIPseq_version=20191203a/MEF_Dox_d7_Etv2_peak_16842\n",
       "  [162008] /panfs/roc/scratch/gongx030/datasets/dataset=Etv2PioneerChIPseq_version=20191203a/MEF_Dox_d1_Etv2_peak_114014\n",
       "  [162009] /panfs/roc/scratch/gongx030/datasets/dataset=Etv2PioneerChIPseq_version=20191203a/MEF_Dox_d1_Etv2_peak_114015\n",
       "  [162010] /panfs/roc/scratch/gongx030/datasets/dataset=Etv2PioneerChIPseq_version=20191203a/MEF_Dox_d1_Etv2_peak_114016\n",
       "               score          source                 group\n",
       "           <numeric>     <character>              <matrix>\n",
       "       [1]  16.22768 MEF_Dox_d2_Etv2   TRUE:TRUE:FALSE:...\n",
       "       [2]   8.36521  EB_Dox_3h_Etv2 FALSE:FALSE:FALSE:...\n",
       "       [3]   9.73217 MEF_Dox_d2_Etv2  FALSE:TRUE:FALSE:...\n",
       "       [4]   3.92081 MEF_Dox_d2_Etv2  FALSE:TRUE:FALSE:...\n",
       "       [5]   12.2973 MEF_Dox_d2_Etv2   TRUE:TRUE:FALSE:...\n",
       "       ...       ...             ...                   ...\n",
       "  [162006] 160.15463 MEF_Dox_d2_Etv2    TRUE:TRUE:TRUE:...\n",
       "  [162007]   8.27339 MEF_Dox_d7_Etv2  FALSE:FALSE:TRUE:...\n",
       "  [162008]   1.53792 MEF_Dox_d1_Etv2  TRUE:FALSE:FALSE:...\n",
       "  [162009]  11.78216 MEF_Dox_d1_Etv2  TRUE:FALSE:FALSE:...\n",
       "  [162010]   5.92807 MEF_Dox_d1_Etv2  TRUE:FALSE:FALSE:...\n",
       "  -------\n",
       "  seqinfo: 21 sequences from an unspecified genome; no seqlengths"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "gr <- Reduce('c', gr_list)\n",
    "gr <- resize(gr, width = 200, fix = 'center')\n",
    "gr0 <- reduce(gr) # a reduced set of Etv2 peak\n",
    "mm <- as.matrix(findOverlaps(gr, gr0))\n",
    "sp <- split(mm[, 1], list(mm[, 2]))\n",
    "j <- unlist(mclapply(sp, function(i) i[which.max(gr$score[i])], mc.cores = 4))\n",
    "gr <- gr[j]\n",
    "G <- do.call('cbind', lapply(1:length(bed_files), function(i) gr %over% gr_list[[i]]))\n",
    "gr$group <- G # group assignment of the union set to each source\n",
    "colnames(gr$group) <- names(bed_files)\n",
    "gr"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "gr_file <- 'data/all_Etv2_peaks.rds'\n",
    "# The file also exist at \n",
    "# https://s3.msi.umn.edu/gongx030/datasets/dataset=Etv2PioneerChIPseq_version=20191203a/all_Etv2_peaks.rds'\n",
    "saveRDS(gr, gr_file)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "R",
   "language": "R",
   "name": "ir"
  },
  "language_info": {
   "codemirror_mode": "r",
   "file_extension": ".r",
   "mimetype": "text/x-r-source",
   "name": "R",
   "pygments_lexer": "r",
   "version": "3.6.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}