Tutorial_4_Snakemake_LS_API.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "99cfa998-06b6-4b42-ae3a-b4e011750d31",
   "metadata": {},
   "source": [
    "# RNA-Seq Analysis using Snakemake and Google Cloud Life Sciences API"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4126cb07-34ee-4780-838f-872015a882b3",
   "metadata": {},
   "source": [
    "## Overview"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f15ea992-faa6-4705-8384-eb5d81f5daff",
   "metadata": {},
   "source": [
    "This short tutorial demonstrates how to run an RNA-Seq workflow using a prokaryotic data set. Steps in the workflow include read trimming, read QC, read mapping, and counting mapped reads per gene to quantitate gene expression. This tutorial uses a popular workflow manager called [Snakemake](https://snakemake.readthedocs.io/en/stable/) run via the [Google Cloud Life Sciences API](https://cloud.google.com/life-sciences/docs/reference/rest). If you completed the other tutorials in this repo, you will see that much of it is a repeat of Tutorial 2, but instead of running Snakemake locally, it uses the Life Sciences API to run in a serverless way. You will need to repeat most of the steps here because this notebook copies the data to a bucket instead of locally.\n",
    "**Make sure you enable the Life Science API before running this tutorial**."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0b6d0785-2d13-476c-b16a-196f74ea277d",
   "metadata": {},
   "source": [
    "### Step 1: Create a new GS Bucket to store input and output files\n",
    "Note that your bucket has to be globally unique, so make sure you don't just copy the example here or it won't work"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4d3dc88f-fa0c-4e7e-972b-055321d3cdbb",
   "metadata": {},
   "outputs": [],
   "source": [
    "#change this bucket name\n",
    "%env BUCKET=cl-testing-bucket"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "17ce680c-4b8c-419c-a6c4-b6caec32d9ba",
   "metadata": {},
   "outputs": [],
   "source": [
    "#will only create the bucket if it doesn't yet exist\n",
    "!gsutil ls gs://$BUCKET >& /dev/null || gsutil mb gs://$BUCKET"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "91292c6d-d5a4-407d-9816-51ca52876fba",
   "metadata": {},
   "outputs": [],
   "source": [
    "#set versioning on the bucket so it can overwrite old files\n",
    "!gsutil versioning set on gs://$BUCKET"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "dd7ab630-955d-43d1-bc43-c7b3e701ed04",
   "metadata": {},
   "source": [
    "### STEP 2: Install mambaforge and snakemake\n",
    "First install mambaforge, then use mamba to install snakemake. Skip this if you have completed the other tutorials."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "682ddf88-e1d9-443f-a423-e1f85ff604a2",
   "metadata": {},
   "outputs": [],
   "source": [
    "! curl -L -O https://github.com/conda-forge/miniforge/releases/latest/download/Mambaforge-$(uname)-$(uname -m).sh\n",
    "! bash Mambaforge-$(uname)-$(uname -m).sh -b -p $HOME/mambaforge"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "de5d0630-1d85-4625-bc04-036aae11ce4a",
   "metadata": {},
   "outputs": [],
   "source": [
    "#add to your path\n",
    "import os\n",
    "os.environ[\"PATH\"] += os.pathsep + os.environ[\"HOME\"]+\"/Mambaforge/bin\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bd5c221b-45ce-47fb-a8e2-29ceee0e296a",
   "metadata": {},
   "outputs": [],
   "source": [
    "#install snakemake\n",
    "! mamba install -y -c conda-forge -c bioconda snakemake"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bb0ce8d5-4b96-4e97-88ed-44e8e85f4fc0",
   "metadata": {},
   "source": [
    "### STEP 3: Copy FASTQ Files\n",
    "In order for this tutorial to run quickly, we will only analyze 50,000 reads from a sample from both sample groups instead of analyzing all the reads from all six samples. These files have been posted on a Google Storage Bucket that we made publicly accessible."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5f3795fd-3e03-476d-9abf-49705a72cc15",
   "metadata": {},
   "outputs": [],
   "source": [
    "! gsutil -m cp -r gs://nigms-sandbox/me-inbre-rnaseq-pipelinev2/data/raw_fastqSub/ gs://$BUCKET/data/raw_fastq"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a2bc4563-951a-45d4-8f01-0accd6b80ea8",
   "metadata": {},
   "source": [
    "Create a fake path to data/fastqc so that snakemake can write files to that bucket path, otherwise the pipeline crashes."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cd2c03dd-2248-4068-8842-ba130f29adc9",
   "metadata": {},
   "outputs": [],
   "source": [
    "! touch blank.txt\n",
    "! gsutil cp blank.txt gs://$BUCKET/data/fastqc/\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ec692c28-f549-43af-bbdf-3c4266fb59ae",
   "metadata": {},
   "source": [
    "### STEP 4: Copy reference files that will be used by Salmon\n",
    "Salmon is a tool that aligns RNA-Seq reads to a set of transcripts rather than the entire genome."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "290238d6-39e0-4575-87e4-880b316ca1f6",
   "metadata": {},
   "outputs": [],
   "source": [
    "!gsutil -m cp -r gs://nigms-sandbox/me-inbre-rnaseq-pipelinev2/data/reference/ gs://$BUCKET/data/\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ea2d17cb-dff6-45d3-9aef-3ec6203508f6",
   "metadata": {},
   "source": [
    "### STEP 5: Copy data file for Trimmomatic"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0be371d5-382e-4a22-a300-2c5249eff825",
   "metadata": {},
   "outputs": [],
   "source": [
    "!gsutil -m cp gs://nigms-sandbox/me-inbre-rnaseq-pipelinev2/config/TruSeq3-PE.fa gs://$BUCKET/TruSeq3-PE.fa"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5ac668db-7851-418e-9b1e-0a2c4abbab6e",
   "metadata": {},
   "source": [
    "### STEP 6: Copy data and config files that will be used in our snakemake environment\n",
    "\n",
    "Next download config files for our snakemake environment."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1dbc460c-50af-4458-8056-c0f6146fff23",
   "metadata": {},
   "outputs": [],
   "source": [
    "!gsutil -m cp -r gs://nigms-sandbox/me-inbre-rnaseq-pipelinev2/envs/ .\n",
    "!gsutil -m cp gs://nigms-sandbox/me-inbre-rnaseq-pipelinev2/config.yaml .\n",
    "!gsutil -m cp gs://nigms-sandbox/me-inbre-rnaseq-pipelinev2/snakefile_ls_api ."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ab3739b6-aa4e-439e-bd73-2ea43be1801b",
   "metadata": {},
   "source": [
    "Add the bucket path to the end of your config file. Since this file was written for running snakemake locally we have to make a few edits to run on the LS API."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "737c0634-171e-489b-8cfb-e93a025cbd01",
   "metadata": {},
   "outputs": [],
   "source": [
    "!echo 'bucket:' >> config.yaml"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4fd4c713-8882-4833-a117-706a4b239374",
   "metadata": {},
   "outputs": [],
   "source": [
    "!echo '   '$BUCKET >>config.yaml"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e45fe428-b9fa-46cb-a69d-2a0e989292e1",
   "metadata": {},
   "source": [
    "Add bucket path to the snakefile"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3bba277a-704a-41ab-8853-7cf324dde727",
   "metadata": {},
   "outputs": [],
   "source": [
    "!sed -i 's/print(SAMPLES)/BUCKET=config[\"bucket\"]/' snakefile_ls_api"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cfc635b4-ed50-4a92-8018-f402fdd563b9",
   "metadata": {},
   "source": [
    "### Step 7: Set up your local environment\n",
    "You need to generate a [service account key](https://cloud.google.com/iam/docs/creating-managing-service-account-keys) for the compute engine default service account to interact with the Life Sciences API using Snakemake. Download the key and copy it to this VM. Then assign the path of the json file to an environment variable."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0e93d6ca-5afe-4f12-9401-be44ce7ac7d9",
   "metadata": {},
   "outputs": [],
   "source": [
    "%env GOOGLE_APPLICATION_CREDENTIALS=cloud_creds.json"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e629d7a5-4ef7-408b-88ff-3c52201879a1",
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "os.environ[\"GOOGLE_APPLICATION_CREDENTIALS\"]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a04afc6c-41cc-46c9-b048-08a7677699c8",
   "metadata": {},
   "source": [
    "Set your project (make sure to replace $Project with your project name)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "10f587c6-9d51-4433-bf13-bec6dabcb9d8",
   "metadata": {},
   "outputs": [],
   "source": [
    "!gcloud config set project $PROJECT"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f7d3b561-7213-4b60-9f8f-c511c6fdc067",
   "metadata": {},
   "source": [
    "Initialize a local git repo"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6d56605b-04bd-48fa-94c0-ba7fffdf08a2",
   "metadata": {},
   "outputs": [],
   "source": [
    "!git init"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f56458ef-e475-4a48-b24b-e903f82eb996",
   "metadata": {},
   "source": [
    "Configure conda"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f72b8fce-be20-4158-8bd5-00b2ad122414",
   "metadata": {},
   "outputs": [],
   "source": [
    "conda config --set channel_priority strict"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ec2c0273-c7f1-4aee-bdf3-43d5773cf2fa",
   "metadata": {},
   "source": [
    "### STEP 8: Run snakemake using the Life Sciences API\n",
    "\n",
    "Aside from the .yaml config files which information about software, dependencies, and versions -- snakemake uses a snakefile which contains information about a workflow.\n",
    "\n",
    "This can be a powerful tool as it allows one to operate and think in terms of workflows instead of individual steps. You should open the snakefile to look at it further. It is composed of 'rules' we have created. Snakefiles work largely based on inputs. For a given input/output, there is an associated 'rule' that runs. Snakefiles may take a while to get the idea of what's going on, but in simplest terms here we take an input of .fastq files, and based on the snakefile rules we created, those fastq files are run through the entire workflow. The rule_all at the top determines which rules are run based on the input files for rule_all (which are outputs from the target rules. Comment out rules you don't want to run. "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d8bf3f71-e394-41c6-9694-b5d4b24cb265",
   "metadata": {},
   "source": [
    "Snakemake requires that you have a service account key to authenticate with the Life Sciences API. This actually is not necessary to use the API from within a notebook, but Snakemake does require it since Snakemake is expecting you to run the command from your own terminal using the SDK. To see all the commands you can run with Snakemake via the Life Sciences API, check out the [docs](https://snakemake.readthedocs.io/en/stable/executor_tutorial/google_lifesciences.html)."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "67c95e75-1619-4694-9411-95edc9f4cee4",
   "metadata": {},
   "source": [
    "Now we can run the Life Sciences APi. You will see that each rule is submitted as a separate job. If the pipeline crashes, the way to troubleshoot is by reading the API logs, or the snakemake rule logs (same info). You can find the Life Sciences API logs by pasting in the gcloud command given in yellow.\n",
    "\n",
    "For example: \n",
    "```\n",
    "gcloud beta lifesciences operations describe <JOB_ID>\n",
    "```\n",
    "Or you can view the logs by finding the path given for logs, and then use gsutil to copy that file locally, or go to the bucket and double click the file. You can get the job ID for the output file in the green section of the rule print out."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bee32318-33df-43b2-98bc-5eb091ceae59",
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "! snakemake --forceall --snakefile snakefile_ls_api --google-lifesciences --default-remote-prefix $BUCKET --use-conda --google-lifesciences-region us-central1 -j 24 --rerun-incomplete --default-resources \"machine_type=n2-standard\""
   ]
  },
  {
   "cell_type": "markdown",
   "id": "df9deb0a-1030-4839-aa16-37c3b32a2c87",
   "metadata": {},
   "source": [
    "### STEP 9: Report the top 10 most highly expressed genes in the samples."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d50f9bd2-dbd2-467f-a9b6-313e63ad304b",
   "metadata": {},
   "source": [
    "Top 10 most highly expressed genes in the wild-type sample. The level of expression is reported in the Transcripts Per Million (`TPM`) and number of reads (`NumReads`) fields:  \n",
    "`Name    Length  EffectiveLength TPM     NumReads`\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0676cbbc-9392-41d1-ab57-e2b4f3cc9aad",
   "metadata": {},
   "outputs": [],
   "source": [
    "!gsutil rm gs://$BUCKET/data/quants/SRR13349122_quant\n",
    "!gsutil rm gs://$BUCKET/data/quants/SRR13349128_quant"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b032ce69-f62d-4f5f-90a3-68c2979d9a85",
   "metadata": {},
   "outputs": [],
   "source": [
    "!gsutil ls gs://$BUCKET/data/quants/*"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c98fd827-6829-400d-af8c-969ad196c3d2",
   "metadata": {},
   "outputs": [],
   "source": [
    "!gsutil cp -r gs://$BUCKET/data/quants/SRR13349122_quant/ .\n",
    "!gsutil cp -r gs://$BUCKET/data/quants/SRR13349128_quant/ ."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7776f671-30a0-4ba8-a9cc-e3434d40cc48",
   "metadata": {},
   "outputs": [],
   "source": [
    "!sort -nrk 5,5 SRR13349122_quant/quant.sf | head -10"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "678efdde-1782-4481-9240-054c34528163",
   "metadata": {},
   "source": [
    "Top 10 most highly expressed genes in the double lysogen sample.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1ceee200-b741-4954-b950-85edec98eb90",
   "metadata": {},
   "outputs": [],
   "source": [
    "!sort -nrk 5,5 SRR13349128_quant/quant.sf | head -10"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "50169f62-e707-4d84-b301-ded51a704130",
   "metadata": {},
   "source": [
    "### STEP 10: Report the expression of a putative acyl-ACP desaturase (BB28_RS16545) that was downregulated in the double lysogen relative to wild-type\n",
    "A acyl-transferase was reported to be downregulated in the double lysogen as shown in the table of the top 20 upregulated and downregulated genes from the paper describing the study.\n",
    "![RNA-Seq workflow](images/table-cushman.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5b3794b0-a477-45fa-aa51-4414d7671441",
   "metadata": {},
   "source": [
    "Use `grep` to report the expression in the wild-type sample. The fields in the Salmon `quant.sf` file are as follows. The level of expression is reported in the Transcripts Per Million (`TPM`) and number of reads (`NumReads`) fields:  \n",
    "`Name    Length  EffectiveLength TPM     NumReads`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a3cb9340-682b-4177-837d-7d803a9775a5",
   "metadata": {},
   "outputs": [],
   "source": [
    "!grep 'BB28_RS16545' SRR13349122_quant/quant.sf"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "79ba6401-261d-43e9-b831-ef76122da623",
   "metadata": {},
   "source": [
    "Use `grep` to report the expression in the double lysogen sample. The fields in the Salmon `quant.sf` file are as follows. The level of expression is reported in the Transcripts Per Million (`TPM`) and number of reads (`NumReads`) fields:  \n",
    "`Name    Length  EffectiveLength TPM     NumReads`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "745ea1c5-79d3-481c-9359-6e0a93b9a286",
   "metadata": {},
   "outputs": [],
   "source": [
    "!grep 'BB28_RS16545' SRR13349128_quant/quant.sf"
   ]
  }
 ],
 "metadata": {},
 "nbformat": 4,
 "nbformat_minor": 5
}