Introduction

In many eukaryotic organisms, such as humans, the genome is tightly packed and organized with the help of nucleosomes (chromatin). A nucleosome is a complex formed by eight histone proteins that is wrapped with ~147bp of DNA. When the DNA is being actively transcribed into RNA, the DNA will be opened and loosened from the nucleosome complex. Many factors, such as the chromatin structure, the position of the nucleosomes, and histone modifications, play an important role in the organization and accessibility of the DNA. Consequently, these factors are also important for the activation and inactivation of genes. Assay for Transposase-Accessible Chromatin using sequencing (ATAC-Seq) is a method to investigate the accessibility of chromatin and thus a method to determine regulatory mechanisms of gene expression. The method can help identify promoter regions and potential enhancers and silencers. A promoter is the DNA region close to the transcription start site (TSS). It contains binding sites for transcription factors that will recruit the RNA polymerase. An enhancer is a DNA region that can be located up to 1 Mb downstream or upstream of the promoter. When transcription factors bind an enhancer and contact a promoter region, the transcription of the gene is increased. In contrast, a silencer decreases or inhibits the gene’s expression. ATAC-Seq has become popular for identifying accessible regions of the genome as it’s easier, faster and requires less cells than alternative techniques, such as FAIRE-Seq and DNase-Seq.

With ATAC-Seq, to find accessible (open) chromatin regions, the genome is treated with a hyperactive derivative of the Tn5 transposase. A transposase can bind to a transposable element, which is a DNA sequence that can change its position (jump) within a genome (read the two links to get a deeper insight). During ATAC-Seq, the modified Tn5 inserts DNA sequences corresponding to truncated Nextera adapters into open regions of the genome and concurrently, the DNA is sheared by the transposase activity. The read library is then prepared for sequencing, including PCR amplification with full Nextera adapters and purification steps. Paired-end reads are recommended for ATAC-Seq for the reasons described here.

In this tutorial we will use data from the study of Buenrostro et al. 2013, the first paper on the ATAC-Seq method. The data is from a human cell line of purified CD4+ T cells, called GM12878. The original dataset had 2 x 200 million reads and would be too big to process in a training session, so we downsampled the original dataset to 200,000 randomly selected reads. We also added about 200,000 reads pairs that will map to chromosome 22 to have a good profile on this chromosome, similar to what you might get with a typical ATAC-Seq sample (2 x 20 million reads in original FASTQ). Furthermore, we want to compare the predicted open chromatin regions to the known binding sites of CTCF, a DNA-binding protein implicated in 3D structure: CTCF. CTCF is known to bind to thousands of sites in the genome and thus it can be used as a positive control for assessing if the ATAC-Seq experiment is good quality. Good ATAC-Seq data would have accessible regions both within and outside of TSS, for example, at some CTCF binding sites. For that reason, we will download binding sites of CTCF identified by ChIP in the same cell line from ENCODE (ENCSR000AKB, dataset ENCFF933NTR).

When working with real data

When you use your own data we suggest you to use this workflow which includes the same steps but is compatible with replicates. If you do not have any control data you can import and edit this workflow, removing all steps with the controls. Controls for the ATAC-Seq procedure are not commonly performed, as discussed here, but could be ATAC-Seq of purified DNA.

Agenda

In this tutorial, we will cover:

1. Introduction
2. Preprocessing
   1. Get Data
   2. Quality Control
   3. Trimming Reads
3. Mapping
   1. Mapping Reads to Reference Genome
4. Filtering Mapped Reads
   1. Filter Uninformative Reads
   2. Filter Duplicate Reads
   3. Check Insert Sizes
5. Peak calling
   1. Call Peaks
6. Visualisation of Coverage
   1. Prepare the Datasets
   2. Create heatmap of coverage at TSS with deepTools
   3. Visualise Regions with pyGenomeTracks
7. Conclusion

Comment: Results may vary

Your results may be slightly different from the ones presented in this tutorial due to differing versions of tools, reference data, external databases, or because of stochastic processes in the algorithms.

Preprocessing

Get Data

We first need to download the sequenced reads (FASTQs) as well as other annotation files. Then, to increase the number of reads that will map to the reference genome (here human genome version 38, GRCh38/hg38), we need to preprocess the reads.

Hands-on: Data upload

1. Create a new history for this tutorial

   Tip: Creating a new history

   Click the new-history icon at the top of the history panel.

   If the new-history is missing:

   1. Click on the galaxy-gear icon (History options) on the top of the history panel
   2. Select the option Create New from the menu

2. Import the files from Zenodo or from the shared data library (GTN - Material -> epigenetics -> ATAC-Seq data analysis):

   https://zenodo.org/record/3862793/files/ENCFF933NTR.bed.gz
   https://zenodo.org/record/3862793/files/SRR891268_chr22_enriched_R1.fastq.gz
   https://zenodo.org/record/3862793/files/SRR891268_chr22_enriched_R2.fastq.gz
   

   Tip: Importing via links

   • Copy the link location

   • Open the Galaxy Upload Manager (galaxy-upload on the top-right of the tool panel)

   • Select Paste/Fetch Data

   • Paste the link into the text field

   • Press Start

   • Close the window

   Tip: Importing data from a data library

   As an alternative to uploading the data from a URL or your computer, the files may also have been made available from a shared data library:

   • Go into Shared data (top panel) then Data libraries
   • Navigate to the correct folder as indicated by your instructor
   • Select the desired files
   • Click on the To History button near the top and select as Datasets from the dropdown menu
   • In the pop-up window, select the history you want to import the files to (or create a new one)
   • Click on Import

