coverage of narrow peaks on genomic region using bedtools

Question

[rthapa]

Hi,

I am trying to estimate the coverage of narrow peaks in genomic region. I have a bed file of both genomic coordinates and narrow peak output. When I use bedtools to find the coverage, all the results are 0 but by mere scanning of the peak region and genomic region, I can see some of the peak regions are inside the genomic region. Is there any suggestion if I am doing anything wrong here?

bedtools coverage -a generegion.bed -b narrowpeak.bed > out.txt

Genomic region

    Chr1    1951    2616    Sobic.001G000100
    Chr1    11180   14899   Sobic.001G000200
    Chr1    23399   24152   Sobic.001G000300
    Chr1    22391   42443   Sobic.001G000400
    Chr1    52891   53594   Sobic.001G000501
    Chr1    53781   63305   Sobic.001G000700
    Chr1    62892   69306   Sobic.001G000800
    Chr1    79159   81636   Sobic.001G000900
    Chr1    81932   83350   Sobic.001G001000
    Chr1    103029  104803  Sobic.001G001066

Narrowpeak region

 Chr01  1629    3099
    Chr01   4084    4355
    Chr01   210952  213035
    Chr01   217139  219527
    Chr01   261167  261473
    Chr01   308533  311133
    Chr01   328267  330498
    Chr01   398709  399366
    Chr01   419666  420161
    Chr01   428019  428806
    Chr01   429432  429909
    Chr01   430970  433767
    Chr01   435036  435796
    Chr01   459030  459301
    Chr01   478976  479992
    Chr01   480778  481057
    Chr01   507360  507912
    Chr01   508132  509148

Output

Chr1    1951    2616    Sobic.001G000100    0   0   665 0.0000000
Chr1    11180   14899   Sobic.001G000200    0   0   3719    0.0000000
Chr1    23399   24152   Sobic.001G000300    0   0   753 0.0000000
Chr1    22391   42443   Sobic.001G000400    0   0   20052   0.0000000
Chr1    52891   53594   Sobic.001G000501    0   0   703 0.0000000
Chr1    53781   63305   Sobic.001G000700    0   0   9524    0.0000000
Chr1    62892   69306   Sobic.001G000800    0   0   6414    0.0000000
Chr1    79159   81636   Sobic.001G000900    0   0   2477    0.0000000
Chr1    81932   83350   Sobic.001G001000    0   0   1418    0.0000000
Chr1    103029  104803  Sobic.001G001066    0   0   1774    0.0000000
Chr1    106847  107376  Sobic.001G001132    0   0   529 0.0000000
Chr1    109968  112306  Sobic.001G001200    0   0   2338    0.0000000
Chr1    112766  116048  Sobic.001G001300    0   0   3282    0.0000000
Chr1    116116  136690  Sobic.001G001400    0   0   20574   0.0000000

Answer 1

[Brice Sarver]

Your first BED has chromosomes in the format ChrX, whereas your second BED has chromosomes in the format Chr0X.

[rthapa]

Thank you! that was a silly mistake. For some of the genes, I get the coverage as 1. But when I checked the start and end position of both peak output file and gene region file, I don't see 100 % coverage. Is the coverage not calculated as peak length/ gene length?

[Brice Sarver]

From the bedtools coverage manual page:

After each interval in A, bedtools coverage will report:

• The number of features in B that overlapped (by at least one base pair) in the A interval.
• The number of bases in A that had non-zero coverage from features in B.
• The length of the entry in A.
• The fraction of bases in A that had non-zero coverage from features in B.

You can also dial-in the type of output you expect to see.