Biostar Beta. Not for public use.
Question: How to compare bigwig tracks of two ATAC libraries.
1
Entering edit mode

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,

ADD COMMENTlink 19 months ago a.rex • 190 • updated 19 months ago ATpoint 17k
4
Entering edit mode

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
ADD COMMENTlink 19 months ago ATpoint 17k
Entering edit mode
1

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

ADD REPLYlink 19 months ago
a.rex
• 190
Entering edit mode
0

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.

ADD REPLYlink 19 months ago
ATpoint
17k
Entering edit mode
0

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

ADD REPLYlink 19 months ago
a.rex
• 190
Entering edit mode
1

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
ADD REPLYlink 19 months ago
ATpoint
17k
Entering edit mode
0

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).

ADD REPLYlink 9 months ago
msimmer92
• 180
Entering edit mode
1

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

ADD REPLYlink 9 months ago
ATpoint
17k

Login before adding your answer.

Similar Posts
Loading Similar Posts
Powered by the version 2.0