ATAC-Seq Pipeline Installation

git clone https://github.com/tobiasrausch/ATACseq.git

cd ATACseq

make all

If one of the above commands fail your operating system probably lacks some build essentials. These are usually pre-installed but if you lack them you need to install these. For instance, for Ubuntu this would require:

apt-get install build-essential g++ git wget unzip

Building promoter regions for QC and downloading motifs

To annotate motifs and estimate TSS enrichments some simple scripts are included in this repository to download these databases.

cd bed/ && Rscript promoter.R && cd ..

cd motif/ && ./downloadMotifs.sh && cd ..

Running the ATAC-Seq analysis pipeline for a single sample

./src/atac.sh <hg38|hg19|mm10> <read1.fq.gz> <read2.fq.gz> <genome.fa> <output prefix>

Plotting the key ATAC-Seq Quality Control metrics

The pipeline produces at various steps JSON QC files (*.json.gz). You can upload and interactively browse these files at https://gear-genomics.embl.de/alfred/. In addition, the pipeline produces a succinct QC file for each sample. If you have multiple output folders (one for each ATAC-Seq sample) you can simply concatenate the QC metrics of each sample.

head -n 1 ./*/*.key.metrics | grep "TssEnrichment" | uniq > summary.tsv

cat ./*/*.key.metrics | grep -v "TssEnrichment" >> summary.tsv

To plot the distribution for all QC parameters.

Rscript R/metrics.R summary.tsv

ATAC-Seq pipeline output files

The ATAC-Seq pipeline produces various output files.

• Bowtie BAM alignment files filtered for duplicates and mitochondrial reads.
• Quality control output files from alfred, samtools, FastQC and cutadapt adapter filter metrics.
• Macs peak calling files and IDR filtered peak lists.
• Succinct browser tracks in bedGraph format and IGV's tdf format.
• Footprint track of nucleosome positions and/or transcription factor bound DNA.
• Homer motif finding results.

Differential peak calling

Merge peaks across samples and create a raw count matrix.

ls ./Sample1/Sample1.peaks ./Sample2/Sample2.peaks ./SampleN/SampleN.peaks > peaks.lst

ls ./Sample1/Sample1.bam ./Sample2/Sample2.bam ./SampleN/SampleN.bam > bams.lst

./src/count.sh hg19 peaks.lst bams.lst <output prefix>

To call differential peaks on a count matrix for TSS peaks, called counts.tss.gz, using DESeq2 we first need to create a file with sample level information (sample.info). For instance, if you have 2 replicates per condition:

echo -e "name\tcondition" > sample.info

zcat counts.tss.gz | head -n 1 | cut -f 5- | tr '\t' '\n' | sed 's/.final$//' | awk '{print $0"\t"int((NR-1)/2);}' >> sample.info

Rscript R/dpeaks.R counts.tss.gz sample.info

Intersecting peaks with annotation tracks

Peaks can of course be intersected with enhancer or conserved element tracks, i.e.:

cd tracks/ && downloadTracks.sh

bedtools intersect -a ./Sample2/Sample2.peaks -b tracks/conserved.bed

Plotting peak density along all chromosomes

There is a basic Rscript available for plotting peak densities.

Rscript R/karyoplot.R input.peaks

Citation

Tobias Rausch, Markus Hsi-Yang Fritz, Jan O Korbel, Vladimir Benes.
Alfred: Interactive multi-sample BAM alignment statistics, feature counting and feature annotation for long- and short-read sequencing.
Bioinformatics. 2018 Dec 6.

B Erarslan, JB Kunz, T Rausch, P Richter-Pechanska et al.
Chromatin accessibility landscape of pediatric T‐lymphoblastic leukemia and human T‐cell precursors
EMBO Mol Med (2020)

License

This ATAC-Seq pipeline is distributed under the BSD 3-Clause license.