lectures/lecture16/lecture16.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Working with Sequence Data using `R` and `Bioconductor`"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`Bioconductor` is a collection ([> 2000](https://www.bioconductor.org/packages/release/BiocViews.html#___Software)) of `R` packages focused on biological data analysis.\n",
    "\n",
    "**Packages we will use today**\n",
    "\n",
    "\n",
    "- [`Biostrings`](https://bioconductor.org/packages/release/bioc/html/Biostrings.html) for working with DNA, RNA, and protein sequences\n",
    "- [`rtracklayer`](https://bioconductor.org/packages/release/bioc/html/rtracklayer.html) for working with sequence annotation files (GFF, BED)\n",
    "- [`plyranges`](https://bioconductor.org/packages/release/bioc/html/plyranges.html) for a `tidyverse`-like interface to Bioconductor data objects. \n",
    "- [`BSgenome.Hsapiens.UCSC.hg38`](https://bioconductor.org/packages/release/data/annotation/html/BSgenome.Hsapiens.UCSC.hg38.html) for human genome sequence data\n",
    "- [`GenomicFeatures`](https://bioconductor.org/packages/release/bioc/html/GenomicFeatures.html) for working with transcript annotations\n",
    "- [`ggbio`](https://www.bioconductor.org/packages/release/bioc/html/ggbio.html) GGplot2-like plotting interface for biological data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Package Installation\n",
    "\n",
    "For installing `Bioconductor` packages using `conda`, I search for the package name + `conda` for finding the relevant `conda` package name and repository.\n",
    "\n",
    "For example, the first link [here](https://www.google.com/search?q=biostrings+conda) is the `conda` package for `biostrings` which I would install as:\n",
    "\n",
    "```bash\n",
    "conda activate tidy # tidy is name of R environment\n",
    "conda install -c bioconda bioconductor-biostrings\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## How to learn further\n",
    "\n",
    "- `Documentation` section of package homepage on `Bioconductor` website, eg. [`plyranges`](https://bioconductor.org/packages/release/bioc/html/plyranges.html)\n",
    "- Bioconductor workshops, eg. [`plyranges`](https://bioconductor.github.io/BiocWorkshops/fluent-genomic-data-analysis-with-plyranges.htmlhttps://bioconductor.org/workshops/plyranges/)\n",
    "- [Q & A](https://support.bioconductor.org/) section of `Bioconductor` website"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Load packages\n",
    "\n",
    "Note: Many functions are named the same across different packages. For example, `select`, `n`, `rename` etc. This can lead to confusing errors.\n",
    "\n",
    "If you are getting errors when you use a function, try loading the packages in a different order or specify the package name explicitly. For example, use `dplyr::rename()`, `plyranges::select()`, etc."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "library(Biostrings)\n",
    "library(rtracklayer)\n",
    "library(plyranges)\n",
    "library(tidyverse)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Working with sequence data using `Biostrings`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "seqs <- readDNAStringSet(\"data/tumor_suppresssors.fasta\") %>%\n",
    "  print()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Reverse complement"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "seqs %>%\n",
    "  reverseComplement()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Translate"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "seqs %>%\n",
    "  translate()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Get a sub-sequence"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "seqs %>%\n",
    "  subseq(1,10)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Calculate codon frequencies"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "seqs %>%\n",
    "  oligonucleotideFrequency(width = 3, step = 3) %>%\n",
    "  as_tibble()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You can also read `fastq` files but not recommended for large files\n",
    "\n",
    "See `ShortRead` package for working with large number of sequence reads"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "readDNAStringSet(\"../lecture08/barcodes_R1.fastq\", format = \"fastq\") %>%\n",
    "  # reverseComplement() %>%\n",
    "  # # convert to a tibble\n",
    "  # as.data.frame() %>%\n",
    "  # rownames_to_column() %>%\n",
    "  # setNames(c(\"name\", \"seq\")) %>%\n",
    "  # # ----------------\n",
    "  # mutate(is_HA = str_detect(seq, \"CCGGATTTGCATATAATGATGCACCAT\")) %>%\n",
    "  # mutate(is_NA = str_detect(seq, \"CACGATAGATAAATAATAGTGCACCAT\")) %>%\n",
    "  # mutate(barcode = substr(seq, nchar(seq) - 16, nchar(seq))) %>%\n",
    "  # # select(-seq)\n",
    "  # group_by(is_HA, is_NA, barcode) %>%\n",
    "  # summarize(n_HA = sum(is_HA), n_NA = sum(is_NA)) %>%\n",
    "  # filter(n_NA > 0) %>%\n",
    "  # arrange(-n_NA) %>%\n",
    "  print()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# How do we get sequence data for a specific gene or transcript?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "genome <- BSgenome.Hsapiens.UCSC.hg38::BSgenome.Hsapiens.UCSC.hg38"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Chromosome names"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "seqnames(genome)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Chromosome sequence"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "chr2 <- genome$chr2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Read annotations for a group of genes\n",
    "\n",
    "These are a subset of the full genome annotations in the [GENCODE](https://www.gencodegenes.org/) project."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "annotations <- import.gff3(\"data/tumor_suppresssors.gff3\") %>%\n",
    "  print()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can convert from `GRanges` to  `tibble` and vice versa."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "annotations %>%\n",
    "  as_tibble() %>%\n",
    "  # distinct(type) %>%\n",
    "  distinct(gene_name) %>%\n",
    "  # GRanges()\n",
    "  print()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Group exons into transcripts\n",
    "\n",
    "We use `split` function with `$` to refer to columns.\n",
    "\n",
    "Each group is list of `GRanges` or a `GRangesList` object.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "transcripts <- annotations %>%\n",
    "  filter(type == \"exon\") %>%\n",
    "  split(.$transcript_id) \n",
    "  # as_tibble() \n",
    "\n",
    "transcripts"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Visualize transcript models for *TP53* gene"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "options(repr.plot.width = 15, repr.plot.height = 10)\n",
    "annotations %>%\n",
    "  filter(type == \"exon\") %>%\n",
    "  filter(gene_name == \"TP53\") %>%\n",
    "  split(.$transcript_id) %>%\n",
    "  ggbio::autoplot() +\n",
    "  theme_minimal()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## How do we extract transcript sequences?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "GenomicFeatures::extractTranscriptSeqs(genome, transcripts)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# In-class exercises\n",
    "\n",
    "(20 minutes)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1. Extract the coding sequence for all the transcripts in the annotation file above. \n",
    "\n",
    "The steps are identical to above except that you need to work with `CDS` features instead of `exon` features.\n",
    "\n",
    "If everything worked correctly, you should at least have the expected start and stop codons at the termini of each coding sequence."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2. Find the transcripts with the longest coding sequence for each gene\n",
    "\n",
    "**Steps**\n",
    "\n",
    "- Convert annotations to a `tibble`\n",
    "- Filter for `CDS` features\n",
    "- Group by gene name and transcript id\n",
    "- Sum the `width` column for each group\n",
    "- Group by gene name and filter for the transcript with the longest `width` column"
   ]
  }
 ],
 "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": "4.0.3"
  },
  "orig_nbformat": 4
 },
 "nbformat": 4,
 "nbformat_minor": 2
}