Bedtools - Understand coverageBed output
neetika nath
Even after going through many discussion to understand coverageBed output, I failed to interpret my output file.
At the moment, I am using bedtools2-2.25 on my ubuntu machine.
After using the following command
/bedtools/bedtools2/bin/coverageBed -hist -abam alignment/output_bowtie.fastq.bam -b bedfiles/Regions.bed > coverage1.txt
I got the following output
chrM 12375 12538 M03129:35:000000000-AD694:1:1112:13735:4599 42 - 12375 12538 0,0,0 1 163, 0,0 163 163 1.0000000
chrM 12382 12538 M03129:35:000000000-AD694:1:1105:17125:4856 42 - 12382 12538 0,0,0 1 156, 0, 0 156 156 1.0000000
chr1 101375 101574 M03129:35:000000000-AD694:1:1101:28243:17684 12 - 101375 101574 0,0,0 1 199,0,0 199 199 1.0000000
chr1 129427 129486 M03129:35:000000000-AD694:1:1103:24032:4006 1 - 129427 129486 0,0,0 1 59,0,0 59 59 1.0000000
chr1 457320 457409 M03129:35:000000000-AD694:1:2118:29453:13272 1 + 457320 457409 0,0,0 1 89,0,0 89 89 1.0000000
chr1 655784 655843 M03129:35:000000000-AD694:1:2108:17677:16517 1 - 655784 655843 0,0,0 1 59, 0,0 59 59 1.0000000
My understanding of the output of column is:
column 1 == chromosome
column 2 == start position
column 3 == stop position
column 4 == ??
column 5 == ??
column 6 == strand
column 7 == ??
column 8 == ??
column 9 == ??
column 10 == ??
column 11 == ??
column 12 == ??
column 13 == ??
column 14 == should be column 12/column 13
It will be highly appreciated for any suggestion what other column means.
In particular, I want to get values of depth and breadth coverage using bedtools so that I can keep a track of information lose during post alignment process. With using earlier version of bedtools coverage I could extract information for breadth of coverage however I could not find the option to extract depth as well. Any suggestion would be helpful.
Thanks in advance!
Luke Goodsell
Hi Neetika,
As a fellow user of Bedtools, I’ll try to help.
With BAM input, coverageBed first converts the input file to bed12 format (similar to bamToBed with the -bed12 option). See https://genome.ucsc.edu/FAQ/FAQformat.html#format1 for the bed12 format spec.
coverageBed adds an extra four fields to the data that you provide. I.e., if you provide 12-column input, it will output 16 columns. Your output actually has 16 columns, but columns 11 and 12 have commas at the end, making it look like they are joined to the next column.
This -hist option changes the final four fields from the coverageBed default, and will cause 1 or more lines of output per entry in your -a file.
column 1 = chromosome
column 2 = start
column 3 = end
column 4 = read name (first field of BAM file)
column 5 = mapping quality (fifth field of BAM file)
column 6 = strand
column 7 = thickStart
column 8 = thickEnd
column 9 = itemRGB (three numbers comma-delimited numbers)
column 10 = blockCount
column 11 = blockSizes (has a trailing comma)
column 12 = blockStarts (has a trailing comma)
column 13 = coverage depth
column 14 = # bases at depth
column 15 = size of A
column 16 = % of A at depth
For the lines that you provided, none of the features has coverage at a depth greater than 0. I.e.: they do not overlap the intervals in your Regions.bed file.
Here is your output with the column numbers:
# 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16
chrM 12375 12538 M03129:35:000000000-AD694:1:1112:13735:4599 42 - 12375 12538 0,0,0 1 163, 0, 0 163 163 1.0000000
chrM 12382 12538 M03129:35:000000000-AD694:1:1105:17125:4856 42 - 12382 12538 0,0,0 1 156, 0, 0 156 156 1.0000000
chr1 101375 101574 M03129:35:000000000-AD694:1:1101:28243:17684 12 - 101375 101574 0,0,0 1 199, 0, 0 199 199 1.0000000
chr1 129427 129486 M03129:35:000000000-AD694:1:1103:24032:4006 1 - 129427 129486 0,0,0 1 59, 0, 0 59 59 1.0000000
chr1 457320 457409 M03129:35:000000000-AD694:1:2118:29453:13272 1 + 457320 457409 0,0,0 1 89, 0, 0 89 89 1.0000000
chr1 655784 655843 M03129:35:000000000-AD694:1:2108:17677:16517 1 - 655784 655843 0,0,0 1 59, 0, 0 59 59 1.0000000
Kind regards,
Luke
neetika nath
Thank you for your reply. This is very useful. Could you also help me to understand bedtools v2.17.0 coverage output. I get lot of information from bedtools manual however I am unable to understand column 5. Here is what I understood:
Col 1 -- chromosome
Col 2 -- start
Col 3 -- Stop
Col 4 -- Name
Col 5 -- ?? (my guess is the depth of the reads or N according to the bedtools manual)
Col 6 -- breadth coverage
Col 7 -- breadth coverage
Col 8 -- %age of col 6/7
1 2 3 4 5 6 7 8
chr1 92733419 92733609 dN 661 190 190 1.0000000
Please correct me if I have misunderstood the column. I am also unable to find a good tutorial for the latest bedtools, could you please direct me to the link.
Thanks in advance!
Luke Goodsell
Hi Neetika,
> Could you also help me to understand bedtools v2.17.0 coverage output.
If possible, I would highly advise you to only use the latest version. Be aware of this notice on the coverageBed documentation page:
> As of version 2.24.0, the coverage tool has changed such that the coverage is computed for the A file, not the B file.
I assume, then, that your Regions.bed file defines the regions d1 and dN and has 4 columns? Columns 5 to 8 will therefore be depth, number of bases at depth, size of B feature, and proportion of B at depth. I don’t know why the depth of coverage is so high, though. Perhaps you could provide your input files (or a sample of them, if they’re large) and the command you run, and I (or someone else) might be able to help?
> I am also unable to find a good tutorial for the latest bedtools, could you please direct me to the link.
The best source of documentation that I know of is the Bedtools documentation site at http://bedtools.readthedocs.org/en/latest/index.html. It’s also good to be familiar with the BED format specification I linked to previously.
Kind regards,
Luke
--
You received this message because you are subscribed to the Google Groups "bedtools-discuss" group.
To unsubscribe from this group and stop receiving emails from it, send an email to
bedtools-discu...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.
apoo...@gmail.com
| compid | start | stop ? | ? | lenth ? | ||
| comp41304_c0_seq1 | 1 | 545 | 9 | 488 | 544 | 0.897059 |
| comp42204_c0_seq1 | 1 | 668 | 24 | 649 | 667 | 0.973014 |
| comp48854_c0_seq1 | 1 | 723 | 63 | 719 | 722 | 0.995845 |
| comp52231_c1_seq1 | 1 | 833 | 138 | 795 | 832 | 0.955529 |
| comp55490_c0_seq2 | 1 | 933 | 174 | 909 | 932 | 0.975322 |
| comp58472_c0_seq1 | 1 | 2014 | 266 | 2013 | 2013 | 1 |
| comp59249_c0_seq1 | 1 | 1708 | 1220 | 1695 | 1707 | 0.99297 |
| comp59249_c0_seq2 | 1 | 848 | 240 | 847 | 847 | 1 |
| comp60037_c1_seq1 | 1 | 1162 | 358 | 1158 | 1161 | 0.997416 |
| comp60182_c1_seq5 | 1 | 1230 | 685 | 1229 | 1229 | 1 |
| comp60182_c1_seq9 | 1 | 562 | 121 | 556 | 561 | 0.991087 |
| comp60261_c1_seq1 | 1 | 1225 | 1357 | 1224 | 1224 | 1 |
| comp60261_c3_seq1 | 1 | 1330 | 1710 | 1325 | 1329 | 0.99699 |
| comp61668_c0_seq1 | 1 | 1863 | 2066 | 1862 | 1862 | 1 |
| comp64573_c1_seq5 | 1 | 2177 | 321 | 2140 | 2176 | 0.983456 |
| comp64793_c1_seq11 | 1 | 1659 | 564 | 1578 | 1658 | 0.951749 |
| comp64793_c1_seq13 | 1 | 1400 | 608 | 1383 | 1399 | 0.988563 |