3. Add a tag called #SRR891268_R1 to the R1 file and a tag called #SRR891268_R2 to the R2 file.

   Tip: Adding a tag

   • Click on the dataset
   • Click on galaxy-tags Edit dataset tags

   • Add a tag starting with #

     Tags starting with # will be automatically propagated to the outputs of tools using this dataset.

   • Check that the tag is appearing below the dataset name

4. Check that the datatype of the 2 FASTQ files is fastqsanger.gz and the peak file (ENCFF933NTR.bed.gz) is encodepeak. If they are not then change the datatype as described below.

   Tip: Changing the datatype

   • Click on the galaxy-pencil pencil icon for the dataset to edit its attributes
   • In the central panel, click on the galaxy-chart-select-data Datatypes tab on the top
   • Select datatypes
     • tip: you can start typing the datatype into the field to filter the dropdown menu
   • Click the Save button

Comment: FASTQ format

If you are not familiar with FASTQ format, see the Quality Control tutorial

Comment: BED / encode narrowPeak format

If you are not familiar with BED format or encode narrowPeak format, see the BED Format

We will visualise regions later in the analysis and obtain the gene information now. We will get information for chromosome 22 genes (names of transcripts and genomic positions) using the UCSC tool.

Hands-on: Obtain Annotation for hg38 genes

1. UCSC Main table browser Tool: ucsc_table_direct1 with the following parameters:
   • “clade”: Mammal
   • “genome”: Human
   • “assembly”: Dec. 2013 (GRCh38/hg38)
   • “group”: Genes and Gene Prediction
   • “track”: All GENCODE V37
   • “table”: Basic
   • “region”: position chr22
   • “output format”: all fields from selected table
   • “Send output to”: Galaxy
2. Click get output

3. Click Send query to Galaxy

   This table contains all the information but is not in a BED format. To transform it into BED format we will cut out the required columns and rearrange:

4. Cut columns from a table Tool: Cut1 with the following parameters:
   • “Cut columns”: c3,c5,c6,c13,c12,c4
   • “Delimited by”: Tab
   • param-file “From”: UCSC Main on Human: wgEncodeGencodeBasicV37 (chr22:1-50,818,468)

5. Check the contents of your file, is this as you expect it to be?

   Question: Expected output

   Our goal here was to convert the data to BED format.

   1. Which columns do you expect in your file? (Tip: read about BED format)
   2. Does your file look like a valid BED format?

   Solution

   1. We expect at least 3 columns, chromosome - start - end, and possibly more optional columns
   2. Your file should look something like this:

      chr22	10736170	10736283	U2	0	-
      chr22	11066417	11068174	CU104787.1	0	+
      chr22	11249808	11249959	5_8S_rRNA	0	-
      chr22	15273854	15273961	U6	0	+
      [..]
      

   • Troubleshooting: Is your second column the Strand column?
     • Make sure you used the correct Cut tool (the one that matches the tool name mentioned in the previous step exactly)
     • There is another tool with (cut) behind the title, we do NOT want to use this tool in this step.
   • Tip: Always check your output files to make sure they match your expectations!

6. Rename galaxy-pencil the dataset as chr22 genes

   Tip: Renaming a dataset

   • Click on the galaxy-pencil pencil icon for the dataset to edit its attributes
   • In the central panel, change the Name field
   • Click the Save button

7. Change galaxy-pencil its datatype to BED

   Tip: Changing the datatype

   • Click on the galaxy-pencil pencil icon for the dataset to edit its attributes
   • In the central panel, click on the galaxy-chart-select-data Datatypes tab on the top
   • Select bed
     • tip: you can start typing the datatype into the field to filter the dropdown menu
   • Click the Save button

8. Click on the galaxy-eye (eye) icon of the file. It should have now column names and they should match the content.

Comment: Gene file

The chr22 genes BED we produced only contains the start, the end, the name, and the strand of each transcript. It does not contain exon information. To be able to have the exon information, you could use a GTF file which can be downloaded from the gencode website but this file would include the information for the whole genome and would slow the analysis.

Quality Control

The first step is to check the quality of the reads and the presence of the Nextera adapters. When we perform ATAC-Seq, we can get DNA fragments of about 40 bp if two adjacent Tn5 transposases cut the DNA Adey et al. 2010. This can be smaller than the sequencing length so we expect to have Nextera adapters at the end of those reads. We can assess the reads with FastQC.

Hands-on: Task description

1. FastQC Tool: toolshed.g2.bx.psu.edu/repos/devteam/fastqc/fastqc/0.72+galaxy1 with the following parameters:
   • “Short read data from your current history”: Choose here either only the SRR891268_R1 file with param-file or use param-files Multiple datasets to choose both SRR891268_R1 and SRR891268_R2.
2. Inspect the web page output of FastQC tool for the SRR891268_R1 sample. Check what adapters are found at the end of the reads.

Question

1. How many reads are in the FASTQ?
2. Which sections have a warning?

Solution

1. There are 285247 reads.

2. The 3 steps below have warnings:

   1. Per base sequence content

      It is well known that the Tn5 has a strong sequence bias at the insertion site. You can read more about it in Green et al. 2012.

   2. Sequence Duplication Levels

      The read library quite often has PCR duplicates that are introduced simply by the PCR itself. We will remove these duplicates later on.

   3. Overrepresented sequences

      One sequence is over represented: you have 306 reads which are exactly the sequence of the Nextera adapter. They correspond to adapters amplified head-to-head. 306 is really low (only 0.1% of reads).

Comment: FastQC Results

This is what you should expect from the Adapter Content section:

The FastQC web page Adapter Content section shows the presence of Nextera Transposase Sequence in the reads. We will remove the adapters with Cutadapt.

