Bug Bulletin: The GenomeLocPArser error in SplitNCigarReads has been fixed; if you encounter it, use the latest nightly build.

DuplicateReadFilter

SystemSystem Posts: 226Administrator admin
edited July 2012 in Tool Bulletin

A new tool has been released!

Check out the documentation at DuplicateReadFilter.

Comments

  • ymwymw Posts: 9Member

    There seem to be no documentation to explain how to use DuplicateReadFilter! Where can I find such documentation?

    Thanks, Chih

  • Geraldine_VdAuweraGeraldine_VdAuwera Posts: 6,192Administrator, GATK Developer admin

    Hi Chih,

    We currently don't have specific documentation for some read filters such as this one, sorry. We will fix that in the near future.

    Basically, all you really need to know about this filter is that it filters out reads that have been marked as duplicates. To use it, you pass it like any other read filter, using -rf DuplicateRead (always omit the "Filter" part in the command line).

    Note that all Locus Walkers (like the Unified Genotyper) apply this filter by default, so for most GATK tools you don't need to specify that you want to use it.

    Geraldine Van der Auwera, PhD

  • drdboehmdrdboehm Posts: 1Member

    Hi Geraldine, first of all, I would like to say thx for the phantastic support on the GATK forum and at the documentation. Mostly in the past it was not neccesary to write an active comment. Any questions and answers so far were included by default. But now, I would like to answer, whether ist is possible to skip a Readfilter by default? At the moment I'm dealing with the human X chromosome and I was wondering, why I do not have a SNV call (neither with UnifiedGenotyper nor HalpotypeCaller) located in the interval chrX:1-2,699,705 in my exome data processed using bwa and the recommended GATK pipeline. I tried several settings using UnifiedGenotyper, try to get so less stringent I ever could... Nevertheless, I saw the reads and the SNV in IGV, but not called in my RAW-SNV. At the end, I became aware, that perhaps I'm dealing with the PAR problem on the X- and Y-chromosome. The questions are: 1. How to discard the DuplicateReadFilter in the analysis when duplicates are marked already? 2. Could be, dealing with the PAR is a more general problem that just switch off the DuplicateReadFilter!? Do you have any recommendations to call reliables datasets on that reagion ?

    Many thanks

    Detlef

  • CarneiroCarneiro Posts: 275Administrator, GATK Developer admin

    I'm not sure what you mean by the PAR problem, do you mean you would like to call ChrX as haploid? If so, you can use the UG with -ploidy 1 and -L X.

    There is no way to deactivate the read filter if it's on by default in the walker. But then, why would you want to use duplicate reads in your analysis?

  • mmajmmaj Posts: 6Member
    edited May 2013

    Hi GATK team,

    The ability to disable "default" readfilters, eg. duplicate read filters for HC, would be very useful for amplicon-sequencing projects! For instance, we're using Agilents Haloplex, which per design requires that duplicate reads are not filtered out. Its not really clear to me why you want to decide what filters people would prefer to apply to their data anyways. Would it be alot of work to add the possibility to discard "default" read filters? Particularly for HC, which is now way too stringent in its filtering of reads from our datasets. Its not a problem in version 2.4.9, only in 2.5.2

    Thanks for a really great toolset, btw!

    Post edited by mmaj on
  • CarneiroCarneiro Posts: 275Administrator, GATK Developer admin

    Why would you mark duplicates on Agilent Haloplex if you don't want them excluded?

  • mmajmmaj Posts: 6Member

    I don't! If Picard MarkDuplicates is applied, we remove like 95% of aligned haloplex reads (a little less with Samtools rmdup) , so clearly, no mark duplicates :) But the default read filters applied to HC in version 2.5.x severely reduces the number of reads used by HC for calling, especially short deletions. HC is ofc run without downsampling. Comparison of the bam file used as input for HC, with the one emitted by HC clearly reveals that alot of - what appears to be - perfectly good reads are disregarded by HC. I don't know which of the default fitlers that are responsible for this. I understand that applying default filters is meant to simplify the use of different GATK tools, but it would still be nice to be able to inactivate these filters as different sequencing projects would most likely require different filters anyways. At the very least, it would give users an opportunity to troubleshoot and understand how individual read filters affect their datasets.

  • CarneiroCarneiro Posts: 275Administrator, GATK Developer admin

    Okay, so the duplicate read filter is not your problem, because if the reads are not marked as duplicates by Picard or Samtools, there is no filtering on them.

    While I do agree with you that it would be nice to have an option to disable the filters for user troubleshooting, I'm afraid this is not a very simple thing to do in a generic way. However, it is very simple for you to edit the source code and remove the read filters from any walker, for local testing.

    The HaplotypeCaller applies the following filters by default:

    • DuplicateReadFilter
    • FailsVendorQualityCheckFilter
    • NotPrimaryAlignmentFilter
    • MalformedReadFilter
    • MappingQualityUnavailableFilter
    • HCMappingQualityFilter
    • UnmappedReadFilter

    If you want to find out which filter is eliminating your reads, you can run PrintReads with -rf for each one of these filters on a small sample of your bam file. If indeed your reads are getting filtered out by one of these filters, you'll see right away.

    Now as a cautionary note, if your reads are getting filtered by these filters, that would be an indication of bad processing or bad sequencing. Most of these filters are basic quality control filters, with the only exception of the HCMappingQuality filter which may be slightly more aggressive than the others (eliminating reads under MQ 20). Good luck!

  • mmajmmaj Posts: 6Member

    Thanks, Carneiro! I'll definitely try the PrintReads suggestion!

  • mmajmmaj Posts: 6Member

    Hi again. I've tested the read filters by PrintReads, and by comparing the AD emitted by HC and UG for specific SNPs, all run on GATK2.6.4. Using PrintReads with each of the default readfilters for HC (as Carneiro list above), does not result in any reads being filtered out for a small bam region encompassing a single SNP.

    Intriguingly, the AD field of the single SNP clearly shows that HC disregard the majority of reads (AD=12,4, DP=16 out of a total of 147 for that particular region), whereas UnifiedGenotyper (UG), which by default applies the same (except one) read filters, includes almost all of these reads for calling that particular SNP (AD=100,44, DP=144 out of 147 reads).

    The HC and UG commands were run with all other parameters default:

    UG: java -jar /gatk264/GenomeAnalysisTK.jar -T UnifiedGenotyper --dbsnp /hg19/vcf/dbsnp_137.hg19.vcf -I Val26.aln.bam -o val26.GATK264.UG.vcf -R /hg19/hg19.fasta -L /interval.files/ROI.bed -nct 8

    HC: java -jar /gatk264/GenomeAnalysisTK.jar -T HaplotypeCaller --dbsnp /hg19/vcf/dbsnp_137.hg19.vcf -I Val26.aln.bam -o val26.GATK264.vcf -R /hg19/hg19.fasta -L /interval.files/ROI.bed -nct 8

    The only difference in the read filters applied by HC and UG is, from what I can see, the use of HCMappingQualityFilter with HC and BadMateFilter with UG. However, neither of these filters resulted in filtering af reads from the test bam file (containing the above mentioned snp +100 bp) when running PrintReads.

    Is there additional walker-specific internal read filters, that could be responsible for this AD discrepancy between UG and HC for calling the same SNP on the same bam? I'd be happy to upload whatever you might need...

    On another note: I'm happy to report that HC - at least in our hands - hasin fact become so sensitive, it cries at the movies... It's so beautiful to witness, it brings water to my eyes :)

  • Geraldine_VdAuweraGeraldine_VdAuwera Posts: 6,192Administrator, GATK Developer admin

    Hi @mmaj, could you upload a small test file so I can reproduce this locally? Instructions are here: http://www.broadinstitute.org/gatk/guide/article?id=1894

    Geraldine Van der Auwera, PhD

  • Geraldine_VdAuweraGeraldine_VdAuwera Posts: 6,192Administrator, GATK Developer admin

    Also,

    On another note: I'm happy to report that HC - at least in our hands - hasin fact become so sensitive, it cries at the movies... It's so beautiful to witness, it brings water to my eyes :)

    Glad you're feeling the HC's emotional high too :))

    Geraldine Van der Auwera, PhD

Sign In or Register to comment.