search

arq5x / bedtools2

bedtools - the swiss army knife for genome arithmetic
MIT License
929 stars 287 forks source link

Empty .bam file from bedtools intersect #81

Closed cwarden45 closed 9 years ago

cwarden45 commented 10 years ago

Hi,

I am running a Duplex-Seq pipeline that requires using bedtools.

All the sequence comes from a relatively small amplicon (90 bp), which is a p53 sequence. I have a sorted, indexed .bam file as the -abam input an a one-line .bed file for the p53 region as the -b input, but the resulting .bam file (which should have the overlapping reads) is empty. I have checked the alignment, and the coordinates are correct. In fact, I also tried using a .bed file specifying the entire chr17, and I still don't get any reads in the resulting .bam file

Can you please help me troubleshoot what might be going on?

Thanks, Charles

arq5x commented 10 years ago

I suspect the chromosome ids in your BAM file do not exactly match owe in your BED file (e.g 17 v. chr17)

On May 21, 2014, at 8:28 PM, cwarden45 notifications@github.com wrote:

Hi,

I am running a Duplex-Seq pipeline that requires using bedtools.

All the sequence comes from a relatively small amplicon (90 bp), which is a p53 sequence. I have a sorted, indexed .bam file as the -abam input an a one-line .bed file for the p53 region as the -b input, but the resulting .bam file (which should have the overlapping reads) is empty. I have checked the alignment, and the coordinates are correct. In fact, I also tried using a .bed file specifying the entire chr17, and I still don't get any reads in the resulting .bam file

Can you please help me troubleshoot what might be going on?

Thanks, Charles

— Reply to this email directly or view it on GitHub.

cwarden45 commented 10 years ago

Thank you for your very prompt response.

I also initially suspected this was the case, but I double-checked that the BED file and reference format (specified in idxstats, for example) were the same (in this case, both were chr17).

Any other suggestions? Please let me know if you need other information. The only thing that I can think of that is kind of unusual is that the paired ends are being analyzed separately for several steps, but I would assume it should appear just like a single-end experiment (with a modified tag for the read name) from the perspective of bedtools.

Thanks, Charles

nkindlon commented 10 years ago

Hi Charles, I'm Neil Kindlon, a developer in Aaron's lab. I wrote the newer versions of intersect, so hopefully I can help with this. Would it be possible to send us your data files, the command you ran, and the version number? Thanks!

cwarden45 commented 10 years ago

Hi Neil,

Sure - I am using bedtools 2.17.0.

I've attached the bed files. I've uploaded the .bam and index in a .zip file to Google Drive, which I believe you can access below.​ 7073.ziphttps://docs.google.com/file/d/0B1xpw6_kQMKubWVvZlM1UmdJRkU/edit?usp=drive_web ​.

I've run the command using bedtools intersect and intersectBed, as follows:

bedtools intersect -abam $sorted_bam1 -b $target_bed > $target_bam1

Please let me know if you need anything else.

Thanks, Charles

On Thu, May 22, 2014 at 8:32 AM, Neil Kindlon notifications@github.comwrote:

Hi Charles, I'm Neil Kindlon, a developer in Aaron's lab. I wrote the newer versions of intersect, so hopefully I can help with this. Would it be possible to send us your data files, the command you ran, and the version number? Thanks!

— Reply to this email directly or view it on GitHubhttps://github.com/arq5x/bedtools2/issues/81#issuecomment-43904173 .

nkindlon commented 10 years ago

Hi Chris, Thanks for the info! I was able to get the bam file, but I didn't see your target bed file in there. Sorry to keep asking for stuff, but could you send that as well, please? I did make a "dummy" bed file, but intersecting your bam file worked with it worked correctly and produced a non-empty bam file, meaning I wasn't able to reproduce the problem with that. Thanks for your patience. -Neil

Neil Kindlon, Genomics Bioinformatics Software Engineer Quinlan Laboratory, http://quinlanlab.org/ Center for Public Health Genomics, University of Virginia School of Medicine nek3d@virginia.edu Work: 434-243-5395 Cell: 203-928-9907