Trimming Reads

To trim the adapters we provide the Nextera adapter sequences to Cutadapt. These adapters are shown in the image below.

The forward and reverse adapters are slightly different. We will also trim low quality bases at the ends of the reads (quality less than 20). We will only keep reads that are at least 20 bases long. We remove short reads (< 20bp) as they are not useful, they will either be thrown out by the mapping or may interfere with our results at the end.

Hands-on: Task description

1. Cutadapt Tool: toolshed.g2.bx.psu.edu/repos/lparsons/cutadapt/cutadapt/1.16.5 with the following parameters:
   • “Single-end or Paired-end reads?”: Paired-end
     • param-file “FASTQ/A file #1”: select SRR891268_R1
     • param-file “FASTQ/A file #2”: select SRR891268_R2
     • In “Read 1 Options”:
       • In “3’ (End) Adapters”:
         • param-repeat “Insert 3’ (End) Adapters”
           • “Source”: Enter custom sequence
             • “Enter custom 3’ adapter name (Optional if Multiple output is ‘No’)”: Nextera R1
             • “Enter custom 3’ adapter sequence”: CTGTCTCTTATACACATCTCCGAGCCCACGAGAC
     • In “Read 2 Options”:
       • In “3’ (End) Adapters”:
         • param-repeat “Insert 3’ (End) Adapters”
           • “Source”: Enter custom sequence
             • “Enter custom 3’ adapter name (Optional)”: Nextera R2
             • “Enter custom 3’ adapter sequence”: CTGTCTCTTATACACATCTGACGCTGCCGACGA
   • In “Filter Options”:
     • “Minimum length”: 20
   • In “Read Modification Options”:
     • “Quality cutoff”: 20
   • In “Output Options”:
     • “Report”: Yes
2. Click on the galaxy-eye (eye) icon of the report and read the first lines.

Comment: Cutadapt Results

You should get similar output to this from Cutadapt:

Question

1. What percentage of reads contain adapters?
2. What percentage of reads are still longer than 20bp after the trimming?

Solution

1. ~14%
2. ~99%

Hands-on: Check Adapter Removal with FastQC

1. FastQC Tool: toolshed.g2.bx.psu.edu/repos/devteam/fastqc/fastqc/0.72+galaxy1 with the following parameters:
   • “Short read data from your current history”: select the output of Cutadapt param-files Multiple datasets to choose both Read 1 Output and Read 2 Output.
2. Click on the galaxy-eye (eye) icon of the report and read the first lines.

Comment: FastQC Results

Now, you should see under Overrepresented sequences that there is no more overrepresented sequences and under Adapter Content that the Nextera adapters are no longer present.

However, you may have noticed that you have a new section with warning: Sequence Length Distribution. This is expected as you trimmed part of the reads.

Mapping

Mapping Reads to Reference Genome

Next we map the trimmed reads to the human reference genome. Here we will use Bowtie2. We will extend the maximum fragment length (distance between read pairs) from 500 to 1000 because we know some valid read pairs are from this fragment length. We will use the --very-sensitive parameter to have more chance to get the best match even if it takes a bit longer to run. We will run the end-to-end mode because we trimmed the adapters so we expect the whole read to map, no clipping of ends is needed. Regarding the genome to choose. The hg38 version of the human genome contains alternate loci. This means that some region of the genome are present both in the canonical chromosome and on its alternate loci. The reads that map to these regions would map twice. To be able to filter reads falling into repetitive regions but keep reads falling into regions present in alternate loci, we will map on the Canonical version of hg38 (only the chromosome with numbers, chrX, chrY, and chrM).

Comment: Dovetailing

We will allow dovetailing of read pairs with Bowtie2. This is because adapters are removed by Cutadapt only when at least 3 bases match the adapter sequence, so it is possible that after trimming a read can contain 1-2 bases of adapter and go beyond it’s mate start site. For example, if the first mate in the read pair is: GCTATGAAGAATAGGGCGAAGGGGCCTGCGGCGTATTCGATGTTGAAGCT and the second mate is CTTCAACATCGAATACGCCGCAGGCCCCTTCGCCCTATTCTTCATAGCCT, where both contain 2 bases of adapter sequence, they will not be trimmed by Cutadapt and will map this way:

<--------------------Mate 1-----------------------
AGCTTCAACATCGAATACGCCGCAGGCCCCTTCGCCCTATTCTTCATAGC
  CTTCAACATCGAATACGCCGCAGGCCCCTTCGCCCTATTCTTCATAGCCT
  ----------------------Mate 2--------------------->

This is what we call dovetailing and we want to consider this pair as a valid concordant alignment.

Hands-on: Mapping reads to reference genome

1. Bowtie2 Tool: toolshed.g2.bx.psu.edu/repos/devteam/bowtie2/bowtie2/2.4.2+galaxy0 with the following parameters:
   • “Is this single or paired library”: Paired-end
     • param-file “FASTQ/A file #1”: select the output of Cutadapt tool “Read 1 Output”
     • param-file “FASTQ/A file #2”: select the output of Cutadapt tool “Read 2 Output”
     • “Do you want to set paired-end options?”: Yes
       • “Set the maximum fragment length for valid paired-end alignments”: 1000
       • “Allow mate dovetailing”: Yes
   • “Will you select a reference genome from your history or use a built-in index?”: Use a built-in genome index
     • “Select reference genome”: Human (Homo sapiens): hg38 Canonical
   • “Set read groups information?”: Do not set
   • “Select analysis mode”: 1: Default setting only
     • “Do you want to use presets?”: Very sensitive end-to-end (--very-sensitive)
   • “Do you want to tweak SAM/BAM Options?”: No
   • “Save the bowtie2 mapping statistics to the history”: Yes
