How to compare bigwig tracks of two ATAC libraries.

Question

[a.rex]

I have two different ATAC-seq libraries that I wish to compare on a genome browser. I have used Macs to generate bedgraph files for each individual library using the command:

macs2 callpeak -t macs2 callpeak -t bamfile --outdir /path/to/ -f BAMPE --keep-dup all --pvalue 1e-2 --call-summits --bdg

I then convert bdg to bigwig.

How can I best normalise the two ATAC-seq libraries that I want to compare by library size? Can I do this in MACS? Or do I do it after I have generated the separate bdg files?

Thank you,

Answer 1

[ATpoint]

I never use these bedGraphs from macs2. They always look kind of weird (a very scientific statement, I know). You can generate a read-per-million normalized bigwig with this command, using the latest deeptools version:

bamCoverage --bam in.bam -o out_normalized.bigwig -bs 1 --normalizeUsing CPM -e

[a.rex]

Presumably as it’s rpm you can then directly compare the two tracks? Also is the bam the one directly following mapping?

[ATpoint]

Well, RPM is good enough for visual inspection. The BAM should be filtered for duplicate reads, e.g. with Picard MarkDuplicates, sambamba markdup or samblaster.

[a.rex]

Ah....I get this error: bamCoverage: error: unrecognized arguments: CPM any clues?

[ATpoint]

Update to the latest version of deeptools. It was introduced recently. Alternatively,do the following (a bit slower but who cares):

## Calculate the normalization factor
SCALE_FACTOR=$(bc <<< "scale=8;1000000/$(samtools idxstats $1 | mawk '{SUM+=$3} END {print SUM}')")

## Get the bedgraph
bedtools genomecov -bg -pc -ibam in.bam -scale ${SCALE_FACTOR} > out.bedGraph

## to bigwig
bedGraphToBigWig out.bedGraph $CHROMSIZES ${BAM%.bam}_CPM.bigwig && rm out_CPM.bigwig

Or, my personal favorit if you are on Linux (not Mac), get mosdepth. This script should give you a normalized bigwig with only one or two cores in a few minutes. Mosdepth is really awesomely fast.

#!/bin/bash

BAM=$1
CHROMSIZES=$2

## 1. Make the depth file:
mosdepth -t 2 ${1%.bam} $1

## 2. Calculate scaling factor:
SCALE_FACTOR=$(bc <<< "scale=8;1000000/$(samtools idxstats $1 | mawk '{SUM+=$3} END {print SUM}')")

## 3. Normalize:
mawk -v SF=${SCALE_FACTOR} 'OFS="\t" {print $1, $2, $3, $4*SF}' <(bgzip -c -d ${BAM%.bam}.per-base.bed.gz) | sort -k1,1 -k2,2n > ${BAM%.bam}_norm.bedGraph.tmp

## 4. to bigwig (cannot read from stdin so far:)
bedGraphToBigWig ${BAM%.bam}_norm.bedGraph.tmp $CHROMSIZES ${BAM%.bam}_CPM.bigwig && rm ${BAM%.bam}_norm.bedGraph.tmp ${BAM%.bam}.per-base* ${BAM%.bam}*mosdepth*

## Use:
./script.sh in.bam chromsizes.txt

[msimmer92]

Does deeptools bamCoverage only accept bams as an input? Because it would be nice to be able to use it directly with the bed files from the MACS2 output (if bigwig tracks with the peak calling results are desired for visualization).

[ATpoint]

No, it doesn't, but bedtools genomecov does.