From: cwarden45 [notifications@github.com] Sent: Thursday, May 22, 2014 2:01 PM To: arq5x/bedtools2 Cc: Kindlon, Neil Edward (nek3d) Subject: Re: [bedtools2] Empty .bam file from bedtools intersect (#81)

Hi Neil,

Sure - I am using bedtools 2.17.0.

I've attached the bed files. I've uploaded the .bam and index in a .zip file to Google Drive, which I believe you can access below.​ 7073.ziphttps://docs.google.com/file/d/0B1xpw6_kQMKubWVvZlM1UmdJRkU/edit?usp=drive_web ​.

I've run the command using bedtools intersect and intersectBed, as follows:

bedtools intersect -abam $sorted_bam1 -b $target_bed > $target_bam1

Please let me know if you need anything else.

Thanks, Charles

On Thu, May 22, 2014 at 8:32 AM, Neil Kindlon notifications@github.comwrote:

Hi Charles, I'm Neil Kindlon, a developer in Aaron's lab. I wrote the newer versions of intersect, so hopefully I can help with this. Would it be possible to send us your data files, the command you ran, and the version number? Thanks!

— Reply to this email directly or view it on GitHubhttps://github.com/arq5x/bedtools2/issues/81#issuecomment-43904173 .

— Reply to this email directly or view it on GitHubhttps://github.com/arq5x/bedtools2/issues/81#issuecomment-43922716.

cwarden45 commented 10 years ago

Hi Neil,

I copied over the code from a Perl wrapper that I am using. The target .bed files are the p53.bed (the actual target) and the chr17.bed file (my dummy bed file). I've reattached them here for convenience.

Also, just to be clear, the resulting .bam file is technically a few KB, but I don't think it actually has any reads (for example, I don't see anything with samtools view).

Thanks, Charles

On Thu, May 22, 2014 at 11:38 AM, Neil Kindlon notifications@github.comwrote:

Hi Chris, Thanks for the info! I was able to get the bam file, but I didn't see your target bed file in there. Sorry to keep asking for stuff, but could you send that as well, please? I did make a "dummy" bed file, but intersecting your bam file worked with it worked correctly and produced a non-empty bam file, meaning I wasn't able to reproduce the problem with that. Thanks for your patience. -Neil

Neil Kindlon, Genomics Bioinformatics Software Engineer Quinlan Laboratory, http://quinlanlab.org/ Center for Public Health Genomics, University of Virginia School of Medicine nek3d@virginia.edu Work: 434-243-5395 Cell: 203-928-9907

From: cwarden45 [notifications@github.com] Sent: Thursday, May 22, 2014 2:01 PM To: arq5x/bedtools2 Cc: Kindlon, Neil Edward (nek3d) Subject: Re: [bedtools2] Empty .bam file from bedtools intersect (#81)

Hi Neil,

Sure - I am using bedtools 2.17.0.

I've attached the bed files. I've uploaded the .bam and index in a .zip file to Google Drive, which I believe you can access below.​ 7073.zip< https://docs.google.com/file/d/0B1xpw6_kQMKubWVvZlM1UmdJRkU/edit?usp=drive_web>

​.

I've run the command using bedtools intersect and intersectBed, as follows:

bedtools intersect -abam $sorted_bam1 -b $target_bed > $target_bam1

Please let me know if you need anything else.

Thanks, Charles

On Thu, May 22, 2014 at 8:32 AM, Neil Kindlon notifications@github.comwrote:

Hi Charles, I'm Neil Kindlon, a developer in Aaron's lab. I wrote the newer versions of intersect, so hopefully I can help with this. Would it be possible to send us your data files, the command you ran, and the version number? Thanks!

— Reply to this email directly or view it on GitHub< https://github.com/arq5x/bedtools2/issues/81#issuecomment-43904173> .

— Reply to this email directly or view it on GitHub< https://github.com/arq5x/bedtools2/issues/81#issuecomment-43922716>.

— Reply to this email directly or view it on GitHubhttps://github.com/arq5x/bedtools2/issues/81#issuecomment-43927005 .

nkindlon commented 10 years ago

Hi Charles, I'm so sorry, I didn't see get the p53.bed or chr17.bed file in that last message. Can you please send them to my e-mail address? It's nek3d@virginia.edu. Thanks, Neil

nkindlon commented 10 years ago

Hey Charles, Thanks for the files. I was able to determine that you were getting blank output because of a somewhat embarrassing bug. Turns out your p53.bed and chr17.bed files didn't end in a newline, and that broke out automated type detection. The best thing for you to do now is to update to a newer version of bedtools. We're releasing 2.20 this weekend. This will get you numerous other improvements and bug fixes since 2.17 as well. If, for some reason, you can't do an update right now, you could also just add newlines at the end of those files after the last characters. Sorry about that.

cwarden45 commented 10 years ago

Hi Niel,

Ok, thanks - at least you can reproduce the error, which is what mattered to me ;)

I will just add the new line character on the .bed file and update the program in the near future.

Thanks, Charles

On Fri, May 23, 2014 at 12:49 PM, Neil Kindlon notifications@github.comwrote:

Hey Charles, Thanks for the files. I was able to determine that you were getting blank output because of a somewhat embarrassing bug. Turns out your p53.bed and chr17.bed files didn't end in a newline, and that broke out automated type detection. The best thing for you to do now is to update to a newer version of bedtools. We're releasing 2.20 this weekend. This will get you numerous other improvements and bug fixes since 2.17 as well. If, for some reason, you can't do an update right now, you could also just add newlines at the end of those files after the last characters. Sorry about that.

— Reply to this email directly or view it on GitHubhttps://github.com/arq5x/bedtools2/issues/81#issuecomment-44052976 .

arq5x commented 9 years ago

Closing as this appears to be resolved.