MNase-seq / fragment size selection

Question

[Constantine]

Hello guys

I have an MNase-seq bam file (~400 million reads) that I'm trying to split into two different fragment sizes

1) 1-80 fragments (for TF-sized nucleosomes)

2) 135 -165 fragments (for mononucleosomes)

To do that I aligned the fastq files with Bowtie2 using the parameters "-I 1 -X 80" and "-I 135 -X 165" and ended up with a 368011030 mapped-read bam file and a 368011030 mapped-read bam file respectively (YES! The same amount of reads)

Later I used the bamCoverage command (--MNase) from Deeptools to normalise the 2 files and subsequently plot Coverage over TSS

The bam file containing the mononucleosomes was properly normalised and Coverage was plotted. On the other hand, the bam containing the TF-sized nucleosomes (~80 bp fragments) returned an error during normalisation

Error: The generated bedGraphFile was empty. Please adjust
your deepTools settings and check your input files.

Is the above approach I used for splitting files appropriate? The fact that I got the exact same amount of reads after alignment is bugging me

Any help would be much appreciated. This is the first time I'm analysing MNase-Seq data and I'm completely lost

Thanks

[Carlo Yague]

Just out of curiosity, what is a TF-sized nucleosome ? Do you mean MNase protection provided by transcription factors binding ?

[Constantine]

Yeap that's it. Sorry, hope it didn't sound confusing the way I addressed it

[Carlo Yague]

No not at all, I'm just not familiar with the concept but it's interesting :)

[rbronste]

Also just curious what you're basing the 80bp on, like any particular papers/studies? Just wanted to take a look, as Im trying to also track down good boundaries to set. Thanks!

[Constantine]

This paper is a great example.

Answer 1

[Devon Ryan]

The --MNase option effectively ignores fragments <130 bases or >200 bases. So passing in your TF-sized-fragments BAM file will result in all of the reads getting ignored, which presumably would result in that error. You can change this behavior by editing deeptools/bamCoverage.py. There's a line near that bottom of the form:

# TODO: this size range could be an input parameter.
if read.is_proper_pair and not read.is_reverse and 130 < abs(read.tlen) < 200:

The next release has a fragment length parameter, though --MNase doesn't use it (maybe I should change that before the 2.3 release this week).

Anyway, change the aforementioned line to this instead:

if read.is_proper_pair and not read.is_reverse:

That'll remove the fragment length filtering. Alternatively, you can just change the bounds:

if read.is_proper_pair and not read.is_reverse and 1 < abs(read.tlen) < 81:

Edit: We've discussed this internally and for the 2.3 release we'll change the code to honour --minFragmentLength and --maxFragmentLength, with defaults of 130 and 200 if nothing else (or 0) is specified. That should allow you to more easily do this without (A) running bowtie2 twice or (B) prefiltering BAM files.

[Constantine]

Thanks a lot man.

[Devon Ryan]

No worries, I've now merged this into the branch that will become version 2.3.