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