EBAII

Supports de cours de l'Ecole de Bioinformatique Aviesan - IFB - Inserm "Initiation au traitement des données de génomique obtenues par séquençage à haut débit"

View the Project on GitHub IFB-ElixirFr/EBAII

ChIP-seq Hands-on Roscoff 2021

  1. Introduction
  2. Downloading ChIP-seq reads from NCBI
  3. Connect to the server and set up your environment
  4. Quality control of the reads and statistics
  5. Mapping the reads with Bowtie
  6. Estimating the number of duplicated reads
  7. ChIP quality controls
  8. Visualizing the data in a genome browser
  9. Peak calling with MACS
  10. Motif analysis
  11. Peak annotation
  12. Bonus: Peak annotation using R
  13. FAQ
  14. References

Introduction

Goal

The aim is to :

Summary

This training gives an introduction to ChIP-seq data analysis, covering the processing steps starting from the reads to the peaks. Among all possible downstream analyses, the practical aspect will focus on motif analyses. A particular emphasis will be put on deciding which downstream analyses to perform depending on the biological question. This training does not cover all methods available today. It does not aim at bringing users to a professional NGS analyst level but provides enough information to allow biologists understand what is DNA sequencing in practice and to communicate with NGS experts for more in-depth needs.

Dataset description

For this training, we will use two datasets:

Downloading ChIP-seq reads from NCBI

Goal: Identify the datasets corresponding to the studied article and retrieve the data (reads as FASTQ files) corresponding to 2 replicates of a condition and the corresponding control.

1 - Obtaining an identifier for a chosen dataset

NGS datasets are (usually) made freely accessible for other scientists, by depositing these datasets into specialized databanks. Sequence Read Archive (SRA) located in USA hosted by NCBI, and its European equivalent European Nucleotide Archive (ENA) located in England hosted by EBI both contains raw reads.

Functional genomic datasets (transcriptomics, genome-wide binding such as ChIP-seq,…) are deposited in the databases Gene Expression Omnibus (GEO) or its European equivalent ArrayExpress.

Within an article of interest, search for a sentence mentioning the deposition of the data in a database. Here, the following sentence can be found at the end of the Materials and Methods section: “All genome-wide data from this publication have been deposited in NCBI’s Gene Expression Omnibus (GSE41195).” We will thus use the GSE41195 identifier to retrieve the dataset from the NCBI GEO (Gene Expression Omnibus) database.

2 - Accessing GSE41195 from GEO

  1. The GEO database hosts processed data files and many details related to the experiments. SRA (Sequence Read Archive) stores the actual raw sequence data.
  2. Search in Google GSE41195. Click on the first link to directly access the correct page on the GEO database. geo
  3. This GEO entry is a mixture of expression analysis (Nimblegen Gene Expression Array), chip-chip and chip-seq. At the bottom of the page, click on the subseries related to the chip-seq datasets. (this subseries has its own identifier: GSE41187). geo2
  4. From this page, we will focus on the experiment FNR IP ChIP-seq Anaerobic A. At the bottom of the page, click on the link “GSM1010219 - FNR IP ChIP-seq Anaerobic A”.
  5. In the new page, go to the bottom to find the SRA identifier. This is the identifier of the raw dataset stored in the SRA database.
    geo3
  6. Copy the identifier SRX189773 (do not click on the link that would take you to the SRA database, see below why)

3 - Downloading FASTQ file from the ENA database

Although direct access to the SRA database at the NCBI is doable, SRA does not store sequences in a FASTQ format. So, in practice, it’s simpler (and quicker!!) to download datasets from the ENA database (European Nucleotide Archive) hosted by EBI (European Bioinformatics Institute) in UK. ENA encompasses the data from SRA.

  1. Go to the EBI website. Paste your SRA identifier (SRX189773) and click on the button “search”. ebi4
  2. Click on the first result. On the next page, there is a link to the FASTQ file. For efficiency, this file has already been downloaded and is available in the “data” folder (SRR576933.fastq.gz)
    ebi5

tip: To download the replicate and control datasets, we should redo the same steps starting from the GEO web page specific to the chip-seq datasets (see step 2.4) and choose FNR IP ChIP-seq Anaerobic B and anaerobic INPUT DNA. Downloaded FASTQ files are available in the data folder (SRR576934.fastq.gz and SRR576938.fastq.gz respectively)

At this point, you have three FASTQ files, two IPs, one control (INPUT).

Connect to the server and set up your environment

1 - Sign in to Jupyterhub and open a Terminal

  1. Go to Jupyterhub
  2. select a job profile: EBAII and click on “Start” jupebai
  3. In the launcher, click on “Terminal” in “Other” section. You should be in your home directory by default. Check it:
    pwd
    

2 - Set up your working environment

  1. Go to your project directory
    cd /shared/projects/<your_project>
    
  2. Create a directory that will contain all results of the upcoming analyses.
    mkdir EBAII2021_chipseq
    
  3. Go to the newly created directory
    cd EBAII2021_chipseq
    
  4. Copy the directory containing data
