Creating coverage files containing paired-end fragments of particular sizes
0
0
Entering edit mode
4.6 years ago
mmmmcandrew ▴ 130

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!

BAM bigwig wig bed paired-end • 2.0k views
0
Entering edit mode

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

0
Entering edit mode

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

0
Entering edit mode

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

0
Entering edit mode

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.

0
Entering edit mode

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.

0
Entering edit mode

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.

0
Entering edit mode

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

0
Entering edit mode

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!

0
Entering edit mode

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?

0
Entering edit mode

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.

0
Entering edit mode

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?

0
Entering edit mode

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.

0
Entering edit mode

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.

0
Entering edit mode

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!