ATACSEQ.ipynb

{
  "nbformat": 4,
  "nbformat_minor": 0,
  "metadata": {
    "colab": {
      "provenance": [],
      "include_colab_link": true
    },
    "kernelspec": {
      "name": "python3",
      "display_name": "Python 3"
    },
    "language_info": {
      "name": "python"
    }
  },
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "view-in-github",
        "colab_type": "text"
      },
      "source": [
        "<a href=\"https://colab.research.google.com/github/vpatel85/cookiecutter-django-rest/blob/master/ATACSEQ.ipynb\" target=\"_parent\"><img src=\"https://colab.research.google.com/assets/colab-badge.svg\" alt=\"Open In Colab\"/></a>"
      ]
    },
    {
      "cell_type": "markdown",
      "source": [
        "#Introduction\n",
        "\n",
        "\n",
        "Over the past decade, advances in Next-Generation Sequencing (NGS) have disrupted life sciences. NGS has opened new horizons in the study of living systems, including their genomes, transcriptomes, chromatin states, and beyond. In recent years, NGS has been used to study gene expression regulation based on chromatin accessibility profiles in a genome-wide manner. The most powerful tool to emerge in this field is the Assay for Transposase-Accessible Chromatin using sequencing – aka, ATAC-seq.\n",
        "\n",
        "Experimentally, ATAC-seq uses a bioengineered hyperactive transposase (Tn5) which selectively integrates its adapter sequences into regions of accessible chromatin. These regions presumably contain genomic elements regulating active transcription of downstream genes, such as enhancers, promoters, transcription factor binding sites, etc. (see Fig. 1a.) After adapter transposition, the captured regions can be amplified and sequenced.\n",
        "\n",
        "Due to its versatility, simplicity, and sensitivity, ATAC-seq has largely superseded other similar techniques such as DNase-seq or MNase-seq (see Fig. 1b, or the original paper by Buenrostro et al, 2015, for additional technical details). In recent years, ATAC-seq protocols have been coupled to microfluidic techniques allowing open chromatin profiling on a single-cell resolution.\n",
        "\n",
        "\n",
        "![](https://drive.google.com/uc?export=view&id=11VdW8Wy2nC_28oBW2V4IqmMqVFxT9DIT)\n",
        "\n",
        "*Fig. 1: The experimental principle of ATAC-seq. (a) The principle of integration of adapter sequences by Tn5 transposase into genomic loci of open chromatin regions; (b) Experimental advantages of ATAC-seq over FAIRE-seq and DNase-seq. Adapted from Buenrostro et al. 2015.*\n",
        "\n",
        "\n",
        "\n",
        "\n",
        "This end-to-end tutorial will guide you through every step of ATAC-seq data analysis. We’ll show you how perform QC and read mapping, peak calling, and differential accessibility analysis, and finally visualize and interpret the results.\n",
        "\n",
        "There are roughly seven steps in our ATAC-seq data analysis. For a simplified schematic representation of this pipeline, see Figure 2 below. \n",
        "\n",
        "\n",
        "![](https://drive.google.com/uc?export=view&id=11WAPv2m4oQmgwHyMym4JaRxOfC2V8_BI)\n",
        "\n",
        "*Fig. 2: The simplified schematic flowchart of the ATAC-seq data analysis pipeline covered in this tutorial. Blue boxes indicate the types of analysis. The text above the boxes indicates the software used for each analysis.*\n",
        "\n",
        "Before we begin, let’s briefly explore this pipeline.\n",
        "First, the data will be downloaded from Sequence Read Archive (SRA) using the online software SRA Explorer. After obtaining the data, we will run a quality control analysis using FastQC and MultiQC (Ewels et al. 2016) software to check and visualize the quality of the data, respectively. Next, we will use Trimmomatic (Bolger, Lohse, and Usadel 2014) software to remove the traces of Nextera adapter sequences that we introduced during the transposase integration process.\n",
        "\n",
        "We now have clean data, and will use BWA (Li and Durbin 2009) software to map the reads against corresponding reference genomes. This will generate alignment files for further ATAC-seq peak calling, which we’ll do with the MACS2 (Zhang et al. 2008) package. Then the obtained data will be visualized with Integrative Genomic Viewer (Robinson et al. 2011). Finally, we will use the R DiffBind package (Ross-Innes et al. 2012) to run quality control and perform the differential chromatin accessibility analysis.\n"
      ],
      "metadata": {
        "id": "EnP7JweTlATq"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "#Tiny installation"
      ],
      "metadata": {
        "id": "shfP5c4OW1na"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "!pip install tiny-cli\n",
        "!pip install pdf2image\n",
        "!apt-get install poppler-utils\n",
        "\n",
        "#importing tiny for the spin up for computing resources to run this analysis\n",
        "import tiny\n",
        "#for inline visualization of paper figures\n",
        "from IPython.display import Image, HTML, display\n",
        "#for inline visualization of relevant table figures\n",
        "import pandas as pd\n",
        "from pdf2image import convert_from_path\n",
        "\n",
        "\n",
        "#define the storage bucket used for accessing the data and output of the experiments.\n",
        "\n"
      ],
      "metadata": {
        "id": "s2gg9W6bXDiQ",
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "outputId": "5d9ebfcb-6fa4-47b6-b430-7dc0e204add1"
      },
      "execution_count": null,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "Looking in indexes: https://pypi.org/simple, https://us-python.pkg.dev/colab-wheels/public/simple/\n",
            "Collecting tiny-cli\n",
            "  Downloading tiny-cli-1.0.32.tar.gz (9.8 kB)\n",
            "  Preparing metadata (setup.py) ... \u001b[?25l\u001b[?25hdone\n",
            "Collecting httpcore[http2]\n",
            "  Downloading httpcore-0.17.0-py3-none-any.whl (70 kB)\n",
            "\u001b[2K     \u001b[90m━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━\u001b[0m \u001b[32m70.6/70.6 KB\u001b[0m \u001b[31m1.4 MB/s\u001b[0m eta \u001b[36m0:00:00\u001b[0m\n",
            "\u001b[?25hCollecting httpx==0.23.1\n",
            "  Downloading httpx-0.23.1-py3-none-any.whl (84 kB)\n",
            "\u001b[2K     \u001b[90m━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━\u001b[0m \u001b[32m85.0/85.0 KB\u001b[0m \u001b[31m2.7 MB/s\u001b[0m eta \u001b[36m0:00:00\u001b[0m\n",
            "\u001b[?25hRequirement already satisfied: tabulate in /usr/local/lib/python3.9/dist-packages (from tiny-cli) (0.8.10)\n",
            "Collecting anytree\n",
            "  Downloading anytree-2.8.0-py2.py3-none-any.whl (41 kB)\n",
            "\u001b[2K     \u001b[90m━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━\u001b[0m \u001b[32m41.7/41.7 KB\u001b[0m \u001b[31m2.6 MB/s\u001b[0m eta \u001b[36m0:00:00\u001b[0m\n",
            "\u001b[?25hCollecting httpcore<0.17.0,>=0.15.0\n",
            "  Downloading httpcore-0.16.3-py3-none-any.whl (69 kB)\n",
            "\u001b[2K     \u001b[90m━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━\u001b[0m \u001b[32m69.6/69.6 KB\u001b[0m \u001b[31m6.1 MB/s\u001b[0m eta \u001b[36m0:00:00\u001b[0m\n",
            "\u001b[?25hRequirement already satisfied: certifi in /usr/local/lib/python3.9/dist-packages (from httpx==0.23.1->tiny-cli) (2022.12.7)\n",
            "Requirement already satisfied: sniffio in /usr/local/lib/python3.9/dist-packages (from httpx==0.23.1->tiny-cli) (1.3.0)\n",
            "Collecting rfc3986[idna2008]<2,>=1.3\n",
            "  Downloading rfc3986-1.5.0-py2.py3-none-any.whl (31 kB)\n",
            "Requirement already satisfied: six>=1.9.0 in /usr/local/lib/python3.9/dist-packages (from anytree->tiny-cli) (1.16.0)\n",
            "Requirement already satisfied: anyio<5.0,>=3.0 in /usr/local/lib/python3.9/dist-packages (from httpcore<0.17.0,>=0.15.0->httpx==0.23.1->tiny-cli) (3.6.2)\n",
            "Collecting h11<0.15,>=0.13\n",
            "  Downloading h11-0.14.0-py3-none-any.whl (58 kB)\n",
            "\u001b[2K     \u001b[90m━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━\u001b[0m \u001b[32m58.3/58.3 KB\u001b[0m \u001b[31m3.7 MB/s\u001b[0m eta \u001b[36m0:00:00\u001b[0m\n",
            "\u001b[?25hCollecting h2<5,>=3\n",
            "  Downloading h2-4.1.0-py3-none-any.whl (57 kB)\n",
            "\u001b[2K     \u001b[90m━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━\u001b[0m \u001b[32m57.5/57.5 KB\u001b[0m \u001b[31m2.7 MB/s\u001b[0m eta \u001b[36m0:00:00\u001b[0m\n",
            "\u001b[?25hRequirement already satisfied: idna>=2.8 in /usr/local/lib/python3.9/dist-packages (from anyio<5.0,>=3.0->httpcore<0.17.0,>=0.15.0->httpx==0.23.1->tiny-cli) (3.4)\n",
            "Collecting hyperframe<7,>=6.0\n",
            "  Downloading hyperframe-6.0.1-py3-none-any.whl (12 kB)\n",
            "Collecting hpack<5,>=4.0\n",
            "  Downloading hpack-4.0.0-py3-none-any.whl (32 kB)\n",
            "Building wheels for collected packages: tiny-cli\n",
            "  Building wheel for tiny-cli (setup.py) ... \u001b[?25l\u001b[?25hdone\n",
            "  Created wheel for tiny-cli: filename=tiny_cli-1.0.32-py3-none-any.whl size=8827 sha256=27c6cb203f5fd6bcb3d1a1a30ac0eebf6f60365abf51dd64156f5a5c73b9aeac\n",
            "  Stored in directory: /root/.cache/pip/wheels/b7/91/8e/ca2474fa4b992a9a573dbef6fe8e55427021fdf59bcd1c469b\n",
            "Successfully built tiny-cli\n",
            "Installing collected packages: rfc3986, hyperframe, hpack, h11, anytree, httpcore, h2, httpx, tiny-cli\n",
            "Successfully installed anytree-2.8.0 h11-0.14.0 h2-4.1.0 hpack-4.0.0 httpcore-0.16.3 httpx-0.23.1 hyperframe-6.0.1 rfc3986-1.5.0 tiny-cli-1.0.32\n",
            "Looking in indexes: https://pypi.org/simple, https://us-python.pkg.dev/colab-wheels/public/simple/\n",
            "Collecting pdf2image\n",
            "  Downloading pdf2image-1.16.3-py3-none-any.whl (11 kB)\n",
            "Requirement already satisfied: pillow in /usr/local/lib/python3.9/dist-packages (from pdf2image) (8.4.0)\n",
            "Installing collected packages: pdf2image\n",
            "Successfully installed pdf2image-1.16.3\n",
            "Reading package lists... Done\n",
            "Building dependency tree       \n",
            "Reading state information... Done\n",
            "The following NEW packages will be installed:\n",
            "  poppler-utils\n",
            "0 upgraded, 1 newly installed, 0 to remove and 24 not upgraded.\n",
            "Need to get 174 kB of archives.\n",
            "After this operation, 754 kB of additional disk space will be used.\n",
            "Get:1 http://archive.ubuntu.com/ubuntu focal-updates/main amd64 poppler-utils amd64 0.86.1-0ubuntu1.1 [174 kB]\n",
            "Fetched 174 kB in 2s (111 kB/s)\n",
            "Selecting previously unselected package poppler-utils.\n",
            "(Reading database ... 122349 files and directories currently installed.)\n",
            "Preparing to unpack .../poppler-utils_0.86.1-0ubuntu1.1_amd64.deb ...\n",
            "Unpacking poppler-utils (0.86.1-0ubuntu1.1) ...\n",
            "Setting up poppler-utils (0.86.1-0ubuntu1.1) ...\n",
            "Processing triggers for man-db (2.9.1-1) ...\n"
          ]
        }
      ]
    },
    {
      "cell_type": "code",
      "source": [
        "auth = tiny.Auth('eyJhbGciOiJIUzI1NiIsInR5cCI6IkpXVCJ9.eyJzdWIiOiIxIiwiYXVkIjpbImZhc3RhcGktdXNlcnM6YXV0aCJdLCJleHAiOjE2ODEyNDQwNTB9.SBS7-8wVO_9lhcb9889WyKqh8Bx8cL4_vqT7-sly8wM')"
      ],
      "metadata": {
        "id": "bJXlihnvMDkd"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "tiny.create_workbench('atac-seq', auth=auth)"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "ivJdmhmLsqkQ",
        "outputId": "58f24a72-727d-49a5-dca0-c2bda4586053"
      },
      "execution_count": null,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "\n",
            "        The atac-seq-20230406175308669217 workbench is now available. \n",
            "\n",
            "        The command workbench.list_files() will return the list of files in your workbench. Note, by default, we've included files for running through an RNA-Seq, ATAC-Seq, and variant calling experiments which are outline here:\n",
            "        \n",
            "        https://docs.tinybio.cloud/docs\n",
            "\n",
            "        The command workbench.run(tool=TOOL_NAME, full_command=COMMAND) will create a super computer (16 cores, 256GB RAM) with the specified tool installed and run the command specified in full_command. \n",
            "\n",
            "        To check the status of your commands for a workbench please run. workbench.logs(). \n",
            "\n",
            "        To upload a file directly from your machine run workbench.upload('file_path_on_your_machine'). Please note, if you're uploading from a colab notebook, you will need to first upload the file to the colab instance and then upload it from that instance. \n",
            "\n",
            "        To upload from a remote machine run workbench.upload_job(method=\"curl/wget\", files=[(\"public_file_url\",\"download_path\")]). \n",
            "\n",
            "        To download a file run the following workbench.download('file_path_on_the_workbench'). This will generate a download URL.\n",
            "    \n"
          ]
        },
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": [
              "Workbench(atac-seq-20230406175308669217)"
            ]
          },
          "metadata": {},
          "execution_count": 3
        }
      ]
    },
    {
      "cell_type": "code",
      "source": [
        "# Use the workbench created above or an existing one\n",
        "workbench = tiny.Workbench('atac-seq-20230406175308669217', auth=auth)"
      ],
      "metadata": {
        "id": "XwTp2X86snM1"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "# Download sample date"
      ],
      "metadata": {
        "id": "myp3ZD0SuMCn"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "files = [\n",
        "    ('ftp://ftp.sra.ebi.ac.uk/vol1/fastq/SRR102/093/SRR10261593/SRR10261593_1.fastq.gz', 'input/SRR10261593_1.fastq.gz'),\n",
        "    ('ftp://ftp.sra.ebi.ac.uk/vol1/fastq/SRR102/093/SRR10261593/SRR10261593_2.fastq.gz', 'input/SRR10261593_2.fastq.gz'),\n",
        "    ('ftp://ftp.sra.ebi.ac.uk/vol1/fastq/SRR102/094/SRR10261594/SRR10261594_1.fastq.gz', 'input/SRR10261594_1.fastq.gz'),\n",
        "    ('ftp://ftp.sra.ebi.ac.uk/vol1/fastq/SRR102/094/SRR10261594/SRR10261594_2.fastq.gz', 'input/SRR10261594_2.fastq.gz'),\n",
        "    ('ftp://ftp.sra.ebi.ac.uk/vol1/fastq/SRR102/096/SRR10261596/SRR10261596_1.fastq.gz', 'input/SRR10261596_1.fastq.gz'),\n",
        "    ('ftp://ftp.sra.ebi.ac.uk/vol1/fastq/SRR102/096/SRR10261596/SRR10261596_2.fastq.gz', 'input/SRR10261596_2.fastq.gz'),\n",
        "    ('ftp://ftp.sra.ebi.ac.uk/vol1/fastq/SRR102/095/SRR10261595/SRR10261595_1.fastq.gz', 'input/SRR10261595_1.fastq.gz'),\n",
        "    ('ftp://ftp.sra.ebi.ac.uk/vol1/fastq/SRR102/095/SRR10261595/SRR10261595_2.fastq.gz', 'input/SRR10261595_2.fastq.gz'),\n",
        "    ('ftp://ftp.sra.ebi.ac.uk/vol1/fastq/SRR102/091/SRR10261591/SRR10261591_1.fastq.gz', 'input/SRR10261591_1.fastq.gz'),\n",
        "    ('ftp://ftp.sra.ebi.ac.uk/vol1/fastq/SRR102/091/SRR10261591/SRR10261591_2.fastq.gz', 'input/SRR10261591_2.fastq.gz'),\n",
        "    ('ftp://ftp.sra.ebi.ac.uk/vol1/fastq/SRR102/092/SRR10261592/SRR10261592_1.fastq.gz', 'input/SRR10261592_1.fastq.gz'),\n",
        "    ('ftp://ftp.sra.ebi.ac.uk/vol1/fastq/SRR102/092/SRR10261592/SRR10261592_2.fastq.gz', 'input/SRR10261592_2.fastq.gz'),\n",
        "]\n",
        "workbench.upload_job(method='curl', files=files)"
      ],
      "metadata": {
        "id": "tbXnWJZ-uK4j"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "#Data Quality Control\n",
        "\n",
        "\n",
        "Now we need to perform quality control of the raw sequencing data to ensure that only high quality reads are used moving forward. For this step, we’ll use FastQC, which analyzes the fastq files (the data format of the raw sequencing data) and reports important statistics and parameters of the data. To run FastQC (v. 0.11.9), first create a directory called fastqc_initial, and then run the fastqc command:"
      ],
      "metadata": {
        "id": "yQQ4B31UFQnj"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "workbench.run(tool='fastqc', full_command='fastqc -t 10 /input/*gz -o /output/fastqc_initial')"
      ],
      "metadata": {
        "id": "5Lg-CtREFXlf",
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "outputId": "530063cc-0de1-46f2-fd76-2341b386e296"
      },
      "execution_count": null,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "+------------------+--------+-----------+----------+-------------------------------------------+---------------------------------------------------+\n",
            "| Job ID           | Tool   | Version   | Status   | Get Logs                                  | Full Command                                      |\n",
            "+==================+========+===========+==========+===========================================+===================================================+\n",
            "| fastqc9f22d8e2af | fastqc | 0.11.8    | Queued   | workbench.jobs('fastqc9f22d8e2af').logs() | fastqc -t 10 /input/*gz -o /output/fastqc_initial |\n",
            "+------------------+--------+-----------+----------+-------------------------------------------+---------------------------------------------------+\n"
          ]
        }
      ]
    },
    {
      "cell_type": "markdown",
      "source": [
        "In this command, the term fastqc.gz means that FastQC will analyze all files in the current directory ending with fastq.gz and store the results in the fastqc_initial directory. The parameter -t 2 instructs FastQC to use two cores of the computer.\n",
        "\n",
        "Because FastQC reports the results for each fastq file separately, we now have 12 html quality reports to review. Rather than inspect these one-by-one (which would be tedious), we’ll use the MultiQC tool* (v. 1.12) to aggregate them. First, navigate to the fastqc_initial folder, then run MultiQC:"
      ],
      "metadata": {
        "id": "7YbqCDVPFbuf"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "workbench.run(tool='multiqc', full_command='multiqc -o /output/fastqc_initial/multiqc /output/fastqc_initial')"
      ],
      "metadata": {
        "id": "UIQ49z0HFir4",
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "outputId": "ffc5f9b8-643d-4fad-d341-c3a6633afc9f"
      },
      "execution_count": null,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "+-------------------+---------+-----------+----------+--------------------------------------------+------------------------------------------------------------------+\n",
            "| Job ID            | Tool    |   Version | Status   | Get Logs                                   | Full Command                                                     |\n",
            "+===================+=========+===========+==========+============================================+==================================================================+\n",
            "| multiqc912fde6de4 | multiqc |      1.14 | Queued   | workbench.jobs('multiqc912fde6de4').logs() | multiqc -o /output/fastqc_initial/multiqc /output/fastqc_initial |\n",
            "+-------------------+---------+-----------+----------+--------------------------------------------+------------------------------------------------------------------+\n"
          ]
        }
      ]
    },
    {
      "cell_type": "markdown",
      "source": [
        "MultiQC will scan the current directory and collect the output files of FastQC, generating a single, beautiful html report of all our samples. (Note that MultiQC can handle output from numerous software types – e.g. mappers, trimmers, etc. – to create informative, interactive, publication-ready reports.) To inspect the aggregated report, open the multiqc_report.html in your web browser. \n",
        "\n"
      ],
      "metadata": {
        "id": "GNKtwJocFlwy"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "This plot shows the overview of our data and its main parameters (Fig. 3). Let's review some of its features:\n",
        "\n",
        "\n",
        "\n",
        "![](https://drive.google.com/uc?export=view&id=11YETLqui5KZMufPHcNO0tCsGQXJ9ZXpf)\n",
        "\n",
        "*Fig. 3: The general statistics of the analyzed fastq samples, showing the sample names, duplication percentages, GC content percentages, and the number of sequencing reads.*\n",
        "\n",
        "\n",
        "First, notice that our percentage of duplicated reads (\"% Dups\") is relatively high. This is to be expected in ATAC-seq data, as the sequenced regions are generally short. This can generate identical reads – especially in cases like ours, where data is high-coverage. These numbers also indicate PCR duplicates, which are technical artifacts and will be removed in later steps of our workflow.\n",
        "\n",
        "Second, our GC content of ~39-40% is expected for S. cerevisiae and S. uvarum, and is a good signal that the data is not contaminated. And finally, the last column indicates the numbers of sequencing reads. As mentioned above, each sample has ~30 million reads per each direction (i.e. 60 million paired-end reads).\n",
        "\n",
        "Let's move on to another part of our MultQC report: the Phred score plot (Fig. 4). Phred score is an important metric of sequencing read quality. It indicates the reliability of each sequenced base (learn more here [link text](https://www.illumina.com/documents/products/technotes/technote_Q-Scores.pdf)), with values between 30 and 40 considered good scores:\n",
        "\n",
        "\n",
        "![](https://drive.google.com/uc?export=view&id=11_E3y9iiGIM7tWomiL5ISNBP0RujowxP)\n",
        "\n",
        "*Fig. 4: Mean Phred-quality scores of each sample*\n",
        "\n"
      ],
      "metadata": {
        "id": "J42J0i4NF0N5"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "The plot above illustrates the quality scores of each base position in all samples (each line corresponds to a sample). Although there is a dip at position 11 in some samples, the overall quality of the ATAC-seq data is very high. That means our data does not need any correction for this metric.\n",
        "\n",
        "Finally, let's review this chart of our adapter content (Fig. 5):\n",
        "\n",
        "![](https://drive.google.com/uc?export=view&id=11b0iOgoqCtIvOPsbdrbCTOyHp2py9v3E)\n",
        "\n",
        "*Fig. 5: Presence of Nextera adapters sequences in all studies samples.*\n",
        "\n",
        "\n",
        "\n",
        "\n",
        "\n",
        "The presence of adapter sequences (which are inserted into the DNA during the experimental sample preparation) is an important parameter of our sequencing data. Adapter sequences can interfere with downstream analyses (for example by reducing the read-mapping rate), and thus should be removed from the data if present. As you can see above, our samples contain trace amounts of Nextera transposase adapter sequence, starting from the middle of the reads. These should be trimmed away in order to improve the quality of the data. (We'll tackle that in our next section.)\n",
        "\n",
        "In summary, our quality control analysis using FastQC and MultiQC has revealed that the ATAC-seq data has high overall quality, but contains trace amounts of adapter sequences which should be removed before further analysis."
      ],
      "metadata": {
        "id": "qHIn-ESYGDl6"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "#Trimming\n",
        "\n",
        "We will use the popular read-trimming software Trimmomatic (v. 0.39) to perform this step. For paired-end sequencing data, Trimmomatic expects two input files (one per sequencing direction), and outputs four different files: forward-paired reads, forward-unpaired reads, reverse-paired reads, and reverse-unpaired reads. The unpaired reads are generated when one read of a pair does not survive the trimming, but the other one does. For the downstream analysis, only paired reads are typically used.\n",
        "\n",
        "The following simple script will run Trimmomatic commands for each sample:\n",
        "\n"
      ],
      "metadata": {
        "id": "GMVJ3eXiGNOF"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "samples=['SRR10261591','SRR10261592','SRR10261593','SRR10261594','SRR10261595','SRR10261596']\n",
        "\n"
      ],
      "metadata": {
        "id": "r50dZUfOHfdz"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "### CHANGE PATH TO ADAPTERS AND OTHET PATHS\n",
        "for sample in samples:\n",
        "  workbench.run(tool='trimmomatic', full_command=f'trimmomatic PE -threads 10 /input/{sample}_1.fastq.gz /input/{sample}_2.fastq.gz /output/trimming/{sample}_tr_1P.fastq.gz /output/trimming/{sample}_tr_1U.fastq.gz /output/trimming/{sample}_tr_2P.fastq.gz /output/trimming/{sample}_tr_2U.fastq.gz ILLUMINACLIP:/input/NexteraPE-PE.fa:2:30:10 MINLEN:30 1>/output/trimming/{sample}_trim_log.txt 2>/output/trimming/{sample}_trim_log_error.txt' )"
      ],
      "metadata": {
        "id": "9d8Oa0B2Hjd1",
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "outputId": "1e429bda-dec7-430f-8687-c7424a9fc455"
      },
      "execution_count": null,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "+-----------------------+-------------+-----------+----------+------------------------------------------------+----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------+\n",
            "| Job ID                | Tool        |   Version | Status   | Get Logs                                       | Full Command                                                                                                                                                                                                                                                                                                                                                                                                               |\n",
            "+=======================+=============+===========+==========+================================================+============================================================================================================================================================================================================================================================================================================================================================================================================================+\n",
            "| trimmomatic6d95477707 | trimmomatic |      0.39 | Queued   | workbench.jobs('trimmomatic6d95477707').logs() | trimmomatic PE -threads 10 /input/SRR10261591_1.fastq.gz /input/SRR10261591_2.fastq.gz /output/trimming/SRR10261591_tr_1P.fastq.gz /output/trimming/SRR10261591_tr_1U.fastq.gz /output/trimming/SRR10261591_tr_2P.fastq.gz /output/trimming/SRR10261591_tr_2U.fastq.gz ILLUMINACLIP:/input/NexteraPE-PE.fa:2:30:10 MINLEN:30 1>/output/trimming/SRR10261591_trim_log.txt 2>/output/trimming/SRR10261591_trim_log_error.txt |\n",
            "+-----------------------+-------------+-----------+----------+------------------------------------------------+----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------+\n",
            "+-----------------------+-------------+-----------+----------+------------------------------------------------+----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------+\n",
            "| Job ID                | Tool        |   Version | Status   | Get Logs                                       | Full Command                                                                                                                                                                                                                                                                                                                                                                                                               |\n",
            "+=======================+=============+===========+==========+================================================+============================================================================================================================================================================================================================================================================================================================================================================================================================+\n",
            "| trimmomatic78066e70f0 | trimmomatic |      0.39 | Queued   | workbench.jobs('trimmomatic78066e70f0').logs() | trimmomatic PE -threads 10 /input/SRR10261592_1.fastq.gz /input/SRR10261592_2.fastq.gz /output/trimming/SRR10261592_tr_1P.fastq.gz /output/trimming/SRR10261592_tr_1U.fastq.gz /output/trimming/SRR10261592_tr_2P.fastq.gz /output/trimming/SRR10261592_tr_2U.fastq.gz ILLUMINACLIP:/input/NexteraPE-PE.fa:2:30:10 MINLEN:30 1>/output/trimming/SRR10261592_trim_log.txt 2>/output/trimming/SRR10261592_trim_log_error.txt |\n",
            "+-----------------------+-------------+-----------+----------+------------------------------------------------+----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------+\n",
            "+-----------------------+-------------+-----------+----------+------------------------------------------------+----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------+\n",
            "| Job ID                | Tool        |   Version | Status   | Get Logs                                       | Full Command                                                                                                                                                                                                                                                                                                                                                                                                               |\n",
            "+=======================+=============+===========+==========+================================================+============================================================================================================================================================================================================================================================================================================================================================================================================================+\n",
            "| trimmomatic6c132c7d29 | trimmomatic |      0.39 | Queued   | workbench.jobs('trimmomatic6c132c7d29').logs() | trimmomatic PE -threads 10 /input/SRR10261593_1.fastq.gz /input/SRR10261593_2.fastq.gz /output/trimming/SRR10261593_tr_1P.fastq.gz /output/trimming/SRR10261593_tr_1U.fastq.gz /output/trimming/SRR10261593_tr_2P.fastq.gz /output/trimming/SRR10261593_tr_2U.fastq.gz ILLUMINACLIP:/input/NexteraPE-PE.fa:2:30:10 MINLEN:30 1>/output/trimming/SRR10261593_trim_log.txt 2>/output/trimming/SRR10261593_trim_log_error.txt |\n",
            "+-----------------------+-------------+-----------+----------+------------------------------------------------+----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------+\n",
            "+-----------------------+-------------+-----------+----------+------------------------------------------------+----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------+\n",
            "| Job ID                | Tool        |   Version | Status   | Get Logs                                       | Full Command                                                                                                                                                                                                                                                                                                                                                                                                               |\n",
            "+=======================+=============+===========+==========+================================================+============================================================================================================================================================================================================================================================================================================================================================================================================================+\n",
            "| trimmomatic56abe5930a | trimmomatic |      0.39 | Queued   | workbench.jobs('trimmomatic56abe5930a').logs() | trimmomatic PE -threads 10 /input/SRR10261594_1.fastq.gz /input/SRR10261594_2.fastq.gz /output/trimming/SRR10261594_tr_1P.fastq.gz /output/trimming/SRR10261594_tr_1U.fastq.gz /output/trimming/SRR10261594_tr_2P.fastq.gz /output/trimming/SRR10261594_tr_2U.fastq.gz ILLUMINACLIP:/input/NexteraPE-PE.fa:2:30:10 MINLEN:30 1>/output/trimming/SRR10261594_trim_log.txt 2>/output/trimming/SRR10261594_trim_log_error.txt |\n",
            "+-----------------------+-------------+-----------+----------+------------------------------------------------+----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------+\n",
            "+-----------------------+-------------+-----------+----------+------------------------------------------------+----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------+\n",
            "| Job ID                | Tool        |   Version | Status   | Get Logs                                       | Full Command                                                                                                                                                                                                                                                                                                                                                                                                               |\n",
            "+=======================+=============+===========+==========+================================================+============================================================================================================================================================================================================================================================================================================================================================================================================================+\n",
            "| trimmomaticadc6244f4e | trimmomatic |      0.39 | Queued   | workbench.jobs('trimmomaticadc6244f4e').logs() | trimmomatic PE -threads 10 /input/SRR10261595_1.fastq.gz /input/SRR10261595_2.fastq.gz /output/trimming/SRR10261595_tr_1P.fastq.gz /output/trimming/SRR10261595_tr_1U.fastq.gz /output/trimming/SRR10261595_tr_2P.fastq.gz /output/trimming/SRR10261595_tr_2U.fastq.gz ILLUMINACLIP:/input/NexteraPE-PE.fa:2:30:10 MINLEN:30 1>/output/trimming/SRR10261595_trim_log.txt 2>/output/trimming/SRR10261595_trim_log_error.txt |\n",
            "+-----------------------+-------------+-----------+----------+------------------------------------------------+----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------+\n",
            "+-----------------------+-------------+-----------+----------+------------------------------------------------+----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------+\n",
            "| Job ID                | Tool        |   Version | Status   | Get Logs                                       | Full Command                                                                                                                                                                                                                                                                                                                                                                                                               |\n",
            "+=======================+=============+===========+==========+================================================+============================================================================================================================================================================================================================================================================================================================================================================================================================+\n",
            "| trimmomaticdf89c2a8a2 | trimmomatic |      0.39 | Queued   | workbench.jobs('trimmomaticdf89c2a8a2').logs() | trimmomatic PE -threads 10 /input/SRR10261596_1.fastq.gz /input/SRR10261596_2.fastq.gz /output/trimming/SRR10261596_tr_1P.fastq.gz /output/trimming/SRR10261596_tr_1U.fastq.gz /output/trimming/SRR10261596_tr_2P.fastq.gz /output/trimming/SRR10261596_tr_2U.fastq.gz ILLUMINACLIP:/input/NexteraPE-PE.fa:2:30:10 MINLEN:30 1>/output/trimming/SRR10261596_trim_log.txt 2>/output/trimming/SRR10261596_trim_log_error.txt |\n",
            "+-----------------------+-------------+-----------+----------+------------------------------------------------+----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------+\n"
          ]
        }
      ]
    },
    {
      "cell_type": "code",
      "source": [
        "workbench.jobs()"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "v0VO1im62rOk",
        "outputId": "839d119d-f51e-4dba-9391-5ea4eb063197"
      },
      "execution_count": null,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "+---------------+--------+-----------+-----------+----------------------------------------+---------------------------------------------------------------------------------------------------------------------------------------------------------------+\n",
            "| Job ID        | Tool   | Version   | Status    | Get Logs                               | Full Command                                                                                                                                                  |\n",
            "+===============+========+===========+===========+========================================+===============================================================================================================================================================+\n",
            "| bwad1a17c6138 | bwa    | 0.7.17    | Succeeded | workbench.jobs('bwad1a17c6138').logs() | bwa index -p /input/ref_gen/SC  /input/ref_gen/S_cerevisiae_masked_rRNA.fasta                                                                                 |\n",
            "+---------------+--------+-----------+-----------+----------------------------------------+---------------------------------------------------------------------------------------------------------------------------------------------------------------+\n",
            "| bwa06161cc243 | bwa    | 0.7.17    | Succeeded | workbench.jobs('bwa06161cc243').logs() | bwa index -p /input/ref_gen/hybrid  /input/ref_gen/hybrid.fasta                                                                                               |\n",
            "+---------------+--------+-----------+-----------+----------------------------------------+---------------------------------------------------------------------------------------------------------------------------------------------------------------+\n",
            "| bwa5a4579fa2f | bwa    | 0.7.17    | Failed    | workbench.jobs('bwa5a4579fa2f').logs() | bwa mem -t 2 /input/ref_gen/SC /output/trimming/SRR10261591_tr_1P.fastq.gz /output/trimming/SRR10261591_tr_2P.fastq.gz > /output/mapping/SRR10261591.sam      |\n",
            "+---------------+--------+-----------+-----------+----------------------------------------+---------------------------------------------------------------------------------------------------------------------------------------------------------------+\n",
            "| bwa669ea9aee4 | bwa    | 0.7.17    | Failed    | workbench.jobs('bwa669ea9aee4').logs() | bwa mem -t 2 /input/ref_gen/SC /output/trimming/SRR10261592_tr_1P.fastq.gz /output/trimming/SRR10261592_tr_2P.fastq.gz > /output/mapping/SRR10261592.sam      |\n",
            "+---------------+--------+-----------+-----------+----------------------------------------+---------------------------------------------------------------------------------------------------------------------------------------------------------------+\n",
            "| bwad46e18eab0 | bwa    | 0.7.17    | Failed    | workbench.jobs('bwad46e18eab0').logs() | bwa mem -t 2 /input/ref_gen/SC /output/trimming/SRR10261593_tr_1P.fastq.gz /output/trimming/SRR10261593_tr_2P.fastq.gz > /output/mapping/SRR10261593.sam      |\n",
            "+---------------+--------+-----------+-----------+----------------------------------------+---------------------------------------------------------------------------------------------------------------------------------------------------------------+\n",
            "| bwa2b3ea1e3fc | bwa    | 0.7.17    | Failed    | workbench.jobs('bwa2b3ea1e3fc').logs() | bwa mem -t 2 /input/ref_gen/hybrid /output/trimming/SRR10261594_tr_1P.fastq.gz /output/trimming/SRR10261594_tr_2P.fastq.gz > /output/mapping/SRR10261594.sam  |\n",
            "+---------------+--------+-----------+-----------+----------------------------------------+---------------------------------------------------------------------------------------------------------------------------------------------------------------+\n",
            "| bwa0c77fba1a4 | bwa    | 0.7.17    | Failed    | workbench.jobs('bwa0c77fba1a4').logs() | bwa mem -t 2 /input/ref_gen/hybrid /output/trimming/SRR10261595_tr_1P.fastq.gz /output/trimming/SRR10261595_tr_2P.fastq.gz > /output/mapping/SRR10261595.sam  |\n",
            "+---------------+--------+-----------+-----------+----------------------------------------+---------------------------------------------------------------------------------------------------------------------------------------------------------------+\n",
            "| bwaccf858c634 | bwa    | 0.7.17    | Failed    | workbench.jobs('bwaccf858c634').logs() | bwa mem -t 2 /input/ref_gen/hybrid /output/trimming/SRR10261596_tr_1P.fastq.gz /output/trimming/SRR10261596_tr_2P.fastq.gz > /output/mapping/SRR10261596.sam  |\n",
            "+---------------+--------+-----------+-----------+----------------------------------------+---------------------------------------------------------------------------------------------------------------------------------------------------------------+\n",
            "| bwac7b579fc56 | bwa    | 0.7.17    | Failed    | workbench.jobs('bwac7b579fc56').logs() | bwa mem -t 2 -o /output/mapping/SRR10261591.sam /input/ref_gen/SC /output/trimming/SRR10261591_tr_1P.fastq.gz /output/trimming/SRR10261591_tr_2P.fastq.gz     |\n",
            "+---------------+--------+-----------+-----------+----------------------------------------+---------------------------------------------------------------------------------------------------------------------------------------------------------------+\n",
            "| bwadabdeb260e | bwa    | 0.7.17    | Failed    | workbench.jobs('bwadabdeb260e').logs() | bwa mem -t 2 -o /output/mapping/SRR10261592.sam /input/ref_gen/SC /output/trimming/SRR10261592_tr_1P.fastq.gz /output/trimming/SRR10261592_tr_2P.fastq.gz     |\n",
            "+---------------+--------+-----------+-----------+----------------------------------------+---------------------------------------------------------------------------------------------------------------------------------------------------------------+\n",
            "| bwa3c7ce6e9b1 | bwa    | 0.7.17    | Failed    | workbench.jobs('bwa3c7ce6e9b1').logs() | bwa mem -t 2 -o /output/mapping/SRR10261593.sam /input/ref_gen/SC /output/trimming/SRR10261593_tr_1P.fastq.gz /output/trimming/SRR10261593_tr_2P.fastq.gz     |\n",
            "+---------------+--------+-----------+-----------+----------------------------------------+---------------------------------------------------------------------------------------------------------------------------------------------------------------+\n",
            "| bwae5e90cc8bf | bwa    | 0.7.17    | Failed    | workbench.jobs('bwae5e90cc8bf').logs() | bwa mem -t 2 -o /output/mapping/SRR10261594.sam /input/ref_gen/hybrid /output/trimming/SRR10261594_tr_1P.fastq.gz /output/trimming/SRR10261594_tr_2P.fastq.gz |\n",
            "+---------------+--------+-----------+-----------+----------------------------------------+---------------------------------------------------------------------------------------------------------------------------------------------------------------+\n",
            "| bwafe502d4704 | bwa    | 0.7.17    | Failed    | workbench.jobs('bwafe502d4704').logs() | bwa mem -t 2 -o /output/mapping/SRR10261595.sam /input/ref_gen/hybrid /output/trimming/SRR10261595_tr_1P.fastq.gz /output/trimming/SRR10261595_tr_2P.fastq.gz |\n",
            "+---------------+--------+-----------+-----------+----------------------------------------+---------------------------------------------------------------------------------------------------------------------------------------------------------------+\n",
            "| bwae1acc0fc17 | bwa    | 0.7.17    | Failed    | workbench.jobs('bwae1acc0fc17').logs() | bwa mem -t 2 -o /output/mapping/SRR10261596.sam /input/ref_gen/hybrid /output/trimming/SRR10261596_tr_1P.fastq.gz /output/trimming/SRR10261596_tr_2P.fastq.gz |\n",
            "+---------------+--------+-----------+-----------+----------------------------------------+---------------------------------------------------------------------------------------------------------------------------------------------------------------+\n",
            "| bwa65622cfc23 | bwa    | 0.7.17    | Succeeded | workbench.jobs('bwa65622cfc23').logs() | bwa mem -t 2 -o /output/mapping/SRR10261591.sam /input/ref_gen/SC /output/trimming/SRR10261591_tr_1P.fastq.gz /output/trimming/SRR10261591_tr_2P.fastq.gz     |\n",
            "+---------------+--------+-----------+-----------+----------------------------------------+---------------------------------------------------------------------------------------------------------------------------------------------------------------+\n",
            "| bwaddc01aaba5 | bwa    | 0.7.17    | Succeeded | workbench.jobs('bwaddc01aaba5').logs() | bwa mem -t 2 -o /output/mapping/SRR10261592.sam /input/ref_gen/SC /output/trimming/SRR10261592_tr_1P.fastq.gz /output/trimming/SRR10261592_tr_2P.fastq.gz     |\n",
            "+---------------+--------+-----------+-----------+----------------------------------------+---------------------------------------------------------------------------------------------------------------------------------------------------------------+\n",
            "| bwa17ee2fbd74 | bwa    | 0.7.17    | Succeeded | workbench.jobs('bwa17ee2fbd74').logs() | bwa mem -t 2 -o /output/mapping/SRR10261593.sam /input/ref_gen/SC /output/trimming/SRR10261593_tr_1P.fastq.gz /output/trimming/SRR10261593_tr_2P.fastq.gz     |\n",
            "+---------------+--------+-----------+-----------+----------------------------------------+---------------------------------------------------------------------------------------------------------------------------------------------------------------+\n",
            "| bwa14113f0940 | bwa    | 0.7.17    | Succeeded | workbench.jobs('bwa14113f0940').logs() | bwa mem -t 2 -o /output/mapping/SRR10261594.sam /input/ref_gen/hybrid /output/trimming/SRR10261594_tr_1P.fastq.gz /output/trimming/SRR10261594_tr_2P.fastq.gz |\n",
            "+---------------+--------+-----------+-----------+----------------------------------------+---------------------------------------------------------------------------------------------------------------------------------------------------------------+\n",
            "| bwa09750580c1 | bwa    | 0.7.17    | Succeeded | workbench.jobs('bwa09750580c1').logs() | bwa mem -t 2 -o /output/mapping/SRR10261595.sam /input/ref_gen/hybrid /output/trimming/SRR10261595_tr_1P.fastq.gz /output/trimming/SRR10261595_tr_2P.fastq.gz |\n",
            "+---------------+--------+-----------+-----------+----------------------------------------+---------------------------------------------------------------------------------------------------------------------------------------------------------------+\n",
            "| bwa0ca620e660 | bwa    | 0.7.17    | Succeeded | workbench.jobs('bwa0ca620e660').logs() | bwa mem -t 2 -o /output/mapping/SRR10261596.sam /input/ref_gen/hybrid /output/trimming/SRR10261596_tr_1P.fastq.gz /output/trimming/SRR10261596_tr_2P.fastq.gz |\n",
            "+---------------+--------+-----------+-----------+----------------------------------------+---------------------------------------------------------------------------------------------------------------------------------------------------------------+\n"
          ]
        }
      ]
    },
    {
      "cell_type": "markdown",
      "source": [
        "Before we continue, let's break down how these commands works. The scripts iterates over the sample names contained in the sample_ids variable  we created earlier, running a Trimmomatic command for each sample. In the output files (which end with tr_1P.fastq.gz), “P” stands for “paired”, and “U” stands for unpaired. The ILLUMINACLIP parameter navigates to the fasta file with Nextera adapters. The parameters 2:30:10 tell Trimmomatic to cut the adapters allowing a maximum of 2 mismatches, a score of 30 for paired reads, and score of 10 for single reads (see the Trimmomatic manual for more details about this score). And the MINLEN option tells Trimmomatic to drop reads shorter than 30 base pairs.\n",
        "\n",
        "\n",
        "After the trimming is finished, we will run FastQC and MultiQC to check the trimmed data. Let's use these commands:\n"
      ],
      "metadata": {
        "id": "aAUqDaxbH8C9"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "workbench.run(tool='fastqc', full_command='fastqc -t 10 /output/trimming/*gz -o /output/trimming/fastqc')"
      ],
      "metadata": {
        "id": "dDdq2Fc-Ijry",
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "outputId": "1a3c2cf4-7cb6-417a-8698-0aa82559fe93"
      },
      "execution_count": null,
      "outputs": [
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": [
              "job-multiqc-20230311-222015356516 running multiqc -o /mnt/gcs/output/trimming/multiqc /mnt/gcs/output/trimming/"
            ]
          },
          "metadata": {},
          "execution_count": 9
        }
      ]
    },
    {
      "cell_type": "code",
      "source": [
        "workbench.run(tool='multiqc', full_command='multiqc -o /output/trimming/multiqc /output/trimming/')"
      ],
      "metadata": {
        "id": "SLD3FKqYwLcw",
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "outputId": "46dbcf34-697e-41d5-ebe7-6b4582946a2b"
      },
      "execution_count": null,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "+-------------------+---------+-----------+----------+--------------------------------------------+-------------------------------------------------------+\n",
            "| Job ID            | Tool    |   Version | Status   | Get Logs                                   | Full Command                                          |\n",
            "+===================+=========+===========+==========+============================================+=======================================================+\n",
            "| multiqc8251bc1636 | multiqc |      1.14 | Queued   | workbench.jobs('multiqc8251bc1636').logs() | multiqc -o /output/trimming/multiqc /output/trimming/ |\n",
            "+-------------------+---------+-----------+----------+--------------------------------------------+-------------------------------------------------------+\n"
          ]
        }
      ]
    },
    {
      "cell_type": "code",
      "source": [
        "workbench.download(\"output/trimming/multiqc/multiqc_report.html\")"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "DduiGl3lxgwb",
        "outputId": "17ea39f2-ab76-4f90-cc65-456f4cca3d2f"
      },
      "execution_count": null,
      "outputs": [
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": [
              "{'download_url': 'https://storage.googleapis.com/atac-seq-20230406175308669217/output/trimming/multiqc/multiqc_report.html?X-Goog-Algorithm=GOOG4-RSA-SHA256&X-Goog-Credential=tiny-api%40tinybio-378903.iam.gserviceaccount.com%2F20230406%2Fauto%2Fstorage%2Fgoog4_request&X-Goog-Date=20230406T210010Z&X-Goog-Expires=3600&X-Goog-SignedHeaders=host&X-Goog-Signature=7aad8fa20b7a3337ff38af8afaf1411409a24bdfbb80b865ede0c11ff3fa5e409f1a154eed0c39054059396cc5670c8aedc02afd0bf6d03d8fa1ff719ef9fd5f4c5866d8fa00064e3fc745a6de3940ddcc3333a09a8b55c21036354e34a12a1dcce0ef2ec6bd10215227833436190e6b40f9df91c30ba517f5c639d7bb1dce350d9c1cf0401671a18e91a93640c88c62d11515c522451c19af0a32b1098b1ceaa2b7be4336386945a69510980df9274ad15b9f1428b79add2e8ce8463a77d7e66a552ab533eeaed4e92cf0b38f786245164314f369cc47f1c4f3343e42b67cf93f37ffd44470d2690b3af7e50ea500d4b7cad44818c0eade02e872ed883c5199'}"
            ]
          },
          "metadata": {},
          "execution_count": 27
        }
      ]
    },
    {
      "cell_type": "markdown",
      "source": [
        "MultiQC will now create another html report – but this time, it will also contain statistics around our Trimmomatic output (Fig. 6):\n",
        "\n",
        "![](https://drive.google.com/uc?export=view&id=11dZ0I0bywQM7MMnTC5GyjnZPGQUleJ_0)\n",
        "\n",
        "*Fig. 6: Summary statistics of trimming results reported by MultiQC.*\n",
        "\n",
        "\n",
        "\n",
        "This report shows that only a small percentage of reads were dropped due to trimming, suggesting that our trimming has gently removed the few unnecessary reads. It also shows that our trimming of adapters was successful and there are no samples found with adapter contamination.\n",
        "\n",
        "In summary, we have just successfully trimmed our ATAC-seq data, removing the unnecessary adapter sequences and remaining short reads while preserving the majority of high-quality data for downstream analysis."
      ],
      "metadata": {
        "id": "38hYkI4HJn0H"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "#Mapping & Related Analyses\n",
        "####Run time: See below for each type of analysis\n",
        "\n",
        "We have arrived at the longest and most computationally demanding step of this workflow: sequencing read mapping to corresponding reference genomes, which will generate the open chromatin profiles of samples originating from the ATAC-seq experiment. This step consists of two parts – reference genome preparations and read mapping itself (with some auxiliary analyses)."
      ],
      "metadata": {
        "id": "GnxmFWfmKLyT"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "###Reference genome preparation\n",
        "\n",
        "\n",
        "In order to map the sequencing data to the reference genomes, our reference genomes must first be obtained, processed and indexed. We aready downloaded the reference genomes of both yeast, so now we will need to process them. \n",
        "\n",
        "\n",
        "One peculiarity of the S. cerevisiae genome is that it contains two clusters of ribosomal (rRNA) genes on chromosome XII (region 450915-469179 and 489349-490611). This can result in problematic alignments during the mapping step, especially in ATAC-seq data analysis. To avoid this, we will mask these regions of chromosome XII using bedtools v. 2.30.0 (Note that this masking step is specific to S. cerevisiae, and may not be necessary for other organisms.).\n",
        "\n",
        "To do this, we first have created a bed file named mask_rRNA.bed containing these regions:\n",
        "\n",
        "\n",
        "\n",
        "```\n",
        "XII    450915    469179\n",
        "XII    489349    490611\n",
        "```\n",
        "\n"
      ],
      "metadata": {
        "id": "fjS1xVmAKe8F"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "Then we run the command below which will create the masked genome file S_cerevisiae_masked_rRNA.fasta."
      ],
      "metadata": {
        "id": "8L9VU3GzLM8f"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "\n",
        "workbench.run(tool='bedtools', full_command='bedtools maskfasta -fi /input/ref_gen/Saccharomyces_cerevisiae.R64-1-1.dna.toplevel.fa -bed /input/ref_gen/mask_rRNA.bed -fo /input/ref_gen/S_cerevisiae_masked_rRNA.fasta')"
      ],
      "metadata": {
        "id": "OM2mcrpVLMUm",
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "outputId": "07cb0a62-cc9e-40d6-d0a1-4275de8fb21e"
      },
      "execution_count": null,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "+--------------------+----------+-----------+----------+---------------------------------------------+----------------------------------------------------------------------------------------------------------------------------------------------------------------------------+\n",
            "| Job ID             | Tool     | Version   | Status   | Get Logs                                    | Full Command                                                                                                                                                               |\n",
            "+====================+==========+===========+==========+=============================================+============================================================================================================================================================================+\n",
            "| bedtools5561087917 | bedtools | 2.30.0    | Queued   | workbench.jobs('bedtools5561087917').logs() | bedtools maskfasta -fi /input/ref_gen/Saccharomyces_cerevisiae.R64-1-1.dna.toplevel.fa -bed /input/ref_gen/mask_rRNA.bed -fo /input/ref_gen/S_cerevisiae_masked_rRNA.fasta |\n",
            "+--------------------+----------+-----------+----------+---------------------------------------------+----------------------------------------------------------------------------------------------------------------------------------------------------------------------------+\n"
          ]
        }
      ]
    },
    {
      "cell_type": "code",
      "source": [
        "workbench.jobs()\n"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "a3Doa10L6NG3",
        "outputId": "8e3478ee-680b-4d28-aa10-78ffa5cd4439"
      },
      "execution_count": null,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "+---------------+--------+-----------+-----------+----------------------------------------+---------------------------------------------------------------------------------------------------------------------------------------------------------------+\n",
            "| Job ID        | Tool   | Version   | Status    | Get Logs                               | Full Command                                                                                                                                                  |\n",
            "+===============+========+===========+===========+========================================+===============================================================================================================================================================+\n",
            "| bwad1a17c6138 | bwa    | 0.7.17    | Succeeded | workbench.jobs('bwad1a17c6138').logs() | bwa index -p /input/ref_gen/SC  /input/ref_gen/S_cerevisiae_masked_rRNA.fasta                                                                                 |\n",
            "+---------------+--------+-----------+-----------+----------------------------------------+---------------------------------------------------------------------------------------------------------------------------------------------------------------+\n",
            "| bwa06161cc243 | bwa    | 0.7.17    | Succeeded | workbench.jobs('bwa06161cc243').logs() | bwa index -p /input/ref_gen/hybrid  /input/ref_gen/hybrid.fasta                                                                                               |\n",
            "+---------------+--------+-----------+-----------+----------------------------------------+---------------------------------------------------------------------------------------------------------------------------------------------------------------+\n",
            "| bwa5a4579fa2f | bwa    | 0.7.17    | Failed    | workbench.jobs('bwa5a4579fa2f').logs() | bwa mem -t 2 /input/ref_gen/SC /output/trimming/SRR10261591_tr_1P.fastq.gz /output/trimming/SRR10261591_tr_2P.fastq.gz > /output/mapping/SRR10261591.sam      |\n",
            "+---------------+--------+-----------+-----------+----------------------------------------+---------------------------------------------------------------------------------------------------------------------------------------------------------------+\n",
            "| bwa669ea9aee4 | bwa    | 0.7.17    | Failed    | workbench.jobs('bwa669ea9aee4').logs() | bwa mem -t 2 /input/ref_gen/SC /output/trimming/SRR10261592_tr_1P.fastq.gz /output/trimming/SRR10261592_tr_2P.fastq.gz > /output/mapping/SRR10261592.sam      |\n",
            "+---------------+--------+-----------+-----------+----------------------------------------+---------------------------------------------------------------------------------------------------------------------------------------------------------------+\n",
            "| bwad46e18eab0 | bwa    | 0.7.17    | Failed    | workbench.jobs('bwad46e18eab0').logs() | bwa mem -t 2 /input/ref_gen/SC /output/trimming/SRR10261593_tr_1P.fastq.gz /output/trimming/SRR10261593_tr_2P.fastq.gz > /output/mapping/SRR10261593.sam      |\n",
            "+---------------+--------+-----------+-----------+----------------------------------------+---------------------------------------------------------------------------------------------------------------------------------------------------------------+\n",
            "| bwa2b3ea1e3fc | bwa    | 0.7.17    | Failed    | workbench.jobs('bwa2b3ea1e3fc').logs() | bwa mem -t 2 /input/ref_gen/hybrid /output/trimming/SRR10261594_tr_1P.fastq.gz /output/trimming/SRR10261594_tr_2P.fastq.gz > /output/mapping/SRR10261594.sam  |\n",
            "+---------------+--------+-----------+-----------+----------------------------------------+---------------------------------------------------------------------------------------------------------------------------------------------------------------+\n",
            "| bwa0c77fba1a4 | bwa    | 0.7.17    | Failed    | workbench.jobs('bwa0c77fba1a4').logs() | bwa mem -t 2 /input/ref_gen/hybrid /output/trimming/SRR10261595_tr_1P.fastq.gz /output/trimming/SRR10261595_tr_2P.fastq.gz > /output/mapping/SRR10261595.sam  |\n",
            "+---------------+--------+-----------+-----------+----------------------------------------+---------------------------------------------------------------------------------------------------------------------------------------------------------------+\n",
            "| bwaccf858c634 | bwa    | 0.7.17    | Failed    | workbench.jobs('bwaccf858c634').logs() | bwa mem -t 2 /input/ref_gen/hybrid /output/trimming/SRR10261596_tr_1P.fastq.gz /output/trimming/SRR10261596_tr_2P.fastq.gz > /output/mapping/SRR10261596.sam  |\n",
            "+---------------+--------+-----------+-----------+----------------------------------------+---------------------------------------------------------------------------------------------------------------------------------------------------------------+\n",
            "| bwac7b579fc56 | bwa    | 0.7.17    | Failed    | workbench.jobs('bwac7b579fc56').logs() | bwa mem -t 2 -o /output/mapping/SRR10261591.sam /input/ref_gen/SC /output/trimming/SRR10261591_tr_1P.fastq.gz /output/trimming/SRR10261591_tr_2P.fastq.gz     |\n",
            "+---------------+--------+-----------+-----------+----------------------------------------+---------------------------------------------------------------------------------------------------------------------------------------------------------------+\n",
            "| bwadabdeb260e | bwa    | 0.7.17    | Failed    | workbench.jobs('bwadabdeb260e').logs() | bwa mem -t 2 -o /output/mapping/SRR10261592.sam /input/ref_gen/SC /output/trimming/SRR10261592_tr_1P.fastq.gz /output/trimming/SRR10261592_tr_2P.fastq.gz     |\n",
            "+---------------+--------+-----------+-----------+----------------------------------------+---------------------------------------------------------------------------------------------------------------------------------------------------------------+\n",
            "| bwa3c7ce6e9b1 | bwa    | 0.7.17    | Failed    | workbench.jobs('bwa3c7ce6e9b1').logs() | bwa mem -t 2 -o /output/mapping/SRR10261593.sam /input/ref_gen/SC /output/trimming/SRR10261593_tr_1P.fastq.gz /output/trimming/SRR10261593_tr_2P.fastq.gz     |\n",
            "+---------------+--------+-----------+-----------+----------------------------------------+---------------------------------------------------------------------------------------------------------------------------------------------------------------+\n",
            "| bwae5e90cc8bf | bwa    | 0.7.17    | Failed    | workbench.jobs('bwae5e90cc8bf').logs() | bwa mem -t 2 -o /output/mapping/SRR10261594.sam /input/ref_gen/hybrid /output/trimming/SRR10261594_tr_1P.fastq.gz /output/trimming/SRR10261594_tr_2P.fastq.gz |\n",
            "+---------------+--------+-----------+-----------+----------------------------------------+---------------------------------------------------------------------------------------------------------------------------------------------------------------+\n",
            "| bwafe502d4704 | bwa    | 0.7.17    | Failed    | workbench.jobs('bwafe502d4704').logs() | bwa mem -t 2 -o /output/mapping/SRR10261595.sam /input/ref_gen/hybrid /output/trimming/SRR10261595_tr_1P.fastq.gz /output/trimming/SRR10261595_tr_2P.fastq.gz |\n",
            "+---------------+--------+-----------+-----------+----------------------------------------+---------------------------------------------------------------------------------------------------------------------------------------------------------------+\n",
            "| bwae1acc0fc17 | bwa    | 0.7.17    | Failed    | workbench.jobs('bwae1acc0fc17').logs() | bwa mem -t 2 -o /output/mapping/SRR10261596.sam /input/ref_gen/hybrid /output/trimming/SRR10261596_tr_1P.fastq.gz /output/trimming/SRR10261596_tr_2P.fastq.gz |\n",
            "+---------------+--------+-----------+-----------+----------------------------------------+---------------------------------------------------------------------------------------------------------------------------------------------------------------+\n",
            "| bwa65622cfc23 | bwa    | 0.7.17    | Succeeded | workbench.jobs('bwa65622cfc23').logs() | bwa mem -t 2 -o /output/mapping/SRR10261591.sam /input/ref_gen/SC /output/trimming/SRR10261591_tr_1P.fastq.gz /output/trimming/SRR10261591_tr_2P.fastq.gz     |\n",
            "+---------------+--------+-----------+-----------+----------------------------------------+---------------------------------------------------------------------------------------------------------------------------------------------------------------+\n",
            "| bwaddc01aaba5 | bwa    | 0.7.17    | Succeeded | workbench.jobs('bwaddc01aaba5').logs() | bwa mem -t 2 -o /output/mapping/SRR10261592.sam /input/ref_gen/SC /output/trimming/SRR10261592_tr_1P.fastq.gz /output/trimming/SRR10261592_tr_2P.fastq.gz     |\n",
            "+---------------+--------+-----------+-----------+----------------------------------------+---------------------------------------------------------------------------------------------------------------------------------------------------------------+\n",
            "| bwa17ee2fbd74 | bwa    | 0.7.17    | Succeeded | workbench.jobs('bwa17ee2fbd74').logs() | bwa mem -t 2 -o /output/mapping/SRR10261593.sam /input/ref_gen/SC /output/trimming/SRR10261593_tr_1P.fastq.gz /output/trimming/SRR10261593_tr_2P.fastq.gz     |\n",
            "+---------------+--------+-----------+-----------+----------------------------------------+---------------------------------------------------------------------------------------------------------------------------------------------------------------+\n",
            "| bwa14113f0940 | bwa    | 0.7.17    | Succeeded | workbench.jobs('bwa14113f0940').logs() | bwa mem -t 2 -o /output/mapping/SRR10261594.sam /input/ref_gen/hybrid /output/trimming/SRR10261594_tr_1P.fastq.gz /output/trimming/SRR10261594_tr_2P.fastq.gz |\n",
            "+---------------+--------+-----------+-----------+----------------------------------------+---------------------------------------------------------------------------------------------------------------------------------------------------------------+\n",
            "| bwa09750580c1 | bwa    | 0.7.17    | Succeeded | workbench.jobs('bwa09750580c1').logs() | bwa mem -t 2 -o /output/mapping/SRR10261595.sam /input/ref_gen/hybrid /output/trimming/SRR10261595_tr_1P.fastq.gz /output/trimming/SRR10261595_tr_2P.fastq.gz |\n",
            "+---------------+--------+-----------+-----------+----------------------------------------+---------------------------------------------------------------------------------------------------------------------------------------------------------------+\n",
            "| bwa0ca620e660 | bwa    | 0.7.17    | Succeeded | workbench.jobs('bwa0ca620e660').logs() | bwa mem -t 2 -o /output/mapping/SRR10261596.sam /input/ref_gen/hybrid /output/trimming/SRR10261596_tr_1P.fastq.gz /output/trimming/SRR10261596_tr_2P.fastq.gz |\n",
            "+---------------+--------+-----------+-----------+----------------------------------------+---------------------------------------------------------------------------------------------------------------------------------------------------------------+\n"
          ]
        }
      ]
    },
    {
      "cell_type": "markdown",
      "source": [
        "For S. uvarum, we don't have to any masking.\n",
        "\n",
        "Because the samples contain data from the S. cerevisiae parent and the hybrid S. cerevisiae x S. uvarum, we will need to artificially create the genome of the hybrid by concatenating the genomes of the two species in one file:"
      ],
      "metadata": {
        "id": "T52kSlzOMV5G"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "#cat S_cerevisiae_masked_rRNA.fasta S_uvarum.fasta > hybrid.fasta\n"
      ],
      "metadata": {
        "id": "6BdcdR_YMeCa"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "The fastq data of S. cerevisiae will be mapped to the parental genome, while the data from the hybrid will be mapped to the hybrid genome. This will separate the reads between the two organisms and further select only the data from S. cerevisiae.\n",
        "\n",
        "Once the final genome files S_cerevisiae_masked_rRNA.fasta and hybrid.fasta are ready, they can be indexed for further read mapping.\n",
        "\n",
        "The read mapping (and thus indexing) will be performed using the popular short-read mapping software BWA v. 0.7.17-r1188 with MEM algorithm. To index the genomes with BWA, lets do:"
      ],
      "metadata": {
        "id": "FUMIF7ImMixn"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "#wa index -p SC  S_cerevisiae_masked_rRNA.fasta\n",
        "#bwa index -p hybrid hybrid.fasta\n",
        "\n",
        "workbench.run(tool='bwa', full_command='bwa index -p /input/ref_gen/SC  /input/ref_gen/S_cerevisiae_masked_rRNA.fasta')\n",
        "workbench.run(tool='bwa', full_command='bwa index -p /input/ref_gen/hybrid  /input/ref_gen/hybrid.fasta')"
      ],
      "metadata": {
        "id": "BCNlX-ljMsP5",
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "outputId": "4cd0932b-638d-4cd6-b66b-898d257e9d09"
      },
      "execution_count": null,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "+---------------+--------+-----------+----------+----------------------------------------+-------------------------------------------------------------------------------+\n",
            "| Job ID        | Tool   | Version   | Status   | Get Logs                               | Full Command                                                                  |\n",
            "+===============+========+===========+==========+========================================+===============================================================================+\n",
            "| bwad1a17c6138 | bwa    | 0.7.17    | Queued   | workbench.jobs('bwad1a17c6138').logs() | bwa index -p /input/ref_gen/SC  /input/ref_gen/S_cerevisiae_masked_rRNA.fasta |\n",
            "+---------------+--------+-----------+----------+----------------------------------------+-------------------------------------------------------------------------------+\n",
            "+---------------+--------+-----------+----------+----------------------------------------+-----------------------------------------------------------------+\n",
            "| Job ID        | Tool   | Version   | Status   | Get Logs                               | Full Command                                                    |\n",
            "+===============+========+===========+==========+========================================+=================================================================+\n",
            "| bwa06161cc243 | bwa    | 0.7.17    | Queued   | workbench.jobs('bwa06161cc243').logs() | bwa index -p /input/ref_gen/hybrid  /input/ref_gen/hybrid.fasta |\n",
            "+---------------+--------+-----------+----------+----------------------------------------+-----------------------------------------------------------------+\n"
          ]
        }
      ]
    },
    {
      "cell_type": "markdown",
      "source": [
        "Now that we have indexed our reference genomes, we can map our fastq data to them.\n",
        "\n",
        "\n"
      ],
      "metadata": {
        "id": "QxSdCTawMv2Y"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "#Read mapping and related analyses\n",
        "\n",
        "Now we will map the reads of the trimmed fastq files to the corresponding reference genomes. We will also sort and index our alignment files, calculate read mapping statistics, mark the duplicated reads, and subset the S. cerevisiae data from the hybrid.\n",
        "\n",
        "For simplicity, the entire analysis is divided into two parts:\n",
        "\n",
        "*  Read mapping, sorting, indexing, and calculation of mapping statistics\n",
        "\n",
        "*   Subsetting of S. cerevisiae from hybrid, marking duplicated reads, and indexing and calculating statistics of the remaining data.\n",
        "\n",
        "\n",
        "\n",
        "For the first part, the following commands will be used:"
      ],
      "metadata": {
        "id": "2etljYsnM5_5"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "SC_samples = ['SRR10261591','SRR10261592','SRR10261593']\n",
        "hybrid_samples = ['SRR10261594','SRR10261595','SRR10261596']\n",
        "\n",
        "\n",
        "for sample_SC in SC_samples:\n",
        "  workbench.run(tool='bwa', full_command=f'bwa mem -t 2 -o /output/mapping/{sample_SC}.sam /input/ref_gen/SC /output/trimming/{sample_SC}_tr_1P.fastq.gz /output/trimming/{sample_SC}_tr_2P.fastq.gz')\n",
        "\n",
        "for sample_hybrid in hybrid_samples:\n",
        "  workbench.run(tool='bwa', full_command=f'bwa mem -t 2 -o /output/mapping/{sample_hybrid}.sam /input/ref_gen/hybrid /output/trimming/{sample_hybrid}_tr_1P.fastq.gz /output/trimming/{sample_hybrid}_tr_2P.fastq.gz')\n",
        "\n"
      ],
      "metadata": {
        "id": "XSClUzzQFitv",
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "outputId": "38f23144-bd47-46c3-a1bc-913cae0ae2ae"
      },
      "execution_count": null,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "+---------------+--------+-----------+----------+----------------------------------------+-----------------------------------------------------------------------------------------------------------------------------------------------------------+\n",
            "| Job ID        | Tool   | Version   | Status   | Get Logs                               | Full Command                                                                                                                                              |\n",
            "+===============+========+===========+==========+========================================+===========================================================================================================================================================+\n",
            "| bwa65622cfc23 | bwa    | 0.7.17    | Queued   | workbench.jobs('bwa65622cfc23').logs() | bwa mem -t 2 -o /output/mapping/SRR10261591.sam /input/ref_gen/SC /output/trimming/SRR10261591_tr_1P.fastq.gz /output/trimming/SRR10261591_tr_2P.fastq.gz |\n",
            "+---------------+--------+-----------+----------+----------------------------------------+-----------------------------------------------------------------------------------------------------------------------------------------------------------+\n",
            "+---------------+--------+-----------+----------+----------------------------------------+-----------------------------------------------------------------------------------------------------------------------------------------------------------+\n",
            "| Job ID        | Tool   | Version   | Status   | Get Logs                               | Full Command                                                                                                                                              |\n",
            "+===============+========+===========+==========+========================================+===========================================================================================================================================================+\n",
            "| bwaddc01aaba5 | bwa    | 0.7.17    | Queued   | workbench.jobs('bwaddc01aaba5').logs() | bwa mem -t 2 -o /output/mapping/SRR10261592.sam /input/ref_gen/SC /output/trimming/SRR10261592_tr_1P.fastq.gz /output/trimming/SRR10261592_tr_2P.fastq.gz |\n",
            "+---------------+--------+-----------+----------+----------------------------------------+-----------------------------------------------------------------------------------------------------------------------------------------------------------+\n",
            "+---------------+--------+-----------+----------+----------------------------------------+-----------------------------------------------------------------------------------------------------------------------------------------------------------+\n",
            "| Job ID        | Tool   | Version   | Status   | Get Logs                               | Full Command                                                                                                                                              |\n",
            "+===============+========+===========+==========+========================================+===========================================================================================================================================================+\n",
            "| bwa17ee2fbd74 | bwa    | 0.7.17    | Queued   | workbench.jobs('bwa17ee2fbd74').logs() | bwa mem -t 2 -o /output/mapping/SRR10261593.sam /input/ref_gen/SC /output/trimming/SRR10261593_tr_1P.fastq.gz /output/trimming/SRR10261593_tr_2P.fastq.gz |\n",
            "+---------------+--------+-----------+----------+----------------------------------------+-----------------------------------------------------------------------------------------------------------------------------------------------------------+\n",
            "+---------------+--------+-----------+----------+----------------------------------------+---------------------------------------------------------------------------------------------------------------------------------------------------------------+\n",
            "| Job ID        | Tool   | Version   | Status   | Get Logs                               | Full Command                                                                                                                                                  |\n",
            "+===============+========+===========+==========+========================================+===============================================================================================================================================================+\n",
            "| bwa14113f0940 | bwa    | 0.7.17    | Queued   | workbench.jobs('bwa14113f0940').logs() | bwa mem -t 2 -o /output/mapping/SRR10261594.sam /input/ref_gen/hybrid /output/trimming/SRR10261594_tr_1P.fastq.gz /output/trimming/SRR10261594_tr_2P.fastq.gz |\n",
            "+---------------+--------+-----------+----------+----------------------------------------+---------------------------------------------------------------------------------------------------------------------------------------------------------------+\n",
            "+---------------+--------+-----------+----------+----------------------------------------+---------------------------------------------------------------------------------------------------------------------------------------------------------------+\n",
            "| Job ID        | Tool   | Version   | Status   | Get Logs                               | Full Command                                                                                                                                                  |\n",
            "+===============+========+===========+==========+========================================+===============================================================================================================================================================+\n",
            "| bwa09750580c1 | bwa    | 0.7.17    | Queued   | workbench.jobs('bwa09750580c1').logs() | bwa mem -t 2 -o /output/mapping/SRR10261595.sam /input/ref_gen/hybrid /output/trimming/SRR10261595_tr_1P.fastq.gz /output/trimming/SRR10261595_tr_2P.fastq.gz |\n",
            "+---------------+--------+-----------+----------+----------------------------------------+---------------------------------------------------------------------------------------------------------------------------------------------------------------+\n",
            "+---------------+--------+-----------+----------+----------------------------------------+---------------------------------------------------------------------------------------------------------------------------------------------------------------+\n",
            "| Job ID        | Tool   | Version   | Status   | Get Logs                               | Full Command                                                                                                                                                  |\n",
            "+===============+========+===========+==========+========================================+===============================================================================================================================================================+\n",
            "| bwa0ca620e660 | bwa    | 0.7.17    | Queued   | workbench.jobs('bwa0ca620e660').logs() | bwa mem -t 2 -o /output/mapping/SRR10261596.sam /input/ref_gen/hybrid /output/trimming/SRR10261596_tr_1P.fastq.gz /output/trimming/SRR10261596_tr_2P.fastq.gz |\n",
            "+---------------+--------+-----------+----------+----------------------------------------+---------------------------------------------------------------------------------------------------------------------------------------------------------------+\n"
          ]
        }
      ]
    },
    {
      "cell_type": "code",
      "source": [
        "for sample in samples:\n",
        "  workbench.run(tool='samtools', full_command=f'samtools view -S -b /output/mapping/{sample}.sam > /output/mapping/{sample}.bam')\n",
        "#workbench.run(tool='samtools', full_command='samtools view -S -b /mnt/gcs/output/mapping/SRR10261592.sam > /mnt/gcs/output/mapping/SRR10261592.bam')\n"
      ],
      "metadata": {
        "id": "T1PZKdHvF8ZV",
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "outputId": "5ad710b3-1ee2-4cd1-c322-41f7697e7c49"
      },
      "execution_count": null,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "+--------------------+----------+-----------+----------+---------------------------------------------+---------------------------------------------------------------------------------------+\n",
            "| Job ID             | Tool     | Version   | Status   | Get Logs                                    | Full Command                                                                          |\n",
            "+====================+==========+===========+==========+=============================================+=======================================================================================+\n",
            "| samtools39ac7adaa8 | samtools | 1.16.1    | Queued   | workbench.jobs('samtools39ac7adaa8').logs() | samtools view -S -b /output/mapping/SRR10261591.sam > /output/mapping/SRR10261591.bam |\n",
            "+--------------------+----------+-----------+----------+---------------------------------------------+---------------------------------------------------------------------------------------+\n",
            "+--------------------+----------+-----------+----------+---------------------------------------------+---------------------------------------------------------------------------------------+\n",
            "| Job ID             | Tool     | Version   | Status   | Get Logs                                    | Full Command                                                                          |\n",
            "+====================+==========+===========+==========+=============================================+=======================================================================================+\n",
            "| samtoolse5fc956a60 | samtools | 1.16.1    | Queued   | workbench.jobs('samtoolse5fc956a60').logs() | samtools view -S -b /output/mapping/SRR10261592.sam > /output/mapping/SRR10261592.bam |\n",
            "+--------------------+----------+-----------+----------+---------------------------------------------+---------------------------------------------------------------------------------------+\n",
            "+--------------------+----------+-----------+----------+---------------------------------------------+---------------------------------------------------------------------------------------+\n",
            "| Job ID             | Tool     | Version   | Status   | Get Logs                                    | Full Command                                                                          |\n",
            "+====================+==========+===========+==========+=============================================+=======================================================================================+\n",
            "| samtoolse3d52c1637 | samtools | 1.16.1    | Queued   | workbench.jobs('samtoolse3d52c1637').logs() | samtools view -S -b /output/mapping/SRR10261593.sam > /output/mapping/SRR10261593.bam |\n",
            "+--------------------+----------+-----------+----------+---------------------------------------------+---------------------------------------------------------------------------------------+\n",
            "+--------------------+----------+-----------+----------+---------------------------------------------+---------------------------------------------------------------------------------------+\n",
            "| Job ID             | Tool     | Version   | Status   | Get Logs                                    | Full Command                                                                          |\n",
            "+====================+==========+===========+==========+=============================================+=======================================================================================+\n",
            "| samtools4f7450d6e6 | samtools | 1.16.1    | Queued   | workbench.jobs('samtools4f7450d6e6').logs() | samtools view -S -b /output/mapping/SRR10261594.sam > /output/mapping/SRR10261594.bam |\n",
            "+--------------------+----------+-----------+----------+---------------------------------------------+---------------------------------------------------------------------------------------+\n",
            "+--------------------+----------+-----------+----------+---------------------------------------------+---------------------------------------------------------------------------------------+\n",
            "| Job ID             | Tool     | Version   | Status   | Get Logs                                    | Full Command                                                                          |\n",
            "+====================+==========+===========+==========+=============================================+=======================================================================================+\n",
            "| samtools5ca2c857db | samtools | 1.16.1    | Queued   | workbench.jobs('samtools5ca2c857db').logs() | samtools view -S -b /output/mapping/SRR10261595.sam > /output/mapping/SRR10261595.bam |\n",
            "+--------------------+----------+-----------+----------+---------------------------------------------+---------------------------------------------------------------------------------------+\n",
            "+--------------------+----------+-----------+----------+---------------------------------------------+---------------------------------------------------------------------------------------+\n",
            "| Job ID             | Tool     | Version   | Status   | Get Logs                                    | Full Command                                                                          |\n",
            "+====================+==========+===========+==========+=============================================+=======================================================================================+\n",
            "| samtoolsf08c6dd83e | samtools | 1.16.1    | Queued   | workbench.jobs('samtoolsf08c6dd83e').logs() | samtools view -S -b /output/mapping/SRR10261596.sam > /output/mapping/SRR10261596.bam |\n",
            "+--------------------+----------+-----------+----------+---------------------------------------------+---------------------------------------------------------------------------------------+\n"
          ]
        }
      ]
    },
    {
      "cell_type": "code",
      "source": [
        "for sample in samples:\n",
        "  workbench.run(tool='samtools', full_command=f'samtools sort -o /output/mapping/{sample}_sorted.bam -@2 /output/mapping/{sample}.bam')\n",
        "#workbench.run(tool='samtools', full_command='samtools sort -o /output/mapping/SRR10261592_sorted.bam -@2 /output/mapping/SRR10261592.bam')\n"
      ],
      "metadata": {
        "id": "zA_CLvkcPBHo",
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "outputId": "aa8b8543-47a1-4ec3-fc97-24b461754d5d"
      },
      "execution_count": null,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "+--------------------+----------+-----------+----------+---------------------------------------------+---------------------------------------------------------------------------------------------+\n",
            "| Job ID             | Tool     | Version   | Status   | Get Logs                                    | Full Command                                                                                |\n",
            "+====================+==========+===========+==========+=============================================+=============================================================================================+\n",
            "| samtools80d8458866 | samtools | 1.16.1    | Queued   | workbench.jobs('samtools80d8458866').logs() | samtools sort -o /output/mapping/SRR10261591_sorted.bam -@2 /output/mapping/SRR10261591.bam |\n",
            "+--------------------+----------+-----------+----------+---------------------------------------------+---------------------------------------------------------------------------------------------+\n",
            "+--------------------+----------+-----------+----------+---------------------------------------------+---------------------------------------------------------------------------------------------+\n",
            "| Job ID             | Tool     | Version   | Status   | Get Logs                                    | Full Command                                                                                |\n",
            "+====================+==========+===========+==========+=============================================+=============================================================================================+\n",
            "| samtoolsb8fc7810be | samtools | 1.16.1    | Queued   | workbench.jobs('samtoolsb8fc7810be').logs() | samtools sort -o /output/mapping/SRR10261592_sorted.bam -@2 /output/mapping/SRR10261592.bam |\n",
            "+--------------------+----------+-----------+----------+---------------------------------------------+---------------------------------------------------------------------------------------------+\n",
            "+--------------------+----------+-----------+----------+---------------------------------------------+---------------------------------------------------------------------------------------------+\n",
            "| Job ID             | Tool     | Version   | Status   | Get Logs                                    | Full Command                                                                                |\n",
            "+====================+==========+===========+==========+=============================================+=============================================================================================+\n",
            "| samtools6526df3386 | samtools | 1.16.1    | Queued   | workbench.jobs('samtools6526df3386').logs() | samtools sort -o /output/mapping/SRR10261593_sorted.bam -@2 /output/mapping/SRR10261593.bam |\n",
            "+--------------------+----------+-----------+----------+---------------------------------------------+---------------------------------------------------------------------------------------------+\n",
            "+--------------------+----------+-----------+----------+---------------------------------------------+---------------------------------------------------------------------------------------------+\n",
            "| Job ID             | Tool     | Version   | Status   | Get Logs                                    | Full Command                                                                                |\n",
            "+====================+==========+===========+==========+=============================================+=============================================================================================+\n",
            "| samtools0dc6f9fdc7 | samtools | 1.16.1    | Queued   | workbench.jobs('samtools0dc6f9fdc7').logs() | samtools sort -o /output/mapping/SRR10261594_sorted.bam -@2 /output/mapping/SRR10261594.bam |\n",
            "+--------------------+----------+-----------+----------+---------------------------------------------+---------------------------------------------------------------------------------------------+\n",
            "+--------------------+----------+-----------+----------+---------------------------------------------+---------------------------------------------------------------------------------------------+\n",
            "| Job ID             | Tool     | Version   | Status   | Get Logs                                    | Full Command                                                                                |\n",
            "+====================+==========+===========+==========+=============================================+=============================================================================================+\n",
            "| samtoolsbeec0c7b3e | samtools | 1.16.1    | Queued   | workbench.jobs('samtoolsbeec0c7b3e').logs() | samtools sort -o /output/mapping/SRR10261595_sorted.bam -@2 /output/mapping/SRR10261595.bam |\n",
            "+--------------------+----------+-----------+----------+---------------------------------------------+---------------------------------------------------------------------------------------------+\n",
            "+--------------------+----------+-----------+----------+---------------------------------------------+---------------------------------------------------------------------------------------------+\n",
            "| Job ID             | Tool     | Version   | Status   | Get Logs                                    | Full Command                                                                                |\n",
            "+====================+==========+===========+==========+=============================================+=============================================================================================+\n",
            "| samtoolsf5e0251cab | samtools | 1.16.1    | Queued   | workbench.jobs('samtoolsf5e0251cab').logs() | samtools sort -o /output/mapping/SRR10261596_sorted.bam -@2 /output/mapping/SRR10261596.bam |\n",
            "+--------------------+----------+-----------+----------+---------------------------------------------+---------------------------------------------------------------------------------------------+\n"
          ]
        }
      ]
    },
    {
      "cell_type": "code",
      "source": [
        "for sample in samples:\n",
        "  workbench.run(tool='samtools', full_command=f'samtools index /output/mapping/{sample}_sorted.bam')\n",
        "#workbench.run(tool='samtools', full_command='samtools index /output/mapping/SRR10261592_sorted.bam')"
      ],
      "metadata": {
        "id": "Je2QI6u8PCpq",
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "outputId": "4b4e7cd8-62e0-4b84-faa9-adad765047d0"
      },
      "execution_count": null,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "+--------------------+----------+-----------+----------+---------------------------------------------+-------------------------------------------------------+\n",
            "| Job ID             | Tool     | Version   | Status   | Get Logs                                    | Full Command                                          |\n",
            "+====================+==========+===========+==========+=============================================+=======================================================+\n",
            "| samtools163b6cfcab | samtools | 1.16.1    | Queued   | workbench.jobs('samtools163b6cfcab').logs() | samtools index /output/mapping/SRR10261591_sorted.bam |\n",
            "+--------------------+----------+-----------+----------+---------------------------------------------+-------------------------------------------------------+\n",
            "+--------------------+----------+-----------+----------+---------------------------------------------+-------------------------------------------------------+\n",
            "| Job ID             | Tool     | Version   | Status   | Get Logs                                    | Full Command                                          |\n",
            "+====================+==========+===========+==========+=============================================+=======================================================+\n",
            "| samtoolsa78b942962 | samtools | 1.16.1    | Queued   | workbench.jobs('samtoolsa78b942962').logs() | samtools index /output/mapping/SRR10261592_sorted.bam |\n",
            "+--------------------+----------+-----------+----------+---------------------------------------------+-------------------------------------------------------+\n",
            "+--------------------+----------+-----------+----------+---------------------------------------------+-------------------------------------------------------+\n",
            "| Job ID             | Tool     | Version   | Status   | Get Logs                                    | Full Command                                          |\n",
            "+====================+==========+===========+==========+=============================================+=======================================================+\n",
            "| samtools2679dcc84f | samtools | 1.16.1    | Queued   | workbench.jobs('samtools2679dcc84f').logs() | samtools index /output/mapping/SRR10261593_sorted.bam |\n",
            "+--------------------+----------+-----------+----------+---------------------------------------------+-------------------------------------------------------+\n",
            "+--------------------+----------+-----------+----------+---------------------------------------------+-------------------------------------------------------+\n",
            "| Job ID             | Tool     | Version   | Status   | Get Logs                                    | Full Command                                          |\n",
            "+====================+==========+===========+==========+=============================================+=======================================================+\n",
            "| samtools5df3bdf53d | samtools | 1.16.1    | Queued   | workbench.jobs('samtools5df3bdf53d').logs() | samtools index /output/mapping/SRR10261594_sorted.bam |\n",
            "+--------------------+----------+-----------+----------+---------------------------------------------+-------------------------------------------------------+\n",
            "+--------------------+----------+-----------+----------+---------------------------------------------+-------------------------------------------------------+\n",
            "| Job ID             | Tool     | Version   | Status   | Get Logs                                    | Full Command                                          |\n",
            "+====================+==========+===========+==========+=============================================+=======================================================+\n",
            "| samtoolsd605db43eb | samtools | 1.16.1    | Queued   | workbench.jobs('samtoolsd605db43eb').logs() | samtools index /output/mapping/SRR10261595_sorted.bam |\n",
            "+--------------------+----------+-----------+----------+---------------------------------------------+-------------------------------------------------------+\n",
            "+--------------------+----------+-----------+----------+---------------------------------------------+-------------------------------------------------------+\n",
            "| Job ID             | Tool     | Version   | Status   | Get Logs                                    | Full Command                                          |\n",
            "+====================+==========+===========+==========+=============================================+=======================================================+\n",
            "| samtoolse42b8485d8 | samtools | 1.16.1    | Queued   | workbench.jobs('samtoolse42b8485d8').logs() | samtools index /output/mapping/SRR10261596_sorted.bam |\n",
            "+--------------------+----------+-----------+----------+---------------------------------------------+-------------------------------------------------------+\n"
          ]
        }
      ]
    },
    {
      "cell_type": "code",
      "source": [
        "for sample in samples:\n",
        "  workbench.run(tool='samtools', full_command=f'samtools flagstat /output/mapping/{sample}_sorted.bam -O /output/mapping/{sample}_sorted.bam_map_stat.txt')\n"
      ],
      "metadata": {
        "id": "n98sxPsPPDy5",
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "outputId": "c932cd28-13ab-4517-fb28-2cf0131dc72b"
      },
      "execution_count": null,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "+--------------------+----------+-----------+----------+---------------------------------------------+-----------------------------------------------------------------------------------------------------------------+\n",
            "| Job ID             | Tool     | Version   | Status   | Get Logs                                    | Full Command                                                                                                    |\n",
            "+====================+==========+===========+==========+=============================================+=================================================================================================================+\n",
            "| samtools8b5d512a4b | samtools | 1.16.1    | Queued   | workbench.jobs('samtools8b5d512a4b').logs() | samtools flagstat /output/mapping/SRR10261591_sorted.bam -O /output/mapping/SRR10261591_sorted.bam_map_stat.txt |\n",
            "+--------------------+----------+-----------+----------+---------------------------------------------+-----------------------------------------------------------------------------------------------------------------+\n",
            "+--------------------+----------+-----------+----------+---------------------------------------------+-----------------------------------------------------------------------------------------------------------------+\n",
            "| Job ID             | Tool     | Version   | Status   | Get Logs                                    | Full Command                                                                                                    |\n",
            "+====================+==========+===========+==========+=============================================+=================================================================================================================+\n",
            "| samtools5f9f73ca88 | samtools | 1.16.1    | Queued   | workbench.jobs('samtools5f9f73ca88').logs() | samtools flagstat /output/mapping/SRR10261592_sorted.bam -O /output/mapping/SRR10261592_sorted.bam_map_stat.txt |\n",
            "+--------------------+----------+-----------+----------+---------------------------------------------+-----------------------------------------------------------------------------------------------------------------+\n",
            "+--------------------+----------+-----------+----------+---------------------------------------------+-----------------------------------------------------------------------------------------------------------------+\n",
            "| Job ID             | Tool     | Version   | Status   | Get Logs                                    | Full Command                                                                                                    |\n",
            "+====================+==========+===========+==========+=============================================+=================================================================================================================+\n",
            "| samtools66b6457818 | samtools | 1.16.1    | Queued   | workbench.jobs('samtools66b6457818').logs() | samtools flagstat /output/mapping/SRR10261593_sorted.bam -O /output/mapping/SRR10261593_sorted.bam_map_stat.txt |\n",
            "+--------------------+----------+-----------+----------+---------------------------------------------+-----------------------------------------------------------------------------------------------------------------+\n",
            "+--------------------+----------+-----------+----------+---------------------------------------------+-----------------------------------------------------------------------------------------------------------------+\n",
            "| Job ID             | Tool     | Version   | Status   | Get Logs                                    | Full Command                                                                                                    |\n",
            "+====================+==========+===========+==========+=============================================+=================================================================================================================+\n",
            "| samtools5bd1a4179a | samtools | 1.16.1    | Queued   | workbench.jobs('samtools5bd1a4179a').logs() | samtools flagstat /output/mapping/SRR10261594_sorted.bam -O /output/mapping/SRR10261594_sorted.bam_map_stat.txt |\n",
            "+--------------------+----------+-----------+----------+---------------------------------------------+-----------------------------------------------------------------------------------------------------------------+\n",
            "+--------------------+----------+-----------+----------+---------------------------------------------+-----------------------------------------------------------------------------------------------------------------+\n",
            "| Job ID             | Tool     | Version   | Status   | Get Logs                                    | Full Command                                                                                                    |\n",
            "+====================+==========+===========+==========+=============================================+=================================================================================================================+\n",
            "| samtools406425553c | samtools | 1.16.1    | Queued   | workbench.jobs('samtools406425553c').logs() | samtools flagstat /output/mapping/SRR10261595_sorted.bam -O /output/mapping/SRR10261595_sorted.bam_map_stat.txt |\n",
            "+--------------------+----------+-----------+----------+---------------------------------------------+-----------------------------------------------------------------------------------------------------------------+\n",
            "+--------------------+----------+-----------+----------+---------------------------------------------+-----------------------------------------------------------------------------------------------------------------+\n",
            "| Job ID             | Tool     | Version   | Status   | Get Logs                                    | Full Command                                                                                                    |\n",
            "+====================+==========+===========+==========+=============================================+=================================================================================================================+\n",
            "| samtoolsf8e4a8bf9c | samtools | 1.16.1    | Queued   | workbench.jobs('samtoolsf8e4a8bf9c').logs() | samtools flagstat /output/mapping/SRR10261596_sorted.bam -O /output/mapping/SRR10261596_sorted.bam_map_stat.txt |\n",
            "+--------------------+----------+-----------+----------+---------------------------------------------+-----------------------------------------------------------------------------------------------------------------+\n"
          ]
        }
      ]
    },
    {
      "cell_type": "code",
      "source": [
        "workbench.jobs()"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "f0kyeyrItCW4",
        "outputId": "baeec961-acae-4f8f-c9d1-459c5bcc39f4"
      },
      "execution_count": null,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "+--------------------+----------+-----------+-----------+---------------------------------------------+-----------------------------------------------------------------------------------------------------------------+\n",
            "| Job ID             | Tool     | Version   | Status    | Get Logs                                    | Full Command                                                                                                    |\n",
            "+====================+==========+===========+===========+=============================================+=================================================================================================================+\n",
            "| samtools163b6cfcab | samtools | 1.16.1    | Succeeded | workbench.jobs('samtools163b6cfcab').logs() | samtools index /output/mapping/SRR10261591_sorted.bam                                                           |\n",
            "+--------------------+----------+-----------+-----------+---------------------------------------------+-----------------------------------------------------------------------------------------------------------------+\n",
            "| samtoolsa78b942962 | samtools | 1.16.1    | Succeeded | workbench.jobs('samtoolsa78b942962').logs() | samtools index /output/mapping/SRR10261592_sorted.bam                                                           |\n",
            "+--------------------+----------+-----------+-----------+---------------------------------------------+-----------------------------------------------------------------------------------------------------------------+\n",
            "| samtools2679dcc84f | samtools | 1.16.1    | Succeeded | workbench.jobs('samtools2679dcc84f').logs() | samtools index /output/mapping/SRR10261593_sorted.bam                                                           |\n",
            "+--------------------+----------+-----------+-----------+---------------------------------------------+-----------------------------------------------------------------------------------------------------------------+\n",
            "| samtools5df3bdf53d | samtools | 1.16.1    | Succeeded | workbench.jobs('samtools5df3bdf53d').logs() | samtools index /output/mapping/SRR10261594_sorted.bam                                                           |\n",
            "+--------------------+----------+-----------+-----------+---------------------------------------------+-----------------------------------------------------------------------------------------------------------------+\n",
            "| samtoolsd605db43eb | samtools | 1.16.1    | Succeeded | workbench.jobs('samtoolsd605db43eb').logs() | samtools index /output/mapping/SRR10261595_sorted.bam                                                           |\n",
            "+--------------------+----------+-----------+-----------+---------------------------------------------+-----------------------------------------------------------------------------------------------------------------+\n",
            "| samtoolse42b8485d8 | samtools | 1.16.1    | Succeeded | workbench.jobs('samtoolse42b8485d8').logs() | samtools index /output/mapping/SRR10261596_sorted.bam                                                           |\n",
            "+--------------------+----------+-----------+-----------+---------------------------------------------+-----------------------------------------------------------------------------------------------------------------+\n",
            "| samtools05740c8d5b | samtools | 1.16.1    | Failed    | workbench.jobs('samtools05740c8d5b').logs() | samtools flagstat /output/mapping/SRR10261591_sorted.bam -o /output/mapping/SRR10261591_sorted.bam_map_stat.txt |\n",
            "+--------------------+----------+-----------+-----------+---------------------------------------------+-----------------------------------------------------------------------------------------------------------------+\n",
            "| samtoolsb7cc382c7e | samtools | 1.16.1    | Failed    | workbench.jobs('samtoolsb7cc382c7e').logs() | samtools flagstat /output/mapping/SRR10261592_sorted.bam -o /output/mapping/SRR10261592_sorted.bam_map_stat.txt |\n",
            "+--------------------+----------+-----------+-----------+---------------------------------------------+-----------------------------------------------------------------------------------------------------------------+\n",
            "| samtools1175195b79 | samtools | 1.16.1    | Failed    | workbench.jobs('samtools1175195b79').logs() | samtools flagstat /output/mapping/SRR10261593_sorted.bam -o /output/mapping/SRR10261593_sorted.bam_map_stat.txt |\n",
            "+--------------------+----------+-----------+-----------+---------------------------------------------+-----------------------------------------------------------------------------------------------------------------+\n",
            "| samtools14d7a52b83 | samtools | 1.16.1    | Failed    | workbench.jobs('samtools14d7a52b83').logs() | samtools flagstat /output/mapping/SRR10261594_sorted.bam -o /output/mapping/SRR10261594_sorted.bam_map_stat.txt |\n",
            "+--------------------+----------+-----------+-----------+---------------------------------------------+-----------------------------------------------------------------------------------------------------------------+\n",
            "| samtools2878c4d2f4 | samtools | 1.16.1    | Failed    | workbench.jobs('samtools2878c4d2f4').logs() | samtools flagstat /output/mapping/SRR10261595_sorted.bam -o /output/mapping/SRR10261595_sorted.bam_map_stat.txt |\n",
            "+--------------------+----------+-----------+-----------+---------------------------------------------+-----------------------------------------------------------------------------------------------------------------+\n",
            "| samtools6756a0c0ca | samtools | 1.16.1    | Failed    | workbench.jobs('samtools6756a0c0ca').logs() | samtools flagstat /output/mapping/SRR10261596_sorted.bam -o /output/mapping/SRR10261596_sorted.bam_map_stat.txt |\n",
            "+--------------------+----------+-----------+-----------+---------------------------------------------+-----------------------------------------------------------------------------------------------------------------+\n",
            "| samtools8b5d512a4b | samtools | 1.16.1    | Succeeded | workbench.jobs('samtools8b5d512a4b').logs() | samtools flagstat /output/mapping/SRR10261591_sorted.bam -O /output/mapping/SRR10261591_sorted.bam_map_stat.txt |\n",
            "+--------------------+----------+-----------+-----------+---------------------------------------------+-----------------------------------------------------------------------------------------------------------------+\n",
            "| samtools5f9f73ca88 | samtools | 1.16.1    | Succeeded | workbench.jobs('samtools5f9f73ca88').logs() | samtools flagstat /output/mapping/SRR10261592_sorted.bam -O /output/mapping/SRR10261592_sorted.bam_map_stat.txt |\n",
            "+--------------------+----------+-----------+-----------+---------------------------------------------+-----------------------------------------------------------------------------------------------------------------+\n",
            "| samtools66b6457818 | samtools | 1.16.1    | Succeeded | workbench.jobs('samtools66b6457818').logs() | samtools flagstat /output/mapping/SRR10261593_sorted.bam -O /output/mapping/SRR10261593_sorted.bam_map_stat.txt |\n",
            "+--------------------+----------+-----------+-----------+---------------------------------------------+-----------------------------------------------------------------------------------------------------------------+\n",
            "| samtools5bd1a4179a | samtools | 1.16.1    | Succeeded | workbench.jobs('samtools5bd1a4179a').logs() | samtools flagstat /output/mapping/SRR10261594_sorted.bam -O /output/mapping/SRR10261594_sorted.bam_map_stat.txt |\n",
            "+--------------------+----------+-----------+-----------+---------------------------------------------+-----------------------------------------------------------------------------------------------------------------+\n",
            "| samtools406425553c | samtools | 1.16.1    | Succeeded | workbench.jobs('samtools406425553c').logs() | samtools flagstat /output/mapping/SRR10261595_sorted.bam -O /output/mapping/SRR10261595_sorted.bam_map_stat.txt |\n",
            "+--------------------+----------+-----------+-----------+---------------------------------------------+-----------------------------------------------------------------------------------------------------------------+\n",
            "| samtoolsf8e4a8bf9c | samtools | 1.16.1    | Succeeded | workbench.jobs('samtoolsf8e4a8bf9c').logs() | samtools flagstat /output/mapping/SRR10261596_sorted.bam -O /output/mapping/SRR10261596_sorted.bam_map_stat.txt |\n",
            "+--------------------+----------+-----------+-----------+---------------------------------------------+-----------------------------------------------------------------------------------------------------------------+\n"
          ]
        }
      ]
    },
    {
      "cell_type": "markdown",
      "source": [
        "The commands will perform read mapping with bwa mem, which will generate the output file in a so-called sam/bam format (learn more about these formats). Then these output files will be sorted (i.e. samtools sort), indexed (i.e. samtools index), and used for calculation of mapping statistics (i.e. samtools flagstat). These three steps will be performed using samtools v. 1.6 software (Li et al. 2009). This sorting and indexing process will allow us to perform various helpful downstream analyses, such as visualization and mapping statistics calculation.\n",
        "\n",
        "\n",
        "\n",
        "When finished, inspect the *map_stats.txt files to view the mapping statistics. You will see that over 80% of reads are mapped and properly paired for all samples, indicating a successful mapping.\n",
        "\n",
        "Before we can proceed in peak-calling our ATAC-seq data, three auxiliary analyses need to be performed:\n",
        "\n",
        "1.   First, because the hybrid samples contain data from both Saccharomyces species, we need to select the S. cerevisiae data in those samples.\n",
        "2.   Second, in ATAC-seq analysis it is generally recommended to remove mitochondrial data before the peak calling analysis. This is because the mitochondrial genome (mtDNA) does not have chromatin and so is “naked,” meaning a high proportion of ATAC-seq data can originate from mtDNA. This can lead to false peaks. In other words, we will select only non-mtDNA S. cerevisiae data from all samples for further analysis.\n",
        "3.   Third, as described above in the PCR Amplification section, ATAC-seq data can contain a high proportion of PCR-duplicated reads. These are technical artifacts from experimental procedures and need to be removed/masked for the downstream analysis.\n",
        "\n",
        "\n",
        "To perform these three analyses, we’ll the following commands:\n"
      ],
      "metadata": {
        "id": "irFYG_-AQOG6"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "for sample in samples:\n",
        "    workbench.run(tool=\"samtools\", full_command=f\"samtools view -bh /output/mapping/{sample}_sorted.bam I II III IV V VI VII VIII IX X XI XII XIII XIV XV XVI > /output/mapping/{sample}_sorted_SC_subset.bam\")\n",
        "    "
      ],
      "metadata": {
        "id": "Wx75FRQpRegU",
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "outputId": "24b2c14b-82c1-49ed-e56b-0a85907966c1"
      },
      "execution_count": null,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "+--------------------+----------+-----------+----------+---------------------------------------------+-------------------------------------------------------------------------------------------------------------------------------------------------------------------+\n",
            "| Job ID             | Tool     | Version   | Status   | Get Logs                                    | Full Command                                                                                                                                                      |\n",
            "+====================+==========+===========+==========+=============================================+===================================================================================================================================================================+\n",
            "| samtools3272c58d91 | samtools | 1.16.1    | Queued   | workbench.jobs('samtools3272c58d91').logs() | samtools view -bh /output/mapping/SRR10261591_sorted.bam I II III IV V VI VII VIII IX X XI XII XIII XIV XV XVI > /output/mapping/SRR10261591_sorted_SC_subset.bam |\n",
            "+--------------------+----------+-----------+----------+---------------------------------------------+-------------------------------------------------------------------------------------------------------------------------------------------------------------------+\n",
            "+--------------------+----------+-----------+----------+---------------------------------------------+-------------------------------------------------------------------------------------------------------------------------------------------------------------------+\n",
            "| Job ID             | Tool     | Version   | Status   | Get Logs                                    | Full Command                                                                                                                                                      |\n",
            "+====================+==========+===========+==========+=============================================+===================================================================================================================================================================+\n",
            "| samtoolsed6745cf4f | samtools | 1.16.1    | Queued   | workbench.jobs('samtoolsed6745cf4f').logs() | samtools view -bh /output/mapping/SRR10261592_sorted.bam I II III IV V VI VII VIII IX X XI XII XIII XIV XV XVI > /output/mapping/SRR10261592_sorted_SC_subset.bam |\n",
            "+--------------------+----------+-----------+----------+---------------------------------------------+-------------------------------------------------------------------------------------------------------------------------------------------------------------------+\n",
            "+--------------------+----------+-----------+----------+---------------------------------------------+-------------------------------------------------------------------------------------------------------------------------------------------------------------------+\n",
            "| Job ID             | Tool     | Version   | Status   | Get Logs                                    | Full Command                                                                                                                                                      |\n",
            "+====================+==========+===========+==========+=============================================+===================================================================================================================================================================+\n",
            "| samtools363c580dda | samtools | 1.16.1    | Queued   | workbench.jobs('samtools363c580dda').logs() | samtools view -bh /output/mapping/SRR10261593_sorted.bam I II III IV V VI VII VIII IX X XI XII XIII XIV XV XVI > /output/mapping/SRR10261593_sorted_SC_subset.bam |\n",
            "+--------------------+----------+-----------+----------+---------------------------------------------+-------------------------------------------------------------------------------------------------------------------------------------------------------------------+\n",
            "+--------------------+----------+-----------+----------+---------------------------------------------+-------------------------------------------------------------------------------------------------------------------------------------------------------------------+\n",
            "| Job ID             | Tool     | Version   | Status   | Get Logs                                    | Full Command                                                                                                                                                      |\n",
            "+====================+==========+===========+==========+=============================================+===================================================================================================================================================================+\n",
            "| samtools19db8479a4 | samtools | 1.16.1    | Queued   | workbench.jobs('samtools19db8479a4').logs() | samtools view -bh /output/mapping/SRR10261594_sorted.bam I II III IV V VI VII VIII IX X XI XII XIII XIV XV XVI > /output/mapping/SRR10261594_sorted_SC_subset.bam |\n",
            "+--------------------+----------+-----------+----------+---------------------------------------------+-------------------------------------------------------------------------------------------------------------------------------------------------------------------+\n",
            "+--------------------+----------+-----------+----------+---------------------------------------------+-------------------------------------------------------------------------------------------------------------------------------------------------------------------+\n",
            "| Job ID             | Tool     | Version   | Status   | Get Logs                                    | Full Command                                                                                                                                                      |\n",
            "+====================+==========+===========+==========+=============================================+===================================================================================================================================================================+\n",
            "| samtoolsead541a4a0 | samtools | 1.16.1    | Queued   | workbench.jobs('samtoolsead541a4a0').logs() | samtools view -bh /output/mapping/SRR10261595_sorted.bam I II III IV V VI VII VIII IX X XI XII XIII XIV XV XVI > /output/mapping/SRR10261595_sorted_SC_subset.bam |\n",
            "+--------------------+----------+-----------+----------+---------------------------------------------+-------------------------------------------------------------------------------------------------------------------------------------------------------------------+\n",
            "+--------------------+----------+-----------+----------+---------------------------------------------+-------------------------------------------------------------------------------------------------------------------------------------------------------------------+\n",
            "| Job ID             | Tool     | Version   | Status   | Get Logs                                    | Full Command                                                                                                                                                      |\n",
            "+====================+==========+===========+==========+=============================================+===================================================================================================================================================================+\n",
            "| samtools5f0de3187c | samtools | 1.16.1    | Queued   | workbench.jobs('samtools5f0de3187c').logs() | samtools view -bh /output/mapping/SRR10261596_sorted.bam I II III IV V VI VII VIII IX X XI XII XIII XIV XV XVI > /output/mapping/SRR10261596_sorted_SC_subset.bam |\n",
            "+--------------------+----------+-----------+----------+---------------------------------------------+-------------------------------------------------------------------------------------------------------------------------------------------------------------------+\n"
          ]
        }
      ]
    },
    {
      "cell_type": "code",
      "source": [
        "for sample in samples:    \n",
        "    workbench.run(tool='picard', full_command=f'picard MarkDuplicates I=/output/mapping/{sample}_sorted_SC_subset.bam O=/output/mapping/{sample}_sorted_SC_subset_dedup.bam M=/output/mapping/{sample}_markdup_metrics.txt')\n"
      ],
      "metadata": {
        "id": "gbQ2qShnCDeA",
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "outputId": "79558d50-403a-4142-d2e8-5d12ea3a675b"
      },
      "execution_count": null,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "+------------------+--------+-----------+----------+-------------------------------------------+-------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------+\n",
            "| Job ID           | Tool   | Version   | Status   | Get Logs                                  | Full Command                                                                                                                                                                        |\n",
            "+==================+========+===========+==========+===========================================+=====================================================================================================================================================================================+\n",
            "| picardf263edacba | picard | 2.27.3    | Queued   | workbench.jobs('picardf263edacba').logs() | picard MarkDuplicates I=/output/mapping/SRR10261591_sorted_SC_subset.bam O=/output/mapping/SRR10261591_sorted_SC_subset_dedup.bam M=/output/mapping/SRR10261591_markdup_metrics.txt |\n",
            "+------------------+--------+-----------+----------+-------------------------------------------+-------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------+\n",
            "+------------------+--------+-----------+----------+-------------------------------------------+-------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------+\n",
            "| Job ID           | Tool   | Version   | Status   | Get Logs                                  | Full Command                                                                                                                                                                        |\n",
            "+==================+========+===========+==========+===========================================+=====================================================================================================================================================================================+\n",
            "| picard07e2dcf5b3 | picard | 2.27.3    | Queued   | workbench.jobs('picard07e2dcf5b3').logs() | picard MarkDuplicates I=/output/mapping/SRR10261592_sorted_SC_subset.bam O=/output/mapping/SRR10261592_sorted_SC_subset_dedup.bam M=/output/mapping/SRR10261592_markdup_metrics.txt |\n",
            "+------------------+--------+-----------+----------+-------------------------------------------+-------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------+\n",
            "+------------------+--------+-----------+----------+-------------------------------------------+-------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------+\n",
            "| Job ID           | Tool   | Version   | Status   | Get Logs                                  | Full Command                                                                                                                                                                        |\n",
            "+==================+========+===========+==========+===========================================+=====================================================================================================================================================================================+\n",
            "| picard6b72f9b849 | picard | 2.27.3    | Queued   | workbench.jobs('picard6b72f9b849').logs() | picard MarkDuplicates I=/output/mapping/SRR10261593_sorted_SC_subset.bam O=/output/mapping/SRR10261593_sorted_SC_subset_dedup.bam M=/output/mapping/SRR10261593_markdup_metrics.txt |\n",
            "+------------------+--------+-----------+----------+-------------------------------------------+-------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------+\n",
            "+------------------+--------+-----------+----------+-------------------------------------------+-------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------+\n",
            "| Job ID           | Tool   | Version   | Status   | Get Logs                                  | Full Command                                                                                                                                                                        |\n",
            "+==================+========+===========+==========+===========================================+=====================================================================================================================================================================================+\n",
            "| picard1c85b0c584 | picard | 2.27.3    | Queued   | workbench.jobs('picard1c85b0c584').logs() | picard MarkDuplicates I=/output/mapping/SRR10261594_sorted_SC_subset.bam O=/output/mapping/SRR10261594_sorted_SC_subset_dedup.bam M=/output/mapping/SRR10261594_markdup_metrics.txt |\n",
            "+------------------+--------+-----------+----------+-------------------------------------------+-------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------+\n",
            "+------------------+--------+-----------+----------+-------------------------------------------+-------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------+\n",
            "| Job ID           | Tool   | Version   | Status   | Get Logs                                  | Full Command                                                                                                                                                                        |\n",
            "+==================+========+===========+==========+===========================================+=====================================================================================================================================================================================+\n",
            "| picard5257f16275 | picard | 2.27.3    | Queued   | workbench.jobs('picard5257f16275').logs() | picard MarkDuplicates I=/output/mapping/SRR10261595_sorted_SC_subset.bam O=/output/mapping/SRR10261595_sorted_SC_subset_dedup.bam M=/output/mapping/SRR10261595_markdup_metrics.txt |\n",
            "+------------------+--------+-----------+----------+-------------------------------------------+-------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------+\n",
            "+------------------+--------+-----------+----------+-------------------------------------------+-------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------+\n",
            "| Job ID           | Tool   | Version   | Status   | Get Logs                                  | Full Command                                                                                                                                                                        |\n",
            "+==================+========+===========+==========+===========================================+=====================================================================================================================================================================================+\n",
            "| picardd95fd25d39 | picard | 2.27.3    | Queued   | workbench.jobs('picardd95fd25d39').logs() | picard MarkDuplicates I=/output/mapping/SRR10261596_sorted_SC_subset.bam O=/output/mapping/SRR10261596_sorted_SC_subset_dedup.bam M=/output/mapping/SRR10261596_markdup_metrics.txt |\n",
            "+------------------+--------+-----------+----------+-------------------------------------------+-------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------+\n"
          ]
        }
      ]
    },
    {
      "cell_type": "code",
      "source": [
        "for sample in samples:    \n",
        "    workbench.run(tool='samtools', full_command=f'samtools index /output/mapping/{sample}_sorted_SC_subset_dedup.bam')\n"
      ],
      "metadata": {
        "id": "IhlATqRgCEqm",
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "outputId": "1368e945-b9d3-4107-d53a-da2b17a13e5d"
      },
      "execution_count": null,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "+--------------------+----------+-----------+----------+---------------------------------------------+-----------------------------------------------------------------------+\n",
            "| Job ID             | Tool     | Version   | Status   | Get Logs                                    | Full Command                                                          |\n",
            "+====================+==========+===========+==========+=============================================+=======================================================================+\n",
            "| samtools1b483dff23 | samtools | 1.16.1    | Queued   | workbench.jobs('samtools1b483dff23').logs() | samtools index /output/mapping/SRR10261591_sorted_SC_subset_dedup.bam |\n",
            "+--------------------+----------+-----------+----------+---------------------------------------------+-----------------------------------------------------------------------+\n",
            "+--------------------+----------+-----------+----------+---------------------------------------------+-----------------------------------------------------------------------+\n",
            "| Job ID             | Tool     | Version   | Status   | Get Logs                                    | Full Command                                                          |\n",
            "+====================+==========+===========+==========+=============================================+=======================================================================+\n",
            "| samtools25d2a9f546 | samtools | 1.16.1    | Queued   | workbench.jobs('samtools25d2a9f546').logs() | samtools index /output/mapping/SRR10261592_sorted_SC_subset_dedup.bam |\n",
            "+--------------------+----------+-----------+----------+---------------------------------------------+-----------------------------------------------------------------------+\n",
            "+--------------------+----------+-----------+----------+---------------------------------------------+-----------------------------------------------------------------------+\n",
            "| Job ID             | Tool     | Version   | Status   | Get Logs                                    | Full Command                                                          |\n",
            "+====================+==========+===========+==========+=============================================+=======================================================================+\n",
            "| samtools3612dccf5d | samtools | 1.16.1    | Queued   | workbench.jobs('samtools3612dccf5d').logs() | samtools index /output/mapping/SRR10261593_sorted_SC_subset_dedup.bam |\n",
            "+--------------------+----------+-----------+----------+---------------------------------------------+-----------------------------------------------------------------------+\n",
            "+--------------------+----------+-----------+----------+---------------------------------------------+-----------------------------------------------------------------------+\n",
            "| Job ID             | Tool     | Version   | Status   | Get Logs                                    | Full Command                                                          |\n",
            "+====================+==========+===========+==========+=============================================+=======================================================================+\n",
            "| samtools753709a242 | samtools | 1.16.1    | Queued   | workbench.jobs('samtools753709a242').logs() | samtools index /output/mapping/SRR10261594_sorted_SC_subset_dedup.bam |\n",
            "+--------------------+----------+-----------+----------+---------------------------------------------+-----------------------------------------------------------------------+\n",
            "+--------------------+----------+-----------+----------+---------------------------------------------+-----------------------------------------------------------------------+\n",
            "| Job ID             | Tool     | Version   | Status   | Get Logs                                    | Full Command                                                          |\n",
            "+====================+==========+===========+==========+=============================================+=======================================================================+\n",
            "| samtoolsa55646d3b2 | samtools | 1.16.1    | Queued   | workbench.jobs('samtoolsa55646d3b2').logs() | samtools index /output/mapping/SRR10261595_sorted_SC_subset_dedup.bam |\n",
            "+--------------------+----------+-----------+----------+---------------------------------------------+-----------------------------------------------------------------------+\n",
            "+--------------------+----------+-----------+----------+---------------------------------------------+-----------------------------------------------------------------------+\n",
            "| Job ID             | Tool     | Version   | Status   | Get Logs                                    | Full Command                                                          |\n",
            "+====================+==========+===========+==========+=============================================+=======================================================================+\n",
            "| samtools900bf1b244 | samtools | 1.16.1    | Queued   | workbench.jobs('samtools900bf1b244').logs() | samtools index /output/mapping/SRR10261596_sorted_SC_subset_dedup.bam |\n",
            "+--------------------+----------+-----------+----------+---------------------------------------------+-----------------------------------------------------------------------+\n"
          ]
        }
      ]
    },
    {
      "cell_type": "code",
      "source": [
        "for sample in samples:    \n",
        "    workbench.run(tool='samtools', full_command=f'samtools flagstat /output/mapping/{sample}_sorted_SC_subset_dedup.bam > /output/mapping/{sample}_sorted_SC_subset_dedup_map_stats.txt')\n"
      ],
      "metadata": {
        "id": "56JFsFURCGAY",
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "outputId": "5ad36bb4-40aa-4eaa-cc24-853cda8379e3"
      },
      "execution_count": null,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "+--------------------+----------+-----------+----------+---------------------------------------------+---------------------------------------------------------------------------------------------------------------------------------------------+\n",
            "| Job ID             | Tool     | Version   | Status   | Get Logs                                    | Full Command                                                                                                                                |\n",
            "+====================+==========+===========+==========+=============================================+=============================================================================================================================================+\n",
            "| samtools1483d8f4a0 | samtools | 1.16.1    | Queued   | workbench.jobs('samtools1483d8f4a0').logs() | samtools flagstat /output/mapping/SRR10261591_sorted_SC_subset_dedup.bam > /output/mapping/SRR10261591_sorted_SC_subset_dedup_map_stats.txt |\n",
            "+--------------------+----------+-----------+----------+---------------------------------------------+---------------------------------------------------------------------------------------------------------------------------------------------+\n",
            "+--------------------+----------+-----------+----------+---------------------------------------------+---------------------------------------------------------------------------------------------------------------------------------------------+\n",
            "| Job ID             | Tool     | Version   | Status   | Get Logs                                    | Full Command                                                                                                                                |\n",
            "+====================+==========+===========+==========+=============================================+=============================================================================================================================================+\n",
            "| samtoolsdafa14c676 | samtools | 1.16.1    | Queued   | workbench.jobs('samtoolsdafa14c676').logs() | samtools flagstat /output/mapping/SRR10261592_sorted_SC_subset_dedup.bam > /output/mapping/SRR10261592_sorted_SC_subset_dedup_map_stats.txt |\n",
            "+--------------------+----------+-----------+----------+---------------------------------------------+---------------------------------------------------------------------------------------------------------------------------------------------+\n",
            "+--------------------+----------+-----------+----------+---------------------------------------------+---------------------------------------------------------------------------------------------------------------------------------------------+\n",
            "| Job ID             | Tool     | Version   | Status   | Get Logs                                    | Full Command                                                                                                                                |\n",
            "+====================+==========+===========+==========+=============================================+=============================================================================================================================================+\n",
            "| samtoolsb416573031 | samtools | 1.16.1    | Queued   | workbench.jobs('samtoolsb416573031').logs() | samtools flagstat /output/mapping/SRR10261593_sorted_SC_subset_dedup.bam > /output/mapping/SRR10261593_sorted_SC_subset_dedup_map_stats.txt |\n",
            "+--------------------+----------+-----------+----------+---------------------------------------------+---------------------------------------------------------------------------------------------------------------------------------------------+\n",
            "+--------------------+----------+-----------+----------+---------------------------------------------+---------------------------------------------------------------------------------------------------------------------------------------------+\n",
            "| Job ID             | Tool     | Version   | Status   | Get Logs                                    | Full Command                                                                                                                                |\n",
            "+====================+==========+===========+==========+=============================================+=============================================================================================================================================+\n",
            "| samtools5d45e60585 | samtools | 1.16.1    | Queued   | workbench.jobs('samtools5d45e60585').logs() | samtools flagstat /output/mapping/SRR10261594_sorted_SC_subset_dedup.bam > /output/mapping/SRR10261594_sorted_SC_subset_dedup_map_stats.txt |\n",
            "+--------------------+----------+-----------+----------+---------------------------------------------+---------------------------------------------------------------------------------------------------------------------------------------------+\n",
            "+--------------------+----------+-----------+----------+---------------------------------------------+---------------------------------------------------------------------------------------------------------------------------------------------+\n",
            "| Job ID             | Tool     | Version   | Status   | Get Logs                                    | Full Command                                                                                                                                |\n",
            "+====================+==========+===========+==========+=============================================+=============================================================================================================================================+\n",
            "| samtoolsce897b1ffd | samtools | 1.16.1    | Queued   | workbench.jobs('samtoolsce897b1ffd').logs() | samtools flagstat /output/mapping/SRR10261595_sorted_SC_subset_dedup.bam > /output/mapping/SRR10261595_sorted_SC_subset_dedup_map_stats.txt |\n",
            "+--------------------+----------+-----------+----------+---------------------------------------------+---------------------------------------------------------------------------------------------------------------------------------------------+\n",
            "+--------------------+----------+-----------+----------+---------------------------------------------+---------------------------------------------------------------------------------------------------------------------------------------------+\n",
            "| Job ID             | Tool     | Version   | Status   | Get Logs                                    | Full Command                                                                                                                                |\n",
            "+====================+==========+===========+==========+=============================================+=============================================================================================================================================+\n",
            "| samtools87ce6ef2b2 | samtools | 1.16.1    | Queued   | workbench.jobs('samtools87ce6ef2b2').logs() | samtools flagstat /output/mapping/SRR10261596_sorted_SC_subset_dedup.bam > /output/mapping/SRR10261596_sorted_SC_subset_dedup_map_stats.txt |\n",
            "+--------------------+----------+-----------+----------+---------------------------------------------+---------------------------------------------------------------------------------------------------------------------------------------------+\n"
          ]
        }
      ]
    },
    {
      "cell_type": "markdown",
      "source": [
        "As in the previous step, this commands will iterate over sample names. The term samtools view selects non-mtDNA chromosomes of S. cerevisiae (chromosome names are denoted with Roman numerals). The term picard MarkDuplicates will tag duplicated reads using Picard MarkDuplicated v. 2.27.3 software. And finally, samtools index and samtools flagstat will respectively index and calculate mapping statistics of the output files.\n",
        "\n",
        "\n",
        "When finished, inspect the *_SC_subset_dedup_map_stats.txt files (i.e mapping statistics after subsetting and deduplication). For S. cerevisiae samples, you will notice that after subsetting, the number of mapped reads remains largely unchanged. This indicates slight mitochondrial contamination. For the hybrid data, after subsetting there are still more than 20 million properly paired mapped reads for S. cerevisiae, which is more than enough for further peak calling analysis (see our previous section [Data Analysis: A Quick Overview](https://tinybio.readme.io/docs/analyze-the-data))."
      ],
      "metadata": {
        "id": "SlVryuTGSDpv"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "#Peak calling\n",
        "\n",
        "Now that we’ve obtained high-quality and high-coverage deduplicated bam files for S. cerevisiae and the hybrid, we can finally perform ATAC-seq peak calling. Peak calling identifies genomic regions which are enriched in aligned reads (to a statistically significant extent) compared to background levels. In the context of ATAC-seq, peaks will ideally correspond to the “naked” nucleosome-free DNA regions attacked by the Tn5 transposase. For this step, we will use MACS2 v. 2.2.7.1 software, which was originally designed for ChIP-seq data, but also works for ATAC-seq.\n",
        "\n",
        "The commands below will generate commands to run MACS2 for all samples:"
      ],
      "metadata": {
        "id": "TGKoxGr8SnaY"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "#while read sample; do\n",
        "#   echo \"macs2 callpeak -t ${sample}_SC_subset_dedup.bam -f BAMPE -n ${sample}_peak -g 12000000 --keep-dup all\"\n",
        "#done<../raw_data/sample_ids.txt \n",
        "\n",
        "\n",
        "for sample in samples:\n",
        "  workbench.run(tool='macs2', full_command=f'macs2 callpeak -t /output/mapping/{sample}_sorted_SC_subset_dedup.bam -f BAMPE -n /output/mapping/{sample}_peak -g 12000000 --keep-dup all')\n",
        "\n",
        "  #workbench.run(tool='macs2', full_command=f'macs2 callpeak -t /mnt/gcs/output/mapping/{sample}_sorted_SC_subset.bam -f BAMPE -n /mnt/gcs/output/mapping/{sample}_peak -g 12000000 --keep-dup all')"
      ],
      "metadata": {
        "id": "EkV4mzxRS0CK",
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "outputId": "955fc76e-e6bf-4437-9c7c-2da1f9915c1d"
      },
      "execution_count": null,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "+-----------------+--------+-----------+----------+------------------------------------------+--------------------------------------------------------------------------------------------------------------------------------------------------+\n",
            "| Job ID          | Tool   | Version   | Status   | Get Logs                                 | Full Command                                                                                                                                     |\n",
            "+=================+========+===========+==========+==========================================+==================================================================================================================================================+\n",
            "| macs2e84080522e | macs2  | 2.2.7.1   | Queued   | workbench.jobs('macs2e84080522e').logs() | macs2 callpeak -t /output/mapping/SRR10261591_sorted_SC_subset_dedup.bam -f BAMPE -n /output/mapping/SRR10261591_peak -g 12000000 --keep-dup all |\n",
            "+-----------------+--------+-----------+----------+------------------------------------------+--------------------------------------------------------------------------------------------------------------------------------------------------+\n",
            "+-----------------+--------+-----------+----------+------------------------------------------+--------------------------------------------------------------------------------------------------------------------------------------------------+\n",
            "| Job ID          | Tool   | Version   | Status   | Get Logs                                 | Full Command                                                                                                                                     |\n",
            "+=================+========+===========+==========+==========================================+==================================================================================================================================================+\n",
            "| macs201a56f821b | macs2  | 2.2.7.1   | Queued   | workbench.jobs('macs201a56f821b').logs() | macs2 callpeak -t /output/mapping/SRR10261592_sorted_SC_subset_dedup.bam -f BAMPE -n /output/mapping/SRR10261592_peak -g 12000000 --keep-dup all |\n",
            "+-----------------+--------+-----------+----------+------------------------------------------+--------------------------------------------------------------------------------------------------------------------------------------------------+\n",
            "+-----------------+--------+-----------+----------+------------------------------------------+--------------------------------------------------------------------------------------------------------------------------------------------------+\n",
            "| Job ID          | Tool   | Version   | Status   | Get Logs                                 | Full Command                                                                                                                                     |\n",
            "+=================+========+===========+==========+==========================================+==================================================================================================================================================+\n",
            "| macs2fbdb269e4c | macs2  | 2.2.7.1   | Queued   | workbench.jobs('macs2fbdb269e4c').logs() | macs2 callpeak -t /output/mapping/SRR10261593_sorted_SC_subset_dedup.bam -f BAMPE -n /output/mapping/SRR10261593_peak -g 12000000 --keep-dup all |\n",
            "+-----------------+--------+-----------+----------+------------------------------------------+--------------------------------------------------------------------------------------------------------------------------------------------------+\n",
            "+-----------------+--------+-----------+----------+------------------------------------------+--------------------------------------------------------------------------------------------------------------------------------------------------+\n",
            "| Job ID          | Tool   | Version   | Status   | Get Logs                                 | Full Command                                                                                                                                     |\n",
            "+=================+========+===========+==========+==========================================+==================================================================================================================================================+\n",
            "| macs2d97578d208 | macs2  | 2.2.7.1   | Queued   | workbench.jobs('macs2d97578d208').logs() | macs2 callpeak -t /output/mapping/SRR10261594_sorted_SC_subset_dedup.bam -f BAMPE -n /output/mapping/SRR10261594_peak -g 12000000 --keep-dup all |\n",
            "+-----------------+--------+-----------+----------+------------------------------------------+--------------------------------------------------------------------------------------------------------------------------------------------------+\n",
            "+-----------------+--------+-----------+----------+------------------------------------------+--------------------------------------------------------------------------------------------------------------------------------------------------+\n",
            "| Job ID          | Tool   | Version   | Status   | Get Logs                                 | Full Command                                                                                                                                     |\n",
            "+=================+========+===========+==========+==========================================+==================================================================================================================================================+\n",
            "| macs22cb5b57017 | macs2  | 2.2.7.1   | Queued   | workbench.jobs('macs22cb5b57017').logs() | macs2 callpeak -t /output/mapping/SRR10261595_sorted_SC_subset_dedup.bam -f BAMPE -n /output/mapping/SRR10261595_peak -g 12000000 --keep-dup all |\n",
            "+-----------------+--------+-----------+----------+------------------------------------------+--------------------------------------------------------------------------------------------------------------------------------------------------+\n",
            "+-----------------+--------+-----------+----------+------------------------------------------+--------------------------------------------------------------------------------------------------------------------------------------------------+\n",
            "| Job ID          | Tool   | Version   | Status   | Get Logs                                 | Full Command                                                                                                                                     |\n",
            "+=================+========+===========+==========+==========================================+==================================================================================================================================================+\n",
            "| macs28aaaa7748d | macs2  | 2.2.7.1   | Queued   | workbench.jobs('macs28aaaa7748d').logs() | macs2 callpeak -t /output/mapping/SRR10261596_sorted_SC_subset_dedup.bam -f BAMPE -n /output/mapping/SRR10261596_peak -g 12000000 --keep-dup all |\n",
            "+-----------------+--------+-----------+----------+------------------------------------------+--------------------------------------------------------------------------------------------------------------------------------------------------+\n"
          ]
        }
      ]
    },
    {
      "cell_type": "code",
      "source": [
        "workbench.jobs()"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "rw4nhEqgPBbt",
        "outputId": "07b00926-6e8d-44d2-8050-5823c72a9815"
      },
      "execution_count": null,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "+--------------------+----------+-----------+-----------+---------------------------------------------+-------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------+\n",
            "| Job ID             | Tool     | Version   | Status    | Get Logs                                    | Full Command                                                                                                                                                                        |\n",
            "+====================+==========+===========+===========+=============================================+=====================================================================================================================================================================================+\n",
            "| samtools163b6cfcab | samtools | 1.16.1    | Succeeded | workbench.jobs('samtools163b6cfcab').logs() | samtools index /output/mapping/SRR10261591_sorted.bam                                                                                                                               |\n",
            "+--------------------+----------+-----------+-----------+---------------------------------------------+-------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------+\n",
            "| samtoolsa78b942962 | samtools | 1.16.1    | Succeeded | workbench.jobs('samtoolsa78b942962').logs() | samtools index /output/mapping/SRR10261592_sorted.bam                                                                                                                               |\n",
            "+--------------------+----------+-----------+-----------+---------------------------------------------+-------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------+\n",
            "| samtools2679dcc84f | samtools | 1.16.1    | Succeeded | workbench.jobs('samtools2679dcc84f').logs() | samtools index /output/mapping/SRR10261593_sorted.bam                                                                                                                               |\n",
            "+--------------------+----------+-----------+-----------+---------------------------------------------+-------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------+\n",
            "| samtools5df3bdf53d | samtools | 1.16.1    | Succeeded | workbench.jobs('samtools5df3bdf53d').logs() | samtools index /output/mapping/SRR10261594_sorted.bam                                                                                                                               |\n",
            "+--------------------+----------+-----------+-----------+---------------------------------------------+-------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------+\n",
            "| samtoolsd605db43eb | samtools | 1.16.1    | Succeeded | workbench.jobs('samtoolsd605db43eb').logs() | samtools index /output/mapping/SRR10261595_sorted.bam                                                                                                                               |\n",
            "+--------------------+----------+-----------+-----------+---------------------------------------------+-------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------+\n",
            "| samtoolse42b8485d8 | samtools | 1.16.1    | Succeeded | workbench.jobs('samtoolse42b8485d8').logs() | samtools index /output/mapping/SRR10261596_sorted.bam                                                                                                                               |\n",
            "+--------------------+----------+-----------+-----------+---------------------------------------------+-------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------+\n",
            "| samtools05740c8d5b | samtools | 1.16.1    | Failed    | workbench.jobs('samtools05740c8d5b').logs() | samtools flagstat /output/mapping/SRR10261591_sorted.bam -o /output/mapping/SRR10261591_sorted.bam_map_stat.txt                                                                     |\n",
            "+--------------------+----------+-----------+-----------+---------------------------------------------+-------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------+\n",
            "| samtoolsb7cc382c7e | samtools | 1.16.1    | Failed    | workbench.jobs('samtoolsb7cc382c7e').logs() | samtools flagstat /output/mapping/SRR10261592_sorted.bam -o /output/mapping/SRR10261592_sorted.bam_map_stat.txt                                                                     |\n",
            "+--------------------+----------+-----------+-----------+---------------------------------------------+-------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------+\n",
            "| samtools1175195b79 | samtools | 1.16.1    | Failed    | workbench.jobs('samtools1175195b79').logs() | samtools flagstat /output/mapping/SRR10261593_sorted.bam -o /output/mapping/SRR10261593_sorted.bam_map_stat.txt                                                                     |\n",
            "+--------------------+----------+-----------+-----------+---------------------------------------------+-------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------+\n",
            "| samtools14d7a52b83 | samtools | 1.16.1    | Failed    | workbench.jobs('samtools14d7a52b83').logs() | samtools flagstat /output/mapping/SRR10261594_sorted.bam -o /output/mapping/SRR10261594_sorted.bam_map_stat.txt                                                                     |\n",
            "+--------------------+----------+-----------+-----------+---------------------------------------------+-------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------+\n",
            "| samtools2878c4d2f4 | samtools | 1.16.1    | Failed    | workbench.jobs('samtools2878c4d2f4').logs() | samtools flagstat /output/mapping/SRR10261595_sorted.bam -o /output/mapping/SRR10261595_sorted.bam_map_stat.txt                                                                     |\n",
            "+--------------------+----------+-----------+-----------+---------------------------------------------+-------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------+\n",
            "| samtools6756a0c0ca | samtools | 1.16.1    | Failed    | workbench.jobs('samtools6756a0c0ca').logs() | samtools flagstat /output/mapping/SRR10261596_sorted.bam -o /output/mapping/SRR10261596_sorted.bam_map_stat.txt                                                                     |\n",
            "+--------------------+----------+-----------+-----------+---------------------------------------------+-------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------+\n",
            "| samtools8b5d512a4b | samtools | 1.16.1    | Succeeded | workbench.jobs('samtools8b5d512a4b').logs() | samtools flagstat /output/mapping/SRR10261591_sorted.bam -O /output/mapping/SRR10261591_sorted.bam_map_stat.txt                                                                     |\n",
            "+--------------------+----------+-----------+-----------+---------------------------------------------+-------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------+\n",
            "| samtools5f9f73ca88 | samtools | 1.16.1    | Succeeded | workbench.jobs('samtools5f9f73ca88').logs() | samtools flagstat /output/mapping/SRR10261592_sorted.bam -O /output/mapping/SRR10261592_sorted.bam_map_stat.txt                                                                     |\n",
            "+--------------------+----------+-----------+-----------+---------------------------------------------+-------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------+\n",
            "| samtools66b6457818 | samtools | 1.16.1    | Succeeded | workbench.jobs('samtools66b6457818').logs() | samtools flagstat /output/mapping/SRR10261593_sorted.bam -O /output/mapping/SRR10261593_sorted.bam_map_stat.txt                                                                     |\n",
            "+--------------------+----------+-----------+-----------+---------------------------------------------+-------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------+\n",
            "| samtools5bd1a4179a | samtools | 1.16.1    | Succeeded | workbench.jobs('samtools5bd1a4179a').logs() | samtools flagstat /output/mapping/SRR10261594_sorted.bam -O /output/mapping/SRR10261594_sorted.bam_map_stat.txt                                                                     |\n",
            "+--------------------+----------+-----------+-----------+---------------------------------------------+-------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------+\n",
            "| samtools406425553c | samtools | 1.16.1    | Succeeded | workbench.jobs('samtools406425553c').logs() | samtools flagstat /output/mapping/SRR10261595_sorted.bam -O /output/mapping/SRR10261595_sorted.bam_map_stat.txt                                                                     |\n",
            "+--------------------+----------+-----------+-----------+---------------------------------------------+-------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------+\n",
            "| samtoolsf8e4a8bf9c | samtools | 1.16.1    | Succeeded | workbench.jobs('samtoolsf8e4a8bf9c').logs() | samtools flagstat /output/mapping/SRR10261596_sorted.bam -O /output/mapping/SRR10261596_sorted.bam_map_stat.txt                                                                     |\n",
            "+--------------------+----------+-----------+-----------+---------------------------------------------+-------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------+\n",
            "| samtools3272c58d91 | samtools | 1.16.1    | Succeeded | workbench.jobs('samtools3272c58d91').logs() | samtools view -bh /output/mapping/SRR10261591_sorted.bam I II III IV V VI VII VIII IX X XI XII XIII XIV XV XVI > /output/mapping/SRR10261591_sorted_SC_subset.bam                   |\n",
            "+--------------------+----------+-----------+-----------+---------------------------------------------+-------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------+\n",
            "| samtoolsed6745cf4f | samtools | 1.16.1    | Succeeded | workbench.jobs('samtoolsed6745cf4f').logs() | samtools view -bh /output/mapping/SRR10261592_sorted.bam I II III IV V VI VII VIII IX X XI XII XIII XIV XV XVI > /output/mapping/SRR10261592_sorted_SC_subset.bam                   |\n",
            "+--------------------+----------+-----------+-----------+---------------------------------------------+-------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------+\n",
            "| samtools363c580dda | samtools | 1.16.1    | Succeeded | workbench.jobs('samtools363c580dda').logs() | samtools view -bh /output/mapping/SRR10261593_sorted.bam I II III IV V VI VII VIII IX X XI XII XIII XIV XV XVI > /output/mapping/SRR10261593_sorted_SC_subset.bam                   |\n",
            "+--------------------+----------+-----------+-----------+---------------------------------------------+-------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------+\n",
            "| samtools19db8479a4 | samtools | 1.16.1    | Succeeded | workbench.jobs('samtools19db8479a4').logs() | samtools view -bh /output/mapping/SRR10261594_sorted.bam I II III IV V VI VII VIII IX X XI XII XIII XIV XV XVI > /output/mapping/SRR10261594_sorted_SC_subset.bam                   |\n",
            "+--------------------+----------+-----------+-----------+---------------------------------------------+-------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------+\n",
            "| samtoolsead541a4a0 | samtools | 1.16.1    | Succeeded | workbench.jobs('samtoolsead541a4a0').logs() | samtools view -bh /output/mapping/SRR10261595_sorted.bam I II III IV V VI VII VIII IX X XI XII XIII XIV XV XVI > /output/mapping/SRR10261595_sorted_SC_subset.bam                   |\n",
            "+--------------------+----------+-----------+-----------+---------------------------------------------+-------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------+\n",
            "| samtools5f0de3187c | samtools | 1.16.1    | Succeeded | workbench.jobs('samtools5f0de3187c').logs() | samtools view -bh /output/mapping/SRR10261596_sorted.bam I II III IV V VI VII VIII IX X XI XII XIII XIV XV XVI > /output/mapping/SRR10261596_sorted_SC_subset.bam                   |\n",
            "+--------------------+----------+-----------+-----------+---------------------------------------------+-------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------+\n",
            "| picardf263edacba   | picard   | 2.27.3    | Succeeded | workbench.jobs('picardf263edacba').logs()   | picard MarkDuplicates I=/output/mapping/SRR10261591_sorted_SC_subset.bam O=/output/mapping/SRR10261591_sorted_SC_subset_dedup.bam M=/output/mapping/SRR10261591_markdup_metrics.txt |\n",
            "+--------------------+----------+-----------+-----------+---------------------------------------------+-------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------+\n",
            "| picard07e2dcf5b3   | picard   | 2.27.3    | Succeeded | workbench.jobs('picard07e2dcf5b3').logs()   | picard MarkDuplicates I=/output/mapping/SRR10261592_sorted_SC_subset.bam O=/output/mapping/SRR10261592_sorted_SC_subset_dedup.bam M=/output/mapping/SRR10261592_markdup_metrics.txt |\n",
            "+--------------------+----------+-----------+-----------+---------------------------------------------+-------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------+\n",
            "| picard6b72f9b849   | picard   | 2.27.3    | Succeeded | workbench.jobs('picard6b72f9b849').logs()   | picard MarkDuplicates I=/output/mapping/SRR10261593_sorted_SC_subset.bam O=/output/mapping/SRR10261593_sorted_SC_subset_dedup.bam M=/output/mapping/SRR10261593_markdup_metrics.txt |\n",
            "+--------------------+----------+-----------+-----------+---------------------------------------------+-------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------+\n",
            "| picard1c85b0c584   | picard   | 2.27.3    | Succeeded | workbench.jobs('picard1c85b0c584').logs()   | picard MarkDuplicates I=/output/mapping/SRR10261594_sorted_SC_subset.bam O=/output/mapping/SRR10261594_sorted_SC_subset_dedup.bam M=/output/mapping/SRR10261594_markdup_metrics.txt |\n",
            "+--------------------+----------+-----------+-----------+---------------------------------------------+-------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------+\n",
            "| picard5257f16275   | picard   | 2.27.3    | Succeeded | workbench.jobs('picard5257f16275').logs()   | picard MarkDuplicates I=/output/mapping/SRR10261595_sorted_SC_subset.bam O=/output/mapping/SRR10261595_sorted_SC_subset_dedup.bam M=/output/mapping/SRR10261595_markdup_metrics.txt |\n",
            "+--------------------+----------+-----------+-----------+---------------------------------------------+-------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------+\n",
            "| picardd95fd25d39   | picard   | 2.27.3    | Succeeded | workbench.jobs('picardd95fd25d39').logs()   | picard MarkDuplicates I=/output/mapping/SRR10261596_sorted_SC_subset.bam O=/output/mapping/SRR10261596_sorted_SC_subset_dedup.bam M=/output/mapping/SRR10261596_markdup_metrics.txt |\n",
            "+--------------------+----------+-----------+-----------+---------------------------------------------+-------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------+\n",
            "| samtools1b483dff23 | samtools | 1.16.1    | Succeeded | workbench.jobs('samtools1b483dff23').logs() | samtools index /output/mapping/SRR10261591_sorted_SC_subset_dedup.bam                                                                                                               |\n",
            "+--------------------+----------+-----------+-----------+---------------------------------------------+-------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------+\n",
            "| samtools25d2a9f546 | samtools | 1.16.1    | Succeeded | workbench.jobs('samtools25d2a9f546').logs() | samtools index /output/mapping/SRR10261592_sorted_SC_subset_dedup.bam                                                                                                               |\n",
            "+--------------------+----------+-----------+-----------+---------------------------------------------+-------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------+\n",
            "| samtools3612dccf5d | samtools | 1.16.1    | Succeeded | workbench.jobs('samtools3612dccf5d').logs() | samtools index /output/mapping/SRR10261593_sorted_SC_subset_dedup.bam                                                                                                               |\n",
            "+--------------------+----------+-----------+-----------+---------------------------------------------+-------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------+\n",
            "| samtools753709a242 | samtools | 1.16.1    | Succeeded | workbench.jobs('samtools753709a242').logs() | samtools index /output/mapping/SRR10261594_sorted_SC_subset_dedup.bam                                                                                                               |\n",
            "+--------------------+----------+-----------+-----------+---------------------------------------------+-------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------+\n",
            "| samtoolsa55646d3b2 | samtools | 1.16.1    | Succeeded | workbench.jobs('samtoolsa55646d3b2').logs() | samtools index /output/mapping/SRR10261595_sorted_SC_subset_dedup.bam                                                                                                               |\n",
            "+--------------------+----------+-----------+-----------+---------------------------------------------+-------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------+\n",
            "| samtools900bf1b244 | samtools | 1.16.1    | Succeeded | workbench.jobs('samtools900bf1b244').logs() | samtools index /output/mapping/SRR10261596_sorted_SC_subset_dedup.bam                                                                                                               |\n",
            "+--------------------+----------+-----------+-----------+---------------------------------------------+-------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------+\n",
            "| samtools1483d8f4a0 | samtools | 1.16.1    | Succeeded | workbench.jobs('samtools1483d8f4a0').logs() | samtools flagstat /output/mapping/SRR10261591_sorted_SC_subset_dedup.bam > /output/mapping/SRR10261591_sorted_SC_subset_dedup_map_stats.txt                                         |\n",
            "+--------------------+----------+-----------+-----------+---------------------------------------------+-------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------+\n",
            "| samtoolsdafa14c676 | samtools | 1.16.1    | Succeeded | workbench.jobs('samtoolsdafa14c676').logs() | samtools flagstat /output/mapping/SRR10261592_sorted_SC_subset_dedup.bam > /output/mapping/SRR10261592_sorted_SC_subset_dedup_map_stats.txt                                         |\n",
            "+--------------------+----------+-----------+-----------+---------------------------------------------+-------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------+\n",
            "| samtoolsb416573031 | samtools | 1.16.1    | Succeeded | workbench.jobs('samtoolsb416573031').logs() | samtools flagstat /output/mapping/SRR10261593_sorted_SC_subset_dedup.bam > /output/mapping/SRR10261593_sorted_SC_subset_dedup_map_stats.txt                                         |\n",
            "+--------------------+----------+-----------+-----------+---------------------------------------------+-------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------+\n",
            "| samtools5d45e60585 | samtools | 1.16.1    | Succeeded | workbench.jobs('samtools5d45e60585').logs() | samtools flagstat /output/mapping/SRR10261594_sorted_SC_subset_dedup.bam > /output/mapping/SRR10261594_sorted_SC_subset_dedup_map_stats.txt                                         |\n",
            "+--------------------+----------+-----------+-----------+---------------------------------------------+-------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------+\n",
            "| samtoolsce897b1ffd | samtools | 1.16.1    | Succeeded | workbench.jobs('samtoolsce897b1ffd').logs() | samtools flagstat /output/mapping/SRR10261595_sorted_SC_subset_dedup.bam > /output/mapping/SRR10261595_sorted_SC_subset_dedup_map_stats.txt                                         |\n",
            "+--------------------+----------+-----------+-----------+---------------------------------------------+-------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------+\n",
            "| samtools87ce6ef2b2 | samtools | 1.16.1    | Succeeded | workbench.jobs('samtools87ce6ef2b2').logs() | samtools flagstat /output/mapping/SRR10261596_sorted_SC_subset_dedup.bam > /output/mapping/SRR10261596_sorted_SC_subset_dedup_map_stats.txt                                         |\n",
            "+--------------------+----------+-----------+-----------+---------------------------------------------+-------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------+\n",
            "| macs2e84080522e    | macs2    | 2.2.7.1   | Running   | workbench.jobs('macs2e84080522e').logs()    | macs2 callpeak -t /output/mapping/SRR10261591_sorted_SC_subset_dedup.bam -f BAMPE -n /output/mapping/SRR10261591_peak -g 12000000 --keep-dup all                                    |\n",
            "+--------------------+----------+-----------+-----------+---------------------------------------------+-------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------+\n",
            "| macs201a56f821b    | macs2    | 2.2.7.1   | Running   | workbench.jobs('macs201a56f821b').logs()    | macs2 callpeak -t /output/mapping/SRR10261592_sorted_SC_subset_dedup.bam -f BAMPE -n /output/mapping/SRR10261592_peak -g 12000000 --keep-dup all                                    |\n",
            "+--------------------+----------+-----------+-----------+---------------------------------------------+-------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------+\n",
            "| macs2fbdb269e4c    | macs2    | 2.2.7.1   | Running   | workbench.jobs('macs2fbdb269e4c').logs()    | macs2 callpeak -t /output/mapping/SRR10261593_sorted_SC_subset_dedup.bam -f BAMPE -n /output/mapping/SRR10261593_peak -g 12000000 --keep-dup all                                    |\n",
            "+--------------------+----------+-----------+-----------+---------------------------------------------+-------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------+\n",
            "| macs2d97578d208    | macs2    | 2.2.7.1   | Succeeded | workbench.jobs('macs2d97578d208').logs()    | macs2 callpeak -t /output/mapping/SRR10261594_sorted_SC_subset_dedup.bam -f BAMPE -n /output/mapping/SRR10261594_peak -g 12000000 --keep-dup all                                    |\n",
            "+--------------------+----------+-----------+-----------+---------------------------------------------+-------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------+\n",
            "| macs22cb5b57017    | macs2    | 2.2.7.1   | Succeeded | workbench.jobs('macs22cb5b57017').logs()    | macs2 callpeak -t /output/mapping/SRR10261595_sorted_SC_subset_dedup.bam -f BAMPE -n /output/mapping/SRR10261595_peak -g 12000000 --keep-dup all                                    |\n",
            "+--------------------+----------+-----------+-----------+---------------------------------------------+-------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------+\n",
            "| macs28aaaa7748d    | macs2    | 2.2.7.1   | Succeeded | workbench.jobs('macs28aaaa7748d').logs()    | macs2 callpeak -t /output/mapping/SRR10261596_sorted_SC_subset_dedup.bam -f BAMPE -n /output/mapping/SRR10261596_peak -g 12000000 --keep-dup all                                    |\n",
            "+--------------------+----------+-----------+-----------+---------------------------------------------+-------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------+\n"
          ]
        }
      ]
    },
    {
      "cell_type": "markdown",
      "source": [
        "Let’s explore some important parameters here. First, -f BAMPE tells MACS2 that the data is paired-end and directs the software to analyze only properly-paired reads, which are abundant in our sample. Second, the option –keep-dup all tells MACS2 that PCR duplicates were handled by other software (picard MarkDuplicates in our case). And third, the option -g 12000000 denotes the effective genome size of S. cerevisiae (~12Mb).\n",
        "\n",
        "\n",
        "Upon completion, MACS2 generates several outputs for each sample:\n",
        "\n",
        "*  peaks.narrowPeak file: A BED6+4 format file which contains the peak locations together with peak summit, p-value, and q-value\n",
        "*  peaks.xls file: A tabular file which contains information about called peaks. Additional information includes pileup and fold enrichment\n",
        "*  summits.bed file: Peak summit locations for every peak.\n",
        "\n",
        "We have now performed an ATAC-seq analysis and obtained the open chromatin peaks (i.e. open chromatin profiles) of S. cerevisiae in its parental and hybrid state. In the next step, we will answer our initial question of whether and to what extent S. cerevisiae changes its chromatin accessibility profiles upon hybridization. To do this, we will need to perform quality control on the obtained peaks and compare them across the two conditions (parent vs hybrid). But first, we will visualize the data generated so far in this tutorial."
      ],
      "metadata": {
        "id": "StURSQcwS4ci"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "\n",
        "#Alignment & Peak Data Visualization\n",
        "\n",
        "Our analyses to this point have yielded an array of abstract (and somewhat cryptic) results, including alignment bam files, peak files in bed format, summit files, etc. To properly integrate and interpret these results, we need data visualization. For this step, we’ll use Integrative Genomics Viewer (IGV) v. 2.4.9. IGV is a powerful and intuitive Java program that makes NGS data clear, compelling, and tangible. IGV can be downloaded on your local machine from [here](https://software.broadinstitute.org/software/igv/download). Additionally, you will need to download the genome and bam files that we will use below.\n",
        "\n",
        "After downloading, open IGV and go to Genome > Load Genomes for File > select the genome fasta file of S. cerevisiae. The chromosomes of the yeast will appear at the top of IGV. Then go to File > Load from file > select the gtf annotation file of S. cerevisiae. The genome annotations will appear underneath the chromosome, showing genes and other annotated elements of the genome. Next, go to File > Load from file > select one of the bam files (for example SRR10261591_SC_subset_dedup.bam).\n",
        "\n",
        "IGV will now create two tracks for this file – one with a coverage line plot (at the top), and a second, larger track visualizing all mapped reads. To see both of them, select one chromosome from the dropdown menu, then zoom in using the + sign at the top-right corner. After zooming in sufficiently, the coverage plot and the reads will appear. The peaks of ATAC-seq should already be visible based on the coverage plot. As expected, these peaks appear in intergenic regions (i.e. between genes), where the regulatory nucleosome-free regions are located.\n",
        "\n",
        "To visualize the peaks found by MACS2, go to File > Load from file > select SRR10261591_peak_peaks.narrowPeak file and in the same way select SRR10261591_peak_summits.bed. After loading both files, IGV will show the peaks’ regions and their summits on the corresponding tracks (tracks in IGV can be dragged and rearranged). It should look like this (Fig. 7):\n",
        "\n",
        "\n",
        "\n",
        "![](https://drive.google.com/uc?export=view&id=11gIG1POPOHCJFefHAwf8BAYfZozXVVEy)\n",
        "\n",
        "\n",
        "*Fig. 7: Screenshot from IGV showing the integrated data for SRR10261591 sample. From the top downward, the tracks are: peaks identified by MACS2, peak summits, coverage line plot based on the bam file, aligned ATAC-seq reads, and genomic annotations of S. cerevisiae.*\n",
        "\n",
        "\n",
        "In this plot, you can see an intergenic ATAC-seq peak (top row, dark blue) and its summit which clearly coincides with the coverage plot of the bam file. This peak likely regulates the expression of both genes AAC3 and RPL19A, harboring some regulatory elements such as TF binding sites or promoter regions.\n",
        "\n",
        "Again, IGV is powerful software that will integrate and clarify your results. It can visualize bam files, genome and genome annotation files, peak files, and more. IGV is a feature-reach software, and we recommend spending more time with it to explore its many functionalities.\n",
        "\n"
      ],
      "metadata": {
        "id": "vqeHusDXXuht"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "#Analysis of chromatin accessibility profile\n",
        "\n",
        "Now that the open chromatin states of the yeast in both conditions have been identified with MACS2, we can compare the two profiles. We will use an R Bioconductor package DiffBind v. 3.4.11 for this step.\n",
        "\n",
        "First, it is necessary to perform the quality control of the open chromatin profiles among biological replicates in both conditions of S. cerevisiae to see whether there are any outlier samples which might bias the downstream analysis.\n",
        "\n",
        "We’ll use the script diffbind_qc.R. Run the script with:"
      ],
      "metadata": {
        "id": "LUZQyrRPaD65"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "workbench.run(tool='r-base', full_command='Rscript /input/rscript/diffbind_qc.R') "
      ],
      "metadata": {
        "id": "qjaryYNLaT-9",
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "outputId": "34f7b946-1cd0-44da-9f12-04e38d1066b6"
      },
      "execution_count": null,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "+------------------+--------+-----------+----------+-------------------------------------------+--------------------------------------+\n",
            "| Job ID           | Tool   | Version   | Status   | Get Logs                                  | Full Command                         |\n",
            "+==================+========+===========+==========+===========================================+======================================+\n",
            "| r-base1d0bc3e9b6 | r-base | 4.2.3     | Queued   | workbench.jobs('r-base1d0bc3e9b6').logs() | Rscript /input/rscript/diffbind_qc.R |\n",
            "+------------------+--------+-----------+----------+-------------------------------------------+--------------------------------------+\n"
          ]
        }
      ]
    },
    {
      "cell_type": "code",
      "source": [
        "workbench.jobs()"
      ],
      "metadata": {
        "id": "V8lktZ6RCilI"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "Let's explore what's happening in the script above. First, it loads the DiffBind library, then loads the previously created peak.xls files into DiffBind. Then it performs quality control analyses among replicates of each condition, creating Venn diagrams showing how many peaks are in common across replicates:\n",
        "\n",
        "\n",
        "```\n",
        "### Parallelize Package installation\n",
        "options(Ncpus = 10)\n",
        "\n",
        "if (!require(\"BiocManager\", quietly = TRUE))\n",
        "  install.packages(\"BiocManager\")\n",
        "BiocManager::install(\"DiffBind\")\n",
        "\n",
        "\n",
        "library(DiffBind)\n",
        "\n",
        "\n",
        "### Set the working directory\n",
        "setwd(\"/output/mapping/\")\n",
        "\n",
        "\n",
        "### Load the peak files from MACS2 to DiffBind\n",
        "SC_vs_hybrid<-dba.peakset(NULL,\n",
        "                          peaks=\"SRR10261591_peak_peaks.xls\",\n",
        "                          peak.caller=\"macs\", sampID=\"SC1\",condition = \"Parent\", replicate=1, bamReads = \"SRR10261591_SC_subset_dedup.bam\")\n",
        "SC_vs_hybrid<-dba.peakset(SC_vs_hybrid,\n",
        "                          peaks=\"SRR10261592_peak_peaks.xls\",\n",
        "                          peak.caller=\"macs\", sampID=\"SC2\",condition = \"Parent\", replicate=2, bamReads = \"SRR10261592_SC_subset_dedup.bam\")\n",
        "SC_vs_hybrid<-dba.peakset(SC_vs_hybrid,\n",
        "                          peaks=\"SRR10261593_peak_peaks.xls\",\n",
        "                          peak.caller=\"macs\", sampID=\"SC3\",condition = \"Parent\", replicate=3, bamReads = \"SRR10261593_SC_subset_dedup.bam\")\n",
        "\n",
        "\n",
        "SC_vs_hybrid<-dba.peakset(SC_vs_hybrid,\n",
        "                          peaks=\"SRR10261594_peak_peaks.xls\",\n",
        "                          peak.caller=\"macs\", sampID=\"hybrid1\",condition = \"Hybrid\", replicate=1, bamReads = \"SRR10261594_SC_subset_dedup.bam\")\n",
        "SC_vs_hybrid<-dba.peakset(SC_vs_hybrid,\n",
        "                          peaks=\"SRR10261595_peak_peaks.xls\",\n",
        "                          peak.caller=\"macs\", sampID=\"hybrid2\",condition = \"Hybrid\", replicate=2, bamReads = \"SRR10261595_SC_subset_dedup.bam\")\n",
        "SC_vs_hybrid<-dba.peakset(SC_vs_hybrid,\n",
        "                          peaks=\"SRR10261596_peak_peaks.xls\",\n",
        "                          peak.caller=\"macs\", sampID=\"hybrid3\",condition = \"Hybrid\", replicate=3, bamReads = \"SRR10261596_SC_subset_dedup.bam\")\n",
        "\n",
        "\n",
        "\n",
        "### Plot Venn diagrams of replicates\n",
        "png('/output/mapping/SC_replicates_atacseq.png', units=\"in\", width=7, heigh=7, res=800)\n",
        "dba.plotVenn(SC_vs_hybrid,SC_vs_hybrid$masks$Parent, main = \"Open chromatic region overlaps in S. cerevisiae replicates\")\n",
        "dev.off()\n",
        "\n",
        "png('/output/mapping/hybrid_SC_replicates_atacseq.png', units=\"in\", width=7, heigh=7, res=800)\n",
        "dba.plotVenn(SC_vs_hybrid,SC_vs_hybrid$masks$Hybrid, main = \"Open chromatic region overlaps in the hybrid replicates\")\n",
        "dev.off()\n",
        "```\n",
        "\n",
        "As you can see in Fig. 8a and 8b (for parental and hybrid respectively), the vast majority of peaks are present between all replicates in both conditions. This indicates a high concordance between biological replicates in this experiment. For example, in the case of the parent, there are 2911 peaks which are in common between all replicates, while there are only 42 peaks which are specific to the replicate 1 (i.e. SC1).\n",
        "\n",
        "\n",
        "![](https://drive.google.com/uc?export=view&id=11l7kaZ9YTjCRUw3EWG_cMqQLw2wBJ_tF)\n",
        "\n",
        "*Fig. 8: Quality control analysis of open chromatin peaks between the replicates in both conditions of S. cerevisiae. (a) Overlap between the peaks in parental samples (“SC” stands for S. cerevisiae); (b) Overlap between the peaks in hybrid samples*\n",
        "\n",
        "\n",
        "DiffBind also created a heatmap plot of Pearson correlations based on previously obtained alignment bam files by calculating the number of reads in each peak across all samples (Fig. 9):\n",
        "\n",
        "```\n",
        "SC_vs_hybrid_counts<-dba.count(SC_vs_hybrid, bParallel = TRUE, score=DBA_SCORE_READS)\n",
        "\n",
        "png('/output/mapping/heatmap.png', units=\"in\", width=7, heigh=7, res=800)\n",
        "plot(SC_vs_hybrid_counts, main=\"Correlation plot of studied samples\")\n",
        "dev.off()\n",
        "```\n",
        "\n",
        "\n",
        "![](https://drive.google.com/uc?export=view&id=11o3F2Q24benkJr4831z5kKth_l26I8z6)\n",
        "\n",
        "*Fig. 9: Heatmap of Pearson correlations between the studies samples.*\n",
        "\n",
        "\n",
        "\n",
        "In the heatmap above, you can see that the second replicate of the hybrid represents an outlier. We will thus remove it from further analysis.\n",
        "\n",
        "In summary, DiffBind has revealed a high overall concordance between replicates of both conditions. It also identified one outlier sample in the hybrid."
      ],
      "metadata": {
        "id": "vfktdlUmahnO"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "###Occupancy and differential chromatin accessibility analysis\n",
        "Having performed quality control between replicates of both conditions, it is now possible to compare the open chromatin profiles of S. cerevisiae between parent and hybrid. We will also use DiffBind for this step.\n",
        "\n",
        "First, the so-called occupancy analysis will qualitatively assess how many open chromatin regions are in common between the two conditions – i.e. it will find the consensus peakset between the conditions. Second, we will perform the so-called differential accessibility analysis – the quantitative analysis showing which peaks within the consensus peakset are differentially more/less accessible due to a condition (in this case being a hybrid).\n",
        "\n",
        "Note: To assess differential accessibility, the DiffBind package relies on DESeq2 software. DESeq2 was originally designed for differential gene expression analysis, but works well for other NGS applications such as ATAC-seq and ChIP-seq.\n",
        "\n",
        "To run both analyses, we will use the script occupancy_and_diff_accessibility.R. "
      ],
      "metadata": {
        "id": "EcJDAFTobx3v"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "workbench.run(tool='r-base', full_command='Rscript /input/rscript/occupancy_and_diff_accessibility.R') "
      ],
      "metadata": {
        "id": "6S1Pt6eVcn4B",
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "outputId": "a6cc4dee-21bb-47f9-d374-f066e1d0c143"
      },
      "execution_count": null,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "+------------------+--------+-----------+----------+-------------------------------------------+-----------------------------------------------------------+\n",
            "| Job ID           | Tool   | Version   | Status   | Get Logs                                  | Full Command                                              |\n",
            "+==================+========+===========+==========+===========================================+===========================================================+\n",
            "| r-base6d3d902c9e | r-base | 4.2.3     | Queued   | workbench.jobs('r-base6d3d902c9e').logs() | Rscript /input/rscript/occupancy_and_diff_accessibility.R |\n",
            "+------------------+--------+-----------+----------+-------------------------------------------+-----------------------------------------------------------+\n"
          ]
        }
      ]
    },
    {
      "cell_type": "code",
      "source": [
        "workbench.jobs('r-basebe2176d99c').__str__()"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "zWzDGyK4L8Cs",
        "outputId": "d222d0c7-54d3-402a-e7a4-f932e380fb55"
      },
      "execution_count": null,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "+------------------+--------+-----------+----------+-------------------------------------------+-----------------------------------------------------------+\n",
            "| Job ID           | Tool   | Version   | Status   | Get Logs                                  | Full Command                                              |\n",
            "+==================+========+===========+==========+===========================================+===========================================================+\n",
            "| r-basebe2176d99c | r-base | 4.2.3     | Failed   | workbench.jobs('r-basebe2176d99c').logs() | Rscript /input/rscript/occupancy_and_diff_accessibility.R |\n",
            "+------------------+--------+-----------+----------+-------------------------------------------+-----------------------------------------------------------+\n"
          ]
        }
      ]
    },
    {
      "cell_type": "markdown",
      "source": [
        "\n",
        "As in our previous quality control analysis, the script first loads the peakset and alignment files to DiffBind. (Note that we've discarded replicate 2 from this analysis because it was an outlier.) DiffBind then performs the occupancy analysis, generating a Venn diagram showing the overlap between the peaksets of the two conditions – i.e. the consensus peakset (Fig. 10). In other words, it reveals the peaks present in both conditions:\n",
        "\n",
        "\n",
        "\n",
        "```\n",
        "if (!require(\"BiocManager\", quietly = TRUE))\n",
        "  install.packages(\"BiocManager\")\n",
        "BiocManager::install(\"DiffBind\")\n",
        "\n",
        "\n",
        "library(DiffBind)\n",
        "### Set the working directory\n",
        "setwd(\"/output/mapping/\")\n",
        "\n",
        "### Load the peak files from MACS2 to DiffBind\n",
        "SC_vs_hybrid<-dba.peakset(NULL,\n",
        "                          peaks=\"SRR10261591_peak_peaks.xls\",\n",
        "                          peak.caller=\"macs\", sampID=\"SC1\",condition = \"Parent\", replicate=1, bamReads = \"SRR10261591_SC_subset_dedup.bam\")\n",
        "SC_vs_hybrid<-dba.peakset(SC_vs_hybrid,\n",
        "                          peaks=\"SRR10261592_peak_peaks.xls\",\n",
        "                          peak.caller=\"macs\", sampID=\"SC2\",condition = \"Parent\", replicate=2, bamReads = \"SRR10261592_SC_subset_dedup.bam\")\n",
        "SC_vs_hybrid<-dba.peakset(SC_vs_hybrid,\n",
        "                          peaks=\"SRR10261593_peak_peaks.xls\",\n",
        "                          peak.caller=\"macs\", sampID=\"SC3\",condition = \"Parent\", replicate=3, bamReads = \"SRR10261593_SC_subset_dedup.bam\")\n",
        "\n",
        "\n",
        "SC_vs_hybrid<-dba.peakset(SC_vs_hybrid,\n",
        "                          peaks=\"SRR10261594_peak_peaks.xls\",\n",
        "                          peak.caller=\"macs\", sampID=\"hybrid1\",condition = \"Hybrid\", replicate=1, bamReads = \"SRR10261594_SC_subset_dedup.bam\")\n",
        "SC_vs_hybrid<-dba.peakset(SC_vs_hybrid,\n",
        "                          peaks=\"SRR10261596_peak_peaks.xls\",\n",
        "                          peak.caller=\"macs\", sampID=\"hybrid3\",condition = \"Hybrid\", replicate=2, bamReads = \"SRR10261596_SC_subset_dedup.bam\")\n",
        "\n",
        "\n",
        "### Occupancy analysis\n",
        "\n",
        "SC_vs_hybrid_occupancy<-dba.peakset(SC_vs_hybrid, consensus = DBA_CONDITION,minOverlap = 0.33)\n",
        "\n",
        "\n",
        "png('/output/mapping/occupancy_consensus_peaks_SC_vs_hybrid.png', units=\"in\", width=7, heigh=7, res=800)\n",
        "dba.plotVenn(SC_vs_hybrid_occupancy,SC_vs_hybrid_occupancy$masks$Consensus, main = \"Overlap of consensus peaks of SC and hybrid\")\n",
        "dev.off()\n",
        "```\n",
        "\n",
        "\n",
        "\n",
        "![](https://drive.google.com/uc?export=view&id=11x6P04PcQgjhBt3ETrF9rPd6WttZnAZD)\n",
        "\n",
        "*Fig. 10: Occupancy analysis and identification of the consensus peaks (n=3060) between the two conditions of S. cerevisiae.*\n",
        "\n",
        "\n",
        "\n",
        "Surprisingly, the majority of peaks are in common between S. cerevisiae in the parental and hybrid state. This indicates that chromatin accessibility profiles are largely preserved during the transition to hybrid state.\n",
        "\n",
        "\n",
        "\n",
        "DiffBind also performs the differential accessibility analysis (sometimes called the affinity analysis), identifying differentially accessible peaks between the parent and hybrid from the consensus peakset. The output of this analysis is a table with the consensus peaks, their coordinates and chromosomes, accessibility values for both conditions, fold-change of these values, and statistical significance of the fold-change. To obtain biologically relevant results, we apply a fold-change filter of more or less than 1.0, as well as a False Discovery Rate (FDR) < 0.05. The resulting table is exported to the diffbind directory (file diff_regions_hybrid_vs_parent.txt).\n",
        "\n",
        "\n",
        "\n",
        "```\n",
        "#### Affinity analysis - differential accessibility analysis\n",
        "counts<-dba.count(SC_vs_hybrid, bParallel = TRUE, score=DBA_SCORE_READS)\n",
        "SC_vs_hybrid_counts<-dba.contrast(counts, categories=DBA_CONDITION,minMembers = 2)\n",
        "\n",
        "\n",
        "SC_vs_hybrid_analyzed<-dba.analyze(SC_vs_hybrid_counts)\n",
        "\n",
        "SC_vs_hybrid_DE_peaks <- dba.report(SC_vs_hybrid_analyzed)\n",
        "SC_vs_hybrid_DE_peaks_DF<-as.data.frame(SC_vs_hybrid_DE_peaks)\n",
        "# Filtering\n",
        "write.table(SC_vs_hybrid_DE_peaks_DF[abs(SC_vs_hybrid_DE_peaks_DF$Fold)>1 & SC_vs_hybrid_DE_peaks_DF$FDR<0.05,], \"/output/mapping/diff_regions_hybrid_vs_parent.txt\", sep = \"\\t\", quote = F)\n",
        "\n",
        "```\n",
        "\n",
        "\n",
        "In our analysis, DiffBind finds only one peak out of more than 3000 consensus peaks which is differentially open in the hybrid state compared to the parent. This strongly supports the hypothesis that open chromatin profiles (and thus gene regulation) are largely preserved in S. cerevisiae upon formation of a hybrid organism. (For scientific implications of this result, see the original paper by [Hovhannisyan et al, 2020](https://pubmed.ncbi.nlm.nih.gov/32457798/).)\n",
        "\n",
        "\n"
      ],
      "metadata": {
        "id": "WHmIPk6fcjKY"
      }
    }
  ]
}