cp -r /shared/home/slegras/EBAII2021_chipseq/data .
  1. Your directory structure should be like this
    /shared/projects/<your_project>/EBAII2021_chipseq
    │
    └───data
    

If you wish, you can check your directory structure:

 tree

Quality control of the reads and statistics

Goal: Get some basic information on the data (read length, number of reads, global quality of datasets)

1 - Generating the FASTQC report

Before you analyze the data, it is crucial to check the quality of the data. We will use the standard tool for checking the quality of data generated on the Illumina platform: FASTQC.

  1. Create a directory named 01-QualityControl in which to output results from fastqc
    mkdir 01-QualityControl
    
  2. Go to the directory you’ve just created
    cd 01-QualityControl
    

Your directory structure should be like this

/shared/projects/<your_project>/EBAII2021_chipseq
│
└───data
│   
└───01-QualityControl <- you should be here
  1. Get FastQC available in your environment
    module add fastqc/0.11.9
    
  2. Check the help page of the program to see its usage and parameters.
fastqc --help
  1. Launch the FASTQC program on the experiment file (SRR576933.fastq.gz)
    • -o: creates all output files in the specified output directory. ‘.’ means current directory.
      fastqc ../data/SRR576933.fastq.gz -o .
      
  2. Wait until the analysis is finished. Check the FastQC result files.
    ls
    

    SRR576933_fastqc.html SRR576933_fastqc.zip

  3. Go to the directory /shared/projects//EBAII2021_chipseq/1-QualityControl in the tree diretory on the left of the jupyterhub panel and double click on SRR576933_fastqc.html to visualize the file. ![fastqc](/EBAII/2021/ebaiin1/chip-seq/images/fastqc.png)

  4. Launch the FASTQC program on the replicate (SRR576934.fastq.gz) and on the control file (SRR576938.fastq.gz)

Analyze the result of the FASTQC program:
How many reads are present in each file ?
What is the read length ?
Is the overall quality good for the three samples ?
Are there any concerns raised by the report ? If so, can you tell where the problem might come from ?

  1. Once you are done with FastQC, unload it
    module rm fastqc/0.11.9
    

Mapping the reads with Bowtie

Goal: Obtain the coordinates of each read to the reference genome.

1 - Choosing a mapping program

There are multiple programs to perform the mapping step. For reads produced by an Illumina machine for ChIP-seq, the currently “standard” programs are BWA (Li and Durbin, 2009) and Bowtie (versions 1 and 2)(Langmead et al., 2009; Langmead and Salzberg, 2012), and STAR (Dobin et al., 2013) is getting popular. We will use Bowtie version 1.2.1.1 for this exercise, as this program remains effective for short reads (< 50bp).

2 - Bowtie

  1. Load Bowtie
    module add bowtie/1.2.3
    
  2. Try out bowtie
    bowtie
    

    This prints the help of the program. However, this is a bit difficult to read ! If you need to know more about the program, it’s easier to directly check out the manual on the website.

  3. bowtie needs the reference genome to align each read on it. This genome needs to be in a specific format (=index) for bowtie to be able to use it. Several pre-built indexes are available for download on bowtie webpage, but our genome is not there. You will need to make this index file.

  4. Create a directory named 02-Mapping in which to output mapping results
    cd ..
    mkdir 02-Mapping
    
  5. Go to the directory you’ve just created
    cd 02-Mapping
    

3 - Prepare the index file

  1. To make the index file, you will need the complete genome, in FASTA format. It has already been downloaded to gain time (Escherichia_coli_K12.fasta in the course folder) (The genome was downloaded from the NCBI).

  2. Create a directory named index in which to output bowtie indexes
    mkdir index
    
  3. Go to the newly created directory
    cd index
    
  4. Try out bowtie-build
    bowtie-build
    
  5. Build the index for bowtie
    ## Creating genome index : provide the path to the genome file and the name to give to the index (Escherichia_coli_K12)
    bowtie-build ../../data/Escherichia_coli_K12.fasta Escherichia_coli_K12
    
  6. Go back to upper directory i.e 02-Mapping
    cd ..
    

4 - Mapping the IP samples

  1. Create a directory named IP to put mapping results for IP
    mkdir bam
    
  2. Go to the newly created directory bam
    cd bam
    

    Your directory structure should be like this:

    /shared/projects/<your_project>/EBAII2021_chipseq
    │
    └───data
    │   
    └───01-QualityControl
    │   
    └───02-Mapping
    |    └───index
    |    └───bam <- you should be here
    
  3. Let’s see the parameters of bowtie before launching the mapping:
    • Escherichia_coli_K12 is the name of our genome index file
    • Number of mismatches for SOAP-like alignment policy (-v): to 2, which will allow two mismatches anywhere in the read, when aligning the read to the genome sequence.
    • Suppress all alignments for a read if more than n reportable alignments exist (-m): to 1, which will exclude the reads that do not map uniquely to the genome.
    • -q indicates the input file is in FASTQ format. SRR576933.fastq is the name of our FASTQ file.
    • -3 will trim x base from the end of the read. As our last position is of low quality, we’ll trim 1 base.
    • -S will output the result in SAM format
    • 2> SRR576933.out will output some statistics about the mapping in the file SRR576933.out
      ## Run alignment
      ## Tip: first type bowtie command line then add quotes around and prefix it with "sbatch --cpus 10 --wrap="
      sbatch --cpus-per-task 10 --wrap="bowtie -p 10 ../../index/Escherichia_coli_K12 ../../../data/SRR576933.fastq.gz -v 2 -m 1 -3 1 -S 2> SRR576933.out > SRR576933.sam"
      

      This should take few minutes as we work with a small genome. For the human genome, we would need either more time and more resources.

