Add bamout read annotations for alignement and callable regions, and supported genotypes #5215

kachulis · 2018-09-24T15:32:58Z

Adds three annotations to reads in bamout produced by HaplotypeCaller and Mutect2:

AR (Alignment Region)- The reference region over which a read is aligned to a haplotype
CR (Callable Region)- The reference region over which the caller is using this read to call variants
SG (Supported Genotypes)- A list of the genotypes supported by this read for variants called in the callable region

davidbenjamin

Small comments

davidbenjamin · 2018-09-27T17:01:05Z

...va/org/broadinstitute/hellbender/tools/walkers/haplotypecaller/AssemblyBasedCallerUtils.java

@@ -35,6 +36,9 @@
 public final class AssemblyBasedCallerUtils {

    static final int REFERENCE_PADDING_FOR_ASSEMBLY = 500;
+    static final String SUPPORTED_GENOTYPES_TAG="SG";


davidbenjamin · 2018-09-27T19:09:09Z

...va/org/broadinstitute/hellbender/tools/walkers/haplotypecaller/AssemblyBasedCallerUtils.java

+     * @param likelihoodsHaplotype ReadLikelihoods containing reads to be annotated along with haplotypes to which these reads have been aligned
+     * @param callableRegion ref region over which caller is using these ReadLikelihoods to call variants
+     */
+    public static void annotateReadLikelihoodsWithRegions(final ReadLikelihoods<Haplotype> likelihoodsHaplotype,


likelihoodsHaplotype should just be likelihoods, IMO

davidbenjamin · 2018-09-27T19:10:34Z

...va/org/broadinstitute/hellbender/tools/walkers/haplotypecaller/AssemblyBasedCallerUtils.java

+    }
+
+    /**
+     * Annotates reads in ReadLiklihoods with genotype supported by each read a particular variant.  If a read already has a


"a particular variant" seems out of place

davidbenjamin · 2018-09-27T19:13:24Z

...va/org/broadinstitute/hellbender/tools/walkers/haplotypecaller/AssemblyBasedCallerUtils.java

+    public static void annotateReadLikelihoodsWithSupportedGenotypes(final VariantContext vc,
+                                                                     final ReadLikelihoods<Allele> likelihoodsAllele) {
+        //assign supported Genotypes to each read
+        final Set<Allele> alleles = new LinkedHashSet<>(vc.getAlleles());


Inline this into next line?

davidbenjamin · 2018-09-27T19:13:48Z

...va/org/broadinstitute/hellbender/tools/walkers/haplotypecaller/AssemblyBasedCallerUtils.java

+        //assign supported Genotypes to each read
+        final Set<Allele> alleles = new LinkedHashSet<>(vc.getAlleles());
+        final Map<Allele, List<Allele>> alleleSubset = alleles.stream().collect(Collectors.toMap(a -> a, Arrays::asList));
+        final ReadLikelihoods<Allele> subsettedLikelihoodsAllele = likelihoodsAllele.marginalize(alleleSubset);


Just call this subsettedLikelihoods

davidbenjamin · 2018-09-27T19:17:09Z

...broadinstitute/hellbender/tools/walkers/haplotypecaller/HaplotypeCallerGenotypingEngine.java

@@ -130,6 +132,11 @@ public CalledHaplotypes assignGenotypeLikelihoods(final List<Haplotype> haplotyp
        final int ploidy = configuration.genotypeArgs.samplePloidy;
        final List<Allele> noCallAlleles = GATKVariantContextUtils.noCallAlleles(ploidy);

+        if (withBamOut) {
+            ////add annotations to reads for alignment regions and calling regions


Just two slashes

davidbenjamin · 2018-09-27T19:17:45Z

...broadinstitute/hellbender/tools/walkers/haplotypecaller/HaplotypeCallerGenotypingEngine.java

@@ -181,6 +188,11 @@ public CalledHaplotypes assignGenotypeLikelihoods(final List<Haplotype> haplotyp
                final VariantContext annotatedCall = makeAnnotatedCall(ref, refLoc, tracker, header, mergedVC, readAlleleLikelihoods, call);
                returnCalls.add( annotatedCall );

+                if (withBamOut) {
+                    //add annotations to reads for which genotypes are supported by each read


This comment is not necessary because the method name is self-documenting

davidbenjamin · 2018-09-27T19:18:15Z

src/main/java/org/broadinstitute/hellbender/tools/walkers/mutect/SomaticGenotypingEngine.java

@@ -103,6 +105,11 @@ public CalledHaplotypes callMutations(
        final Set<Haplotype> calledHaplotypes = new HashSet<>();
        final List<VariantContext> returnCalls = new ArrayList<>();

+        if(withBamOut){
+            ////add annotations to reads for alignment regions and calling regions


Just two slashes

davidbenjamin · 2018-09-27T19:18:34Z

src/main/java/org/broadinstitute/hellbender/tools/walkers/mutect/SomaticGenotypingEngine.java

@@ -172,6 +179,10 @@ public CalledHaplotypes callMutations(
            final ReadLikelihoods<Allele> trimmedLikelihoods = log10Likelihoods.marginalize(trimmedToUntrimmedAlleleMap);

            final VariantContext annotatedCall =  annotationEngine.annotateContext(trimmedCall, featureContext, referenceContext, trimmedLikelihoods, a -> true);
+            if(withBamOut) {
+                //add annotations to reads for which genotypes are supported by each read


This comment is not necessary

davidbenjamin · 2018-09-27T19:19:15Z

...roadinstitute/hellbender/tools/walkers/haplotypecaller/AssemblyBasedCallerUtilsUnitTest.java

+        final int start = 13763;
+        final SimpleInterval loc = new SimpleInterval(contig, start, start + hap1String.length());
+        final SimpleInterval loc2 = new SimpleInterval(contig, loc.getStart() + 4, loc.getEnd() - 4);
+        List<Integer> snpIndecies = Arrays.asList(null, 12, 23);


codecov-io · 2018-09-28T15:05:04Z

Codecov Report

Merging #5215 into master will increase coverage by 0.01%.
The diff coverage is 94.91%.

@@             Coverage Diff              @@
##             master    #5215      +/-   ##
============================================
+ Coverage     86.76%   86.77%   +0.01%     
- Complexity    29944    30000      +56     
============================================
  Files          1838     1838              
  Lines        138820   139051     +231     
  Branches      15285    15329      +44     
============================================
+ Hits         120448   120667     +219     
- Misses        12804    12810       +6     
- Partials       5568     5574       +6

Impacted Files	Coverage Δ	Complexity Δ
...walkers/haplotypecaller/HaplotypeCallerEngine.java	`78.12% <100%> (+0.15%)`	`74 <0> (ø)`	⬇️
...hellbender/tools/walkers/mutect/Mutect2Engine.java	`91.08% <100%> (ø)`	`59 <0> (ø)`	⬇️
...kers/haplotypecaller/AssemblyBasedCallerUtils.java	`76.42% <100%> (+5.71%)`	`32 <6> (+6)`	⬆️
...plotypecaller/HaplotypeCallerGenotypingEngine.java	`79.87% <80%> (ø)`	`39 <0> (+2)`	⬆️
.../tools/walkers/mutect/SomaticGenotypingEngine.java	`89.67% <80%> (-0.28%)`	`73 <0> (+2)`
...hellbender/utils/haplotype/HaplotypeBAMWriter.java	`96.87% <85.71%> (-1.46%)`	`15 <1> (+1)`
...lotypecaller/AssemblyBasedCallerUtilsUnitTest.java	`95.77% <95.28%> (-4.23%)`	`45 <43> (+43)`
...utils/smithwaterman/SmithWatermanIntelAligner.java	`50% <0%> (-30%)`	`1% <0%> (-2%)`
...ithwaterman/SmithWatermanIntelAlignerUnitTest.java	`60% <0%> (ø)`	`2% <0%> (ø)`	⬇️
...te/hellbender/utils/genotyper/ReadLikelihoods.java	`90.14% <0%> (+0.4%)`	`143% <0%> (ø)`	⬇️
... and 1 more

Continue to review full report at Codecov.

Legend - Click here to learn more
Δ = absolute <relative> (impact), ø = not affected, ? = missing data
Powered by Codecov. Last update 8103bde...7d53fb9. Read the comment docs.

ldgauthier · 2018-10-01T13:50:12Z

...va/org/broadinstitute/hellbender/tools/walkers/haplotypecaller/AssemblyBasedCallerUtils.java

+        for (ReadLikelihoods<Allele>.BestAllele bestAllele : bestAlleles) {
+            String attribute = "";
+            GATKRead read = bestAllele.read;
+            Allele allele = bestAllele.allele;


It's not obvious to me why these are supported genotypes. What is the value of the tag supposed to be?

I think you're right. The value of the tag is really the alleles of called variants that each read supports. For example, if a read overlaps two called variants in the active region at 1:10000 and 1:10010, and supports the reference allele at the first and the first alt allele at the second, the annotation will be 1:10000=0, 1:10010=1. Essentially this is the same information encoded by allele depth in a vcf. So I think it is more correct to call these supported alleles as opposed to supported genotypes.

ldgauthier · 2018-10-01T13:50:15Z

...va/org/broadinstitute/hellbender/tools/walkers/haplotypecaller/AssemblyBasedCallerUtils.java

+    }
+
+    /**
+     * Annotates reads in ReadLiklihoods with genotype supported by each read a particular variant.  If a read already has a


"Annotates reads in ReadLikelihoods with genotype supported by each read at a particular variant." ?

kachulis · 2018-10-04T14:38:54Z

back to you when you get a chance to take another look, @davidbenjamin and @ldgauthier

…edCallerUtils

…upported alleles (broadinstitute#5215) Add bamout read annotations for alignment and callable regions, and supported alleles (broadinstitute#5215)

droazen requested review from ldgauthier and davidbenjamin September 26, 2018 22:01

droazen assigned ldgauthier and davidbenjamin Sep 26, 2018

davidbenjamin requested changes Sep 27, 2018

View reviewed changes

ldgauthier reviewed Oct 1, 2018

View reviewed changes

davidbenjamin approved these changes Oct 11, 2018

View reviewed changes

davidbenjamin removed their assignment Oct 11, 2018

kachulis added 11 commits October 12, 2018 13:10

moving active region labeling to HaplotupeBAMWriter implementations

3ca7095

stuff

cd0cf0f

Adding bamout region and supported genotype annotation to AssemblyBas…

b4eae17

…edCallerUtils

adding bamout region and supported genotype annotations to Mutect2

71a74b8

tests for bamout region and supported haplotype annotations

ae72b67

removing extra line break in HaplotypeCallerGenotypingEngine

85e286b

aesthetics

16f8a2b

refactoring building of reads in AssemblyBasedCallerUtilsUnitTest

e5c84d9

cleaning things up and adding comments

c7b168d

changes requested by David

7d74917

changing supported genotypes to supported alleles in bamout annotation

7d53fb9

kachulis force-pushed the ck_bamout_annotate_reads_active_regions branch from 93c7a9e to 7d53fb9 Compare October 12, 2018 17:20

kachulis merged commit 17d17fd into master Oct 12, 2018

kachulis deleted the ck_bamout_annotate_reads_active_regions branch October 12, 2018 19:32

kachulis mentioned this pull request Nov 15, 2018

Changing SUPPORTED_ALLELES_TAG from SA to XA in bamout #5418

Merged

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Add bamout read annotations for alignement and callable regions, and supported genotypes #5215

Add bamout read annotations for alignement and callable regions, and supported genotypes #5215

kachulis commented Sep 24, 2018

davidbenjamin left a comment

davidbenjamin Sep 27, 2018

davidbenjamin Sep 27, 2018

davidbenjamin Sep 27, 2018

davidbenjamin Sep 27, 2018

davidbenjamin Sep 27, 2018

davidbenjamin Sep 27, 2018

davidbenjamin Sep 27, 2018

davidbenjamin Sep 27, 2018

davidbenjamin Sep 27, 2018

davidbenjamin Sep 27, 2018

codecov-io commented Sep 28, 2018 •

edited

Loading

ldgauthier Oct 1, 2018

kachulis Oct 1, 2018

ldgauthier Oct 1, 2018

kachulis commented Oct 4, 2018

Add bamout read annotations for alignement and callable regions, and supported genotypes #5215

Add bamout read annotations for alignement and callable regions, and supported genotypes #5215

Conversation

kachulis commented Sep 24, 2018

davidbenjamin 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

codecov-io commented Sep 28, 2018 • edited Loading

Codecov Report

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

kachulis commented Oct 4, 2018

codecov-io commented Sep 28, 2018 •

edited

Loading