ATAC-Seq data analysis

Overview

question Questions
  • Which DNA regions are accessible in the human lymphoblastoid cell line GM12878?

  • How to analyse and visualise ATAC-Seq data?

objectives Objectives
  • Apply appropriate analysis and quality control steps for ATAC-Seq

  • Generate a heatmap of transcription start site accessibility

  • Visualise peaks for specific regions

requirements Requirements

time Time estimation: 3 hours

Supporting Materials

last_modification Last modification:

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.

ATAC-Seq
Figure 1: Buenrostro et al. 2013 Nat Methods

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. Preprocessing
    1. Get Data
    2. Quality Control
    3. Trimming Reads
  2. Mapping
    1. Mapping Reads to Reference Genome
  3. Filtering Mapped Reads
    1. Filter Uninformative Reads
    2. Filter Duplicate Reads
    3. Check Insert Sizes
  4. Peak calling
    1. Call Peaks
  5. Visualisation of Coverage
    1. Prepare the Datasets
    2. Create heatmap of genes
    3. Visualise Regions with pyGenomeTracks
  6. Conclusion

comment Note: 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 Hands-on: Data upload

  1. Create a new history for this tutorial

    tip 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 and ENCODE or from the shared data library

    https://zenodo.org/record/3270536/files/SRR891268_R1.fastq.gz
    https://zenodo.org/record/3270536/files/SRR891268_R2.fastq.gz
    https://www.encodeproject.org/files/ENCFF933NTR/@@download/ENCFF933NTR.bed.gz
    
    • 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

    By default, Galaxy uses the URL as the name, so rename the files with a more useful name.

    tip 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

    • Find the correct folder (ask 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 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 BED file is bed. If they are not then change the datatype as described below.

    tip 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
    • Click the Change datatype button

comment FASTQ format

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

comment BED format

If you are not familiar with BED 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 Hands-on: Obtain Annotation for hg38 genes

  1. UCSC Main tool with the following parameters:
    • “clade”: Mammal
    • “genome”: Human
    • “assembly”: Dec. 2013 (GRCh38/hg38)
    • “group”: Genes and Gene Prediction
    • “track”: All GENCODE V31
    • “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:

  1. Cut columns from a table tool with the following parameters:
    • param-text “Cut columns”: c3,c5,c6,c13,c12,c4
    • param-text “Delimited by”: Tab
    • param-file “From”: UCSC Main on Human: wgEncodeGencodeBasicV31 (chr22:1-50,818,468)
  2. Rename the dataset as chr22 genes

    tip 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
  3. Change its datatype to BED

    tip 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
    • Click the Change datatype button

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 Hands-on: Task description

  1. FastQC tool with the default 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 Questions

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

    solution 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

    Nextera adapter sequences are observable in the Adapter Content section.

comment FastQC Results

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

FastQC screenshot of the Adapter Content section
Figure 2: FastQC screenshot on 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.

Nextera library with the sequence of adapters
Figure 3: Nextera library with the sequence of adapters

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 Hands-on: Task description

  1. Cutadapt tool 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
              • param-text “Enter custom 3’ adapter name (Optional if Multiple output is ‘No’)”: Nextera R1
              • param-text “Enter custom 3’ adapter sequence”: CTGTCTCTTATACACATCTCCGAGCCCACGAGAC
      • In “Read 2 Options”:
        • In “3’ (End) Adapters”:
          • param-repeat “Insert 3’ (End) Adapters”
            • “Source”: Enter custom sequence
              • param-text “Enter custom 3’ adapter name (Optional)”: Nextera R2
              • param-text “Enter custom 3’ adapter sequence”: CTGTCTCTTATACACATCTGACGCTGCCGACGA
    • In “Filter Options”:
      • param-text “Minimum length”: 20
    • In “Read Modification Options”:
      • param-text “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:

Summary of cutadapt
Figure 4: Summary of cutadapt

question Questions

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

solution Solution

  1. ~14%
  2. ~99%

hands_on Hands-on: Check Adapter Removal with FastQC

  1. FastQC tool with the default 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

If we run FastQC again we should see under Adapter Content that the Nextera adapters are no longer present.

FastQC screenshot on the adapter content section after cutadapt
Figure 5: FastQC screenshot on the adapter content section after cutadapt

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.

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 Hands-on: Mapping reads to reference genome

  1. Bowtie2 tool 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 Dec. 2013 (GRCh38/hg38 (hg38)
    • “Select analysis mode”: 1: Default setting only
      • “Do you want to use presets?”: Very sensitive end-to-end (--very-sensitive)
    • “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:

Mapping statistics of bowtie2
Figure 6: Mapping statistics of bowtie2

question Questions

What percentage of read pairs mapped concordantly?

solution Solution

54.07+43.63=97.7%

comment 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 Hands-on: Filtering of uninformative reads

  1. Filter BAM datasets on a variety of attributes tool with the following parameters:
    • param-file “BAM dataset(s) to filter”: Select the output of Bowtie2 tool “alignments”
    • In “Condition”:
      • “1: Condition”
        • In “Filter”:
          • “1: 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 Questions

  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 Solution

  1. The original BAM file is 28 MB, the filtered one is 14.8 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.

tip 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 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 sorting):

Samtools idxstats result
Figure 7: Samtools idxstats result

There are 221 000 reads which map to chrM and 170 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 Hands-on: Remove duplicates

  1. MarkDuplicates tool 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 Comment: Default of MarkDuplicates tool

    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:

Metrics of MarkDuplicates
Figure 8: Metrics of MarkDuplicates

tip Tip: Formatting the MarkDuplicate metrics for readability

  1. Select tool 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 they are not then change the datatype as described above.
  3. Transpose tool:
    • param-file “Select lines from”: Select the output of Select tool
Metrics of MarkDuplicates
Figure 9: Metrics of MarkDuplicates

question Questions

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

solution Solution

  1. 133284
  2. 3549

Check Insert Sizes

We will check the insert sizes with Picard CollectInsertSizeMetrics. 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 Hands-on: Plot the distribution of fragment sizes.

  1. CollectInsertSizeMetrics tool with the following parameters:
    • param-file “Select SAM/BAM dataset or dataset collection”: Select the output of MarkDuplicates tool “BAM output”
    • “Load reference genome from”: Local cache
      • “Using reference genome”: Human Dec. 2013 (GRCh38/hg38) (hg38)
  2. Click on the galaxy-eye (eye) icon of the upper one of the 2 outputs (the pdf file).

comment CollectInsertSizeMetrics Results

This is what you get from CollectInsertSizeMetrics:

Fragment size distribution
Figure 10: Fragment size distribution

question Questions

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

solution 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:

Fragment size distribution of a purified DNA
Figure 11: Fragment size distribution of a purified DNA

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

Fragment size distribution of a failed ATAC-Seq
Figure 12: Fragment size distribution of a failed ATAC-Seq
Fragment size distribution of another failed ATAC-Seq
Figure 13: Fragment size distribution of another very noisy ATAC-Seq

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

Fragment size distribution of a good ATAC-Seq
Figure 14: Fragment size distribution of a good ATAC-Seq

comment Comment on FR and RF

FR stands for forward reverse orientation of the read pairs, meaning, your reads are oriented as -> <- so the first read is on the forward and the second on the reverse strand. RF stands for reverse forward oriented, i.e., <- ->. It really depends on your experiment, how your reads are oriented and if the orientation plays a role.

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. We will use Genrich. It is very important at this point that we center the reads on the 5’ extremity (read start site) as this is where Tn5 cut. You want your peaks around the nucleosomes and not directly on the nucleosome:

Scheme of ATAC-Seq reads relative to nucleosomes
Figure 15: Scheme of ATAC-Seq reads relative to nucleosomes

comment Comment on Tn5 insertion

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

Nextera Library Construction
Figure 16: Nextera Library Construction

This means that to have the read start site reflect 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.

If we only assess the coverage of the start sites 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 100bp (50 bp in each direction) to assess coverage.

hands_on Hands-on: Identifying enriched genomic regions

  1. Genrich tool with the following parameters:
    • “Are you pooling Treatment Files?”: No
    • param-file “Treatment File(s)”: Select the output of MarkDuplicates tool
    • “Do you have a Control File?”: No
    • “Filter Options”:
      • “Remove PCR duplicates”: Yes
    • “ATAC Options”:
      • “Use ATAC-seq mode.”: Yes
      • “Expand cut sites.”: 100
    • “Output Options”:
      • Bedgraph-ish Pileups”: Yes

Visualisation of Coverage

Prepare the Datasets

Thanks to Genrich we now have a coverage file which represents the coverage of the read start sites extended 50 bp to each side. The output of Genrich is a BedGraph-ish pileup (6 columns text format with a comment line and a header). We will first need to convert it to a bedgraph format (4 columns text format with no header) to be able to visualise it. 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.

Convert BedGraph-ish pileup to bigWig

hands_on Hands-on: Convert bedgraph-ish pileup to bigWig.

  1. Text reformatting with awk tool with the following parameters:
    • param-file “File to process”: Select the output of Genrich tool Bedgraph Pileups”.
    • “AWK Program”: NR>=3 {print $1,$2,$3,$4}

    comment Comment: From BedGraph-ish pileup to bedgraph

    The awk program will read each line of the output of Genrich tool, when the number of the line is greater or equal to 3 (NR>=3), it will write the first 4 columns (print $1,$2,$3,$4) into a new file.

  2. Wig/BedGraph-to-bigWig tool with the following parameters:
    • param-file “Convert”: Select the output of Text reformatting with awk tool Bedgraph”.
    • “Converter settings to use”: Default

Sort CTCF Peaks

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. The pyGenomeTracks tool needs all BED files sorted, thus we sort the CTCF peaks.

hands_on Hands-on: Sort the BED files

  1. bedtools SortBED order the intervals tool with the following parameters:
    • param-file “Sort the following BED/bedGraph/GFF/VCF file”: ENCFF933NTR.bed.gz

Convert the Genrich peaks to BED

At the moment, pyGenomeTracks does not deal with the datatype encodepeak which is a special bed. So we need to change the datatype of the output of Genrich from encodepeak to bed.

tip 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
  • Click the Change datatype button

Create heatmap of genes

You might also be interested in 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).

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 Hands-on: Generate the matrix

  1. computeMatrix tool with the following parameters:
    • In “Select regions”:
        1. “Select regions”
          • param-file “Regions to plot”: Select the dataset hg38_Gencode_V28_chr22_geneName.bed
    • “Sample order matters”: No
      • param-file “Score file”: Select the output of Wig/BedGraph-to-bigWig tool.
    • “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
      • “Labels for the samples (each bigwig)”: ATAC-Seq

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.

