SEQanswers

Go Back   SEQanswers > Bioinformatics > Bioinformatics



Similar Threads
Thread Thread Starter Forum Replies Last Post
Bedtools fastafrombed output lwhitmore Bioinformatics 5 06-20-2014 10:17 AM
BEDtools intersect output is BED instead of BAM syfo Bioinformatics 1 12-18-2012 05:26 AM
bedtools cmccabe General 2 10-31-2012 12:54 PM
intersectBed (BEDtools) generating empty output file palc Bioinformatics 1 08-28-2012 09:36 AM
samtools vs. bedtools bair Bioinformatics 2 02-13-2012 08:57 AM

Reply
 
Thread Tools
Old 09-17-2014, 09:45 AM   #1
cmccabe
Senior Member
 
Location: chicago

Join Date: Jul 2012
Posts: 354
Default bedtools output

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

Code:
 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
Attached Files
File Type: txt output.bam.hist.txt (7.6 KB, 9 views)
cmccabe is offline   Reply With Quote
Old 09-17-2014, 10:03 AM   #2
dpryan
Devon Ryan
 
Location: Freiburg, Germany

Join Date: Jul 2011
Posts: 3,480
Default

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.
dpryan is offline   Reply With Quote
Old 09-17-2014, 01:39 PM   #3
cmccabe
Senior Member
 
Location: chicago

Join Date: Jul 2012
Posts: 354
Default

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.
cmccabe is offline   Reply With Quote
Old 09-17-2014, 01:44 PM   #4
dpryan
Devon Ryan
 
Location: Freiburg, Germany

Join Date: Jul 2011
Posts: 3,480
Default

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.
dpryan is offline   Reply With Quote
Old 09-18-2014, 07:19 AM   #5
cmccabe
Senior Member
 
Location: chicago

Join Date: Jul 2012
Posts: 354
Default

Code:
 coverageBed -hist -abam -d IonXpress_001.bam -b LCH_IDT > output.bam.hist.txt
? Thanks.
cmccabe is offline   Reply With Quote
Old 09-18-2014, 07:23 AM   #6
dpryan
Devon Ryan
 
Location: Freiburg, Germany

Join Date: Jul 2011
Posts: 3,480
Default

I expect you need -d before -abam, otherwise it'll look for a BAM file called "-d".
Code:
coverageBed -hist -d -abam IonXpress_001.bam -b LCH_IDT > output.bam.hist.txt
dpryan is offline   Reply With Quote
Old 09-18-2014, 07:54 AM   #7
cmccabe
Senior Member
 
Location: chicago

Join Date: Jul 2012
Posts: 354
Default

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
Attached Files
File Type: zip output.bam.hist.zip (7.7 KB, 3 views)
cmccabe is offline   Reply With Quote
Old 09-18-2014, 08:26 AM   #8
dpryan
Devon Ryan
 
Location: Freiburg, Germany

Join Date: Jul 2011
Posts: 3,480
Default

Maybe, but you can also just use awk.

Code:
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:
Code:
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
dpryan is offline   Reply With Quote
Old 06-18-2015, 10:17 AM   #9
cmccabe
Senior Member
 
Location: chicago

Join Date: Jul 2012
Posts: 354
Default

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 .

Code:
chr12:112884064-112884217	158.209150      36x
chr12:112888106-112888331	220.533333      40x
chr12:112890983-112891206	228.286996      42x
cmccabe is offline   Reply With Quote
Reply

Thread Tools

Posting Rules
You may not post new threads
You may not post replies
You may not post attachments
You may not edit your posts

BB code is On
Smilies are On
[IMG] code is On
HTML code is Off




All times are GMT -8. The time now is 06:29 PM.


Powered by vBulletin® Version 3.8.9
Copyright ©2000 - 2021, vBulletin Solutions, Inc.
Single Sign On provided by vBSSO