2. Click on the galaxy-eye (eye) icon of the mapping stats.

Comment: Bowtie2 Results

You should get similar results to this from Bowtie2:

Question

What percentage of read pairs mapped concordantly?

Solution

54.8+42.87=97.67%

Comment: On the number of uniquely mapped.

You might be surprised by the number of uniquely mapped compared to the number of multi-mapped reads (reads mapping to more than one location in the genome). One of the reasons is that we have used the parameter --very-sensitive. Bowtie2 considers a read as multi-mapped even if the second hit has a much lower quality than the first one. Another reason is that we have reads that map to the mitochondrial genome. The mitochondrial genome has a lot of regions with similar sequence.

Filtering Mapped Reads

Filter Uninformative Reads

We apply some filters to the reads after the mapping. ATAC-Seq datasets can have a lot of reads that map to the mitchondrial genome because it is nucleosome-free and thus very accessible to Tn5 insertion. The mitchondrial genome is uninteresting for ATAC-Seq so we remove these reads. We also remove reads with low mapping quality and reads that are not properly paired.

Hands-on: Filtering of uninformative reads

1. Filter BAM datasets on a variety of attributes Tool: toolshed.g2.bx.psu.edu/repos/devteam/bamtools_filter/bamFilter/2.4.1 with the following parameters:
   • param-file “BAM dataset(s) to filter”: Select the output of Bowtie2 tool “alignments”
   • In “Condition”:
     • param-repeat “Insert Condition”
       • In “Filter”:
         • param-repeat “Insert Filter”
           • “Select BAM property to filter on”: mapQuality
             • “Filter on read mapping quality (phred scale)”: >=30
         • param-repeat “Insert Filter”
           • “Select BAM property to filter on”: isProperPair
             • “Select properly paired reads”: Yes
         • param-repeat “Insert Filter”
           • “Select BAM property to filter on”: reference
             • “Filter on the reference name for the read”: !chrM
   • “Would you like to set rules?”: No
2. Click on the input and the output BAM files of the filtering step. Check the size of the files.

Question

1. Based on the file size, what proportion of alignments was removed (approximately)?
2. Which parameter should be modified if you are interested in repetitive regions?

Solution

1. The original BAM file is 28.1 MB, the filtered one is 15.2 MB. Approximately half of the alignments were removed.

2. You should modify the mapQuality criteria and decrease the threshold.

High numbers of mitochondrial reads can be a problem in ATAC-Seq. Some ATAC-Seq samples have been reported to be 80% mitochondrial reads and so wet-lab methods have been developed to deal with this issue Corces et al. 2017 and Litzenburger et al. 2017. It can be a useful QC to assess the number of mitochondrial reads. However, it does not predict the quality of the rest of the data. It is just that sequencing reads have been wasted.

Tip: Getting the number of mitochondrial reads

To get the number of reads that mapped to the mitochondrial genome (chrM) you can run Samtools idxstats Tool: toolshed.g2.bx.psu.edu/repos/devteam/samtools_idxstats/samtools_idxstats/2.0.3 on the output of Bowtie2 tool “alignments”. The columns of the output are: chromosome name, chromosome length, number of reads mapping to the chromosome, number of unaligned mate whose mate is mapping to the chromosome. The first 2 lines of the result would be (after using Sort Tool: toolshed.g2.bx.psu.edu/repos/bgruening/text_processing/tp_sort_header_tool/1.1.1 ):

There are 220 000 reads which map to chrM and 165 000 which map to chr22.

Filter Duplicate Reads

Because of the PCR amplification, there might be read duplicates (different reads mapping to exactly the same genomic region) from overamplification of some regions. As the Tn5 insertion is random within an accessible region, we do not expect to see fragments with the same coordinates. We consider such fragments to be PCR duplicates. We will remove them with Picard MarkDuplicates.

Hands-on: Remove duplicates

1. MarkDuplicates Tool: toolshed.g2.bx.psu.edu/repos/devteam/picard/picard_MarkDuplicates/2.18.2.2 with the following parameters:
   • param-file “Select SAM/BAM dataset or dataset collection”: Select the output of Filter tool “BAM”
   • “If true do not write duplicates to the output file instead of writing them with appropriate flags set”: Yes

   Comment: Defaults of MarkDuplicates

   By default, the tool will only “Mark” the duplicates. This means that it will change the Flag of the duplicated reads to enable them to be filtered afterwards. We use the parameter “If true do not write duplicates to the output file instead of writing them with appropriate flags set” to directly remove the duplicates.

2. Click on the galaxy-eye (eye) icon of the MarkDuplicate metrics.

Comment: MarkDuplicates Results

You should get similar output to this from MarkDuplicates:

Tip: Formatting the MarkDuplicate metrics for readability

1. Select lines that match an expression Tool: Grep1 with the following parameters:
   • param-file “Select lines from”: Select the output of MarkDuplicates tool
   • “that: Matching
   • “the pattern: (Library|LIBRARY)
2. Check that the datatype is tabular. If not, change it. <blockquote class="tip"> <div id="tip-changing-the-datatype-2" class="box-title"> Tip: Changing the datatype</div> <ul> <li>Click on the galaxy-pencil pencil icon for the dataset to edit its attributes</li> <li>In the central panel, click on the galaxy-chart-select-data Datatypes tab on the top</li> <li>Select tabular <ul> <li>tip: you can start typing the datatype into the field to filter the dropdown menu</li> </ul> </li> <li>Click the Save button</li></ul></blockquote>
3. Transpose rows/columns in a tabular file Tool: toolshed.g2.bx.psu.edu/repos/iuc/datamash_transpose/datamash_transpose/1.1.0 :
   • param-file “Select lines from”: Select the output of Select tool

