tutorial_4-methylseq2.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "bbe8280a-cc0b-41b1-a18b-3f43131ddb58",
   "metadata": {},
   "source": [
    "# The WGBS data analysis tutorial 4 - \n",
    "# run nf-core/methylseq using Google Batch"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "80fcb1bb-572d-41d2-906e-c06813203765",
   "metadata": {},
   "source": [
    "For real-world datasets, the sequence file sizes are usually too large to process using a single virtual machine (Vertex AI notebook), or take a long time. In this tutorial, we will show how to run a nf-core/methyseq pipeline to process WGBS data using the Google Batch.  \n",
    "\n",
    "<img src=\"images/notebook4_2.png\" width=\"900\" />\n",
    "\n",
    "The [**Google Batch**](https://www.nextflow.io/docs/latest/google.html#cloud-batch) is a managed computing service that allows the execution of containerized workloads in the Google Cloud Platform infrastructure. It provides a simple way to execute a series of Compute Engine containers on Google Cloud. The most common use case when using Google Batch is to run an existing tool or custom script that reads and writes files, typically to and from Cloud Storage. Nextflow provides built-in support for Google Batch, which allows the seamless deployment of a Nextflow pipeline in the cloud, offloading the process executions through the Google Cloud service. Google Batch can run independently over hundreds or thousands of these files. \n",
    "\n",
    "\n",
    "There were several steps before we can submit a job to Google Batch through Nextflow. And these steps (listed below) will be covered in this tutorial:\n",
    "\n",
    "- [Create a Nexflow Service Account](#CNSA) -- if does not exist\n",
    "- [Create a Notebook with Service Account Permissions](#Create-a-Notebook-with-Service-Account-Permissions)\n",
    "\n",
    "Now, you can open a new Vertex AI notebook with a Nextflow service account, and:\n",
    "\n",
    "- [Install Nextflow, and Create a Config File for Google Batch](#INCC) \n",
    "- [Download and Test nf-core/methylseq using Google Batch](#TEST)\n",
    "- [An Example of a Real-World Dataset](#REAL)\n",
    "- [Configuration of a Full-scale Dataset](#FULL) -- Troubleshooting\n",
    "\n",
    "## Create a Nextflow Service Account<a name=\"CNSA\"></a>\n",
    "**Most of what is in this section is unnecessary if you are using an NIH Cloud Lab Project/Account**. If you are a Cloud Lab user, feel free to skip ahead to `Install Nextflow, and create a config file for Google Batch` without creating a Service Account. Everything should run fine without the Service Account. You will still need to enable the APIs.\n",
    "`\n",
    "\n",
    "Before creating a new service account, please check if there is already a Nextflow service account available (`Menu` > `IAM & Admin` > `Service Accounts`). There is no need to create a new one if there is one that already exists. If not, follow the steps below to create one.\n",
    "\n",
    "#### Enable APIs  \n",
    "- Enable the Batch, Compute Engine, Cloud Logging, and Cloud Storage APIs by searching each of the GCP products and clicking **`ENABLE`** button (for the whole project, should have already been done in the beginning README.md)\n",
    "\n",
    "#### Create a Nextflow service account  \n",
    "- Click the main navigation menu, go to **IAM & Admin** click **Service Accounts**\n",
    "- Select **+ CREATE SERVICE ACOUNT**\n",
    "- Type in 'nextflow-service-account' as the service account name and press **`DONE`**\n",
    "\n",
    "<img src=\"images/4_create_service_account.png\" width=\"800\">\n",
    "\n",
    "#### Add roles to the service account:  \n",
    "- On the **IAM & Admin** menu click **IAM** then click `edit` next to the Nextflow service account just created\n",
    "- Add the following roles and click **`SAVE`**:  \n",
    "    - Service Account User\n",
    "    - Batch Agent Reporter\n",
    "    - Storage Admin\n",
    "    - Storage Object Admin\n",
    "    - Batch Job Editor\n",
    "\n",
    "<img src=\"images/4_create_service_account_roles.png\" width=\"800\"> \n",
    "\n",
    "## Create a Notebook with Service Account Permissions\n",
    "\n",
    "When creating a notebook you can edit the permissions to utilize the Nextflow service account.  \n",
    "- Using the 'IAM & Admin' menu on the left click 'Service Accounts' (if you aren't there already) locate your Nextflow service account and copy the entire email name\n",
    "- Edit the Permissions section by **unclicking** 'Use Compute Engine default service account' and enter your service account email.\n",
    "- then click 'Create'\n",
    "\n",
    "<img src=\"images/4_create_notebook.png\" width=\"800\">\n",
    "\n",
    "## Create a storage bucket if you haven't\n",
    "\n",
    "1. In the navigation menu `(≡)`, select `Cloud Storage` and then __Create bucket__.\n",
    "2. Enter a name for your bucket. You will reference this name when you need to transfer the output results from the GCP or running the nf-core/methylseq pipeline. You can also upload your own dataset to the bucket to use in GCP. (**NOTE**: Do not use underscores (_) in your bucket name. Use hyphens (-) instead.) \n",
    "3. Select __Region__ for the __Location type__ and select the __Location__ for your bucket.\n",
    "4. Select __Standard__ for the default storage class.\n",
    "5. Select __Uniform__ for the Access control.\n",
    "6. Select __Create__.\n",
    "7. Once the bucket is created, you will be redirected to the Bucket details page.\n",
    "8. Select Permissions, then + Add.\n",
    "9. Copy the email address of the Nextflow service account into New principals.\n",
    "10. Select the following roles:\n",
    "    - Storage Admin\n",
    "    - Storage Legacy Bucket Owner\n",
    "    - Storage Legacy Object Owner\n",
    "    - Storage Object Creator\n",
    "11. If you have a service account that need to access the bucket, repeat step 9 to enter the service account email, and step 10 to select the following roles: \n",
    "    - Storage Admin\n",
    "    - Storage Object Admin\n",
    "\n",
    "<div class=\"alert alert-block alert-danger\">\n",
    "    <i class=\"fa fa-exclamation-circle\" aria-hidden=\"true\"></i>\n",
    "    <b>Alert: </b>  Please <b>do not create a service key</b> if instructed by any tutorial. API keys are generally not considered secure; they are typically accessible to clients, making it easy for someone to steal an API key. Once the key is stolen, it has no expiration, so it may be used indefinitely, unless the project owner revokes or regenerates the key.\n",
    "</div>\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9513b1d1-ff29-451c-b916-23fb2ab42840",
   "metadata": {},
   "source": [
    "## Now open the JupyterLab with the Nextflow service account and download the tutorials from the repository as shown in the README.md before. \n",
    "\n",
    "Using command: `! git clone https://github.com/NIGMS/DNA-Methylation-Sequencing-Analysis-with-WGBS.git` and open **tutorial4_methylseq2.ipynb** \n",
    "\n",
    "**Note**: if your notebook has the Nextflow service account permission in the beginning, then you don't need to create a notebook and re-download the notebooks (.ipynb).\n",
    "\n",
    "---"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "26432290-b8e1-4fc6-882c-c23b9dede729",
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "os.environ[\"PATH\"] += os.pathsep + os.environ[\"HOME\"]+\"/mambaforge/bin\"\n",
    "\n",
    "! mamba install -c bioconda nextflow -y"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "32040353-96da-49c8-8f1e-b4e8213ce2c3",
   "metadata": {},
   "source": [
    "### Create and modify your own config file to include a 'gcb' profile block\n",
    "\n",
    "The config file allows Nextflow to utilize executors like Google Batch. Below is an example config file to run a Nextflow job using Google Batch:  \n",
    "```bash\n",
    "profiles{\n",
    "  gcb{\n",
    "      google.project = '<PROJECT_ID>'\n",
    "      google.location = 'us-central1'\n",
    "      google.region  = 'us-central1'\n",
    "      \n",
    "      process.executor = 'google-batch'\n",
    "      process.machineType = 'c2-standard-30'\n",
    "      \n",
    "      workDir = 'gs://BUCKET_NAME/work'\n",
    "      params.outdir = 'gs://BUCKET_NAME/result'\n",
    "     }\n",
    "}\n",
    "```  \n",
    "There are some fields that you need to define or pay attention to:\n",
    "- **Your project ID**. Not the project name, but the project ID. It can be found when you click the project name in the menu at the top of the home page: <img src=\"images/4_project_ID.png\" width=\"800\" />\n",
    "- **Executor**. To run the job using Google Batch, the executor must be defined here using: `process.executor = 'google-batch'`\n",
    "- **Region**. Make sure that your region is a region included in Google Batch. A comprehensive list is available [here](https://cloud.google.com/batch/docs/locations). \n",
    "- **Machine type**. Specify the machine type you would like to use, ensuring that there is enough memory and CPUs for the workflow. Google Cloud provides different machine types within several machine families that you can choose from to create a virtual machine (VM) instance with the resources you need. For example, the machine type we defined here is `c2-standard-30` that belongs the compute-optimized machine family, which has the highest performance per core on Compute Engine and optimized for compute-intensive workloads.`c2` is the machine series, and `30` is the CPU thread number it has. `c2-standard-30` also has a memory of 120GB and you can attach up to 3TB of local storage to these VMs for applications that require higher storage performance. For more information about the machine types please visit https://cloud.google.com/compute/docs/machine-types.\n",
    "    - If not defined, Google Batch will automatically use 1 CPU\n",
    "    - Something to consider is that `c2-standard-30 machine` type is compute intensive and a little more expensive than the `e2` or `n1` machine types. You can use a cheaper option if running time is not the first priority.\n",
    "- **Data storage**. For a full-scale dataset, make sure you create the bucket ahead of time and a directory in your specified bucket to store the input, output, and intermediate files. The `workDir` define the working directory for intermediate file to store. You can also define the input files and output path using the parameter `params.input` and `params.outdir` to specify your working directory bucket and output directory bucket\n",
    "    - If not defined, the work directory and output directory with be in your local notebook directory named `work`, and `results`. This is risky, since the intermediate and final outputs can be too large to store in the notebook instance.\n",
    "\n",
    "If many parameters need to be specified, you can write the config file using **scopes**. Configuration settings can be organized in different scopes by dot prefixing the property names with a scope identifier or grouping the properties in the same scope using the curly brackets notation:\n",
    "```bash\n",
    "profiles{\n",
    "  gcb{\n",
    "      workDir = 'gs://BUCKET_NAME/work'\n",
    "      process {\n",
    "          executor = 'google-batch'\n",
    "          machineType = 'c2-standard-30'\n",
    "      }\n",
    "      google {\n",
    "          location = 'us-central1'\n",
    "          region  = 'us-central1'\n",
    "          project = '<PROJECT_ID>'\n",
    "      }\n",
    "      params {\n",
    "          outdir = 'gs://BUCKET_NAME/output'\n",
    "          input = 'gs://BUCKET_NAME/*_R{1,2}.fastq.gz'\n",
    "          max_memory = 120.GB\n",
    "          max_cpus = 30\n",
    "          max_time = 24.h\n",
    "      }   \n",
    "}\n",
    "```\n",
    "\n",
    "__Note:__ Best practices are to make sure your working directory (`workDir`) and output directory (`outdir`) are **different**! Google Batch creates temporary files in the working directory within your bucket that do take up space. So once your pipeline has completed successfully, feel free to delete the temporary files.\n",
    "\n",
    "An example of the config file is `docs/test_LS.config`, we will use this config file to run the test profile (a very small dataset). "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a36934f9-846a-43fe-8897-762cbc59abd8",
   "metadata": {},
   "source": [
    "## Create a Google Cloud Storage Bucket\n",
    "\n",
    "You can create a customized bucket to use to store the output from the pipeline. Bucket names must be **globally unique** across all Google Cloud projects, including those outside of your organization.\n",
    "\n",
    "\n",
    "<b>Note:</b> We will use bucket name \"dna-methyl\" as an example in the following steps, but you need to create your own bucket to store the data. Please <b>replace</b> the \"dna-methyl\" below with your own bucket name\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1656a114-d2ca-476b-ace8-f0fee93ecb03",
   "metadata": {},
   "source": [
    "## Download and **Test** nf-core/methylseq using Google Batch<a name=\"TEST\"/>\n",
    "\n",
    "The `test` profile (`-profile test`) uses a small dataset allowing you to ensure the workflow works with your config file without long run times. Ensure you include:\n",
    "- Version of the nf-core tool [-r]\n",
    "- Location of the config file [-c]\n",
    "\n",
    "<div class=\"alert alert-block alert-success\">\n",
    "    <i class=\"fa fa-hand-paper-o\" aria-hidden=\"true\"></i>\n",
    "<b>Note:</b> Please <b>replace</b> the \"dna-methyl\" located in docs/test_LS.config with your own bucket name and save it. You can double-click the file to open and edit it, or edit it using terminal.\n",
    "</div>\n",
    "\n",
    "`docs/test_LS.config`:\n",
    "```\n",
    "profiles{\n",
    "  gcb{\n",
    "      workDir = 'gs://dna-methyl/test/work'\n",
    "      process {\n",
    "          executor = 'google-batch'\n",
    "          machineType = 'e2-standard-4'\n",
    "      }\n",
    "      google {\n",
    "          location = 'us-central1'\n",
    "          region  = 'us-central1'\n",
    "          project = '<PROJECT_ID>'\n",
    "      }\n",
    "      params {\n",
    "          outdir = 'gs://dna-methyl/test/results'\n",
    "      } \n",
    "     }\n",
    "}\n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c24aeed5-5f93-4af1-8c37-e9789aabaf58",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create the output directory for this tutorial\n",
    "! mkdir Tutorial_4"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "68f683b5-62a2-44f6-92c8-f27d46863124",
   "metadata": {},
   "outputs": [],
   "source": [
    "! rm -rf Tutorial_4/test\n",
    "!nextflow self-update\n",
    "!nextflow run nf-core/methylseq -r 2.4.0 -profile test,gcb -c docs/test_LS.config"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c4600c1c-e6c2-4d00-928b-3c097bcbc564",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Remove the remote trace file diretory\n",
    "! rm -rf gs:"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "599af349-877b-4f47-b2ab-fd3bc8f944be",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-block alert-warning\">\n",
    "    <i class=\"fa fa-pencil\" aria-hidden=\"true\"></i>\n",
    "<b>Note</b>: The <code>preseq</code> process may failed but ignored in the pipeline. This won't affect the output results. The preseq package is aimed at predicting and estimating the complexity of a genomic sequencing library, equivalent to predicting and estimating the number of redundant reads from a given sequencing depth and how many will be expected from additional sequencing using an initial sequencing experiment.\n",
    "    </div>\n",
    "\n",
    "This nf-core/methyseq test profile takes about 20 minutes to finish. When compared with the test profile running time (about 3 minutes) from Tutorial 3, we can see that there is extra time needed for Nextflow to talk to Google Batch, Cloud Storage, and VMs. It is not worthwhile for a small dataset, but this time difference can be ignored when running large datasets that need more computational resources. In other words, Google Batch works well for coarse-grained workloads i.e. long-running jobs. It’s not suggested to use this feature for pipelines spawning many short-lived tasks."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1a01c9ed-e2b7-48af-b5ef-1cd27a92e0c3",
   "metadata": {},
   "source": [
    "## An Example of a Real World Dataset<a name=\"REAL\"/>\n",
    "\n",
    "1. Install SRA-tools\n",
    "2. Download the Data\n",
    "3. Modify the config file\n",
    "4. Run the job"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9da4fb16-a862-480e-9002-8ee9861202db",
   "metadata": {},
   "source": [
    "#### Install SRA-tools\n",
    "\n",
    "The **[SRA Toolkit](https://github.com/ncbi/sra-tools/wiki)** and SDK from NCBI is a collection of tools and libraries for using data in the INSDC Sequence Read Archives. Here we use `mamba` to install `sra-tools`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1314272a-5980-4c46-9217-838de16313d4",
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "os.environ[\"PATH\"] += os.pathsep + os.environ[\"HOME\"]+\"/mambaforge/bin\"\n",
    "\n",
    "! mamba install -c bioconda \"sra-tools > 2.11\" -y"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0db87bc0-f745-4aba-aecd-5c4606238594",
   "metadata": {},
   "source": [
    "#### Download the Data\n",
    "\n",
    "The data was from [Molaro, Antoine, et al. Cell 146.6 (2011): 1029-1041](https://www.sciencedirect.com/science/article/pii/S0092867411009421) and [Laurent, Louise, et al. \" Genome research 20.3 (2010): 320-331](https://genome.cshlp.org/content/20/3/320.full). During germ cell and preimplantation development, mammalian cells undergo nearly complete reprogramming of DNA methylation patterns. The studies profiled the methylomes of human and chimp sperm as a basis for comparison to methylation patterns of embryonic stem cells (ESCs).   \n",
    "<img src=\"images/4_data_graph.jpg\" width=\"300\" />\n",
    "\n",
    "We use one sample from human sperm and one sample from ESCs as examples to demonstrate the workflow here.\n",
    "\n",
    "Use `fasterq-dump` to download data from SRA using accession numbers. The data will be store at `Tutorial_4/sra_download`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "04e74cb8-1ffc-484f-8b26-4c11703cbc66",
   "metadata": {},
   "outputs": [],
   "source": [
    "! fasterq-dump --threads 4 --progress SRR306435 SRR033942 -O Tutorial_4/sra_download"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "633c0fd0-0309-4451-96ae-ea4156a5b734",
   "metadata": {},
   "source": [
    "Remove the temporary output directory from running `fasterq-dump`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3ab5a465-a3a7-4968-84c4-f6d7bed30dee",
   "metadata": {},
   "outputs": [],
   "source": [
    "! rm -rf fasterq.tmp.*"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6152126a-9c1a-4eff-a8ea-e51f45bb9a44",
   "metadata": {},
   "source": [
    "Compress the files."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8efc1fcd-7496-4186-9b36-94678283e877",
   "metadata": {},
   "outputs": [],
   "source": [
    "#faster way to compress the files\n",
    "!pigz Tutorial_4/sra_download/SRR*"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1775f2b8-83a4-4e96-aa6b-64108377f22e",
   "metadata": {},
   "source": [
    "#### Create a samplesheet (located in Tutorial_3) to provide all sample information\n",
    "\n",
    "**Format:**    \n",
    "sample, fastq1, fastq2    \n",
    "sample1,sample1_R1.fastq,sample1_R2.fastq    \n",
    "control1,control1_R1.fastq,control1  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "27799ac8-6fd2-4c2e-abb6-af9c81bbef57",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Pandas DataFrame by lists of dicts.\n",
    "import pandas as pd\n",
    " \n",
    "# Initialize data to lists.\n",
    "samples = [{'sample': 'SRR033942', 'fastq_1': 'Tutorial_4/sra_download/SRR033942_1.fastq.gz', 'fastq_2': 'Tutorial_4/sra_download/SRR033942_2.fastq.gz'},\n",
    "        {'sample': 'SRR306435', 'fastq_1': 'Tutorial_4/sra_download/SRR306435_1.fastq.gz', 'fastq_2': 'Tutorial_4/sra_download/SRR306435_2.fastq.gz'}\n",
    "       ]\n",
    " \n",
    "# Creates DataFrame.\n",
    "df2 = pd.DataFrame(samples)\n",
    " \n",
    "# Print the data\n",
    "df2"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1e7b35d6-c886-410f-a040-371d26c4c3e0",
   "metadata": {},
   "source": [
    "Export dataframe to CSV file."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fa9bf795-bdd9-486e-ae8e-2ae4fa6c25df",
   "metadata": {},
   "outputs": [],
   "source": [
    "df2.to_csv('Tutorial_4/samplesheet.csv', index=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c436deda-9341-40db-9af1-b760478cee08",
   "metadata": {},
   "source": [
    "#### Create/Modify the config file \n",
    "\n",
    "As mentioned [previously](#Create-and-modify-your-own-config-file-to-include-a-'gcb'-profile-block), we need to modify the config file for the methylseq to run in Google Batch. The config file is located at `docs/human_sperm.config`. In this example, we set the working and output directory in the GCP Cloud Storage bucket `Tutorial_4/methyseq_sperm`. You need to change the path to your own result bucket. The input file is the sample sheet in the directory `Tutorial_4`, that we just created.\n",
    "\n",
    "```bash\n",
    "profiles{\n",
    "  gcb{\n",
    "      process.executor = 'google-batch'\n",
    "      google.location = 'us-central1'\n",
    "      google.region  = 'us-central1'\n",
    "      google.project = '<PROJECT_ID>'\n",
    "      workDir = \"gs://nosi-hawaii-dna-27fa/methyseq_sperm/work\"\n",
    "      params.outdir = \"gs://methyl/methyseq_sperm/results\"\n",
    "      process.machineType = 'c2-standard-16'\n",
    "     }\n",
    "}\n",
    "```\n",
    "\n",
    "<div class=\"alert alert-block alert-success\">\n",
    "    <i class=\"fa fa-hand-paper-o\" aria-hidden=\"true\"></i>\n",
    "<b>Again:</b> Please <b>replace</b> the \"dna-methyl\" in the docs/human_sperm.config with your own bucket name.\n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b000a077-b203-4b31-949b-39fb4b56fc0d",
   "metadata": {},
   "source": [
    "#### Run methylseq using Google Batch\n",
    "\n",
    "If not defined in the config file, you can always use command line parameters:\n",
    "- `-r` pipeline version\n",
    "- `-profile` profile to use ('gcb' was defined in docs/human_sperm.config)\n",
    "- `-c` config file\n",
    "- `--genome` the reference genome to use. Here we use human assembly GRCh38 as the reference genome\n",
    "- `--clip_r1` instructs Trim Galore to remove certein number of bps from the 5' end of read 1\n",
    "- `--tracedir` defines a local diretory to save the pipeline information\n",
    "\n",
    "There will be some pipeline information saved to the default `results` directory. So please make sure the directory is empty before running the pipeline."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5f15f673-f29e-4d75-a6b4-5801ca4a1dc0",
   "metadata": {},
   "outputs": [],
   "source": [
    "! rm -rf Tutorial_4/methyseq_sperm\n",
    "\n",
    "! nextflow run nf-core/methylseq \\\n",
    "    -profile gcb \\\n",
    "    -r 2.6.0 \\\n",
    "    -c docs/human_sperm.config  \\\n",
    "    --input 'Tutorial_4/samplesheet.csv' \\\n",
    "    --genome GRCh38 \\\n",
    "    --clip_r1 2 \\\n",
    "    --tracedir 'Tutorial_4/methyseq_sperm/pipeline_info' \\\n",
    "    -resume\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2695dbca-104a-4bab-b0bc-393c0b3d547d",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Remove the remote trace file diretory\n",
    "! rm -rf gs:"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "02012eb4-cc48-49b8-a7df-5c270b776734",
   "metadata": {},
   "source": [
    "#### Check to see if files are in your output directory bucket\n",
    "\n",
    "The output files should be saved in your bucket's methylseq_sperm/results directory. You can list the results directory to see the file structures. You can also copy the files to your local directory to view them. For example, the MultiQC report file is located at `gs://dna-methyl/methyseq_sperm/results/MultiQC/multiqc_report.html`. Let's copy and view it using the commands below:\n",
    "\n",
    "<div class=\"alert alert-block alert-success\">\n",
    "    <i class=\"fa fa-hand-paper-o\" aria-hidden=\"true\"></i>\n",
    "<b>Note:</b> Please <b>replace</b> the \"dna-methyl\" in commands below with your own bucket name.\n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "093cf8e2-eb12-46fd-9fd3-ca0ac5e35a3d",
   "metadata": {},
   "outputs": [],
   "source": [
    "# List the output files/directories in the results folder\n",
    "! gsutil ls gs://dna-methyl/methyseq_sperm/results\n",
    "\n",
    "# Copy the multiQC output multiqc_report.html to local notebook:\n",
    "! gsutil cp -r gs://dna-methyl/methyseq_sperm/results/multiqc/bismark/multiqc_report.html .\n",
    "\n",
    "# View the MultiQC output HTML file:\n",
    "from IPython.display import IFrame\n",
    "IFrame(src='multiqc_report.html', width=900, height=600)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c1ad5eda-2318-4ec0-a353-0dac3955bbf7",
   "metadata": {},
   "source": [
    "There are two files (`execution_timeline.html` and `execution_report.html`) about the pipeline running information will be saved in the null/pipeline_info directory locally in the notebook, which can provide detailed information about the running time for each process and the their resource usages. This can provide more insights for potential optimizations. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1e0a760c-8805-4354-a2e6-f4a317d3f7ed",
   "metadata": {},
   "outputs": [],
   "source": [
    "#use the file names listed in here\n",
    "!ls null/pipeline_info"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "79af3cc4-caf1-4d4f-8133-d307dbaa2916",
   "metadata": {},
   "outputs": [],
   "source": [
    "from IPython.display import IFrame\n",
    "IFrame(src='null/pipeline_info/execution_timeline.html', width=800, height=600)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ac49deb8-3b8d-4554-9d0c-6a6a9a9aa772",
   "metadata": {},
   "outputs": [],
   "source": [
    "from IPython.display import IFrame\n",
    "IFrame(src='null/pipeline_info/execution_report.html', width=800, height=600)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b4052257-cebd-4c6d-b942-f03a40269540",
   "metadata": {},
   "source": [
    "## <a name=\"FULL\" />Configuration of a Full-scale Dataset - Troubleshooting\n",
    "\n",
    "For a full-scale WGBS study, the sequencing data size can range from several hundred GBs to several TBs. For example, the data we downloaded in this tutorial: GSE30340 and GSE19418, both have many runs with the size add up to several hundred GBs. Given the large data files, the storage and memory can become an issue when running the pipeline as instructed in this tutorial. \n",
    "\n",
    "#### Download the data\n",
    "\n",
    "There are several options that we can use to download the data:\n",
    "1. Download the data in a notebook. You need to make sure that the disk size you assigned to the notebook is enough for the data that you want to download. Also, when you use `prefetch` from SRA toolkit, there is a default maximum download-size of 20G; you will need to increase that limit.  \n",
    "2. Cloud Data Delivery Service. SRA has created a cloud data delivery service to deliver the source files and other file types from NCBI cold storage buckets to individual data consumers' buckets in AWS and GCP. This service is provided for both public and authorized access (dbGaP) data. [More detailed information here](https://www.ncbi.nlm.nih.gov/sra/docs/data-delivery/).\n",
    "3. Upload to the storage bucket directly. You can upload the data to the GCP storage bucket directly from your local computer, HPC, or service server using the `gsutil` tool. [More detailed information here](https://cloud.google.com/storage/docs/discover-object-storage-gsutil). \n",
    "\n",
    "#### Troubleshooting the nf-core pipeline\n",
    "\n",
    "If the nf-core pipeline does not complete successfully, you can refer to the [troubleshooting](https://nf-co.re/usage/troubleshooting) page that nf-core provided for more information. For our tutorial here, the most likely reasons that the pipeline fails are:\n",
    "- service account is not set up correctly\n",
    "- file paths are not correct\n",
    "- memory or storage issues for large dataset\n",
    "\n",
    "If you have a command exit status of 104, 134, 137, 139, 143, 247, the most probable cause is an \"out of memory\" issue. To solve the memory issue, you need to increase the memory limit in the configuration file for the process that fails. For example:   \n",
    "``` bash\n",
    "profiles {\n",
    "    gcb {\n",
    "      process {\n",
    "        withName: qualimap { \n",
    "              machineType = 'c2-standard-16'\n",
    "              cpus = 16\n",
    "              memory = 64.GB\n",
    "        }\n",
    "      }\n",
    "    }\n",
    "}\n",
    "```\n",
    "\n",
    "In GCP, the memory is also limited by the [machine type](https://cloud.google.com/compute/docs/machine-types) you select to run the process. For example, if you choose `c2-standard-8` then the memory is limited to 32GB. You can change the machine types to increase the memory. There are [memory-optimized machine families](https://cloud.google.com/compute/docs/memory-optimized-machines) (m1, m2) that you can use for workloads that require higher memory-to-vCPU ratios.  \n",
    "\n",
    "#### Optimize nf-core/methylseq configuration\n",
    "\n",
    "The nf-core/methylseq workflow contains multiple processes, and the requirements of computational and memory resources for each process vary a lot. For better performance or billing purposes, you can change the configuration for each process. You can check the default settings for each process at the pipeline's [base.config](https://github.com/nf-core/methylseq/blob/master/conf/base.config) file. \n",
    "\n",
    "As an example of running a 12 sample WGBS data, [docs/optimization_example.config](docs/optimization_example.config) was the config file that finish processing these 24 fastq files (pair-end, averge size 325M reads per fastq file, ) using Google Batch in less than 30 hours."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "71539336-703c-4c71-ad0d-b6aabae94189",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-block alert-danger\">\n",
    "    <i class=\"fa fa-exclamation-circle\" aria-hidden=\"true\"></i>\n",
    "    <b>Don't forget:</b> after finish running the notebook, stop the notebook in Vertex AI Workbench to avoid cost accumulation.\n",
    "</div>"
   ]
  }
 ],
 "metadata": {},
 "nbformat": 4,
 "nbformat_minor": 5
}