*// Edit to make the post more clear (Mapping done via Bowtie2). My problem is that when counting Exome Coverage via coverageBed gives different results than via genomeCoverageBed. So I'm not sure if I'm doing something wrong, or which of the 2 methods is correct.
1) My first step is to build an .bed file of my Illumina Paired-End reads, returning the positions that only fall in targeted exon regions. I'm doing that via intersectBed -a [data.bed] -b [illuminaexonregions.bed].
2) My next step is to calculate the coverage of my new datafile via coverageBed -a [newdata.bed] -b [illuminaexonregions.bed]. I calculated some statistics:
Number of exons 214126 with a total length of 45326818
Number of matched nucleotides 10993449.0
Nucleotides/Length*100 24.253740909 % Coverage.
3) The next step was to calculate the coverage of my new datafile via genomeCoverageBed -i [newdata.bed] -g [genome.txt] -d awk '$3>0 {print $1"\t"$2"\t"$3}'. I calculated some statistics:
Number of exons 214126 with a total length of 45326818
Number of matched nucleotides 10576907.0
Nucleotides/Length*100 23.3347661863 % Coverage.
Somehow there's a difference in matched nucleotides, which I can't explain. What am I doing wrong?
3 answers
If you already have a .bed file with your targeted regions (in this case, the exome), I think the easiest solution to to use CalculateHsMetrics in Picard (which will also provide additional statistics beyond just on-target coverage)
http://picard.sourceforge.net/command-line-overview.shtml#CalculateHsMetrics
Picard has been essentially folded in to the GATK. The new documentation is here: https://broadinstitute.github.io/picard/command-line-overview.html
Assuming you are using the same files for both commands (from your examples it wasn't clear) it looks like you are only dealing with a few bases. IIRC genomeCoverageBed only counts the bases on reads that lie within the targeted regions. IntersectBed returns all reads that intersect with the targeted regions in the bedfile, and then you are counting all of those bases. So it would include the bases that lie within intronic regions, etc as long as part of the read overlaps the exons.
The first step is to intersectBed the derived data against the Illumina Targeted Exome file. In that case all the intronic regions are filtered out right? I thought only the regions that fall in the exon positions are returned, and in that case no intronic regions should be in my new file?
When you run intersectBed it gives you all reads that intersect that bed file. Many of the reads though can span the intron/exon boundary. Intersectbed doesn't just give you the bases that intersect, at least I don't think so. I would use HSMetrics as suggested above, or one of the GATK tools for calculating coverage as you can get much more detailed metrics.
Additional information;
intersectBed -a chromFA.bed -b nexterarapidcapture_exome_targetedregions_v1.2.bed > new_exome_regions_intersectBed.bed
coverageBed -a new_exome_regions_intersectBed.bed -b nexterarapidcapture_exome_targetedregions_v1.2.bed > new_exome_regions_coverageBed.bed
genomeCoverageBed -i new_exome_regions_intersectBed.bed -g genome.txt -d | awk '$3>0 {print $1"\t"$2"\t"$3}' > new_exome_regions_genomeCoverageBed.bed
genomeCoverageBed
Number of exons 214126 with a total length of 45326818
Number of bases 10576907.0
Bases/Length*100 23.3347661863
coverageBed
Number of exons 214126 with a total length of 45326818
Number of bases 10993449.0
Bases/Length*100 24.253740909
Log in to answer this question.
-
3'tsRNA alignment
written by aranyak111 •
-
How has StringTie caculated the transcript coverage?
written by ddzhangzz
-
Use Bedtools to extract stranded read depth from BAM files
written by gtho123
-
Strategy for RBP sequence logo alignments to other RNA sequences
written by syntonic.c •
-
Snp Density
written by int11ap1
-
How Much Of The Genome Will Remain Un-Sequenced At A Given Coverage?
written by kerpalha •
-
Getting The Average Coverage From The Coverage Counts At Each Depth.
written by Jordan
-
Intersectbed/Coveragebed -Split Purify Exon?
written by Puriney
-
Coveragebed, Depth/Breadth Of Coverage
written by Newbionf •
-
What Metrics Are Best To Describe The "Coverage" Of Rna-Seq Data?
written by Bio_X2Y