Question

1. How many pairs were in the input?
2. How many pairs are duplicates?

Solution

1. 135813
2. 3584

Once again, if you have a high number of replicates it does not mean that your data are not good, it just means that you sequenced too much compared to the diversity of the library you generated. Consequently, libraries with a high portion of duplicates should not be resequenced as this would not increase the amount of data.

Check Insert Sizes

We will check the insert sizes with Paired-end histogram of insert size frequency. The insert size is the distance between the R1 and R2 read pairs. This tells us the size of the DNA fragment the read pairs came from. The fragment length distribution of a sample gives a very good indication of the quality of the ATAC-Seq.

Hands-on: Plot the distribution of fragment sizes.

1. Paired-end histogram Tool: toolshed.g2.bx.psu.edu/repos/iuc/pe_histogram/pe_histogram/1.0.1 with the following parameters:
   • param-file “BAM file”: Select the output of MarkDuplicates tool “BAM output”
   • “Lower bp limit (optional)”: 0
   • “Upper bp limit (optional)”: 1000
2. Click on the galaxy-eye (eye) icon of the lower one of the 2 outputs (the png file).

Comment: CollectInsertSizeMetrics Results

This is what you get from CollectInsertSizeMetrics:

Question

Could you guess what the peaks at approximately 50bp, 200bp, 400bp and 600bp correspond to?

Solution

The first peak (50bp) corresponds to where the Tn5 transposase inserted into nucleosome-free regions. The second peak (a bit less than 200bp) corresponds to where Tn5 inserted around a single nucleosome. The third one (around 400bp) is where Tn5 inserted around two adjacent nucleosomes and the fourth one (around 600bp) is where Tn5 inserted around three adjacent nucleosomes.

This fragment size distribution is a good indication if your experiment worked or not. In absence of chromatin (without nucleosome), this is the profile you would get:

Here are examples of Fragment size distributions of ATAC-Seq which were very noisy:

A final example of a Fragment size distribution of a very good ATAC-Seq, even if we cannot see the third nucleosome “peak”.

Peak calling

Call Peaks

We have now finished the data preprocessing. Next, in order to find regions corresponding to potential open chromatin regions, we want to identify regions where reads have piled up (peaks) greater than the background read coverage. The tools which are currently used are Genrich and MACS2. MACS2 is more widely used. Genrich has a mode dedicated to ATAC-Seq but is still not published and the more reads you have, the less peaks you get (see the issue here). That’s why we will not use Genrich in this tutorial.

At this step, two approaches exists:

• The first one is to select only paired whose fragment length is below 100bp corresponding to nucleosome-free regions and to use a peak calling like you would do for a ChIP-seq, joining signal between mates. The disadvantages of this approach is that you can only use it if you have paired-end data and you will miss small open regions where only one Tn5 bound.
• The second one chosen here is to use all reads to be more exhaustive. In this approach, it is very important to re-center the signal of each reads on the 5’ extremity (read start site) as this is where Tn5 cuts. Indeed, you want your peaks around the nucleosomes and not directly on the nucleosome:

Comment: on Tn5 insertion

When Tn5 cuts an accessible chromatin locus it inserts adapters separated by 9bp (Kia et al. 2017):

This means in order to have the read start site reflecting the centre of where Tn5 bound, the reads on the positive strand should be shifted 4 bp to the right and reads on the negative strands should be shifted 5 bp to the left as in Buenrostro et al. 2013. Genrich can apply these shifts when ATAC-seq mode is selected. In most cases, we do not have 9bp resolution so we don’t take it into account but if you are interested in the footprint, this is important.

If we only assess the coverage of the 5’ extremity of the reads, the data would be too sparse and it would be impossible to call peaks. Thus, we will extend the start sites of the reads by 200bp (100bp in each direction) to assess coverage.

Using MACS2

We convert the BAM file to BED format because when we set the extension size in MACS2, it will only consider one read of the pair while here we would like to use the information from both.

Hands-on: Convert the BAM to BED

1. bedtools BAM to BED converter Tool: toolshed.g2.bx.psu.edu/repos/iuc/bedtools/bedtools_bamtobed/2.30.0 with the following parameters:
   • param-file “Convert the following BAM file to BED”: Select the output of MarkDuplicates tool

We call peaks with MACS2. In order to get the coverage centered on the 5’ extended 100bp each side we will use --shift -100 and --extend 200:

Hands-on: Call peaks with MACS2

1. MACS2 callpeak Tool: toolshed.g2.bx.psu.edu/repos/iuc/macs2/macs2_callpeak/2.1.1.20160309.6 with the following parameters:
   • “Are you pooling Treatment Files?”: No
     • param-file Select the output of bedtools BAM to BED converter tool
   • “Do you have a Control File?”: No
   • “Format of Input Files”: Single-end BED
   • “Effective genome size”: H. sapiens (2.7e9)
   • “Build Model”: Do not build the shifting model (--nomodel)
     • “Set extension size”: 200
     • “Set shift size”: -100. It needs to be - half the extension size to be centered on the 5’.
   • “Additional Outputs”:
     • Check Peaks as tabular file (compatible with MultiQC)
     • Check Peak summits
     • Check Scores in bedGraph files
   • In “Advanced Options”:
     • “Composite broad regions”: No broad regions
       • “Use a more sophisticated signal processing approach to find subpeak summits in each enriched peak region”: Yes
     • “How many duplicate tags at the exact same location are allowed?”: all

   Comment: Why keeping all duplicates is important

   We previously removed duplicates using MarkDuplicates tool using paired-end information. If two pairs had identical R1 but different R2, we knew it was not a PCR duplicate. Because we converted the BAM to BED we lost the pair information. If we keep the default (removing duplicates) one of the 2 identical R1 would be filtered out as duplicate.

