coverage of narrow peaks on genomic region using bedtools
1
0
Entering edit mode
2.3 years ago
rthapa ▴ 50

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

ChIP-Seq • 555 views
1
Entering edit mode
2.3 years ago
Brice Sarver ★ 3.7k

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

0
Entering edit mode

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?

0
Entering edit mode

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.