BEDTools coverage is giving empty output when computing coverage of bam vs bed

Question

[DrDoom]

Hi,

I have an exome comprised of 21 samples. I'm trying to follow this in order to plot the coverage of multiple samples in a single plot, using the bed file of the exome capture kit. However I'm stuck at the beginning because I get an empty output / zero bytes file from bedtools coverage:

I first tried bedtools coverage -hist -abam mybam.bam -b exome_capture.bed | grep ^all > mybam.bam.hist.all.txt

Then read that the commands had evolved (am currently using the latest: bedtools 2.29.2) so instead used bedtools coverage -hist -a exome_capture.bed -b mybam.bam | grep ^all > mybam.bam.hist.all.txt

Both gave an empty output.

I finally tried on Galaxy with samtools bedcov and the output was also empty.

So I'm guessing that something must be wrong with either the bam or bed, causing the tools to not recognize the overlaps, but I'm just beginning and I can't see what is wrong with them...

BAM file :

A00461:58:H7L77DRXX:1:2126:15374:1031   99  chr1    65356   0   101M    =65521  266 CCAAGGGAAAACTTGTGAGACTATAAAAGTTAGTCTCAGTACACAAAGCTCAGACTGGCTATTCCCAGATCTCTTCAGGTACATCTAGTCCATTCATAAAG   AB>BBGHFGGGHHEJDJGIFHHDDDGGGIDEDIDIHICIDDHCHCGGIHHICIGHHJHHHDDEIGGCIG=IHIHEICIHDDHCDIHDIDIGCDEICB@BBD   XA:Z:chr19,+106940,101M,0;chr15,-102465404,101M,0;  MC:Z:101M   BD:Z:PPSPSTPSQFFPQPQNNQPPRQPOONFFPQPNQQRPPQRRQROOOPGQSSQRSRSRSSSSQPPQSOSSRSSQQQQQRSSSRROSTRRSSTTTSRSTTPOGQ  MD:Z:101    RG:Z:group1 BI:Z:QQRPUSQSSKKSTRUSRSTRTTPSQQKKTURQTUUTSTUUQURRRRKUVUSUUSUUVUUUQTSVUQTUSUUUTUSVUUSTRURUUUQTUTSRSRTUVPPJS  NM:i:0  MQ:i:0  AS:i:202    XS:i:202    ms:i:3628
A00461:58:H7L77DRXX:1:2142:29713:15969  163 chr1    65358   0   101M    =65580  323 AAGGGAAAACTTGTGAGACTATAAAAGTTAGTCTCAGTACACAAAGCTCAGACTGGCTATTCCCAGATCTCTTCAGGTACATCAAGTCCATTCATAAAGGG   ABACEEFFFGG/ICIGHFG47CDFFGHC/C@DI4HC@CCGCGBFF@GHIBIGHHJGHHDDEIGGCHFDHHI4EICIHDDHCDI)GIDIGB,EHCCCEB@BB   XA:Z:chr15,-102465402,101M,1;chr19,+106942,101M,1;  MC:Z:101M   BD:Z:PPQSPSQGFPQOPNNQPPQPONNMEFPQPNQQRPPQRQQQNNNOFPRRPQSRSRSSSSQPPQSOSSRSSQQQQQRSSSRRORTRQRSTTTSQSTTRQGQSO  MD:Z:83T17  RG:Z:group1 BI:Z:QQSQPRRJJQRPSQPQRPRRNQOPJJSTQPSTTSRSSTPSPQPPJSTSRSSRSTUTTTPSRUTPSTRTTTSTRUTTRSQTQTTTQTUUSQRPSSSQRJSQO  NM:i:1  MQ:i:0  AS:i:192    XS:i:192    ms:i:3689
A00461:58:H7L77DRXX:2:1155:5810:19163   1187    chr1    65358   0   101M    =65580  323 AAGGGAAAACTTGTGAGACTATAAAAGTTAGTCTCAGTACACAAAGCTCAGACTGGATATTCCCAGATCTCTTCAGGTACATCTAGTCCATTCATAAAGGG   A=A?EE+FF-CDICIGHFGGD,DFFGHCEC/DI4/C/C&GCGBFF@GHIBIGHBJG*DD,EI>GCHFD>HIHE/CIHDDHCD/4DADIGBD//C,CEB0BA   XA:Z:chr19,+106942,101M,1;chr15,-102465402,101M,1;  MC:Z:101M   BD:Z:PPQSPSQGFPQOPNNQPPQPONNMEFPQPNQQRPPQRQQQNNNOFPRRPQSRSRSSSSPPPQSOSSRSSQQQQQRSSSRRORTQRSSTTTSQSTTRQGQSO  MD:Z:56C44  RG:Z:group1 BI:Z:QQSQPRRJJQRPSQPQRPRRNQOPJJSTQPSTTSRSSTPSPQPPJSTSRSSRSTUTSTQSRUTPSTRTTTSTRUTTRSQTQTTTPTUUSQRPSSSQRJSQO  NM:i:1  MQ:i:0  AS:i:192    XS:i:192    ms:i:3640
A00461:58:H7L77DRXX:2:1176:25681:25238  163 chr1    65367   0   101M    =65496  230 CTTGTGAGACTATAAAAGTTAGTCTCAGTACACAAAGCTCAGACTGGCTATTCCCAGATCTCTTCAGGTACATCTAGTCCATTCATAAAGGGCTTTTAATT   AC@EAHFHFGGCCCFGFHCDDHDHGIBHDCHCHBFGHGGHCHFGGIGHHCDEIGGBIGDIHIHEIBHHCDHCDIHDIDIGCDEICDDGGHHHGHDDC?A>?   XA:Z:chr19,+106951,101M,0;chr15,-102465393,101M,0;  MC:Z:101M   BD:Z:PPQRQPSRQRQONNMEEOPOMPPQOPQRQQQNNNOFPRRPQRQRQRRRRPPPQSOSSRSSQQQQQRSSSRRORSQQRRSSSRQQSSQPHRTPUUSJJOOQP  MD:Z:101    RG:Z:group1 BI:Z:QQQTSQRSQRRNQOOIIRSPORSTSRSSTPSPQPPJSTSRSSQSSTSSSORRTTPSTRTTTSTRUTTRSQTQTTTPTUUTSTRUTTQQJRPNRSQMNPPRQ  NM:i:0  MQ:i:0  AS:i:202    XS:i:202    ms:i:3677
A00461:58:H7L77DRXX:1:2236:10321:5588   163 chr1    65384   0   101M    =65556  273 GTTAGTCTCAGTACACAAAGCTCAGACTGGCTATTCCCAGATCTCTTCAGGTACATCTAGTCCATTCATAAAGGGCTTTTAATTAACCAAGTGGTTTACTA   A?@?FBHGHBHCCGBHBFFHHGIBHGGGJGHHDCDIFFB@GCHGHGDICHHDDHCCIHDIDIGCDDHCCDG:IHHHHEEEDGDEDGHGCFADIHCDC?BB>   XA:Z:chr15,-102465376,101M,0;chr19,+106968,101M,0;  MC:Z:101M   BD:Z:PPQOSRSQPQRPPPMMMNEOQQOPQQRQRRRRPOOPRNRRQRRPPPPPQRSSRRORSQQRRSSSRPQRRPOGQSOSSQHHOORPPPRSTQRSRUUSJORRQ  MD:Z:101    RG:Z:group1 BI:Z:QQQPTTTSRRRSOROPOOIRSRQSSQSSTSSSORQTSORSQSSSRSQTSSQSPTQTTTPTUUTSTRUTTQQKTRPTTRMMQQSRQQSSRORSPSRRNPSSO  NM:i:0  MQ:i:0  AS:i:202    XS:i:202    ms:i:3725