Visualisation of Coverage

Prepare the Datasets

Extract CTCF peaks on chr22 in intergenic regions

As our training dataset is focused on chromosome 22 we will only use the CTCF peaks from chr22. We expect to have ATAC-seq coverage at TSS but only good ATAC-seq have coverage on intergenic CTCF. Indeed, the CTCF protein is able to position nucleosomes and creates a region depleted of nucleosome of around 120bp Fu et al. 2008. This is smaller than the 200bp nucleosome-free region around TSS and also probably not present in all cells. Thus it is more difficult to get enrichment. In order to get the list of intergenic CTCF peaks of chr22, we will first select the peaks on chr22 and then exclude the one which overlap with genes.

Hands-on: Select CTCF peaks from chr22 in intergenic regions:

1. Filter data on any column using simple expressions Tool: Filter1 with the following parameters:
   • param-file “Filter”: Select the first dataset: ENCFF933NTR.bed.gz
   • “With following condition”: c1=='chr22'
2. bedtools Intersect intervals find overlapping intervals in various ways Tool: toolshed.g2.bx.psu.edu/repos/iuc/bedtools/bedtools_intersectbed/2.30.0 with the following parameters:
   • param-file “File A to intersect with B”: Select the output of Filter data on any column using simple expressions tool
   • “Combined or separate output files”: One output file per 'input B' file
     • param-file “File B to intersect with A”: Select the dataset chr22 genes
   • “What should be written to the output file?”: Write the original entry in A for each overlap (-wa)
   • “Required overlap”: Default: 1bp
   • “Report only those alignments that do not overlap with file(s) B”: Yes
3. Rename the datasets intergenic CTCF peaks chr22.

Convert bedgraph from MACS2 to bigwig

The bedgraph format is easily readable for human but it can be very large and visualising a specific region is quite slow. We will change it to bigwig format which is a binary format, so we can visualise any region of the genome very quickly.

Tip: Speed-up the bedgraph to bigwig conversion

In this tutorial we focus on chr22, thus we could restrict our bedgraph to chr22 before doing the conversion. This will both speed-up the conversion and decrease the amount of memory needed.

To do so, we will use Filter data on any column using simple expressions Tool: Filter1 with the following parameters:

• param-file “Filter”: Select the output of MACS2 tool (Bedgraph Treatment).
• “With following condition”: c1=='chr22'

This will decrease by half the size of the file. In the next step, choose the output of Filter instead of the output of MACS2.

Hands-on: Convert bedgraphs to bigWig.

1. Wig/BedGraph-to-bigWig Tool: wig_to_bigWig with the following parameters:
   • param-file “Convert”: Select the output of MACS2 tool (Bedgraph Treatment).
   • “Converter settings to use”: Default
2. Rename the datasets MACS2 bigwig.

Create heatmap of coverage at TSS with deepTools

You might be interested in checking the coverage on specific regions. For this, you can compute a heatmap. We will use the deepTools plotHeatmap. As an example, we will here make a heatmap centered on the transcription start sites (TSS) and another one centered on intergenic CTCF peaks. First, on the TSS:

Generate computeMatrix

The input of plotHeatmap is a matrix in a hdf5 format. To generate it we use the tool computeMatrix that will evaluate the coverage at each locus we are interested in.

Hands-on: Generate the matrix

1. computeMatrix Tool: toolshed.g2.bx.psu.edu/repos/bgruening/deeptools_compute_matrix/deeptools_compute_matrix/3.3.2.0.0 with the following parameters:
   • In “Select regions”:
     • param-repeat “Insert Select regions”
       • param-file “Regions to plot”: Select the dataset chr22 genes
   • “Sample order matters”: No
     • param-file “Score file”: Select the output of Wig/BedGraph-to-bigWig tool that should be named MACS2 bigwig.
   • “computeMatrix has two main output options”: reference-point
   • “The reference point for the plotting”: beginning of region (e.g. TSS)
   • “Show advanced output settings”: no
   • “Show advanced options”: yes
     • “Convert missing values to 0?”: Yes

Plot with plotHeatmap

We will now generate a heatmap. Each line will be a transcript. The coverage will be summarized with a color code from red (no coverage) to blue (maximum coverage). All TSS will be aligned in the middle of the figure and only the 2 kb around the TSS will be displayed. Another plot, on top of the heatmap, will show the mean signal at the TSS. There will be one heatmap per bigwig.

Hands-on: Generate the heatmap

1. plotHeatmap Tool: toolshed.g2.bx.psu.edu/repos/bgruening/deeptools_plot_heatmap/deeptools_plot_heatmap/3.3.2.0.1 with the following parameters:
   • param-file “Matrix file from the computeMatrix tool”: Select the output of computeMatrix tool.
   • “Show advanced output settings”: no
   • “Show advanced options”: no

Comment: plotHeatmap Results

This is what you get from plotHeatmap:

Question

1. Is the coverage symmetric?
2. What is the mean value in genes?

Solution

1. No, it is higher on the left which is expected as usually the promoter of active genes is accessible.
2. Around 5.5.

Now we will repeat the procedure for CTCF peaks of chr22 in intergenic regions:

Hands-on: Generate the matrix

