1 package org.broadinstitute.hellbender.tools.walkers.annotator;
2 
3 import com.google.common.collect.ImmutableMap;
4 import com.google.common.collect.TreeMultiset;
5 import htsjdk.variant.variantcontext.Allele;
6 import htsjdk.variant.variantcontext.VariantContext;
7 import htsjdk.variant.vcf.VCFHeaderLineCount;
8 import htsjdk.variant.vcf.VCFHeaderLineType;
9 import htsjdk.variant.vcf.VCFInfoHeaderLine;
10 import org.broadinstitute.hellbender.engine.ReferenceContext;
11 import org.broadinstitute.hellbender.utils.QualityUtils;
12 import org.broadinstitute.hellbender.utils.Utils;
13 import org.broadinstitute.hellbender.utils.genotyper.AlleleLikelihoods;
14 import org.broadinstitute.hellbender.utils.read.GATKRead;
15 
16 import java.util.*;
17 import java.util.stream.Collectors;
18 
19 public class BaseQualityHistogram extends InfoFieldAnnotation {
20 
21     public static final String KEY = "BQHIST";
22 
annotate(final ReferenceContext ref, final VariantContext vc, final AlleleLikelihoods<GATKRead, Allele> likelihoods)23     public Map<String, Object> annotate(final ReferenceContext ref,
24                                         final VariantContext vc,
25                                         final AlleleLikelihoods<GATKRead, Allele> likelihoods) {
26         Utils.nonNull(vc);
27         if ( likelihoods == null ) {
28             return Collections.emptyMap();
29         }
30 
31         final Map<Allele, TreeMultiset<Integer>> values = likelihoods.alleles().stream()
32                 .collect(Collectors.toMap(a -> a, a -> TreeMultiset.create()));
33 
34         Utils.stream(likelihoods.bestAllelesBreakingTies())
35                 .filter(ba -> ba.isInformative() && isUsableRead(ba.evidence))
36                 .forEach(ba -> BaseQuality.getBaseQuality(ba.evidence, vc).ifPresent(v -> values.get(ba.allele).add(v)));
37 
38 
39         final List<Integer> distinctBaseQualities = likelihoods.alleles().stream()
40                 .flatMap(a -> values.get(a).stream())
41                 .distinct()
42                 .sorted()
43                 .collect(Collectors.toList());
44 
45         final List<Integer> output = new ArrayList<>();
46 
47         for (final int qual : distinctBaseQualities) {
48             output.add(qual);
49             likelihoods.alleles().forEach(allele -> output.add(values.get(allele).count(qual)));
50         }
51 
52 
53         return ImmutableMap.of(KEY, output);
54     }
55 
56     @Override
getDescriptions()57     public List<VCFInfoHeaderLine> getDescriptions() {
58         return Arrays.asList(new VCFInfoHeaderLine(KEY, VCFHeaderLineCount.A, VCFHeaderLineType.Integer,
59                 "Base quality counts for each allele represented sparsely as alternating entries of qualities and counts for each allele." +
60                 "For example [10,1,0,20,0,1] means one ref base with quality 10 and one alt base with quality 20."));
61     }
62 
63     @Override
getKeyNames()64     public List<String> getKeyNames() { return Arrays.asList(KEY); }
65 
isUsableRead(final GATKRead read)66     private static boolean isUsableRead(final GATKRead read) {
67         return read.getMappingQuality() != 0 && read.getMappingQuality() != QualityUtils.MAPPING_QUALITY_UNAVAILABLE;
68     }
69 }
70