lectures/lecture15/Lecture15_GenomicData.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "b5480d1d-59d0-4c46-94de-88da8f93e305",
   "metadata": {
    "tags": []
   },
   "source": [
    "# `GenomicRanges` Object to Store Genomic Data\n",
    "\n",
    "Genomic data is often described using chromosomes and coordinates. A locus can be a single base position or a region that includes a start and end coordinate. In R, there is a Bioconductor package called `GenomicRanges` that stores this in a convenient structure for efficient querying using routine operations. `GRanges` object class is in which genomic data will be stored. We will demonstrate the most common operation, `findOverlaps`, to determine intersecting positions or regions in the genome. See https://bioconductor.org/packages/release/bioc/vignettes/GenomicRanges/inst/doc/GenomicRangesIntroduction.html\n",
    "\n",
    "An additional package called `plyranges` provides convenient syntax similar to that used in `tidyverse` to manipulate and apply operations on `Granges` objects. See https://www.bioconductor.org/packages/devel/bioc/vignettes/plyranges/inst/doc/an-introduction.html\n",
    "\n",
    "In this tutorial, we will work with The Cancer Genome Atlas (TCGA) data for primary breast cancer patient samples. Specifically, these are segmentation data used for copy number alteration analysis. See Lecture 16: Slide 47.\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9a5d67aa-4a23-41da-a1ac-a9425afa9de8",
   "metadata": {},
   "source": [
    "## 0. Load the `GenomicRanges` Bioconductor package"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "6f1c4d84-3e24-48a0-8faa-69ecf66a6cad",
   "metadata": {},
   "outputs": [],
   "source": [
    "#Sys.unsetenv(\"DISPLAY\")\n",
    "\n",
    "suppressPackageStartupMessages({\n",
    "    #library(tidyverse)\n",
    "    library(GenomicRanges)\n",
    "    library(plyranges)\n",
    "    #library(VariantAnnotation)\n",
    "    #library(Rsamtools)\n",
    "})"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "43262063-f461-4338-a941-2039a9c49daa",
   "metadata": {},
   "source": [
    "## 1. Create a GRanges object.\n",
    "A `GRanges` object must contain an attribute called `seqnames` to represent chromosomes and `ranges` attribute to represent the `start` and `end` coordinates. The range is 1-index-based (as opposed to 0-index), The `start` and `end` can be the same value if it is a single base-pair.  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "74547f46-b65d-4240-8abb-a9556d8e3aee",
   "metadata": {},
   "outputs": [],
   "source": [
    "myGRange <- GRanges(seqnames = \"17\",\n",
    "                    ranges = IRanges(start = 37844393, end = 37844393))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fab8315c-f5df-45c7-be32-4e2966e6d900",
   "metadata": {},
   "source": [
    "Alternatively, using `plyranges` we can use familiar syntax to create the same `GRanges` object."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "7ff1c62e-af45-463e-891c-bcbe87163178",
   "metadata": {},
   "outputs": [],
   "source": [
    "myGRange <- data.frame(seqnames = \"17\", start = 37844393, end = 37844393) %>% as_granges()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "23852a67-dd82-4160-8510-c8ed4ecbb6e8",
   "metadata": {
    "tags": []
   },
   "source": [
    "## 2. Load Genomic Data From A File\n",
    "There are numerous text file formats for representing genomic data and some of these were discussed in Lecture 16. Here, we will show you that a `GRanges` can be easily created from any text file that contains delimited columns specifying genomic coorindates.\n",
    "\n",
    "### 2.1 SEG format\n",
    "SEGment Data (http://software.broadinstitute.org/software/igv/SEG) format is tab-delimited and a flexible way to define any genomic data.\n",
    "\n",
    "There are 4 required columns:\n",
    "\n",
    "  1. Name\n",
    "  2. Chromosome\n",
    "  3. Start Coordinate\n",
    "  4. End Coordinate\n",
    "\n",
    "This is similar to the BED file format but with the additional requirement for *Name* as the first column.\n",
    "\n",
    "### a. Load the SEG file containing the segments into a `data.frame` object."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "3e81b5fb-30e5-4b1f-876d-72b20c8dc4be",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "'/fh/fast/henikoff_s/user/jgreene/projects/sciCnT/code'"
      ],
      "text/latex": [
       "'/fh/fast/henikoff\\_s/user/jgreene/projects/sciCnT/code'"
      ],
      "text/markdown": [
       "'/fh/fast/henikoff_s/user/jgreene/projects/sciCnT/code'"
      ],
      "text/plain": [
       "[1] \"/fh/fast/henikoff_s/user/jgreene/projects/sciCnT/code\""
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/html": [
       "'/fh/fast/henikoff_s/user/jgreene/projects/TFCB/data'"
      ],
      "text/latex": [
       "'/fh/fast/henikoff\\_s/user/jgreene/projects/TFCB/data'"
      ],
      "text/markdown": [
       "'/fh/fast/henikoff_s/user/jgreene/projects/TFCB/data'"
      ],
      "text/plain": [
       "[1] \"/fh/fast/henikoff_s/user/jgreene/projects/TFCB/data\""
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "#old wdir\n",
    "getwd()\n",
    "\n",
    "#new wdir\n",
    "data_dir = '/fh/fast/henikoff_s/user/jgreene/projects/TFCB/data'\n",
    "setwd(data_dir)\n",
    "getwd()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "bd6af1e4-ad17-4f09-81ca-af9cfedc6fcb",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<style>\n",
       ".list-inline {list-style: none; margin:0; padding: 0}\n",
       ".list-inline>li {display: inline-block}\n",
       ".list-inline>li:not(:last-child)::after {content: \"\\00b7\"; padding: 0 .5ex}\n",
       "</style>\n",
       "<ol class=list-inline><li>'BRCA_IDC_cfDNA.bam'</li><li>'BRCA_IDC_cfDNA.bam.bai'</li><li>'BRCA.genome_wide_snp_6_broad_Level_3_scna.seg'</li><li>'drive-download-20221115T203016Z-003.zip'</li><li>'GIAB_highconf_v.3.3.2.vcf.gz'</li><li>'GIAB_highconf_v.3.3.2.vcf.gz.tbi'</li><li>'hs37d5-002.fa'</li><li>'hs37d5.fa.gz.fai'</li></ol>\n"
      ],
      "text/latex": [
       "\\begin{enumerate*}\n",
       "\\item 'BRCA\\_IDC\\_cfDNA.bam'\n",
       "\\item 'BRCA\\_IDC\\_cfDNA.bam.bai'\n",
       "\\item 'BRCA.genome\\_wide\\_snp\\_6\\_broad\\_Level\\_3\\_scna.seg'\n",
       "\\item 'drive-download-20221115T203016Z-003.zip'\n",
       "\\item 'GIAB\\_highconf\\_v.3.3.2.vcf.gz'\n",
       "\\item 'GIAB\\_highconf\\_v.3.3.2.vcf.gz.tbi'\n",
       "\\item 'hs37d5-002.fa'\n",
       "\\item 'hs37d5.fa.gz.fai'\n",
       "\\end{enumerate*}\n"
      ],
      "text/markdown": [
       "1. 'BRCA_IDC_cfDNA.bam'\n",
       "2. 'BRCA_IDC_cfDNA.bam.bai'\n",
       "3. 'BRCA.genome_wide_snp_6_broad_Level_3_scna.seg'\n",
       "4. 'drive-download-20221115T203016Z-003.zip'\n",
       "5. 'GIAB_highconf_v.3.3.2.vcf.gz'\n",
       "6. 'GIAB_highconf_v.3.3.2.vcf.gz.tbi'\n",
       "7. 'hs37d5-002.fa'\n",
       "8. 'hs37d5.fa.gz.fai'\n",
       "\n",
       "\n"
      ],
      "text/plain": [
       "[1] \"BRCA_IDC_cfDNA.bam\"                           \n",
       "[2] \"BRCA_IDC_cfDNA.bam.bai\"                       \n",
       "[3] \"BRCA.genome_wide_snp_6_broad_Level_3_scna.seg\"\n",
       "[4] \"drive-download-20221115T203016Z-003.zip\"      \n",
       "[5] \"GIAB_highconf_v.3.3.2.vcf.gz\"                 \n",
       "[6] \"GIAB_highconf_v.3.3.2.vcf.gz.tbi\"             \n",
       "[7] \"hs37d5-002.fa\"                                \n",
       "[8] \"hs37d5.fa.gz.fai\"                             "
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "list.files()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "199d73fb-1c0b-40a1-839d-e9c564241d97",
   "metadata": {},
   "outputs": [],
   "source": [
    "segs <- read.table(\"BRCA.genome_wide_snp_6_broad_Level_3_scna.seg\", header = TRUE)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e1842cab-a987-43bf-bf6c-e5c4798ea978",
   "metadata": {},
   "source": [
    "Small processing of this file to correct a few legacy hacks. We need to change chromosome 23 to chromosome X."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "99bdf8fe-023a-4be3-b9f9-036881b85f2d",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "'data.frame':\t284458 obs. of  6 variables:\n",
      " $ Sample      : chr  \"TCGA-3C-AAAU-10A-01D-A41E-01\" \"TCGA-3C-AAAU-10A-01D-A41E-01\" \"TCGA-3C-AAAU-10A-01D-A41E-01\" \"TCGA-3C-AAAU-10A-01D-A41E-01\" ...\n",
      " $ Chromosome  : int  1 1 1 1 1 1 1 1 1 2 ...\n",
      " $ Start       : int  3218610 95676511 95680124 167057495 167059760 181603120 181610685 201474400 201475220 484222 ...\n",
      " $ End         : int  95674710 95676518 167057183 167059336 181602002 181609567 201473647 201474544 247813706 51515129 ...\n",
      " $ Num_Probes  : int  53225 2 24886 3 9213 6 12002 2 29781 30300 ...\n",
      " $ Segment_Mean: num  0.0055 -1.6636 0.0053 -1.0999 -0.0008 ...\n"
     ]
    }
   ],
   "source": [
    "str(segs) # show the class type for each column\n",
    "mode(segs$Chromosome) <- \"character\" # change the class of the chromosome to character\n",
    "segs[segs$Chromosome == 23, \"Chromosome\"] <- \"X\""
   ]
  },
  {
   "cell_type": "markdown",
   "id": "806965e2-7ee9-4b2d-9c6d-0b9804e4ad9c",
   "metadata": {},
   "source": [
    "### b. Convert the `data.frame` object into a `GRanges`. \n",
    "You can use the `as()` function, as long as the 3 required columns are present. It is also flexible how the columns are named. For example, the column can be `Start`, `start`, `Chr`, `chr`, `Chromosome`, `End`, `Stop`, etc."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "80b67f8f-186a-41e5-9f65-86cefae245d3",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "GRanges object with 284458 ranges and 3 metadata columns:\n",
       "           seqnames              ranges strand |                 Sample\n",
       "              <Rle>           <IRanges>  <Rle> |            <character>\n",
       "       [1]        1    3218610-95674710      * | TCGA-3C-AAAU-10A-01D..\n",
       "       [2]        1   95676511-95676518      * | TCGA-3C-AAAU-10A-01D..\n",
       "       [3]        1  95680124-167057183      * | TCGA-3C-AAAU-10A-01D..\n",
       "       [4]        1 167057495-167059336      * | TCGA-3C-AAAU-10A-01D..\n",
       "       [5]        1 167059760-181602002      * | TCGA-3C-AAAU-10A-01D..\n",
       "       ...      ...                 ...    ... .                    ...\n",
       "  [284454]       19     284018-58878226      * | TCGA-Z7-A8R6-01A-11D..\n",
       "  [284455]       20     455764-62219837      * | TCGA-Z7-A8R6-01A-11D..\n",
       "  [284456]       21   15347621-47678774      * | TCGA-Z7-A8R6-01A-11D..\n",
       "  [284457]       22   17423930-49331012      * | TCGA-Z7-A8R6-01A-11D..\n",
       "  [284458]        X   3157107-154905589      * | TCGA-Z7-A8R6-01A-11D..\n",
       "           Num_Probes Segment_Mean\n",
       "            <integer>    <numeric>\n",
       "       [1]      53225       0.0055\n",
       "       [2]          2      -1.6636\n",
       "       [3]      24886       0.0053\n",
       "       [4]          3      -1.0999\n",
       "       [5]       9213      -0.0008\n",
       "       ...        ...          ...\n",
       "  [284454]      23950      -0.1170\n",
       "  [284455]      37283       0.3435\n",
       "  [284456]      20582      -0.1117\n",
       "  [284457]      16927      -0.1231\n",
       "  [284458]      63797       0.0014\n",
       "  -------\n",
       "  seqinfo: 23 sequences from an unspecified genome; no seqlengths"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "segs.gr <- as(segs, \"GRanges\")\n",
    "segs.gr"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1d8740b2-11b4-4d2b-bff2-5ecdca69d1b1",
   "metadata": {},
   "source": [
    "Alternatively, using `plyranges`. Here, we need to rename the columns: `Chromosome`->`seqnames`, `Start`->`start`, `End`->`end`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "2378a7c9-2948-413c-9ade-f43ab3ee4ba1",
   "metadata": {},
   "outputs": [],
   "source": [
    "colnames(segs)[2:4] <- c(\"seqnames\", \"start\", \"end\")\n",
    "segs.gr <- segs %>% as_granges()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3137703e-9d8b-4d4f-9edb-f20aa6690681",
   "metadata": {},
   "source": [
    "## 3. Operations and features of GenomicRanges\n",
    "Some of the most useful features of `GRanges` object is the fast and easy methods for determining overlaps between sets of ranges. Here, we will describe examples using some of the common functions.\n",
    "\n",
    "### 3.1 Tiling the genome\n",
    "Often we would like to *find* or *count* events overlapping regions in the genome. In an unbiased fashion, we could do this genome-wide by dividing the genome into tiles/windows/bins. \n",
    "\n",
    "We will use the `tileGenome()` for this task, which requires three arguments: length of the chromosomes, number of tiles and the size of each tile.\n",
    "\n",
    "### a. We need the lengths of the chromosomes in the human genome.\n",
    "We need to load human genome information for build `hg19`. Since there are non-standard chromosomes, we only want to keep the standard chromosomes using `keepStandardChromosomes`. Then, since our `segs` data uses `NCBI` chromosome naming convention (i.e. `1` instead of `chr1`), we need set the `seqlevelStyle`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "eeea07c7-5f7f-48a0-9e36-09ae99bfd3f1",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Warning message in (function (seqlevels, genome, new_style) :\n",
      "“cannot switch some of hg19's seqlevels from UCSC to NCBI style”\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "Seqinfo object with 25 sequences (1 circular) from 2 genomes (GRCh37.p13, hg19):\n",
       "  seqnames seqlengths isCircular     genome\n",
       "  1         249250621      FALSE GRCh37.p13\n",
       "  2         243199373      FALSE GRCh37.p13\n",
       "  3         198022430      FALSE GRCh37.p13\n",
       "  4         191154276      FALSE GRCh37.p13\n",
       "  5         180915260      FALSE GRCh37.p13\n",
       "  ...             ...        ...        ...\n",
       "  21         48129895      FALSE GRCh37.p13\n",
       "  22         51304566      FALSE GRCh37.p13\n",
       "  X         155270560      FALSE GRCh37.p13\n",
       "  Y          59373566      FALSE GRCh37.p13\n",
       "  chrM          16571       TRUE       hg19"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "seqinfo <- Seqinfo(genome = \"hg19\")\n",
    "seqinfo <- keepStandardChromosomes(seqinfo) \n",
    "seqlevelsStyle(seqinfo) <- \"NCBI\"\n",
    "seqinfo"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "30fb6344-ff40-4e84-8be7-5751d3f65dcc",
   "metadata": {},
   "source": [
    "### b. Split the genome into 500kb tiles or windows."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "fbcbd245-3ff9-4159-9e13-43f718b12508",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "GRanges object with 6207 ranges and 0 metadata columns:\n",
       "         seqnames            ranges strand\n",
       "            <Rle>         <IRanges>  <Rle>\n",
       "     [1]        1          1-500000      *\n",
       "     [2]        1    500001-1000000      *\n",
       "     [3]        1   1000001-1500000      *\n",
       "     [4]        1   1500001-2000000      *\n",
       "     [5]        1   2000001-2500000      *\n",
       "     ...      ...               ...    ...\n",
       "  [6203]        Y 57500001-58000000      *\n",
       "  [6204]        Y 58000001-58500000      *\n",
       "  [6205]        Y 58500001-59000000      *\n",
       "  [6206]        Y 59000001-59373566      *\n",
       "  [6207]     chrM           1-16571      *\n",
       "  -------\n",
       "  seqinfo: 25 sequences from an unspecified genome"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "slen <- seqlengths(seqinfo) # get the length of the chromosomes\n",
    "tileWidth <- 500000 # tile size of 500kb\n",
    "tiles <- tileGenome(seqlengths = slen, tilewidth = tileWidth,\n",
    "                    cut.last.tile.in.chrom = TRUE)\n",
    "tiles"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "007db130-94ed-4923-b032-ca4039c07465",
   "metadata": {},
   "source": [
    "### 3.2 Finding overlap of ranges\n",
    "One of the most useful features of `GenomicRanges` is to simply identify the ranges that overlap between two `GRanges` objects. The `findOverlaps` function is a basic method in the `GRanges` class for finding the overlaps of the elements that overlap between two `GRanges`. The argmuents `query` for your main `tiles.subset` and `subject` for the `segs.gr`. The `type` argument describes the type of overlap, such as `any`, `within`, `start`, `end`, `equal`, and there are additional arguments for criteria for overlap such as `minoverlap` size.\n",
    "\n",
    "For this example, let's find which copy number alteration segments from `segs.gr` overlap in *any* way with our ranges in `tiles.subset` (`17:35500000-37000000`). \n",
    "\n",
    "### a. Find the tiled ranges for chromosome `17`, starting `35500000` and ending `37000000`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "c7bfd88b-6686-40e6-93c2-686b97b429e5",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "GRanges object with 3 ranges and 0 metadata columns:\n",
       "      seqnames            ranges strand\n",
       "         <Rle>         <IRanges>  <Rle>\n",
       "  [1]       17 35500001-36000000      *\n",
       "  [2]       17 36000001-36500000      *\n",
       "  [3]       17 36500001-37000000      *\n",
       "  -------\n",
       "  seqinfo: 25 sequences from an unspecified genome"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "tiles.subset <- tiles[seqnames(tiles) == \"17\" & start(tiles) >= 35500000 & end(tiles) <= 37000000]\n",
    "tiles.subset"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "98195ecf-c9f2-4cfd-a333-02589399dafd",
   "metadata": {},
   "source": [
    "Alternatively, using `plyranges` and the `filter` verb."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "7059ca7b-70fd-4df9-960b-1081389a159c",
   "metadata": {},
   "outputs": [],
   "source": [
    "tiles.subset <- tiles %>% filter(seqnames == \"17\" & start >= 35500000 & end <= 37000000)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "120d0dd4-09d6-43cb-ab60-2f69155c01d3",
   "metadata": {},
   "source": [
    "### b. Find the overlap between `segs.gr` and `tiles.subset`.\n",
    "Let's find the segments in `segs.gr` (`query`) that overlap our `tiles.subset` (`subject`).\n",
    "`plyranges` provides convenient functions that can bypass having to deal with hits/indices and returns the overlapped regions. Here, we use the function `find_overlaps`. This function will return all of the ranges in the `query` that overlap with the `subject`.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "f138cb5e-8ba6-4dbe-83c2-ff86e47cf1cc",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "GRanges object with 2 ranges and 3 metadata columns:\n",
       "      seqnames          ranges strand |                 Sample Num_Probes\n",
       "         <Rle>       <IRanges>  <Rle> |            <character>  <integer>\n",
       "  [1]       17 987221-73296953      * | TCGA-3C-AAAU-10A-01D..      33859\n",
       "  [2]       17 987221-73296953      * | TCGA-3C-AAAU-10A-01D..      33859\n",
       "      Segment_Mean\n",
       "         <numeric>\n",
       "  [1]       0.0088\n",
       "  [2]       0.0088\n",
       "  -------\n",
       "  seqinfo: 23 sequences from an unspecified genome; no seqlengths"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "segs.overlap <- find_overlaps(segs.gr, tiles.subset)  # arguments: find_overlaps(query, subject)\n",
    "segs.overlap[1:2] # show first 2 segments from segs.gr that overlapped tiles.subset"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a6564984-bbf4-4cc7-a69e-f8708efb7f93",
   "metadata": {},
   "source": [
    "## Exercise 1: \n",
    "\n",
    "### a. Create a range for `11:69400000-69500000`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "e8dff18b-fab5-4ade-9b2f-e9de1fd45a07",
   "metadata": {},
   "outputs": [],
   "source": [
    "# GRanges()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6ea3b2fc-276a-4718-bb24-917d20a1ac93",
   "metadata": {},
   "source": [
    "### b. Find overlap between `11:69400000-69500000` and `segs.gr`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "25beff24-9abc-4e15-a000-5db6f7a32de6",
   "metadata": {},
   "outputs": [],
   "source": [
    "# find_overlaps() from plyranges"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6fde5c7e-c2ae-48ba-b2f9-daf7297ca0c9",
   "metadata": {},
   "source": [
    "### c. What is the `Segment_Mean` for the 2nd segment that overlaps `11:69400000-69500000`?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "id": "f76a6a13-5142-4500-8fbf-c1a55a1ce335",
   "metadata": {},
   "outputs": [],
   "source": [
    "# index the 2nd segment in the result to (b)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "R 4.1.2",
   "language": "R",
   "name": "ir412"
  },
  "language_info": {
   "codemirror_mode": "r",
   "file_extension": ".r",
   "mimetype": "text/x-r-source",
   "name": "R",
   "pygments_lexer": "r",
   "version": "4.1.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}