1. computeMatrix Tool: toolshed.g2.bx.psu.edu/repos/bgruening/deeptools_compute_matrix/deeptools_compute_matrix/3.3.2.0.0 with the following parameters:
   • In “Select regions”:
     • param-repeat “Insert Select regions”
       • param-file “Regions to plot”: Select the dataset intergenic CTCF peaks chr22
   • “Sample order matters”: No
     • param-file “Score file”: Select the output of Wig/BedGraph-to-bigWig tool that should be named MACS2 bigwig.
   • “Would you like custom sample labels?”: No, use sample names in the history
   • “computeMatrix has two main output options”: reference-point
     • “The reference point for the plotting”: center of region
   • “Show advanced output settings”: no
   • “Show advanced options”: yes
     • “Convert missing values to 0?”: Yes
2. plotHeatmap Tool: toolshed.g2.bx.psu.edu/repos/bgruening/deeptools_plot_heatmap/deeptools_plot_heatmap/3.3.2.0.1 with the following parameters:
   • param-file “Matrix file from the computeMatrix tool”: Select the output of computeMatrix tool.
   • “Show advanced output settings”: no
   • “Show advanced options”: yes
     • In “Colormap to use for each sample”:
       • param-repeat “Insert Colormap to use for each sample”
         • “Color map to use for the heatmap”: Blues # Or what you want
     • “The x-axis label”: distance from peak center (bp)
     • “The y-axis label for the top panel”: CTCF peaks
     • “Reference point label”: peak center
     • “Labels for the regions plotted in the heatmap”: CTCF_peaks
     • “Did you compute the matrix with more than one groups of regions?”: Yes, I used multiple groups of regions

Comment: plotHeatmap Results

This is what you get from plotHeatmap, this is much more symetric:

Visualise Regions with pyGenomeTracks

In order to visualise a specific region (e.g. the gene RAC2), we can either use a genome browser like IGV or UCSC browser, or use pyGenomeTracks to make publishable figures. We will use pyGenomeTracks.

Hands-on: Task description

1. pyGenomeTracks Tool: toolshed.g2.bx.psu.edu/repos/iuc/pygenometracks/pygenomeTracks/3.6 with the following parameters:
   • “Region of the genome to limit the operation”: chr22:37,193,000-37,252,000
   • In “Include tracks in your plot”:
     • param-repeat “Insert Include tracks in your plot”
       • “Choose style of the track”: Bigwig track
         • “Plot title”: Coverage from MACS2 (extended +/-100bp)
         • param-file “Track file(s) bigwig format”: Select the output of Wig/BedGraph-to-bigWig tool called MACS2 bigwig.
         • “Color of track”: Select the color of your choice
         • “Minimum value”: 0
         • “height”: 5
         • “Show visualization of data range”: Yes
     • param-repeat “Insert Include tracks in your plot”
       • “Choose style of the track”: NarrowPeak track
         • “Plot title”: Peaks from MACS2 (extended +/-100bp)
         • param-file “Track file(s) encodepeak or bed format”: Select the output of MACS2 tool (narrow Peaks).
         • “Color of track”: Select the color of your choice
         • “display to use”: box: Draw a box
         • “Plot labels (name, p-val, q-val)”: No
     • param-repeat “Insert Include tracks in your plot”
       • “Choose style of the track”: Gene track / Bed track
         • “Plot title”: Genes
         • param-file “Track file(s) bed or gtf format”: chr22 genes
         • “Color of track”: Select the color of your choice
         • “height”: 5
         • “Plot labels”: yes
           • “Put all labels inside the plotted region”: Yes
           • “Allow to put labels in the right margin”: Yes
     • param-repeat “Insert Include tracks in your plot”
       • “Choose style of the track”: NarrowPeak track
         • “Plot title”: CTCF peaks
         • param-file “Track file(s) encodepeak or bed format”: Select the first dataset: ENCFF933NTR.bed.gz
         • “Color of track”: Select the color of your choice
         • “display to use”: box: Draw a box
         • “Plot labels (name, p-val, q-val)”: No
     • param-repeat “Insert Include tracks in your plot”
       • “Choose style of the track”: X-axis
2. Click on the galaxy-eye (eye) icon of the output.

Comment: pyGenomeTracks Results

You should get similar to results to this from pyGenomeTracks:

Question

In the ATAC-Seq sample in this selected region we see four peaks detected by MACS2.

1. How many TSS are accessible in the sample in the displayed region?
2. How many CTCF binding loci are accessible?
3. Can you spot peaks with no TSS and no CTCF peak?

Solution

1. In total, we can see 3 TSS for 6 transcripts for 2 genes. The TSS of RAC2 corresponds to an ATAC-Seq peak whereas there is no significant coverage on both TSS of SSTR3.

2. Only the first peak on the left overlaps with a CTCF binding site.

3. Amongst the 4 peaks in this region, the 2 in the middle do not correspond to CTCF peaks or TSS.

As CTCF creates accessible regions, a region containing a peak with no corresponding CTCF peak or TSS could be a putative enhancer. In the pyGenomeTracks plot we see a region like this located in the intron of a gene and another one between genes. However, it is impossible to guess from the position which would be the gene controlled by this region. And of course, more analyses are needed to assess if it is a real enhancer, for example, histone ChIP-seq, 3D structure, transgenic assay, etc.

Conclusion

