GATK variant calling on uniquely mapped reads?

[pravee1216]

Hello All,

Does GATK use uniquely mapped or reads that satisfy the parameter criteria for calling variants by default? How do I instruct GATK to use only uniquely mapped reads to variant calling?

Any comments or experience on this

Thanks

Raj

[neha]

Hi ,

I think if you want to call variants only on uniquely mapped reads you can first remove duplicates using picardtools and then call variants using GATK.

Also you can do per base score recalibration with GATK and then call variants so that there is more confidence in calling variants.

Just sharing my thoughts.

Neha

[pravee1216]

I don't think so, because duplicate removal doesn't mean that reads mapped to multiple locations are eliminated. Both are two different concepts. why I am saying this because even if you dedup reads, still you would see reads with mapQ=0 (means reads mapped to more than one location in reference). I think, applying a mapQ cut-off would *also* be required to get uniquely mapped reads. What do you say?

[brofallon]

I believe by default that reads with zero mapping quality are ignored by the UnifiedGenotyper. If you're curious, you could filter out such reads using the PrintReads tool and the MappingQualityZero filter (see the first example here: http://www.broadinstitute.org/gsa/ga...adsWalker.html

[neha]

I don't think so, because duplicate removal doesn't mean that reads mapped to multiple locations are eliminated. Both are two different concepts. why I am saying this because even if you dedup reads, still you would see reads with mapQ=0 (means reads mapped to more than one location in reference). I think, applying a mapQ cut-off would *also* be required to get uniquely mapped reads. What do you say?

Maping quality does not tell whether reads are uniquely mapped or duplicate. It signifies the hypothesis of having confidence that read mapping to particular genomic position is correct. Higher the mapping quality value more is the confidence.

[jflowers]

Maping quality does not tell whether reads are uniquely mapped or duplicate. It signifies the hypothesis of having confidence that read mapping to particular genomic position is correct. Higher the mapping quality value more is the confidence.

That is not necessarily true. If you have aligned your data with BWA then the recommended way to get "unique" reads (by Heng Li) is to filter on map quality. See the samtools FAQ under the section "I want to get `unique' alignments from SAM/BAM":

http://sourceforge.net/apps/mediawiki/samtools/index.php?title=SAM_FAQ

In cases where reads are reliable (ie map "uniquely") mapq will be set to > 0 in BWA. I suspect this is specific to the aligner.