Benchmarking

For benchmarking and validation purposes, we processed previously published ATAC-seq data, quantified peak coverages, and generated scaled coverage tracks. As comparisons, we used BEDtools for performance benchmarking and validation of raw read counts, and coverage track generation benchmarks with IGVTools and Deeptools (bamCoverage).

Input dataset

For benchmarking purposes, data from ATAC-seq was used:

"SLFN11 Blocks Stressed Replication Forks Independently of ATR. Mol Cell, 2018"

Pubmed and GEO page.

Performance and validation of peak quantification

In the six processed ATAC-seq samples, a total of 32,819 peaks were identified (transposase-accessible chromatin regions).

Performance

Since BAMscale is multithreaded, we quantified peaks using 1, 4 and 8 threads. BEDtools' multicov is a single threaded program, and was run once. BAMscale is >2-fold quicker than BEDtools using 4 threads, with similar performance using one thread. It is important to note, that BEDtools generates raw read counts only, while BAMscale outputs three different normalized tables as well (FPKM, TPM and library size). We suggest using 4 threads, since higher thread counts have minimal impact on performance.

Validation

Correlation of raw read counts between BEDTools and BAMscale was above 0.999 in all cases.

Performance and validation of scaled coverage track generation

Next we compared performance of coverage track generation to IGVTools and Deeptools bamCoverage programs. IGVTools can only generate unscaled tracks, with no read filter functionality. Deeptools "bamCoverage" on the other hand is capable of filtering reads and output scaled coverage tracks.

Performance

BAMscale using a single thread outperforms IGVtools, and outperforms deeptools (using 4 threads) by ~7-fold using a single thread. We see some performance increase by adding additional threads to BAMscale, but >4 threads is not suggested.

Comparison

Generated results from the three methods output similar results, with some minor differences in track height, due to differences in binning technique, and scaling.

Pre-processing ATAC-seq data

1) Aligning fastq files

Raw sequencing reads were downloaded with the fasterq-dump program sratoolkit, and aligned with the BWA mem aligner using default alignment settings (but multithreaded) to the hg19 human genome. Reads aligned with BWA were directly converted from SAM format to BAM format using samtools.

bwa mem -t 8 <hg19.idx> <PAIR1.fq> <PAIR2.fq> \
      | samtools view -Sbh \
      > out.unsorted.bam

2) Sorting reads and marking duplicate reads (alignments)

The aligned and unsorted reads were then sorted by samtools. Sorted reads were then indexed, duplicate reads were marked using picard-tools, and indexed one final time.

samtools sort -o out.sorted.bam out.unsorted.bam

samtools index out.sorted.bam

java -Xmx8g -jar picartools.jar MarkDuplicates \
      I=out.sorted.bam \
      O=out.sorted.dedup.bam \
      M=out.sorted.dedup.metrics

samtools out.sorted.dedup.bam

3) Peak calling

Peaks were called using the MACS2 peak caller. Since the downloaded data was paired-end, the "BAMPE" flag was set.

macs2 callpeak -t out.sorted.dedup.bam --outdir outdir \
      -f BAMPE -g hs --cutoff-analysis -q 0.01 \

4) Generating Uniform peak sets

Called peaks files were merged to create one BED file for all samples (using the unix cat command). Since many peak positions are common across samples, peaks were first sorted based on their position, then merged into a unique set of peaks using bedtools.

cat *.narrowPeak \
      | bedtools sort \
      | bedtools merge \
      > merged_peaks.bed   

5) Quantifying peaks

Peaks were quantified using the BAMscale "cov" function. The output of this step are tables with: raw read conts, FPKM, TPM and library size scaled peak scores. For comparison, we calculated raw read counts using the BEDTools "multicov" program.

./BAMscale cov --bed merged_peaks.bed \
      --bam SAMPLE1.sorted.dedup.bam \
      --bam SAMPLE1.sorted.dedup.bam \
      --bam SAMPLEn.sorted.dedup.bam

bedtools multicov -bed merged_peaks.bed \
      -p \
      -bams SAMPLE1.sorted.dedup.bam SAMPLE2.sorted.dedup.bam SAMPLE3.sorted.dedup.bam \
      > Raw_coverages.bedtools.txt

6) Generating scaled coverage tracks with BAMscale

Scaled coverage tracks for each BAM file were generated with the BAMscale "scale" function. The default settings were used, meaning the bin coverages were scaled to the number of aligned bases divided by the gene length. As comparison, we generated coverage tracks with IGVtools (using the hg19 chromSize file) and deeptools bamcoverage

./BAMscale scale --bam SAMPLE1.sorted.dedup.bam

igvtools count -z 5 SAMPLE1.sorted.dedup.bam SAMPLE1.sorted.dedup.bam.tdf hg19.chrom.sizes

bamCoverage --ignoreDuplicates --binSize 5 --effectiveGenomeSize 2864785220 \
    -o SAMPLE1.sorted.dedup.bam.bw \
    --normalizeUsing RPGC \
    -b SAMPLE1.sorted.dedup.bam