Bowtie output is a SAM file. The SAM format corresponds to large text files, that can be compressed (“zipped”) into a BAM format. The BAM files takes up to 4 time less disk space and are usually sorted and indexed for fast access to the data it contains. The index of a given .bam file is named .bam.bai or .bai file. Some tools require to have the index of the bam file to process it.

  1. Sort the sam file and create a bam file using samtools
    • -b: output BAM
      ## First load samtools
      module add samtools/1.10
      ## Then run samtools
      samtools sort SRR576933.sam | samtools view -b > SRR576933.bam
      
  2. Create an index for the bam file
    samtools index SRR576933.bam
    
  3. Compress the .sam file (you could also delete the file)
    gzip SRR576933.sam
    

Analyze the result of the mapped reads:
Open the file SRR576933.out (for example using the ` less ` command), which contains some statistics about the mapping. How many reads were mapped? How many multi-mapped reads were originally present in the sample? To quit less press ‘q’

  1. Once it’s done, unload the tools you used
    module rm samtools/1.10 bowtie/1.2.3
    

5 - Map the second replicate and the control

  1. Repeat the steps above (in 4 - Mapping the experiment) for the files SRR576934.fastq.gz and SRR576938.fastq.gz in the directory named “bam” within the directory 02-Mapping.

Analyze the result of the mapped reads:
Open the file SRR576938.out. How many reads were mapped?

Estimating the number of duplicated reads

Goal: Duplicated reads i.e reads mapped at the same positions in the genome are present in ChIP-seq results. They can arise for several reasons including a biased amplification during the PCR step of the library prep, DNA fragment coming from repetitive elements of the genome, sequencing saturation or the same clusters read several times on the flowcell (i.e optical duplicates). As analyzing ChIP-Seq data consist at some point in detecting signal enrichment, we can not keep duplicated reads for subsequent analysis. So let’s detect them using Picard

  1. Go to the directory with alignment files
    cd /shared/projects/<your_project>/EBAII2021_chipseq/02-Mapping/bam
    
  2. Run Picard markDuplicates to mark duplicated reads (= reads mapping at the exact same location on the genome)
## Load picard
module add picard/2.22.0

## Run picard
picard MarkDuplicates \

CREATE_INDEX=true \
INPUT=SRR576933.bam \
OUTPUT=Marked_SRR576933.bam \
METRICS_FILE=metrics \
VALIDATION_STRINGENCY=STRICT

picard MarkDuplicates \
CREATE_INDEX=true \
INPUT=SRR576938.bam \
OUTPUT=Marked_SRR576938.bam \
METRICS_FILE=metrics \
VALIDATION_STRINGENCY=STRICT

To determine the number of duplicated reads marked by Picard, we can run the samtools flagstat command:

## Add samtools to your environment
module add samtools/1.10
## run samtools
samtools flagstat Marked_SRR576933.bam

