ATACseq-all4samples_v3.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# ATAC-Seq analysis example\n",
    "### Original authors\n",
    "Délara Sabéran-Djoneidi, Anne-Laure Schang, Kate Wooley-Allen, Sascha Ott\n",
    "### Update\n",
    "Alix Silvert, Magali Hennion\n",
    "## Intro and general description\n",
    "\n",
    "This document aims to be a detailed example of an ATAC-Seq analysis used in a research context. The goal here is to help you draft your own experiments from this particular base, by walking you through the steps we took and by pointing out the parts where the specificity of your experimental set will require some tinkering.\n",
    "\n",
    "## Clone the git repository\n",
    "\n",
    "The documents related to this tutorial are available at https://github.com/parisepigenetics/ATAC-seq. The rest of this description assumes that you have cloned this repository (using `git clone`) and that you are in `ATAC-seq` folder. The `pwd` command allows you to check where you are. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pwd"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## if you're not in the ATAC-seq folder, use cd to enter it.\n",
    "## cd PATHto/ATAC-seq\n",
    "cd ATAC-seq"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Load Conda environment\n",
    "\n",
    "In order to be sure to have the appropriate tools to run this analysis, we provide a Conda environment summarizing them. In order to use it, please download [Conda](https://docs.conda.io/projects/conda/en/latest/user-guide/install/download.html) and run the following line.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "conda env create -f conda_environment.yml"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "All the remaining of the analysis must be run in the ATAC-seq Conda environment. You have to activate the environment. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "conda activate ATAC-Seq"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Download test dataset\n",
    "\n",
    "To test this tutorial, you can download a small test dataset. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "wget https://raw.githubusercontent.com/nf-core/test-datasets/atacseq/testdata/SRR182215{3,4,7,8}_1.fastq.gz https://raw.githubusercontent.com/nf-core/test-datasets/atacseq/testdata/SRR182215{3,4,7,8}_2.fastq.gz -P data/raw_reads/"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ls data/raw_reads"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And the corresponding reference files. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "wget https://github.com/nf-core/test-datasets/raw/atacseq/reference/genome.fa -P reference"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Initial quality control\n",
    "\n",
    "In order to see how much the pre-processing improves the data, it is a good practice to look at various statistics on the raw dataset. FastQC provides a good summary of those."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mkdir -p output/stats/raw_fastQC/\n",
    "fastqc -o output/stats/raw_fastQC/ data/raw_reads/SRR182215{3,4,7,8}_1.fastq.gz data/raw_reads/SRR182215{3,4,7,8}_2.fastq.gz "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If you have multiple samples, multiQC is a good way to look at the global quality over your samples."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "multiqc output/stats/raw_fastQC -o output/stats/raw_fastQC"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Trimming\n",
    "\n",
    "The fastq files are trimmed using Trimmomatic v0.39 to remove any adapter sequences in the reads caused by read through associated with DNA fragments shorter in size than the read length being sequenced.\n",
    "\n",
    "This particular step need to be adapted depending on your exact data and need. We found these options to be the most efficient in the case of ATAC-Seq. If you want to know more we encourage you to read the [documentation](http://www.usadellab.org/cms/uploads/supplementary/Trimmomatic/TrimmomaticManual_V0.32.pdf) of the Trimmomatic tool in order to find the best possible use for your case.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mkdir -p output/trimmed_files\n",
    "samples=\"SRR1822153 SRR1822154 SRR1822157 SRR1822158\"\n",
    "for sample in $samples; do\n",
    "   trimmomatic PE -threads 2 data/raw_reads/${sample}_1.fastq.gz data/raw_reads/${sample}_2.fastq.gz \\\n",
    "   output/trimmed_files/${sample}_1_trimmed.fastq output/trimmed_files/${sample}_1_trimmed_unpaired.fa \\\n",
    "   output/trimmed_files/${sample}_2_trimmed.fastq output/trimmed_files/${sample}_2_trimmed_unpaired.fa \\\n",
    "   ILLUMINACLIP:data/adapters/NexteraPE-PE.fa:2:30:10:1:true TRAILING:3 SLIDINGWINDOW:4:15;\n",
    "done"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Alignement\n",
    "\n",
    "This step must also be adapted to your particular experiment, namely the genomes used may vary and the aligner used as well as the alignment options can be optimized.\n",
    "\n",
    "In our case, alignment is performed by bowtie2.4.4. \n",
    "\n",
    "### Bowtie2 Index\n",
    "\n",
    "In order to align reads to genome(s) and particular cellular compartments (i.e. mitochondria), it is necessary to build a bowtie index for each of them and potential confounders. In the end, only the reads mapping solely to the target of interest will be saved. In our case, the mitochondrial hits will be removed.\n",
    "\n",
    "Bowtie2 necessitate and index to be built before the alignment can be done. By default, the exact method used is chosen depending on the machine doing the computation, please read the [documentation of bowtie2-build](http://bowtie-bio.sourceforge.net/bowtie2/manual.shtml#the-bowtie2-build-indexer) for further information, and be careful when repeating a study on different machines.\n",
    "\n",
    "If you already have a bowtie index built for this type of alignment, this step can be skipped.\n",
    "\n",
    "In this example, we start from the fasta of a complete genome that we split into chromosomes using a small bash script."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "bash scripts/split_reference.sh reference/genome.fa"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For bigger genome, it might be more efficient to download the fasta already split by chromosomes, for example at https://hgdownload.soe.ucsc.edu/downloads.html. Choose your organism and then click on `Sequence data by chromosome`. \n",
    "\n",
    "Then we build one index for the mitochondrie, and one index for the other chromosomes. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cd reference\n",
    "bowtie2-build genome.fa_MT yeast_MT > MT_index.log\n",
    "bowtie2-build genome.fa_I,genome.fa_II,genome.fa_III,genome.fa_IV,genome.fa_V,genome.fa_VI,genome.fa_VII,genome.fa_VIII,genome.fa_IX,genome.fa_X,genome.fa_XI,genome.fa_XII,genome.fa_XIII,genome.fa_XIV,genome.fa_XV,genome.fa_XVI yeast_NC > NC_index.log\n",
    "cd .. # to go back to ATAC-seq folder. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Alignment\n",
    "\n",
    "#### Removal of Mitochondrial DNA\n",
    "\n",
    "Let's breakdown the various options used:\n",
    "\n",
    "* `-p 4` We are using 4 threads for alignment. This option can also be accessed via `--threads`.\n",
    "* `-X 2000` The maximum fragment length accepted between paired end. This number should be adapted depending on your sequencing protocol.\n",
    "* `--very-sensitive` Preset option corresponding to `-D 20 -R 3 -N 0 -L 20 -i S,1,0.50`.\n",
    "* `-x` Selected bowtie index.\n",
    "*  `-1 and -2` The input fastq files in the case of paired-end alignment\n",
    "*  `-S` File in which to write the output SAM file of the aligned reads. This file is not useful but printing it allows us to use less memory.\n",
    "*  `--un-conc` We only keep pairs that DIDN'T align on mitochondrial genome.\n",
    "\n",
    "for more informations and option, you can run `bowtie2 --help`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mkdir -p tmp\n",
    "for sample in $samples; do\n",
    "  echo ${sample} >> output/stats/Bowtie2_output_MT.txt;\n",
    "  bowtie2 -p 4 -X 2000 --very-sensitive -x reference/yeast_MT \\\n",
    "  -1 output/trimmed_files/${sample}_1_trimmed.fastq \\\n",
    "  -2 output/trimmed_files/${sample}_2_trimmed.fastq \\\n",
    "  -S tmp/${sample}_MT.sam --un-conc-gz tmp/${sample}_notMT.fastq.gz \\\n",
    "  2>> output/stats/Bowtie2_output_MT.txt;\n",
    "done"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`output/stats/Bowtie2_output_MT.txt` contains information about the mapping including the % of mapped reads on the mitochondrial genome (about 1% for these test samples). \n",
    "The file containing the alignment on the mitochondria won't be used, so we remove it. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rm tmp/SRR182215{3,4,7,8}_MT.sam"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Alignment on chromosomal genome\n",
    "\n",
    "We now align the remaining reads (stored in `tmp` folder) on the nuclear genome, in a very similar manner.  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ls tmp"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "mkdir -p output/aligned_files\n",
    "for sample in $samples; do\n",
    "   echo ${sample} >> output/stats/Bowtie2_output_notMT.txt;\n",
    "   bowtie2 -p 4 -X 2000 --very-sensitive -x reference/yeast_NC \\\n",
    "   -1 tmp/${sample}_notMT.fastq.1.gz \\\n",
    "   -2 tmp/${sample}_notMT.fastq.2.gz \\\n",
    "   -S output/aligned_files/${sample}_notMT_aligned.sam \\\n",
    "   2>> output/stats/Bowtie2_output_notMT.txt;\n",
    "done"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`output/stats/Bowtie2_output_notMT.txt` contains information about the mapping including the % of mapped reads on the nuclear genome (about 95% in this test dataset).  \n",
    "The alignment SAM files are stored in `output/aligned_files`. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ls output/aligned_files"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can now remove temporary files. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rm -fr tmp/"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We now sort the sam file, convert it to bam, and then index it."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "for sample in $samples; do\n",
    "   echo ${sample};\n",
    "   samtools sort output/aligned_files/${sample}_notMT_aligned.sam \\\n",
    "   -o output/aligned_files/${sample}_notMT_aligned_sorted.bam;\n",
    "   samtools index output/aligned_files/${sample}_notMT_aligned_sorted.bam;\n",
    "done"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The sorted BAM files and the indexes are in `output/aligned_files`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ls output/aligned_files"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Quality Controls and Filters\n",
    "\n",
    "The various outputs of this part, located in `output/stats/` must be studied to see whether something has gone wrong with the alignment or previous steps. What to do when you have an unexpected behavior will very much depend on your experiment.\n",
    "\n",
    "### SAMtools flags\n",
    "SAMtools flags give several information about a read, such as whether it is paired or not. Here we check the distribution of those flags within our sample.\n",
    "\n",
    "The [wikipedia page](https://en.wikipedia.org/wiki/SAM_(file_format)) about SAM files does explicit the meaning of each SAMtools flag."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for sample in $samples; do\n",
    "   echo ${sample};\n",
    "   samtools flagstat output/aligned_files/${sample}_notMT_aligned_sorted.bam \\\n",
    "   > output/stats/${sample}_aligned_samflags.txt;\n",
    "   done"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Samtools MAPQ score\n",
    "The MAPQ score is related to the probability that a particular read is misaligned. The higher the MAPQ score is, the lower the risk that the read is misaligned is. In our experiment we used a particular threshold of 22, corresponding to a probability of misalignment of 10<sup>-12</sup>. \n",
    "Let's create the MAPQ score distribution using SAMtools for reads that have the flag that corresponds to 2 (each segment properly aligned according to aligner)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for sample in $samples; do\n",
    "   echo ${sample};\n",
    "   samtools view -f2 output/aligned_files/${sample}_notMT_aligned_sorted.bam | \\\n",
    "   awk '{print $5}' | sort -n | uniq -c | sed 's/^ *//g'  \\\n",
    "   > output/stats/${sample}_aligned_aligned_mapq.txt;\n",
    "done"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Using Bowtie2, the reads that mapped several times have a MAPQ score of 1. The reads with MAPQ > 22 are uniquely mapped with good quality. We will now filter the reads to keep only these ones.   \n",
    "\n",
    "### Filter base on the quality scores threshold you choose"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for sample in $samples; do\n",
    "   echo ${sample};\n",
    "   samtools view -f2 -q22 -b output/aligned_files/${sample}_notMT_aligned_sorted.bam \\\n",
    "   > output/aligned_files/${sample}_notMT_aligned_filtered.bam;\n",
    "   done"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Template length\n",
    "\n",
    "Let's observe the Template Length. This is the distance between the mapped end of the DNA fragment and the mapped start of the DNA fragment, inclusively. If the mapping is correct, it corresponds to the size of the original DNA fragment. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for sample in $samples; do\n",
    "   echo ${sample};\n",
    "   samtools view output/aligned_files/${sample}_notMT_aligned_filtered.bam | \\\n",
    "   cut -f 9| sed 's/^-//' | sort -n | uniq -c |sed  's/^ *//g' \\\n",
    "   > output/stats/${sample}_aligned_length.txt;\n",
    "done"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### fastQC - general quality check"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mkdir output/stats/fastQC/\n",
    "fastqc -o output/stats/fastQC/ output/aligned_files/SRR182215{3,4,7,8}_notMT_aligned_filtered.bam"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Same as before, if you are treating several files at once, multiQC can be a good way to check whether the quality vary depending on each sample/condition, which might affect the results later on."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "multiqc output/stats/fastQC/ -o output/stats/fastQC/"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Peak calling\n",
    "\n",
    "To call the ATAC-Seq peaks, we use MACS2. It might be more efficient to pool the reads from similar samples. You can do it by giving several bam files when calling MACS2. If you have several conditions that you want to compare, keep them separated and create one peak file per condition.\n",
    "\n",
    "This part might need to be adapted to your particular case, so let's break down the various arguments : \n",
    "* `-t output/aligned_filtered/SRR182215{3,4,7,8}_notMT_aligned_filtered.bam` name of the file (or files) on which to do the peak calling. If you input several files, they will be pooled to give a single callpeak file.\n",
    "* `-f BAM` Type of the output file desired\n",
    "* `-g 1.2e7` Effective length of the genome (S.cerevisiae here). This has to be adapted to your genome (some are precompiled : hs: 2.7e9; mm: 1.87e9; ce: 9e7; dm: 1.2e8)\n",
    "* `-q 0.05` minimum q-value for report of the peaks\n",
    "* `-n output/ATAC_Seq_peaks/SRR182215{3,4,7,8}` base name of the output, several files containing a variety results will be generated\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here we create one files with the peaks from all the samples. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mkdir -p output/ATAC_Seq_peaks\n",
    "macs2 callpeak -t output/aligned_files/SRR182215{3,4,7,8}_notMT_aligned_filtered.bam \\\n",
    "-f BAMPE -g 1.2e7 -q 0.05 -n output/ATAC_Seq_peaks/All_samples"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Mask regions\n",
    "\n",
    "Some peaks might be in regions known to create artefactual peaks. It is important to remove the corresponding peaks from the results. \n",
    "\n",
    "The blacklists can be found [here](http://mitra.stanford.edu/kundaje/akundaje/release/blacklists/) or from [ENCODE](https://github.com/Boyle-Lab/Blacklist/tree/master/lists) and should be adapted to your project. In yeast, there is no such region and we only remove rDNA peaks.  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "bedtools subtract -A -a output/ATAC_Seq_peaks/All_samples_summits.bed \\\n",
    "-b data/blacklist/sacCer3_rDNA.bed > output/ATAC_Seq_peaks/All_samples_summits_BL.bed"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "bedtools subtract -A -a output/ATAC_Seq_peaks/All_samples_peaks.narrowPeak \\\n",
    "-b data/blacklist/sacCer3_rDNA.bed > output/ATAC_Seq_peaks/All_samples_peaks_BL.bed"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Peak coverage per sample\n",
    "\n",
    "We use bedtools coverage to get the peak coverage of every sample. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "mkdir output/peak_coverage/\n",
    "for sample in $samples; do\n",
    "   echo ${sample};\n",
    "   bedtools coverage -a output/ATAC_Seq_peaks/All_samples_peaks_BL.bed \\\n",
    "   -b output/aligned_files/${sample}_notMT_aligned_filtered.bam \\\n",
    "   > output/peak_coverage/${sample}_peaks.bed;\n",
    "   done"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Those outputs can be used to run differential analysis between 2 conditions with DEseq2 or edgeR. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "head output/peak_coverage/SRR1822153_peaks.bed"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Output description\n",
    "\n",
    "Columns in the peaks.bed files: \n",
    "1. chr \n",
    "2. start \n",
    "3. end \n",
    "4. peak_ID \n",
    "5. MACS_score  \n",
    "6. MACS_strand \n",
    "7. MACS_fold_enrichment \n",
    "8. MACS_-log10(pval) \n",
    "9. MACS_-log10(qval) \n",
    "10. MACS_summit (0-based offset from chrom start)\n",
    "11. reads_in_peak \n",
    "12. peak_non-zero_bases \n",
    "13. peak_length \n",
    "14. peak_fraction_covered\n",
    "\n",
    "Columns 1 to 10 are from MACS: https://hbctraining.github.io/Intro-to-ChIPseq/lessons/05_peak_calling_macs.html. \n",
    "\n",
    "Columns 11 to 14 are from bedtools coverage: https://bedtools.readthedocs.io/en/latest/content/tools/coverage.html\n",
    "\n",
    "11. The number of features in B that overlapped (by at least one base pair) the A interval.\n",
    "12. The number of bases in A that had non-zero coverage from features in B.\n",
    "13. The length of the entry in A.\n",
    "14. The fraction of bases in A that had non-zero coverage from features in B.\n",
    "\n",
    "## To go further\n",
    "You can now run the R Jupyter Notebook `Differential_accessibility.ipynb` to compare different conditions. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Bash",
   "language": "bash",
   "name": "bash"
  },
  "language_info": {
   "codemirror_mode": "shell",
   "file_extension": ".sh",
   "mimetype": "text/x-sh",
   "name": "bash"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}