Allele Frequency QC Metric for Arrays #6039

skwalker · 2019-07-16T14:51:32Z

This extends Variant Eval to compare AFs between variants in binned AF buckets based on Thousand Genomes VCF, between the expected AF from Thousand Genomes and the seen one in the actual VCF, to be used as a QC metric for our arrays pipeline.

codecov · 2019-07-16T15:32:56Z

Codecov Report

❗ No coverage uploaded for pull request base (master@f499656). Click here to learn what that means.
The diff coverage is 81.513%.

@@            Coverage Diff             @@
##             master     #6039   +/-   ##
==========================================
  Coverage          ?   87.572%           
  Complexity        ?     36801           
==========================================
  Files             ?      2044           
  Lines             ?    166417           
  Branches          ?     19264           
==========================================
  Hits              ?    145735           
  Misses            ?     14365           
  Partials          ?      6317

Impacted Files	Coverage Δ	Complexity Δ
...itute/hellbender/utils/report/GATKReportTable.java	`70.522% <100%> (ø)`	`68 <1> (?)`
...lbender/tools/walkers/varianteval/VariantEval.java	`90.769% <100%> (ø)`	`144 <3> (?)`
.../varianteval/AlleleFrequencyQCIntegrationTest.java	`100% <100%> (ø)`	`3 <3> (?)`
...nder/metrics/analysis/AlleleFrequencyQCMetric.java	`100% <100%> (ø)`	`1 <1> (?)`
...s/varianteval/stratifications/AlleleFrequency.java	`53.333% <60.87%> (ø)`	`6 <4> (?)`
...kers/varianteval/evaluators/PVariantEvaluator.java	`78.571% <78.571%> (ø)`	`9 <9> (?)`
...r/tools/walkers/varianteval/AlleleFrequencyQC.java	`86.538% <86.538%> (ø)`	`10 <10> (?)`

takutosato

Great work Sarah. This is good stuff.

takutosato · 2019-07-19T19:26:06Z

src/main/java/org/broadinstitute/hellbender/tools/walkers/varianteval/AlleleFrequencyQC.java

+
+    @Argument(shortName = "pvalue-threshold",
+            doc = "Threshold to cut off the pvalue.")
+    protected Double threshold = 0.05;


use the primitive double unless you have a reason to use the object Double

takutosato · 2019-07-19T19:34:41Z

src/main/java/org/broadinstitute/hellbender/tools/walkers/varianteval/AlleleFrequencyQC.java

+
+        super.onTraversalSuccess();
+
+        GATKReportTable table= new GATKReport(outFile).getTable(MODULES_TO_USE.get(0));


all these variables declared here should be final

takutosato · 2019-07-29T15:31:18Z

src/main/java/org/broadinstitute/hellbender/tools/walkers/varianteval/AlleleFrequencyQC.java

+            doc = "File to save the results from variant eval for debugging", optional = true)
+    protected File debugFile;
+
+    private String R_SCRIPT = "plotAlleleFrequencyQC.R";


R_SCRIPT -> rScript.

Reserve the "all caps and underscore" case for static final constants

I think this should actually be a static final constant

takutosato · 2019-07-29T16:28:52Z

...va/org/broadinstitute/hellbender/tools/walkers/varianteval/evaluators/PVariantEvaluator.java

+    public int totalCalledSites;
+    @DataPoint(description = "Number of called heterozygous sites;", format = "%d")
+    public int totalHetSites;
+    @DataPoint(description = "Number of called heterozygous sites;", format = "%d")


heterozygous -> hom var

takutosato · 2019-07-29T16:29:11Z

...va/org/broadinstitute/hellbender/tools/walkers/varianteval/evaluators/PVariantEvaluator.java

+    public int totalHetSites;
+    @DataPoint(description = "Number of called heterozygous sites;", format = "%d")
+    public int totalHomVarSites;
+    @DataPoint(description = "Number of called heterozygous sites;", format = "%d")


heterozygous -> hom ref

takutosato · 2019-07-29T18:45:30Z

...org/broadinstitute/hellbender/tools/walkers/varianteval/stratifications/AlleleFrequency.java

