In this tutorial, we are going to show how to use HINT for identification of footprints using open chromatin assay and differential footprint analysis.

DNase-seq footprinting

Download here. Execute the following commands in order to perform footprints identification from histone modification data:

cd HINT_DNaseTest 
rgt-hint --dnase-footprints --output-location=./ --output-prefix=test DNase.bam DNasePeaks.bed

If you want to use bias correction version of HINT:

rgt-hint --dnase-footprints --bias-correction --output-location=./ --output-prefix=test DNase.bam DNasePeaks.bed

Histone footprinting

Download here. Execute the following commands in order to perform footprints identification from histone modification data:

cd HINT_HistoneTest
rgt-hint --histone-footprints --output-location=./ --output-prefix=test histone.bam histonePeaks.bed

DNase+Histone footprinting

Download here. Execute the following commands in order to perform footprints identification from DNase-seq combining with histone modification data:

cd HINT_DNaseHistoneTest
rgt-hint --dnase-histone-footprints --output-location=./ --output-prefix=test DNase.bam H3K4me1.bam H3K4me3.bam regions.bed

ATAC-seq footprinting and differential analysis

we will analyze ATAC-seq data from LSK cells (equivalent to MPP cells), B cells and T CD4 cells obtained from Lara-Astiaso et al 2014. We have performed low-level analysis steps including read alignment and peaks calling in chromosome 1, which can be found here and the script of commands used togenerate these files for B cell ATAC-seq experiments is found here.

1. First, go to the HINT_ATAC directory and generate an output folder for the result files:

mkdir output 

2. Execute the following commands to call footprints on B, LSK and CD4 cells for chromosome 1:

rgt-hint --atac-footprints --organism=mm10 input/B_ATAC_chr1.bam input/B_ATACPeaks_chr1.bed --output-location=output --output-prefix=B_ATAC_chr1_footprints 

rgt-hint --atac-footprints --organism=mm10 input/LSK_ATAC_chr1.bam input/LSK_ATACPeaks_chr1.bed --output-location=output --output-prefix=LSK_ATAC_chr1_footprints 

rgt-hint --atac-footprints --organism=mm10 input/CD4_ATAC_chr1.bam input/CD4_ATACPeaks_chr1.bed --output-location=output --output-prefix=CD4_ATAC_chr1_footprints 

This will generate an output file, i.e ./output/B_ATAC_chr1_footprints.bed, containing the genomic locations of the footprints. HINT also produces a file with ending “.info”, which has general statistics from the analysis as no. of footprints, total number of reads and so on.

We can use IGV to vizualise ATAC-seq signals and footprint predictions in particular loci. First, we can use a special HINT command to generate genomic profiles (bigWig files).

rgt-hint --print-signal --bc-signal --bigWig --organism=mm10 --reads-file=input/B_ATAC_chr1.bam --regions-file=input/B_ATACPeaks_chr1.bed --output-location=output --output-prefix=B_ATAC_chr1

rgt-hint --print-signal --bc-signal --bigWig --organism=mm10 --reads-file=input/CD4_ATAC_chr1.bam --regions-file=input/CD4_ATACPeaks_chr1.bed --output-location=output --output-prefix=CD4_ATAC_chr1

rgt-hint --print-signal --bc-signal --bigWig --organism=mm10 --reads-file=input/LSK_ATAC_chr1.bam --regions-file=input/LSK_ATACPeaks_chr1.bed --output-location=output --output-prefix=LSK_ATAC_chr1

This bigwig file contains number of ATAC-seq (or DNAse-seq) reads at each genomic position as estimated by HINT after signal normalization and cleveage bias correction. This is therefore more accurate than simply looking a coverage profiles of a bam file.

Open all bigwig and footprint files (bed) generated above in IGV. Remember to set the genome version to mm10 beforehand. You can also enrich the data by opening bam files of histone modifications as H3K4me3, H3K4me1 and H3K27ac (provided in result folder). Check for example the genomic profiles around the gene Zp70 in the below, which is part of T cell receptor. We observe that this gene has several open chromatin regions and high levels of histone levels in the start of the gene, which are specific to CD4 T cells. There is also some open chromatin level around exon 3, which is supported by H3K4me1 modifications. This potential enhancer region is also present in distinct degrees in B and LSK cells.

3. An important question when doing footprint analysis is to evaluate which TF motifs overlap with footprints and evaluate the ATAC-seq profiles around these motifs. RGT suite also offers a tool for finding motif predicted binding sites (mpbs). For example, we analyze here motifs from factors SPI1 and ELK4, which were discussed in Lara-Astiaso et al. 2014 to be associated respectively associated to LSK and CD4 cells.

Execute the following commands to do motif matching inside footprints for chromosome 1:

rgt-motifanalysis --matching --organism=mm10 --output-location=output --use-only-motifs=input/motifs.txt output/B_ATAC_chr1_footprints.bed

rgt-motifanalysis --matching --organism=mm10 --output-location=output --use-only-motifs=input/motifs.txt output/CD4_ATAC_chr1_footprints.bed

rgt-motifanalysis --matching --organism=mm10 --output-location=output --use-only-motifs=input/motifs.txt output/LSK_ATAC_chr1_footprints.bed

The file input/motifs.txt contains a list of JASPAR motif ids to be used in the analysis. Ignoring this option will search for all JASPAR motifs. The above commands will generate bed files (i.e. LSK_ATAC_footprints_mpbs.bed) containing mpbs overlapping with distinct footprint regions. The 4th column contains the motif name and the 5th column the bitscore of the motif match. If you open these files in IGV, you will observe that a ELK4 binding site near the 3rd exon of the gene Zp70. Check here for more details about motif matching.

4. Finally, we use HINT to generate average ATAC-seq profiles around binding sites of particular TF. This analysis allows us to inspect the chromatin accessibility and the underlying sequence. Moreover, by comparing the cut profiles from two ATAC-seq libraries (i.s. LSK vs T CD4 cells ), we can get insights on changes in binding in two cells. For this, execute the following commands:

mkdir output/LSK_B
rgt-hint --diff-footprints --organism=mm10 --mpbs-file=result/LSK_B_ATAC_footprints_mpbs.bed --reads-file1=input/LSK_ATAC.bam --reads-file2=input/B_ATAC.bam --output-location=output/LSK_B --output-prefix=LSK_B

mkdir output/LSK_CD4
rgt-hint --diff-footprints --organism=mm10 --mpbs-file=result/LSK_CD4_ATAC_footprints_mpbs.bed --reads-file1=input/LSK_ATAC.bam --reads-file2=input/CD4_ATAC.bam --output-location=output/LSK_CD4 --output-prefix=LSK_CD4

The above commands will generate eps files with a ATAC-seq profile for each of the motifs founds in the provided mpbs bed files as shown in the below. Let’s check the profiles in the comparison LSK and CD4, you will see that ELK4 has higher number of ATAC-seq counts in CD4 cells, while SPI1 has more ATAC-seq in LSK cells. Higher ATAC counts indicates higher activity of the factor in that particular cell. This fits with the results discussed in Lara-Astiaso that SPI1 are more relevant/active in LSK, while ELK4 in CD4 cells.