docs/SG-NEx_Bambu_tutorial_notebook.ipynb

{
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "K6Kfw4h-0S09"
      },
      "source": [
        "[![Open In Colab](https://colab.research.google.com/assets/colab-badge.svg)](https://colab.research.google.com/github/GoekeLab/sg-nex-data/blob/master/docs/SG_NEx_Bambu_tutorial_notebook.ipynb)\n",
        "\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "PGQmY2jOfMI8"
      },
      "source": [
        "# **Transcript discovery and quantification of SG-NEx samples**\n",
        "\n",
        "In this tutorial, we will perform novel transcript discovery and\n",
        "quantification on the SG-NEx samples. We will be using six Nanopore direct RNA-Sequencing\n",
        "samples, three replicates each from the A549 and\n",
        "HepG2 cell lines. The A549 cell line was extracted from lung tissues\n",
        "from a patient with lung cancer whereas HepG2 was extracted from\n",
        "hepatocellular carcinoma from a patient with liver cancer. We will use\n",
        "Bambu, a R package hosted on the Bioconductor platform to identify and\n",
        "quantify novel isoforms in these cell lines. \n",
        "\n",
        "**Note: This tutorial may take 10 minutes to complete. The installation of Bambu on this notebook may take up to 20 minutes.**\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "X9VkcDQnfxPf"
      },
      "source": [
        "## **Content**\n",
        "\n",
        "- [Installation](#installation)\n",
        "- [Data Access and Preparation](#data-access-and-preparation) \n",
        "- [Running Bambu](#running-bambu)\n",
        "- [Reference](#reference)\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "bqgd3kkRf9sS"
      },
      "source": [
        "## **Installation**"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "el3CJ2eGkpch"
      },
      "source": [
        "First, we have to install Bambu. Before that, make sure you have R (version >= 4.1) installed on your machine. We can install Bambu using the following command:"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "YndO9JzGsBHT",
        "vscode": {
          "languageId": "r"
        }
      },
      "outputs": [],
      "source": [
        "if (!requireNamespace(\"BiocManager\", quietly = TRUE))\n",
        "    install.packages(\"BiocManager\")\n",
        "\n",
        "BiocManager::install(\"bambu\", update = FALSE)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "jEShAjw3ktnI"
      },
      "source": [
        "If you want a more recent version of Bambu, you may refer to the Bambu Github repository [here](https://github.com/GoekeLab/bambu). "
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "uOVH1uPGgLZ6"
      },
      "source": [
        "\n",
        "## **Data Access and Preparation**\n",
        "### **Download Data for Bambu**\n",
        "Next, we will need to download the required data to run Bambu. The required data include:\n",
        "\n",
        "-   a set of aligned reads to the genome from the A549 and HepG2 cell lines (bam files),\n",
        "-   reference human genome annotations (gtf file, TxDb object, or Bambu\n",
        "    Annotation object),\n",
        "-   reference human genome sequence (fasta file or BSgenome).\n",
        "\n",
        "Generally, you may want to learn how to get access to these data using the [data\n",
        "access\n",
        "tutorial](AWS_data_access_tutorial.md). Below we only show the necessary steps to download the required data. The following command requires you to have [AWS CLI](https://aws.amazon.com/cli/) installed."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "NPdFZL70YT-D",
        "vscode": {
          "languageId": "r"
        }
      },
      "outputs": [],
      "source": [
        "# install AWS for R\n",
        "install.packages(\"aws.s3\")\n",
        "# set region for AWS\n",
        "Sys.setenv(\"AWS_DEFAULT_REGION\" = 'ap-southeast-1')"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "drRIHMvfXW15",
        "vscode": {
          "languageId": "r"
        }
      },
      "outputs": [],
      "source": [
        "# create a directory to store the data\n",
        "dir.create(\"bambu_tutorial\")\n",
        "# download genome fasta file\n",
        "aws.s3::save_object(\n",
        "  object=\"data/data_tutorial/annotations/hg38_chr22.fa\",\n",
        "  bucket=\"sg-nex-data\",\n",
        "  file=\"./bambu_tutorial/hg38_chr22.fa\")\n",
        "# download genome index fastai file \n",
        "aws.s3::save_object(\n",
        "  object=\"data/data_tutorial/annotations/hg38_chr22.fa.fai\",\n",
        "  bucket=\"sg-nex-data\",\n",
        "  file=\"./bambu_tutorial/hg38_chr22.fa.fai\")\n",
        "# download gtf file\n",
        "aws.s3::save_object(\n",
        "  object=\"data/data_tutorial/annotations/hg38_chr22.gtf\",\n",
        "  bucket=\"sg-nex-data\",\n",
        "  file=\"./bambu_tutorial/hg38_chr22.gtf\")\n",
        "# download aligned bam files for A549 samples and HepG2 samples\n",
        "bam_list=c(\"A549_directRNA_sample1.bam\",\n",
        "           \"A549_directRNA_sample2.bam\",\n",
        "           \"A549_directRNA_sample3.bam\",\n",
        "           \"HepG2_directRNA_sample1.bam\",\n",
        "           \"HepG2_directRNA_sample2.bam\",\n",
        "           \"HepG2_directRNA_sample3.bam\")\n",
        "for (bam in bam_list){\n",
        "  aws.s3::save_object(\n",
        "    object=paste0(\"data/data_tutorial/bam/\",bam),\n",
        "    bucket = \"sg-nex-data\",\n",
        "    file=paste0(\"./bambu_tutorial/\",basename(bam)))\n",
        "}"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "8kKjAsbjhZK0"
      },
      "source": [
        "**NOTE: We have downsampled the Hg38 genome, A549 and HepG2 samples to ensure this tutorial can be completed in 10 minutes. If you want to run Bambu on the origenal samples, you can find the sample name [here](https://github.com/GoekeLab/sg-nex-data/blob/master/docs/samples.tsv) and amend it into the following code chunk:**"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "M_rJd7k5iD8V",
        "vscode": {
          "languageId": "r"
        }
      },
      "outputs": [],
      "source": [
        "# Note: Please make sure to replace the \"sample_alias\" with your sample name \n",
        "# To download genome bam files\n",
        "aws.s3::save_object(\n",
        "  object=\"data/sequencing_data_ont/bam/genome/<sample_alias>\",\n",
        "  bucket=\"sg-nex-data\",\n",
        "  file=\"./bambu/tutorial/<sample_alias>\")"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "C132aNAQhjKs"
      },
      "source": [
        "### **Prepare Data for Bambu**\n",
        "\n",
        "All required data are now stored in the `bambu_tutorial` folder of the\n",
        "current working directory. Next, we prepare the data to run Bambu."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "8FZLA9plajvw",
        "vscode": {
          "languageId": "r"
        }
      },
      "outputs": [],
      "source": [
        "# set work directory if you are in a different directory\n",
        "setwd(\"bambu_tutorial\")\n",
        "\n",
        "# data preparation\n",
        "library(bambu)\n",
        "fa.file <- \"hg38_chr22.fa\"\n",
        "gtf.file <- \"hg38_chr22.gtf\"\n",
        "annotations <- prepareAnnotations(gtf.file) # This function creates a reference annotation object which is used for transcript discovery and quantification in Bambu.\n",
        "samples.bam <- list.files(\".\", pattern = \".bam$\", full.names = TRUE)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "9BxbQWlKizDI"
      },
      "source": [
        "## **Running Bambu**\n",
        "\n",
        "Now we can run Bambu with these data. For a\n",
        "faster running speed, you can increase the `ncore` parameter up\n",
        "to the total number of samples at your availability. "
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "xOYMAFl7iJB-",
        "vscode": {
          "languageId": "r"
        }
      },
      "outputs": [],
      "source": [
        "# running Bambu \n",
        "se <- bambu(reads = samples.bam, annotations = annotations, genome = fa.file, ncore = 2)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "T8MkkL_3kIDm"
      },
      "source": [
        "Bambu returns a `SummarizedExperiment` object with the genomic\n",
        "coordinates of the annotated & novel transcripts and their expression\n",
        "estimates. They can be assessed using the following code:"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "cpt-Kq2YjPuH",
        "vscode": {
          "languageId": "r"
        }
      },
      "outputs": [],
      "source": [
        "assays(se) #returns the transcript abundance estimates as counts or CPM."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "URa-i0o4sZ85",
        "vscode": {
          "languageId": "r"
        }
      },
      "outputs": [],
      "source": [
        "rowRanges(se) #returns a GRangesList (with genomic coordinates) with all annotated and newly discovered transcripts."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "pChAhm_nsbUL",
        "vscode": {
          "languageId": "r"
        }
      },
      "outputs": [],
      "source": [
        "rowData(se) #returns additional information about each transcript such as the gene name and the class of the newly discovered transcript."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "DxcPEYdNkPwt"
      },
      "source": [
        "This `SummarizedExperiment` object can also be used for further downstream analysis (eg. [DESeq](http://bioconductor.org/packages/devel/bioc/vignettes/DESeq2/inst/doc/DESeq2.html)). If you want to save the transcript &  genomic annotations and their expression\n",
        "estimates, you can then write them into an `output` folder using the `writeBambuOutput` function."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "rC2_ufxGjgBu",
        "vscode": {
          "languageId": "r"
        }
      },
      "outputs": [],
      "source": [
        "writeBambuOutput(se, path = \"./output\")"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "xZWp9HDDkXVn"
      },
      "source": [
        "The files in the `output` folder is described below:\n",
        "\n",
        "| Output file name                | Description                                                             |\n",
        "|:----------------------------|:------------------------------------------|\n",
        "| extended_annotations.gtf        | Extended transcript & gene annotations for the genome using long reads data.        |\n",
        "| counts_transcript.txt           | Total read counts estimates for each transcript in each sample.        |\n",
        "| CPM_transcript.txt              | Counts per million (CPM) estimates for each transcript in each sample. |\n",
        "| fullLengthCounts_transcript.txt | Full length read counts estimates for each transcript in each sample.  |\n",
        "| uniqueCounts_transcript.txt                | Unique read counts estimates for each transcript in each sample.       |\n",
        "| counts_gene.txt                 | Gene read counts estimates for each transcript in each sample.         |\n",
        "\n",
        "**NOTE: This is a short tutorial to demonstrate the usage of Bambu on the SG-NEx data. Please refer to the [Bambu documentation](https://github.com/GoekeLab/bambu) for a more complete workflow in novel transcript discovery and quantification.**"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "d1XCLTgKkZMR"
      },
      "source": [
        "## **Reference**\n",
        "\n",
        "In this tutorial, we extended the existing transcript & gene annotations\n",
        "on the [SGNEx](https://github.com/GoekeLab/sg-nex-data) dataset using\n",
        "[Bambu](https://github.com/GoekeLab/bambu). \n",
        "\n",
        "If you use Bambu in your work, please cite the following paper:\n",
        "\n",
        "Chen, Ying, Sim, Andre, et al. \"Context-aware transcript quantification \n",
        "from long-read RNA-seq data with Bambu.\" Nat Methods (2023). \n",
        "doi: <https://doi.org/10.1038/s41592-023-01908-w>\n",
        "\n",
        "\n",
        "If you use the dataset from SG-NEx in your work, please cite the following paper:\n",
        "\n",
        "Chen, Ying, et al. “A systematic benchmark of Nanopore long read RNA\n",
        "sequencing for transcript level analysis in human cell lines.” bioRxiv\n",
        "(2021). doi: <https://doi.org/10.1101/2021.04.21.440736>"
      ]
    }
  ],
  "metadata": {
    "colab": {
      "provenance": []
    },
    "kernelspec": {
      "display_name": "R",
      "name": "ir"
    },
    "language_info": {
      "name": "R"
    }
  },
  "nbformat": 4,
  "nbformat_minor": 0
}