hands_on Hands-on: Generate the heatmap

  1. plotHeatmap tool 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:

plotHeatmap output
Figure 17: plotHeatmap output

question Questions

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

solution Solution

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

Visualise Regions with pyGenomeTracks

hands_on Hands-on: Task description

  1. pyGenomeTracks tool 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”:
      • “1. Include tracks in your plot”
        • “Choose style of the track”: Bigwig track
          • “Plot title”: Coverage from Genrich (extended +/-50bp)
          • param-file “Track file bigwig format”: Select the output of Wig/BedGraph-to-bigWig tool.
          • “Color of track”: Select the color of your choice
          • “height”: 5
          • “Show visualization of data range”: Yes
          • “Include spacer at the end of the track”: 0.5
      • param-repeat “Insert Include tracks in your plot” - “Plot title”: Peaks from Genrich (extended +/-50bp) - param-file “Track file bed format”: Select the output of Genrich tool (the one you converted from encodepeak to bed). - “Color of track”: Select the color of your choice - “height”: 3 - “Plot labels”: No - “Include spacer at the end of the track”: 0.5
      • param-repeat “Insert Include tracks in your plot”
        • “Choose style of the track”: Gene track / Bed track
          • “Plot title”: Genes
          • param-file “Track file bed format”: chr22 genes
          • “Color of track”: Select the color of your choice
          • “height”: 5
          • “Include spacer at the end of the track”: 0.5
      • param-repeat “Insert Include tracks in your plot”
        • “Choose style of the track”: Gene track / Bed track
          • “Plot title”: CTCF peaks
          • param-file “Track file bed format”: Select the dataset bedtools SortBED of ENCFF933NTR.bed.gz
          • “Color of track”: Select the color of your choice
          • “Plot labels”: No
          • “Include spacer at the end of the track”: 0.5
    • “Configure x-axis”: Yes
      • “Where to place the x-axis”: Bottom
  2. Click on the galaxy-eye (eye) icon of the output.