In this training you have learned the general principles of ATAC-Seq data analysis. ATAC-Seq is a method to investigate the chromatin accessibility and the genome is treated with a transposase (enzyme) called Tn5. It marks open chromatin regions by cutting and inserting adapters for sequencing. The training material gave you an insight into how to quality control the data. You should look for low quality bases, adapter contamination, correct insert size and PCR duplicates (duplication level). We showed you how to remove adapters and PCR duplicates, if FastQC, shows a warning in these areas. We mapped the reads with Bowtie2, filtered our reads for properly paired, good quality and reads that do not map to the mitochondrial genome. We found open chromatin regions with MACS2, a tool to find regions of genomic enrichment (peaks). We investigated the read coverage around TSS with the help of computeMatrix and plotHeatmap. Last but not least, we visualised the peaks and other informative tracks, such as CTCF binding regions and hg38 genes, with the help of pyGenomeTracks. At the end, we found open chromatin regions that did not overlap with CTCF sites or TSS, which could be potential putative enhancer regions detected by the ATAC-Seq experiment.

Key points

• ATAC-Seq can be used to identify accessible gene promoters and enhancers

• Several filters are applied to the reads, such as removing those mapped to mitochondria

• Fragment distribution can help determine whether an ATAC-Seq experiment has worked well

Frequently Asked Questions

Have questions about this tutorial? Check out the tutorial FAQ page or the FAQ page for the Epigenetics topic to see if your question is listed there. If not, please ask your question on the GTN Gitter Channel or the Galaxy Help Forum

References

1. Fu, Y., M. Sinha, C. L. Peterson, and Z. Weng, 2008 The insulator binding protein CTCF positions 20 nucleosomes around its binding sites across the human genome. PLoS genetics 4: e1000138. 10.1371/journal.pgen.1000138
2. Adey, A., H. G. Morrison, A. (no last name), X. Xun, J. O. Kitzman et al., 2010 Rapid, low-input, low-bias construction of shotgun fragment libraries by high-density in vitro transposition. Genome Biology 11: R119. 10.1186/gb-2010-11-12-r119
3. Green, B., C. Bouchier, C. Fairhead, N. L. Craig, and B. P. Cormack, 2012 Insertion site preference of Mu, Tn5, and Tn7 transposons. Mobile DNA 3: 3. 10.1186/1759-8753-3-3
4. Buenrostro, J. D., P. G. Giresi, L. C. Zaba, H. Y. Chang, and W. J. Greenleaf, 2013 Transposition of native chromatin for fast and sensitive epigenomic profiling of open chromatin, DNA-binding proteins and nucleosome position. Nature Methods 10: 1213–1218. 10.1038/nmeth.2688
5. Corces, M. R., A. E. Trevino, E. G. Hamilton, P. G. Greenside, N. A. Sinnott-Armstrong et al., 2017 An improved ATAC-seq protocol reduces background and enables interrogation of frozen tissues. Nature Methods 14: 959–962. 10.1038/nmeth.4396
6. Kia, A., C. Gloeckner, T. Osothprarop, N. Gormley, E. Bomati et al., 2017 Improved genome sequencing using an engineered transposase. BMC Biotechnology 17: 10.1186/s12896-016-0326-1
7. Litzenburger, U. M., J. D. Buenrostro, B. Wu, Y. Shen, N. C. Sheffield et al., 2017 Single-cell epigenomic variability reveals functional cancer heterogeneity. Genome Biology 18: 10.1186/s13059-016-1133-7

Feedback

Did you use this material as an instructor? Feel free to give us feedback on how it went.
Did you use this material as a learner or student? Click the form below to leave feedback.

Citing this Tutorial

1. Lucille Delisle, Maria Doyle, Florian Heyl, 2022 ATAC-Seq data analysis (Galaxy Training Materials). https://training.galaxyproject.org/training-material/topics/epigenetics/tutorials/atac-seq/tutorial.html Online; accessed TODAY
2. Batut et al., 2018 Community-Driven Data Analysis Training for Biology Cell Systems 10.1016/j.cels.2018.05.012

BibTeX

@misc{epigenetics-atac-seq,
author = "Lucille Delisle and Maria Doyle and Florian Heyl",
title = "ATAC-Seq data analysis (Galaxy Training Materials)",
year = "2022",
month = "10",
day = "18"
url = "\url{https://training.galaxyproject.org/training-material/topics/epigenetics/tutorials/atac-seq/tutorial.html}",
note = "[Online; accessed TODAY]"
}
@article{Batut_2018,
    doi = {10.1016/j.cels.2018.05.012},
    url = {https://doi.org/10.1016%2Fj.cels.2018.05.012},
    year = 2018,
    month = {jun},
    publisher = {Elsevier {BV}},
    volume = {6},
    number = {6},
    pages = {752--758.e1},
    author = {B{\'{e}}r{\'{e}}nice Batut and Saskia Hiltemann and Andrea Bagnacani and Dannon Baker and Vivek Bhardwaj and Clemens Blank and Anthony Bretaudeau and Loraine Brillet-Gu{\'{e}}guen and Martin {\v{C}}ech and John Chilton and Dave Clements and Olivia Doppelt-Azeroual and Anika Erxleben and Mallory Ann Freeberg and Simon Gladman and Youri Hoogstrate and Hans-Rudolf Hotz and Torsten Houwaart and Pratik Jagtap and Delphine Larivi{\`{e}}re and Gildas Le Corguill{\'{e}} and Thomas Manke and Fabien Mareuil and Fidel Ram{\'{\i}}rez and Devon Ryan and Florian Christoph Sigloch and Nicola Soranzo and Joachim Wolff and Pavankumar Videm and Markus Wolfien and Aisanjiang Wubuli and Dilmurat Yusuf and James Taylor and Rolf Backofen and Anton Nekrutenko and Björn Grüning},
    title = {Community-Driven Data Analysis Training for Biology},
    journal = {Cell Systems}
}
                   

Congratulations on successfully completing this tutorial!