+
+    private int getLogitBucket(Double alleleFraction) {
+        // calculate logit value
+        Float score = (float)( -10 * Math.log10((alleleFraction/(1-alleleFraction))));


If it's not too much work, I would remove the factor of -10 and rescale the bin width and limits. Logit by definition doesn't have the factor of -10, and you shouldn't take phred scale of something that's not a probability i.e. p/(1-p). It's an unnecessary layer that could potentially mislead the user.

Ok so I talked about this with Yossi, and I think we actually want to keep it like this; first because it means a bin width of 1 works really well, and he thought that the idea of phred-scale was actually beneficial. So I think I will keep it like this.

OK. And this way seems to have binned the variants nicely so I guess there's no reason to change.

takutosato · 2019-07-29T18:46:49Z

...org/broadinstitute/hellbender/tools/walkers/varianteval/stratifications/AlleleFrequency.java

+                }
+                break;
+            case LOGARITHMIC:
+                for (int a = -logLimit; a <= logLimit; a += 1) {


extract a variable private static int LOGIT_BIN_WIDTH = 1;

takutosato · 2019-07-29T19:02:27Z

...va/org/broadinstitute/hellbender/tools/walkers/varianteval/evaluators/PVariantEvaluator.java

+import org.broadinstitute.hellbender.tools.walkers.varianteval.util.Analysis;
+import org.broadinstitute.hellbender.tools.walkers.varianteval.util.DataPoint;
+
+@Analysis(description = "Computes different estimates of theta based on variant sites and genotypes")


What's theta?

I think it'd be good to document the method succinctly; i.e. divide sites by AF in 1000G, then within each bin add up the alt alleles and divide by 2*num_sites, and do the goodness of fit test to get the p-value.

I think that documentation should be for AllleFrequencyQC; this is just the evaluator being used to calculate average var allele fraction (P --> why its PVariantEvaluator) and count the number of hets/hom vars/ hom refs; definitely updating description but not sure about a new name unless you have a better idea ?

Yea, docs should go wherever you think is appropriate, I just happened to think that when I was reviewing this part.

As for the name, maybe something like AggregateAltFraction

How about VariantAFEvaluator ?
I think there is decent documentation in AlleleFrequencyQC.java -- let me know if you think it's unclear?

I think the first sentence in AlleleFreqeuncyQC.java docs "This tool uses Variant Eval..." could be improved. What do you think?

takutosato · 2019-07-29T19:02:52Z

...va/org/broadinstitute/hellbender/tools/walkers/varianteval/evaluators/PVariantEvaluator.java

+
+@Analysis(description = "Computes different estimates of theta based on variant sites and genotypes")
+
+public class PVariantEvaluator extends VariantEvaluator {


Could we pick a more descriptive name than PVariantEvaluator?

takutosato · 2019-07-29T19:31:25Z

src/main/java/org/broadinstitute/hellbender/tools/walkers/varianteval/AlleleFrequencyQC.java

+        return null;
+    }
+
+    // This creates a modified chi squared statistic that allows for a constant variance (1%) rather than scaling


I'm going to have to think about this.

takutosato · 2019-07-30T16:11:02Z

...org/broadinstitute/hellbender/tools/walkers/varianteval/stratifications/AlleleFrequency.java

        }
    }

    public List<Object> getRelevantStates(ReferenceContext referenceContext, ReadsContext readsContext, FeatureContext featureContext, VariantContext comp, String compName, VariantContext eval, String evalName, String sampleName, String FamilyName) {
        if (eval != null) {
            try {
-                return Collections.singletonList((Object)String.format("%.3f", (5.0 * MathUtils.roundToNDecimalPlaces(eval.getAttributeAsDouble("AF", 0.0) / 5.0, 3))));
+                Double alleleFraction  = Collections.max(eval.getAttributeAsDoubleList("AF", 0.0));


alleleFraction -> alleleFrequency to minimize confusion (and also because the class is called AlleleFrequency)

...va/org/broadinstitute/hellbender/tools/walkers/varianteval/evaluators/PVariantEvaluator.java

takutosato

Rewrite the doc as you see fit and it's ready to go

lbergelson

@skwalker Looks good. Retag the docker image and push it to docker hub, then update the pr and make sure tests pass with the new image. 👍 from me when that's green.

lbergelson · 2019-09-18T18:46:40Z

scripts/docker/gatkbase/build_docker_base.sh

@@ -6,7 +6,7 @@ set -e

 REPO=broadinstitute
 PROJECT=gatk
-VERSION=2.0.3
+VERSION=2.1.0


Update to 2.2.0 to avoid conflict with the existing wierd 2.1.0

lbergelson · 2019-09-18T18:47:08Z

Dockerfile

@@ -1,5 +1,5 @@
 # Using OpenJDK 8
-FROM broadinstitute/gatk:gatkbase-2.0.3
+FROM skwalker/gatk:gatkbase-2.1.0


Retag as broadinstitute/gatk:gatkbase-2.2.0 and push to dockerhub.

lbergelson · 2019-09-18T18:52:44Z

Oh, be sure to mention that you updated the base image and included dplyr in the commit message.

skwalker requested a review from takutosato July 16, 2019 14:51

takutosato self-assigned this Jul 29, 2019

takutosato requested changes Jul 29, 2019

View reviewed changes

takutosato reviewed Jul 30, 2019

View reviewed changes

...va/org/broadinstitute/hellbender/tools/walkers/varianteval/evaluators/PVariantEvaluator.java Outdated Show resolved Hide resolved

takutosato approved these changes Sep 6, 2019

View reviewed changes

skwalker added 10 commits September 17, 2019 15:26

allele frequency qc with sampling from thousand genomes

42d7219

updates to af qc metric for arrays

dabde31

adding updates to gatk report table

3680e81

minor tweaks

015a493

adding dplyr dependencies

44761d5

trying my docker base to get integration test to pass hopefully

24fdeba

try 2

681e876

responding to takuto's comments

2147bf6

updated doc

41a1cc4

docker changes

a317675

skwalker force-pushed the skwalker_af_arrays_qc branch from 00c8ee2 to a317675 Compare September 17, 2019 20:38

lbergelson reviewed Sep 18, 2019

View reviewed changes

updating docker base image that includes dplyr

4243512

skwalker merged commit 895843f into master Sep 19, 2019

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Allele Frequency QC Metric for Arrays #6039

Allele Frequency QC Metric for Arrays #6039

skwalker commented Jul 16, 2019

codecov bot commented Jul 16, 2019 •

edited

Loading

takutosato left a comment

takutosato Jul 19, 2019

takutosato Jul 19, 2019

takutosato Jul 29, 2019

skwalker Sep 5, 2019

takutosato Sep 6, 2019

takutosato Jul 29, 2019

takutosato Jul 29, 2019

takutosato Jul 29, 2019

skwalker Sep 5, 2019

takutosato Sep 6, 2019

takutosato Jul 29, 2019

takutosato Jul 29, 2019

takutosato Jul 29, 2019

skwalker Jul 30, 2019

takutosato Jul 30, 2019

skwalker Sep 5, 2019

takutosato Sep 6, 2019

takutosato Jul 29, 2019

takutosato Jul 29, 2019

takutosato Jul 30, 2019

takutosato left a comment

lbergelson left a comment

lbergelson Sep 18, 2019

lbergelson Sep 18, 2019

lbergelson commented Sep 18, 2019


		super.onTraversalSuccess();

		GATKReportTable table= new GATKReport(outFile).getTable(MODULES_TO_USE.get(0));


		@Analysis(description = "Computes different estimates of theta based on variant sites and genotypes")

		public class PVariantEvaluator extends VariantEvaluator {

Allele Frequency QC Metric for Arrays #6039

Allele Frequency QC Metric for Arrays #6039

Conversation

skwalker commented Jul 16, 2019

codecov bot commented Jul 16, 2019 • edited Loading

Codecov Report

takutosato left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

takutosato left a comment

Choose a reason for hiding this comment

lbergelson left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

lbergelson commented Sep 18, 2019

codecov bot commented Jul 16, 2019 •

edited

Loading