Creating coverage files containing paired-end fragments of particular sizes

Question

[mmmmcandrew]

Hi all- This is somewhat similar to this question, but my question goes a little deeper than this. My goal is to create coverage files (bigwig) that contain ONLY paired-end fragments below a certain basepair size threshold from BAM files. I thought that using deepTools bamCoverage would allow me to do this relatively quickly and painlessly by using the --maxFragmentLength argument. So, I generated bigwig files that should contain only paired-end fragments at or below certain size thresholds, then used UCSC bigWigToWig to generate wig files. I would then like to convert these to bed files. I have done this routinely for other bigWig files generated with bamCoverage. In the past, my wig files had the format:

fixedStep chrom=chr1 start=3000001 step=1 span=1

0.17

0.17

0.17

0.17

etc.

However, the most recent wig files I generated have the format:

bedGraph section chr1:0-3151026

chr1 0 3000006 0

chr1 3000006 3000016 0.065786

And the subsequent conversion to bed fails. My question is two-fold: 1) is this the best way of doing this? and 2) do I simply need to use the script from the quoted post, or does this have something to do with an indexing problem or something I'm unaware of? Note that in my previous conversions to bed using BedOps convert2bed, I always used the [--zero-indexed] argument to ensure that these files would be compatible with other files I am using. Thanks!

[Devon Ryan]

Which version of deepTools did you use? What you're describing should work correctly.

[mmmmcandrew]

I believe I'm using v. 2.4.0. That's what _version.py says anyway.

[Devon Ryan]

I suspect this was fixed in 2.4.1, but you might as well upgrade to 2.5 anyway.

[mmmmcandrew]

Okie doke, I've upgraded to 2.5 and I'm going to regenerate the bigwig files, then move through the process of creating wig and bed. It's curious, though, that I've used the deeptools 2.4 version and not had this issue before. I'll let you know.

[Devon Ryan]

If that doesn't work then I can try to fix this if you make the BAM file available and let me know the value you're using for the filtering.

[mmmmcandrew]

It seems that the new bigwig files are still having the same issue when converted to wig, and the wig header is the same as above. Do you have a suggestion for where to make the file available? The file is 28 GB. I have been trying to use bamCoverage with --maxFragmentLength of either 100 bp, 115 bp, or 135 bp.

[Devon Ryan]

You can create an account on deeptools.ie-freiburg.mpg.de and then FTP the file there.

[mmmmcandrew]

OK, after a couple of false starts, I have created an account and used FTP via the Unix terminal to transfer the BAM file to my home directory. It is called "rMm.bam". Let me know if you have any trouble. Thanks for looking into this!

[Devon Ryan]

I found your BAM file but I'm not able to reproduce what you're seeing. I tried a simple bamCoverage -b rMm.bam -o rMm.bw --maxFragmentLength 100 -p 20 and the output bigWig file looks appropriate. What exact command were you using?

[mmmmcandrew]

I was using the command

bamCoverage -b rMm.bam -o rMm_100bp.bw -of bigwig -bs 1 --maxFragmentLength 100

I only varied it when changing the base pairs of the fragment lengths.

[Devon Ryan]

Sorry, I still can't reproduce the issue. I ran what you posted and then in python:

import pyBigWig
bw = pyBigWig.open("rMm.bw")
bw.intervals("chr1", 3000000, 3000700)

The output was:

((0, 3000628, 0.0), (3000628, 3000651, 1.0), (3000651, 3000678, 2.0), (3000678, 3000693, 1.0), (3000693, 3000695, 2.0), (3000695, 3000699, 3.0), (3000699, 3000728, 2.0))

So it has appropriate regions and the bin size is 1. Shall I just send you the bigWig file?

[mmmmcandrew]

Hmmm.... OK, I am getting the same output after running the same command with deepTools 2.5. But the issue with the wig file is still present. If I use "less" to look at the wig file generated after conversion with UCSC tools, it looks like:

bedGraph section chr1:0-3133866

chr1 0 3000628 0

chr1 3000628 3000651 1

The subsequent wig2bed operation then fails, making me believe it is a problem with the wig file generated.

[Devon Ryan]

I wonder if something is going wrong when converted to wiggle format. Have a look at the original bigWig in IGV and see if it looks more reasonable.

[mmmmcandrew]

It looks like the conversion didn't work because the bigwig file contained multiple regions with a score of 0. When I removed regions with a score of 0, the conversion to wiggle format worked perfectly. Thanks for your help!