comment pyGenomeTracks Results

You should get similar to results to this from pyGenomeTracks:

pyGenomeTracks output
Figure 18: pyGenomeTracks output

Unfortunately, Genrich does not work very well with our small training dataset (every covered region is called a peak). This is because most of the data is on chr22 whereas the background model was built on the whole genome. When the pipeline described here was run on 20 million of pairs from the original dataset, this is the output of pyGenomeTracks:

pyGenomeTracks output for 20 million of pairs on the whole genome
Figure 19: pyGenomeTracks output for 20 million of pairs on the whole genome

question Questions

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

  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 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 binds so ubiquitously and by itself can displace the nucleosome creating 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 Genrich, 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.

ATAC workflow
Figure 20: ATAC workflow

keypoints 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

References

  1. 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
  2. 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
  3. 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
  4. 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
  5. 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
  6. 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

Citing this Tutorial

  1. Lucille Delisle, Maria Doyle, Florian Heyl, ATAC-Seq data analysis (Galaxy Training Materials). /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

details BibTeX

@misc{epigenetics-atac-seq,
    author = "Lucille Delisle and Maria Doyle and Florian Heyl",
    title = "ATAC-Seq data analysis (Galaxy Training Materials)",
    year = "",
    month = "",
    day = ""
    url = "\url{/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 Congratulations on successfully completing this tutorial!



Did you use this material as an instructor? Feel free to give us feedback on how it went.