bedtools output

[cmccabe]

I'm not sure what the output is (there are no headers) and can't find any documentation. Thanks.

 coverageBed -hist -abam IonXpress_001.bam -b LCH_IDT.bed | grep ^all > output.bam.hist.txt

output.bam.hist.txt (attached)

LCH_IDT.bed
chr12 112884064 112884217
chr12 112888106 112888331
chr12 112890983 112891206
chr12 112892352 112892499
chr12 112893738 112893882
chr12 112910732 112910859
chr12 112915439 112915549
chr12 112915645 112915834
chr12 112919862 112920024
chr12 112924263 112924448
chr12 112926231 112926329
chr12 112926812 112926994
chr12 112939932 112940075
chr12 112942483 112942583

[dpryan]

See here.

The second column is the depth, the third is the number of bases with that depth, the fourth is the total number of bases in your bed file (2188 bases) and the last column the percentage with that coverage.

[cmccabe]

So I am trying to understand the output.

There are 14 input regions in the bed and ~300 rows in the output. Is there a way to identify the reads at a baspair for each bait?

For example, of those 300 or so rows how many reads belong to the 1 bait in the bed, how many belong to the second, etc...

So according to the output there are 127 bases with 0 reads? Thanks.

I am assuming the all in the output is beacause of ^all.

[dpryan]

Try removing the grep for "all", you may also want the -d option, depending on your goals. Yes, there are 127 bases with no coverage.

[cmccabe]

 coverageBed -hist -abam -d IonXpress_001.bam -b LCH_IDT > output.bam.hist.txt

? Thanks.

[dpryan]

I expect you need -d before -abam, otherwise it'll look for a BAM file called "-d".

coverageBed -hist -d -abam IonXpress_001.bam -b LCH_IDT > output.bam.hist.txt

[cmccabe]

Attached is the out put of the command. Below are the first 3 lines, if I am reading it correctly it is saying that in bp 1 of chr12:112884064-112884217 there are 141 reads, in bp 2 of chr12:112884064-112884217 there are 140 reads. In bp 3 of chr12:112884064-112884217 there are 141 reads. Is there a command that will output the average depth per the 14 targets in the bed file? Thanks.

chr12 112884064 112884217 chr12:112884064-112884217
1 141
chr12 112884064 112884217 chr12:112884064-112884217
2 140
chr12 112884064 112884217 chr12:112884064-112884217
3 141

[dpryan]

Maybe, but you can also just use awk.

awk '{if(len==0){last=$4;total=$6;len=1;getline}if($4!=last){printf("%s\t%f\n", last, total/len);last=$4;total=$6;len=1}else{total+=$6;len+=1}}END{printf("%s\t%f\n", last, total/len)}' output.bam.hist.txt

BTW, you have an incorrect line ending in the middle of each line of that file. Have you been using windows again?

Anyway, here's the output of the awk command:

chr12:112884064-112884217	158.209150
chr12:112888106-112888331	220.533333
chr12:112890983-112891206	228.286996
chr12:112892352-112892499	182.102041
chr12:112893738-112893882	202.076389
chr12:112910732-112910859	0.000000
chr12:112915439-112915549	82.590909
chr12:112915645-112915834	217.269841
chr12:112919862-112920024	279.586420
chr12:112924263-112924448	452.535135
chr12:112926231-112926329	162.234694
chr12:112926812-112926994	189.131868
chr12:112939932-112940075	291.020979
chr12:112942483-112942583	160.510000

[cmccabe]

Can the awk be modified to include the average coverage of a bait as well as the amount of reads? I am aligning to hg19. Thank you .

chr12:112884064-112884217	158.209150      36x
chr12:112888106-112888331	220.533333      40x
chr12:112890983-112891206	228.286996      42x