GATK - DepthOfCoverage

[giverny]

Hi

I try to use DepthOfCovareage. My command line is

java -jar GenomeAnalysisTK.jar
-T DepthOfCoverage
-L 5:125921513-125923713
-o test
-I input.bam
-R hg18.fa

For the firsts lines I obtain these results :
Locus Total_Depth Average_Depth_sample Depth_for_N/A
5:125921513 253 253.00 253
5:125921514 285 285.00 285
5:125921515 292 292.00 292
5:125921516 310 310.00 310
5:125921517 319 319.00 319
5:125921518 320 320.00 320
(...)

but when I'm looking for the reads aligned at the position 125921513 (for example ... first line of my output file from GATK) via this command line :

samtools view input.bam | cut -f 3-4 | grep 125921513 | wc -l
I obtain 7214 data (and all of them are within chromosome 5!)

When I do the same position 125921514, I obtain 41 reads.

I'd like to understand the difference between my interpretation of the bam file and the interpretation done by GATK. Which one is good? If it's GATK what about my command line output? What does GATK use?

THanks a lot for your time

[Heisman]

Are you sure nothing is wrong with your files prior to running DepthOfCoverage? I feel it's unlikely to have 253x coverage for one base and then have it jump to 285x coverage for the next base, although I don't know how the data was generated.

[giverny]

Hi Heisman.
It seems my command line was wrong. Someone told me that "
"Your command is finding all records in the SAM file with a start position of 125921513. (The depth at that location is a combination of the reads that start there, and those that start prior to that location, yet span it).

In addition, there may be a number of pieces of information about each record such that it does not count as 'depth' at a particular site. For instance, unmapped reads can still have a start location (from the SAM spec: "Bit 0x4 [of the flag field] is the only reliable place to tell whether the segment is unmapped. If 0x4 is set, no
assumptions can be made about RNAME, POS, CIGAR, MAPQ..."). Similarly, the records could be marked as optical duplicates, or as having failed vendor quality checks; in both cases they should not be counted as a part of coverage. "
so it make sense I think...
Thanks for your time!