The tool is at: http://bedtools.readthedocs.org/en/latest/content/tools/coverage.html
I want to get base coverage for specific regions (I'll provide a bed file). However, the output is not in a standard BedGraph format that I can export into R for visualisation.
Q: How do I draw density plot of coverage over the regions from the outputs generated by the coverage tool?
Q: Alternatively, how to convert the non-standard outputs to a proper BedGraph? For example, the Sushi R package takes a bedgraph input file.
3 answers
You should also consider bedtools genomecov, with the -bg option (bg being bedgraph)
http://bedtools.readthedocs.org/en/latest/content/tools/genomecov.html
Edit: After reading comments, and after some experimenting, what I think you want to do is:
bedtools intersect -a myBam.bam -b myRegions.bed > intersected.bam
bedtools genomecov -trackline -bg -ibam intersected.bam > intersected.bedgraph
Bedtools intersect will give only bam reads that overlap your regions (in myRegions.bed)
Then genomecov will draw a bedgraph file (coverage) only for the read depth for your filtered regions.
•
written
by
Joseph Pearson
Fair enough, I didn't catch the subtlety in your question, but James Ashmore did.
I think this might be what you want to do - genomecov and then pull out the regions you want using bedtools intersect. But kind of hard to tell.
Get your output into a standard five- or more columned BED file and then use awk:
$ awk '{ print $1"\t"$2"\t"$3"\t"$5 }' foo.bed > foo.bedgraph
@ Alex: How can I convert my .sam file to standard .bed file? I have tried sam2bed < input.sam > output.bed it gives me output file but the column 5 always contains the same value (which is 99). According to UCSC bed format this is score column and to me this look strange that for every region the score is same. Then by using your mentioned command I obtain the .bedgraph file which I later use for ChromImpute. Problem is that when ChromeImpute convert these raw signal into certain sized bin signal (25 bp) then for every region I get the 0 signal.. I have tried many ways (using Pyisco, make wig files, genomeCoverageBed etc..) but I am unable to solve this problem. Can you help me with this one? Thanks.
To generate bedGraph output from Bedtools coverage command you need to specify the -counts flag, for example:
bedtools coverage -a sample.bed -b sample.bam -counts > sample.bedgraph
Strictly speaking each row of a bedGraph file contains the chromosome name, the start position, the end position and then some value, which in your case would be the coverage, therefore you will want to keep only these fields, so if your had a bed file with six columns, you would want to keep only the first 3 columns (chrom, start, end) and the last column (coverage):
bedtools coverage -a sample.bed -b sample.bam -counts | cut -f1,2,3,7 > sample.bedgraph
It's not really clear to me from the original question that this is what they want. This would give one value for each genome region. I suspect they want "wiggly" values for each genome region. But could be wrong.
Log in to answer this question.
-
How to obtain zero-based coordinates read depth using bedtools coverage for a specific region?
written by vladimiralan •
-
find coverage of specific exomes from a bam file
written by lait
-
Is BED file generated by bametobed sorted ?
written by pacman
-
Include values of zero in TSS metaplot?
written by goodez
-
Retrieving specific gene location
written by akylzhanov.r •
-
BEDtools coverage output
written by novice
-
how to merge overlapping intervals in a bedgraph file
written by Fidel
-
bedtools genomecov .genome file usage
written by kaston •
-
Bed File Of Mapq Sliding Window On A Bam File?
written by 14134125465346445
-
Plotting Read Densities And Gene Tracks In R?
written by user