BED file: (Agilent SureSelect DNA - SureSelectXT Human All Exon V5 and clinical - Covered bed file extended by 100bp on either side)

chr1    65409   65725   -
chr1    65731   66073   -
chr1    69381   69700   ref|OR4F5,ref|NM_001005484,ens|ENST00000335137,ccds|CCDS30547.1
chr1    721281  721619  ens|ENST00000593022,ens|ENST00000586288,ens|ENST00000587530,ens|ENST00000585768,ens|ENST00000590848,ens|ENST00000585745,ens|ENST00000591440,ens|ENST00000586928,ens|ENST00000588951,ens|ENST00000591702,ens|ENST00000429505,ens|ENST00000589531,ens|ENST00000358533,mRNA|AK125248,mRNA|AK290103
chr1    721430  721906  ens|ENST00000593022,ens|ENST00000586288,ens|ENST00000587530,ens|ENST00000585768,ens|ENST00000590848,ens|ENST00000585745,ens|ENST00000591440,ens|ENST00000586928,ens|ENST00000588951,ens|ENST00000591702,ens|ENST00000429505,ens|ENST00000589531,ens|ENST00000358533,mRNA|AK125248,mRNA|AK290103
chr1    721751  722042  ens|ENST00000593022,ens|ENST00000586288,ens|ENST00000587530,ens|ENST00000585768,ens|ENST00000590848,ens|ENST00000585745,ens|ENST00000591440,ens|ENST00000586928,ens|ENST00000588951,ens|ENST00000591702,ens|ENST00000429505,ens|ENST00000589531,ens|ENST00000358533,mRNA|AK125248,mRNA|AK290103
chr1    752816  753135  ref|FAM87B,ref|NR_103536,ens|ENST00000326734,ens|ENST00000435300,mRNA|AK097327
chr1    761995  762375  ref|LINC00115,ref|NR_024321,ens|ENST00000473798,ens|ENST00000536430,mRNA|AK026292,mRNA|KJ901126,mRNA|BC017762
chr1    762180  762514  ref|LINC00115,ref|NR_024321,ens|ENST00000473798,ens|ENST00000536430,mRNA|AK026292,mRNA|KJ901126,mRNA|BC017762

Thanks for your help!

[Asaf]

What happens without grep ^all?

[DrDoom]

I've tried a "classic" bedtools coverage -hist on our institute's Galaxy instance on the off chance that something might be wrong with my version of bedtools and the output was also an empty file. So without the grep.

You don't see anything wrong with the structure of the bam or the bed?

Answer 1

[DrDoom]

Files and command were ok, problem was actually coming from the lack of memory as I was computing locally and not on a cluster. Even the much larger memory capacity of our Galaxy server wasn't enough (bam file was around 12GB).

I had to use the -sorted option as per the bedtools documentation:

If you are trying to compute coverage for very large files and are having trouble with excessive memory usage, please presort your data by chromosome and then by start position (e.g., sort -k1,1 -k2,2n in.bed > in.sorted.bed for BED files) and then use the -sorted option. This invokes a memory-efficient algorithm designed for large files.

As my files were already sorted, I just had to modify my original line:

bedtools coverage -hist -sorted -a exome_capture.bed -b mybam.bam | grep ^all > mybam.bam.hist.all.txt

I also increased the swap size just to be sure, and everything worked perfectly!

Too bad there was no error message or warning about the lack of memory...