Test-drive the GATK tools and Best Practices pipelines on Terra
Check out this blog post to learn how you can get started with GATK and try out the pipelines in preconfigured workspaces (with a user-friendly interface!) without having to install anything.
FindCoveredIntervals -minBQ -minMQ -cov
shahr2
Member ✭
Hello Admin's,
When running FindCoveredIntervals for GATK version 3.3-0:
java -Xmx20g -jar GenomeAnalysisTK.jar -T FindCoveredIntervals -R Homo_sapiens_assembly19.fasta -I s_C_000036_T001_d.bam -minBQ 20 -minMQ 20 -cov 20 -o s_C_000036_T001_d.list -rf FailsVendorQualityCheck -rf BadMate -rf UnmappedRead -rf BadCigar
When we open the interval in IGV (file attached) you can observe that the coverage at the start and end of the interval is lower than the set threshold of 20.
If you don't mind explaining as to why that would be the case and how can we get intervals that are satisfy the above criteria that will be great.
Best,
Ronak
Tagged:
Answers
Hi Ronak (@shahr2),
My best guess is that coverage is some average coverage or smoothed, e.g. gaussian, coverage over a sliding window. This sliding window can range from the minimum region size to the maximum region size:
Would it make sense that the covered region includes at least one site with 20 coverage within a given window size?
If you need a more literal coverage profile, then I recommend IGV's Tools>igvtools>Count feature, which produces a binary TDF format coverage file. You'll see where you can adjust the window size. You can then convert the binary TDF to bedGraph format with
igvtools tdftobedgraph.Thank you for your reply on this @shlee.
Our interest lie not in finding the coverage profile but in finding all the interval that support that minimum coverage.
With regards to your question:
Yes in the sliding approach where you would have at least one site with 20X coverage.
But here we are getting this padding around that interval which my concern. But your explanation for taking an average coverage make sense in this case.
Not sure if using the igvtools will work directly but thanks for the suggestion.
@shahr2,
We have a more sophisticated coverage tool, [DepthOfCoverage] (https://software.broadinstitute.org/gatk/documentation/tooldocs/current/org_broadinstitute_gatk_tools_walkers_coverage_DepthOfCoverage.php) that can emit depth of coverage at each base. An example of how to interpret the results from this tool is in this thread. Good luck.
Hi @shahr2,
One of our team points out that perhaps FindCoveredIntervals includes duplicate reads in its count that your IGV setting is filtering. Can you check your IGV View>Preferences>Alignments to see if you are filtering duplicate reads as well as any other type of reads that the tool may count towards coverage?
Yes our IGV does filter duplicate reads when it displays.
But I am not sure if that statement is correct of duplicate reads being included, as in the usage it seems that filter is applied by default as it states in the Read filter section
https://software.broadinstitute.org/gatk/documentation/tooldocs/current/org_broadinstitute_gatk_tools_walkers_diagnostics_FindCoveredIntervals.php
If in the older version of GATK this filter was not applied then you might be true.
Thanks for your help.
Best,
Ronak
Hmm, what you want seems like something that should be simple to do. Do you need us to look into this more? If so, we have to wait for members of our team to get back from vacation/workshops. It might be a while but we'd be happy to follow up. Let me know.
Thanks @shlee, i agree that it should be simple to do, but want to understand the correct behaviour of FindCoveredIntervals. We actually have already done what we need but yes it will be great once people are back to have better explanation.
Issue · Github
by shlee
@shahr2, I've asked someone to look into this and they should get back to you.
@shahr2
Hi Ronak,
We did not forget about your question. I will get back to you as soon as the team reaches a conclusion.
-Sheila
@shahr2 This turned out to be a dumb bug: the value passed in by
-covis incorrectly used as the minimum base quality threshold instead of the minimum coverage threshold. We'll get this fixed. Thanks for reporting it!Thanks @Geraldine_VdAuwera, so was the default value still used, which is set to 20
Right now the argument values are getting flipped, so it looks like if you specify any parameters at all the tool will do the wrong thing. If you run with all default values I think you should get the correct behavior (for those settings).
Thanks, will be looking forward for the fix in the next release.
This has been fixed in development; if you're in a rush the fix will be available in the nightly builds (see Downloads) starting tomorrow.