Go back to working home directory (i.e /shared/projects//EBAII2021_chipseq/)

## Unload picard and samtools
module rm samtools/1.10 picard/2.22.0
## If you are in 02-Mapping/bam
cd ../..

ChIP quality controls

Goal: The first exercise aims at plotting the Lorenz curve to assess the quality of the chIP. The second exercise aims at calculating the NSC and RSC ENCODE quality metrics. These metrics allow to classify the datasets (after mapping, contrary to FASTQC that works on raw reads) in regards to the NSC and RSC values observed in the ENCODE datasets (see ENCODE guidelines)

1 - Plot the Lorenz curve with Deeptools

  1. Create a directory named 03-ChIPQualityControls in which to put mapping results for IP
    mkdir 03-ChIPQualityControls
    
  2. Go to the newly created directory
    cd 03-ChIPQualityControls
    
  3. Run Deeptools plotFingerprint (Fidel et al, NAR, 2016) to draw the Lorenz curve

Look at the result files fingerprint.png. What do you think of it?

Go back to working home directory (i.e /shared/projects/training/<login>/EBAII2021_chipseq)

## Unload deepTools
module rm deeptools/3.2.0
## If you are in 03-ChIPQualityControls
cd ..

Visualizing the data in a genome browser

Goal: View the data in their genomic context, to check whether the IP worked

1 - Choosing a genome browser

There are several options for genome browsers, divided between the local browsers (need to install the program, eg. IGV) and the online web browsers (eg. UCSC genome browser, Ensembl). We often use both types, depending on the aim and the localization of the data. If the data are on your computer, to prevent data transfer, it’s easier to visualize the data locally (IGV). Note that if you’re working on a non-model organism, the local viewer will be the only choice. If the aim is to share the results with your collaborators, view many tracks in the context of many existing annotations, then the online genome browsers are more suitable.

2 - Viewing the raw alignment data in IGV

  1. Download the following files from the server onto your computer

Browse around in the genome. Do you see peaks?
Browse into IGV. Go to the following genes: pepT (geneID:b1127), ycfP (geneID:b1108)

However, looking at BAM file as such does not allow to directly compare the two samples as data are not normalized. Let’s generate normalized data for visualization.

3 - Viewing scaled data

bamCoverage from deepTools generates BigWigs out of BAM files

  1. Try it out
    ## Load deeptools in your environment
    module add deeptools/3.2.0
    ## run bamCoverage
    bamCoverage --help
    
  2. Create a directory named 04-Visualization to store bamCoverage outputs
    mkdir 04-Visualization
    
  3. Go to the newly created directory
    cd 04-Visualization
    

Your directory structure should be like this:

/shared/projects/<your_project>/EBAII2021_chipseq
│
└───data
│   
└───01-QualityControl
│   
└───02-Mapping
|    └───index
|    └───bam
│   
└───03-ChIPQualityControls
│   
└───04-Visualization <- you should be in this folder
  1. Generate a scaled bigwig file on the IP with bamCoverage
  1. Do it for the replicate and the control.
  2. Download the three bigwig files you have just generated

Go back to the genes we looked at earlier: pepT, ycfP. Look at the shape of the signal.
Keep IGV opened.

Go back to working home directory (i.e /shared/projects//EBAII2021_chipseq)

## If you are in 04-Visualization
cd ..

Peak calling with MACS2

Goal: Define the peaks, i.e. the region with a high density of reads, where the studied factor was bound

1 - Choosing a peak-calling program

There are multiple programs to perform the peak-calling step. Some are more directed towards histone marks (broad peaks) while others are specific to narrow peaks (transcription factors). Here we will use the callpeak function of MACS2 (version 2.1.1.20160309) because it’s known to produce generally good results, and it is well-maintained by the developer.

2 - Calling the peaks

  1. Create a directory named 05-PeakCalling and one directory named replicates within to store peaks coordinates.
    mkdir 05-PeakCalling
    mkdir 05-PeakCalling/replicates
    
  2. Go to the newly created directory replicates
    cd 05-PeakCalling/replicates
    
  3. Try out MACS2
    ## Load macs2 in your environment
    module add macs2/2.2.7.1
    macs2 callpeak
    

    This prints the help of the program.

  4. Let’s see the parameters of MACS before launching the mapping:
  1. In a new directory called pool, run macs2 for the pooled replicates A and B by giving both bam files as input treatment files (-t).
# If you are in 05-PeakCalling
mkdir pool
cd pool

# Run macs2 for pooled replicates
macs2 callpeak -t ../../02-Mapping/bam/SRR576933.bam ../../02-Mapping/bam/SRR576934.bam \
-c ../../02-Mapping/bam/SRR576938.bam --format BAM \
--gsize 4639675 --name 'FNR_Anaerobic_pool' --bw 400 \
--fix-bimodal -p 1e-2 &> pool_MACS.out

3 - Analyzing the MACS results

Look at the files that were created by MACS. Which files contains which information ?
How many peaks were detected by MACS ?

At this point, you should have two BED files containing the peak coordinates, one for each replicate.

4 - Combine replicate together (IDR)

In order to take advantage of having biological replicates, we will create a combine set of peaks based on the reproducibility of each individual replicate peak calling. We will use the Irreproducible Discovery Rate (IDR) algorithm.

  1. Create a new directory to store the peak coordinates resulting after idr analysis
    ## If you are in 05-PeakCalling
    mkdir idr
    cd idr
    

    Your directory structure should be like this:

    /shared/projects/<your_project>/EBAII2021_chipseq
    │
    └───data
    │   
    └───01-QualityControl
    │   
    └───02-Mapping
    |    └───index
    |    └───bam
    │   
    └───03-ChIPQualityControls
    │   
    └───04-Visualization
    |
    └───05-PeakCalling
    |    └───replicates
    |    └───pool
    |    └───idr <- you should be in this folder
    
  2. Load the module idr and have a look at its parameters
    ## Load idr in your environment
    module add idr/2.0.4.2
    idr --help
    
  1. Run idr
    idr --samples ../replicates/FNR_Anaerobic_A_peaks.narrowPeak ../replicates/FNR_Anaerobic_B_peaks.narrowPeak \
    --peak-list ../pool/FNR_Anaerobic_pool_peaks.narrowPeak \
    --input-file-type narrowPeak --output-file FNR_anaerobic_idr_peaks.bed \
    --plot
    
  2. Remove IDR and MACS2 from your environment and go back to working home directory (i.e /shared/projects//EBAII2021_chipseq) ```bash module rm macs2/2.1.1.20160309 module rm idr/2.0.4.2

If you are in 05-PeakCalling/idr

cd ../..


### 5 - Visualize peaks into IGV

1. Download the following BED files from the server into your computer to visualise in IGV :

* 05-PeakCalling/repA/FNR_Anaerobic_A_peaks.bed
* 05-PeakCalling/repB/FNR_Anaerobic_B_peaks.bed
* 05-PeakCalling/idr/FNR_anaerobic_idr_peaks.bed

**Go back again to the genes we looked at earlier: pepT, ycfP. Do you see peaks?**


## Motif analysis <a name="motif"></a>
**Goal**: Define binding motif(s) for the ChIPed transcription factor and identify potential cofactors

### 1 - Retrieve the peak sequences corresponding to the peak coordinate file (BED)

For the motif analysis, you first need to extract the sequences corresponding to the peaks. There are several ways to do this (as usual...). If you work on a UCSC-supported organism, the easiest is to use RSAT fetch-sequences or Galaxy. Here, we will use Bedtools, as we have the genome of interest on our computer (Escherichia_coli_K12.fasta).
1. Create a directory named **06-MotifAnalysis** to store data needed for motif analysis
```bash
mkdir 06-MotifAnalysis
  1. Go to the newly created directory
    cd 06-MotifAnalysis
    

Your directory structure should be like this:

/shared/projects/<your_project>/EBAII2021_chipseq
│
└───data
│   
└───01-QualityControl
│   
└───02-Mapping
|    └───index
|    └───bam
│   
└───03-ChIPQualityControls
│   
└───04-Visualization
│   
└───05-PeakCalling
│   
└───06-MotifAnalysis <- you should be in this folder
  1. Extract peak sequence in fasta format ```bash

    First load samtools

    module add samtools/1.10

    Create an index of the genome fasta file

    samtools faidx ../data/Escherichia_coli_K12.fasta

First load bedtools

module add bedtools/2.29.2

Extract fasta sequence from genomic coordinate of peaks

bedtools getfasta -fi ../data/Escherichia_coli_K12.fasta
-bed ../05-PeakCalling/idr/FNR_anaerobic_idr_peaks.bed -fo FNR_anaerobic_idr_peaks.fa

4. Download the file FNR_anaerobic_idr_peaks.fa on your computer

### 2 - Motif discovery with RSAT
1. Open a connection to a Regulatory Sequence Analysis Tools server. You can choose between various website mirrors.

  * Teaching Server  (recommended for this training) [http://pedagogix-tagc.univ-mrs.fr/rsat/](http://pedagogix-tagc.univ-mrs.fr/rsat/)
2. In the left menu, click on **NGS ChIP-seq** and then click on **peak-motifs**. A new page opens, with a form
3. The default peak-motifs web form only displays the essential options. There are only two mandatory parameters.

  * The **title box**, which you will set as **FNR Anaerobic** . The **sequences**, that you will **upload from your computer**, by clicking on the button Choose file, and select the file **FNR_anaerobic_idr_peaks.fa** from your computer.
4. We will now modify some of the advanced options in order to fine-tune the analysis according to your data set.

  * Open the "Reduce peak sequences" title, and make sure the **Cut peak sequences: +/- ** option is set to **0** (we wish to analyze our full dataset)
  * Open the “Motif Discovery parameters” title, and check the **oligomer sizes 6 and 7** (but not 8). Check "Discover over-represented spaced word pairs **[dyad-analysis]**"
  * Under “Compare discovered motifs with databases”, **remove** "JASPAR core vertebrates" and **add RegulonDB prokaryotes** (2015_08) as the studied organism is the bacteria E. coli.
5. You can indicate your email address in order to receive notification of the task submission and completion. This is particularly useful because the full analysis may take some time for very large datasets.
6. Click “**GO**”. As soon as the query has been launched, you should receive an email indicating confirming the task submission, and providing a link to the future result page.
7. The Web page also displays a link, You can already click on this link. The report will be progressively updated during the processing of the workflow.

### 3 - OPTIONAL : Motif discovery with RSAT (short peaks)
1. Restrict the dataset to the summit of the peaks +/- 100bp using bedtools slop. Using bedtools slop to extend genomic coordinates allow not to go beyond chromosome boundaries as the user give the size of chromosomes as input (see fai file).
```bash
bedtools slop -b 100 -i ../05-PeakCalling/replicates/FNR_Anaerobic_A_summits.bed -g ../data/Escherichia_coli_K12.fasta.fai > FNR_Anaerobic_A_summits+-100.bed
  1. Extract the sequences for this BED file ```bash

    Extract fasta sequence from genomic coordinate of peaks

    bedtools getfasta -fi ../data/Escherichia_coli_K12.fasta -bed FNR_Anaerobic_A_summits+-100.bed -fo FNR_Anaerobic_A_summits+-100.fa

Compress the genome file as we won’t need it anymore

gzip ../data/Escherichia_coli_K12.fasta

3. Run RSAT peak-motifs with the same options, but choosing as input file this new dataset (FNR_Anaerobic_A_summits+-100.fa)
and setting the title box to **FNR Anaerobic A summit +/-100bp**

## Peak annotation <a name="annotation"></a>

**Goals**: Associate ChIP-seq peaks to genomic features, draw metagenes, identify closest genes and run ontology analyses

1. Create a directory named **07-PeakAnnotation**
```bash
mkdir 07-PeakAnnotation
  1. Go to the newly created directory
    cd 07-PeakAnnotation
    

1-Associate peaks to closest genes

annotatePeaks.pl from the Homer suite associates peaks with nearby genes.

  1. We will need the annotation file data/Escherichia_coli_K_12_MG1655.annotation.fixed.gtf.gz and the genome file data/Escherichia_coli_K12.fasta.gz. First start by uncompress the files ```bash

    Uncompress annotation file

    gunzip ../data/Escherichia_coli_K_12_MG1655.annotation.fixed.gtf.gz

Uncompress genome file

gunzip ../data/Escherichia_coli_K12.fasta.gz

2. Create a file suitable for annotatePeaks.pl.
```bash
cut -f1-5 ../05-PeakCalling/idr/FNR_anaerobic_idr_peaks.bed | \
  awk -F "\t" '{print $0"\t+"}'  > FNR_anaerobic_idr_peaks.bed
  1. Try annotatePeaks.pl ```bash

    First load bedtools

    module add homer/4.10

run Homer annotatePeaks

annotatePeaks.pl –help

Let's see the parameters:

annotatePeaks.pl peak/BEDfile genome > outputfile
	User defined annotation files (default is UCSC refGene annotation):
		annotatePeaks.pl accepts GTF (gene transfer formatted) files to annotate positions relative
		to custom annotations, such as those from de novo transcript discovery or Gencode.

		-gtf <gtf format file> (Use -gff and -gff3 if appropriate, but GTF is better)


4. Annotation peaks with nearby genes with Homer
```bash
annotatePeaks.pl \
  FNR_anaerobic_idr_peaks.bed \
  ../data/Escherichia_coli_K12.fasta \
  -gtf ../data/Escherichia_coli_K_12_MG1655.annotation.fixed.gtf \
  > FNR_anaerobic_idr_annotated_peaks.tsv
  1. Open an Rstudio tab in Jupyterhub launchrstudio

  2. Add gene symbol annotation using R ```R

set working directory

setwd(“/shared/projects//EBAII2021_chipseq/07-PeakAnnotation")

read the file with peaks annotated with homer

data are loaded into a data frame

sep=”\t”: this is a tab separated file

h=T: there is a line with headers (ie. column names)

d <- read.table(“FNR_anaerobic_idr_annotated_peaks.tsv”, sep=”\t”, h=T)

Load a 2-columns files which contains in the first column gene IDs

and in the second column gene symbols

data are loaded into a data frame

h=F: there is no header line

gene.symbol <- read.table(“../data/Escherichia_coli_K_12_MG1655.annotation.tsv.gz”, h=F)

Merge the 2 data frames based on a common field

by.x gives the columns name in which the common field is for the d data frame

by.y gives the columns name in which the common field is for the gene.symbol data frame

d contains several columns with no information. We select only interesting columns

d.annot <- merge(d[,c(1,2,3,4,5,6,8,10,11)], gene.symbol, by.x=”Nearest.PromoterID”, by.y=”V1”)

Change column names of the resulting data frame

colnames(d.annot)[2] <- “PeakID” # name the 2d column of the new file “PeakID” colnames(d.annot)[dim(d.annot)[2]] <- “Gene.Symbol”

output the merged data frame to a file named “FNR_anaerobic_idr_final_peaks_annotation.tsv”

col.names=T: output column names

row.names=F: don’t output row names

sep=”\t”: table fields are separated by tabs

quote=F: don’t put quote around text.

write.table(d.annot, “FNR_anaerobic_idr_final_peaks_annotation.tsv”, col.names=T, row.names=F, sep=”\t”, quote=F)


**What information is listed in each column of the file?**

**In which column number is the official gene symbol of the nearest gene?**

**What are all the possible gene types?**

7. Retrieve the list of closest genes

```bash
tail -n 2 FNR_anaerobic_idr_final_peaks_annotation.tsv | awk '{print $11}'
  1. Count occurences of each associated feature type
    # sort | uniq -c to list and count occurences of each item
    tail -n +2 FNR_anaerobic_idr_final_peaks_annotation.tsv | awk '{print $8}' | sort | uniq -c
    

How many peaks are associated with promoters ?

tail -n +2 FNR_anaerobic_idr_final_peaks_annotation.tsv | awk '{if ($8=="promoter-TSS") print $11}'

Is the number of genes in your file consistent with the previous reply?

tail -n +2 FNR_anaerobic_idr_final_peaks_annotation.tsv | awk '{if ($8=="promoter-TSS") print $11}' | wc -l
tail -n +2 FNR_anaerobic_idr_final_peaks_annotation.tsv | awk '{if ($8=="promoter-TSS") print $11}' \
> FNR_anaerobic_idr_final_peaks_annotation_officialGeneSymbols.tsv
  1. Compress back the annotation file
    ## Compress annotation file
    gzip ../data/Escherichia_coli_K_12_MG1655.annotation.fixed.gtf
    

Go back to working home directory (i.e /shared/projects/training/<login>/EBAII2021_chipseq)

## If you are in 07-PeakAnnotation
cd ..

2- Search for Biological Processes, Molecular Functions or Cellular Compartments enrichment

This gene list can then be used with Gene Ontology search tools such as Database for Annotation, Visualization and Integrated Discovery (DAVID) or Ingenuity Pathway Analysis (IPA).

Input your gene list (filename : FNR_anaerobic_idr_final_peaks_annotation_officialGeneSymbols.tsv) on the DAVID website: https://david.ncifcrf.gov/

Are there biological processes enriched in the list of genes associated to the peaks?

Are these genes enriched in some KEGG map?

Bonus: Annotation of ChIP-peaks using R tools

In this part, we will use a different set of peaks obtained using a peak caller from a set of p300 ChIP-seq experiments in different mouse embryonic tissues (midbrain, forebrain and limb).

1 - Obtain the bed files from GEO

From now on, we will work locally on your personal machine.

  1. We will download the already called peak files in bed format from GEO. Create a new folder and go in it.
    cd /shared/projects/<your_project>/EBAII2021_chipseq
    mkdir 07-PeakAnnotation-bonus
    cd 07-PeakAnnotation-bonus
    
  2. Search for the dataset GSE13845 either using Google or from the front page of GEO
  3. On the description page, find the three GSM files, and click on each of then
  4. On each page, select and download the GSMxxxxx_p300_peaks.txt.gz file to the newly created folder (where xxxxx represents the GSM number) You should now have downloaded 3 files:

    GSM348064_p300_peaks.txt.gz (Forebrain) GSM348065_p300_peaks.txt.gz (Midbrain) GSM348066_p300_peaks.txt.gz (limb)

Beware: Make sure to check which genome version was used to call the peaks (remember: this is mouse data!)

2 - Performing a first evaluation of peak sets using R

Now, we will use RStudio to perform the rest of the analysis in R. For the analysis, we will need some R/Bioconductor libraries

  1. Go to Rstudio
    # load the required libraries
    library(RColorBrewer)
    library(ChIPseeker)
    library(TxDb.Mmusculus.UCSC.mm9.knownGene)
    library(org.Mm.eg.db)
    # define the annotation of the mouse genome
    txdb = TxDb.Mmusculus.UCSC.mm9.knownGene
    # define colors
    col = brewer.pal(9,'Set1')
    
  2. read the peak files for the three datasets:
# set the working directory to the folder in which the peaks are stored
setwd("/shared/projects/<your_project>/EBAII2021_chipseq/07-PeakAnnotation-bonus")

# read the peaks for each dataset
peaks.forebrain = readPeakFile('GSM348064_p300_peaks.txt.gz')
peaks.midbrain = readPeakFile('GSM348065_p300_peaks.txt.gz')
peaks.limb = readPeakFile('GSM348066_p300_peaks.txt.gz')

# create a list containing all the peak sets
all.peaks = list(forebrain=peaks.forebrain,
midbrain=peaks.midbrain,
limb=peaks.limb)

The peaks are stored as GenomicRanges object; this is an R format which look like the bed format, but is optimized in terms of memory requirements and speed of execution.

We can start by computing some basic statistics on the peak sets.

How many peaks?

# check the number of peaks for the forebrain dataset
length(peaks.forebrain)

# compute the number of peaks for all datasets using the list object
sapply(all.peaks,length)

# display this as a barplot
barplot(sapply(all.peaks,length),col=col)

How large are these peaks?

# statistics on the peak length for forebrain
summary(width(peaks.forebrain))

# size distribution of the peaks
peaks.width = lapply(all.peaks,width)
lapply(peaks.width,summary)

# boxplot of the sizes
boxplot(peaks.width,col=col)

What is the score of these peaks?

Can you adapt the previous code to display a boxplot of the peak score distribution for the Forebrain peak set (column Maximum.Peak.Height)?

Where are the peaks located?

We can now display the genomic distribution of the peaks along the chromosomes, including the peak scores, using the covplot function from ChIPSeeker:

# genome wide distribution
covplot(peaks.forebrain, weightCol="Maximum.Peak.Height")

Exercice: use the option “lower” in covplot to display only the peaks with a score (Max.Peak.Height) above 10

How does the signal look like at TSS?

In addition to the genome wide plot, we can check if there is a tendency for the peaks to be located close to gene promoters.

# define gene promoters
promoter = getPromoters(TxDb=txdb, upstream=5000, downstream=5000)

# compute the density of peaks within the promoter regions
tagMatrix = getTagMatrix(peaks.limb, windows=promoter)

# plot the density
tagHeatmap(tagMatrix, xlim=c(-5000, 5000), color="red")

3 - Functional annotation of the peaks

We can now assign the peaks to the closest genes and genomic compartments (introns, exons, promoters, distal regions, etc…) This is done using the function annotatePeak which compares the peak files with the annotation file of the mouse genome. This function returns a complex object which contains all this information.

peakAnno.forebrain = annotatePeak(peaks.forebrain, tssRegion=c(-3000, 3000), TxDb=txdb, annoDb="org.Mm.eg.db")
peakAnno.midbrain = annotatePeak(peaks.midbrain, tssRegion=c(-3000, 3000), TxDb=txdb, annoDb="org.Mm.eg.db")
peakAnno.limb = annotatePeak(peaks.limb, tssRegion=c(-3000, 3000), TxDb=txdb, annoDb="org.Mm.eg.db")

genomic localization

We can now analyze more in details the localization of the peaks (introns, exons, promoters, distal regions,…)

# distribution of genomic compartments for forebrain peaks
plotAnnoPie(peakAnno.forebrain)

# for all the peaks
plotAnnoBar(list(forebrain=peakAnno.forebrain, midbrain=peakAnno.midbrain,limb=peakAnno.limb))

Question: do you see differences between the three peak sets?

functional annotation

An important step in ChIP-seq analysis is to interpret genes that are located close to the ChIP peaks. Hence, we need to

  1. assign genes to peaks
  2. compute functional enrichments of the target genes.

Beware: By doing so, we assume that the target gene of the peak is always the closest one. Hi-C/4C analysis have shown that in higher eukaryotes, this is not always the case. However, in the absence of data on the real target gene of ChIP-peaks, we can work with this approximation.

We will compute the enrichment of the Gene Ontology “Biological Process” categories in the set of putative target genes.

# load the library
library(clusterProfiler)

# define the list of all mouse genes as a universe for the enrichment analysis
universe = mappedkeys(org.Mm.egACCNUM)

## extract the gene IDs of the forebrain target genes
genes.forebrain = peakAnno.forebrain@anno$geneId
ego.forebrain = enrichGO(gene          = genes.forebrain,
                universe      = universe,
                OrgDb         = org.Mm.eg.db,
                ont           = "BP",
                pAdjustMethod = "BH",
                pvalueCutoff  = 0.01,
                qvalueCutoff  = 0.05,
        readable      = TRUE)

# display the results as barplots        
barplot(ego.forebrain,showCategory=10)

Question: do you see an enrichment of the expected categories? What does the x-axis mean? What does the color mean?

Exercise: redo this analysis for the limb dataset and check if the enriched categories make sense.

FAQ

How to extract peaks from the supplementary data of a publication ?

The processed peaks (BED file) is sometimes available on the GEO website, or in supplementary data. Unfortunately, most of the time, the peak coordinates are embedded into supplementary tables and thus not usable “as is”. This is the case for the studied article. To be able to use these peaks (visualize them in a genome browser, compare them with the peaks found with another program, perform downstream analyses…), you will need to (re)-create a BED file from the information available. Here, Table S5 provides the coordinates of the summit of the peaks. The coordinates are for the same assembly as we used.

  1. copy/paste the first column into a new file, and save it as retained_peaks.txt
  2. use a PERL command (or awk if you know this language) to create a BED-formatted file. As we need start and end coordinates, we will arbitrarily take +/-50bp around the summit.
    perl -lane 'print "gi|49175990|ref|NC_000913.2|\t".($F[0]-50)."\t".($F[0]+50)."\t" ' retained_peaks.txt > retained_peaks.bed
    
  3. The BED file looks like this:
    gi 49175990 ref NC_000913.2 120 220
    gi 49175990 ref NC_000913.2 20536 20636
    gi 49175990 ref NC_000913.2 29565 29665
    gi 49175990 ref NC_000913.2 34015 34115
  4. Depending on the available information, the command will be different.

How to obtain the annotation (=Gene) GTF file for IGV?

Annotation files can be found on genome websites, NCBI FTP server, Ensembl, … However, IGV required GFF format, or BED format, which are often not directly available. Here, I downloaded the annotation from the UCSC Table browser as “Escherichia_coli_K_12_MG1655.annotation.gtf”. Then, I changed the “chr” to the name of our genome with the following PERL command:

perl -pe 's/^chr/gi\|49175990\|ref\|NC_000913.2\|/' Escherichia_coli_K_12_MG1655.annotation.gtf > Escherichia_coli_K_12_MG1655.annotation.fixed.gtf

This file will work directly in IGV

References

Dobin, A., Davis, C.A., Schlesinger, F., Drenkow, J., Zaleski, C., Jha, S., Batut, P., Chaisson, M., and Gingeras, T.R. (2013). STAR: ultrafast universal RNA-seq aligner. Bioinformatics 29, 15–21. Langmead, B., and Salzberg, S.L. (2012). Fast gapped-read alignment with Bowtie 2. Nat. Methods 9, 357–359. Langmead, B., Trapnell, C., Pop, M., and Salzberg, S.L. (2009). Ultrafast and memory-efficient alignment of short DNA sequences to the human genome. Genome Biol. 10, R25. Li, H., and Durbin, R. (2009). Fast and accurate short read alignment with Burrows–Wheeler transform. Bioinformatics 25, 1754–1760. Ramírez, F., Ryan, D.P., Grüning, B., Bhardwaj, V., Kilpert, F., Richter, A.S., Heyne, S., Dündar, F., and Manke, T. (2016). deepTools2: a next generation web server for deep-sequencing data analysis. Nucleic Acids Res. 44, W160–W165.