diff --git a/src/main/java/org/broad/igv/sam/AlignmentCounts.java b/src/main/java/org/broad/igv/sam/AlignmentCounts.java index 98ce5a69a2..f0079e7895 100644 --- a/src/main/java/org/broad/igv/sam/AlignmentCounts.java +++ b/src/main/java/org/broad/igv/sam/AlignmentCounts.java @@ -38,6 +38,10 @@ public interface AlignmentCounts extends Feature { int getTotalCount(int pos); + public int getTotalPositiveCount(int pos); + + public int getTotalNegativeCount(int pos); + int getTotalQuality(int pos); int getCount(int pos, byte b); diff --git a/src/main/java/org/broad/igv/sam/AlignmentDataManager.java b/src/main/java/org/broad/igv/sam/AlignmentDataManager.java index 6412852b53..222fe3a55f 100644 --- a/src/main/java/org/broad/igv/sam/AlignmentDataManager.java +++ b/src/main/java/org/broad/igv/sam/AlignmentDataManager.java @@ -36,6 +36,7 @@ import org.broad.igv.feature.genome.Genome; import org.broad.igv.prefs.IGVPreferences; import org.broad.igv.prefs.PreferencesManager; +import org.broad.igv.sam.mods.BaseModificationKey; import org.broad.igv.sam.mods.BaseModificationSet; import org.broad.igv.sam.reader.AlignmentReader; import org.broad.igv.sam.reader.AlignmentReaderFactory; @@ -69,7 +70,9 @@ public class AlignmentDataManager implements IGVEventObserver { private Map peStats; private SpliceJunctionHelper.LoadOptions loadOptions; - private Set allBaseModifications = new HashSet<>(); + private Set allBaseModificationKeys = new HashSet<>(); + + private Set simplexBaseModfications = new HashSet<>(); private Range currentlyLoading; @@ -198,8 +201,8 @@ public Map getPEStats() { /** * Return all base modfications seen in loaded alignments */ - public Set getAllBaseModifications() { - return allBaseModifications; + public Set getAllBaseModificationKeys() { + return allBaseModificationKeys; } public boolean isPairedEnd() { @@ -407,14 +410,28 @@ AlignmentInterval loadInterval(String chr, int start, int end, AlignmentTrack.Re } private void updateBaseModfications(List alignments) { + for(Alignment a : alignments) { List bmSets = a.getBaseModificationSets(); if(bmSets != null) { for(BaseModificationSet bms : bmSets) { - allBaseModifications.add(bms.getModification()); + allBaseModificationKeys.add(BaseModificationKey.getKey(bms.getBase(), bms.getStrand(), bms.getModification())); } } } + // Search for simplex modifications (single strand read, e.g. C+m with no G-m) + Map tmp = new HashMap<>(); + for(BaseModificationKey key : allBaseModificationKeys) { + if(tmp.containsKey(key.getModification())) { + tmp.remove(key); + } else { + tmp.put(key.getModification(), key); + } + } + for(Map.Entry entries : tmp.entrySet()) { + simplexBaseModfications.add(entries.getKey()); + simplexBaseModfications.add("NONE_" + entries.getValue().getCanonicalBase()); + } } public AlignmentTrack.ExperimentType inferType() { @@ -573,6 +590,10 @@ public Collection getLoadedIntervals() { return intervalCache; } + public Set getSimplexBaseModifications() { + return simplexBaseModfications; + } + public static class DownsampleOptions { private boolean downsample; private int sampleWindowSize; diff --git a/src/main/java/org/broad/igv/sam/AlignmentRenderer.java b/src/main/java/org/broad/igv/sam/AlignmentRenderer.java index 2cd70ad052..54c5929d38 100644 --- a/src/main/java/org/broad/igv/sam/AlignmentRenderer.java +++ b/src/main/java/org/broad/igv/sam/AlignmentRenderer.java @@ -1299,7 +1299,7 @@ private Color getAlignmentColor(Alignment alignment, AlignmentTrack track) { case BISULFITE: case BASE_MODIFICATION: - case BASE_MODIFICATION_C: + case BASE_MODIFICATION_2COLOR: case SMRT_SUBREAD_IPD: case SMRT_SUBREAD_PW: case SMRT_CCS_FWD_IPD: diff --git a/src/main/java/org/broad/igv/sam/AlignmentTrack.java b/src/main/java/org/broad/igv/sam/AlignmentTrack.java index 11c076b034..80e15337f2 100644 --- a/src/main/java/org/broad/igv/sam/AlignmentTrack.java +++ b/src/main/java/org/broad/igv/sam/AlignmentTrack.java @@ -41,6 +41,7 @@ import org.broad.igv.prefs.IGVPreferences; import org.broad.igv.prefs.PreferencesManager; import org.broad.igv.renderer.GraphicUtils; +import org.broad.igv.sam.mods.BaseModficationFilter; import org.broad.igv.session.Persistable; import org.broad.igv.track.*; import org.broad.igv.ui.FontManager; @@ -101,8 +102,7 @@ public enum ColorOption { LINK_STRAND, YC_TAG, BASE_MODIFICATION, - BASE_MODIFICATION_5MC, - BASE_MODIFICATION_C, + BASE_MODIFICATION_2COLOR, SMRT_SUBREAD_IPD, SMRT_SUBREAD_PW, SMRT_CCS_FWD_IPD, @@ -111,7 +111,7 @@ public enum ColorOption { SMRT_CCS_REV_PW; public boolean isBaseMod() { - return this == BASE_MODIFICATION || this == BASE_MODIFICATION_C; + return this == BASE_MODIFICATION || this == BASE_MODIFICATION_2COLOR; } public boolean isSMRTKinetics() { @@ -1330,7 +1330,7 @@ public static class RenderOptions implements Cloneable, Persistable { private Boolean hideSmallIndels; private Integer smallIndelThreshold; - private String basemodFilter; + private BaseModficationFilter basemodFilter; BisulfiteContext bisulfiteContext = BisulfiteContext.CG; Map peStats; @@ -1602,11 +1602,11 @@ public int getSmallIndelThreshold() { return smallIndelThreshold == null ? getPreferences().getAsInt(SAM_SMALL_INDEL_BP_THRESHOLD) : smallIndelThreshold; } - public String getBasemodFilter() { + public BaseModficationFilter getBasemodFilter() { return basemodFilter; } - public void setBasemodFilter(String basemodFilter) { + public void setBasemodFilter(BaseModficationFilter basemodFilter) { this.basemodFilter = basemodFilter; } @@ -1707,7 +1707,7 @@ public void marshalXML(Document document, Element element) { element.setAttribute("showInsertionMarkers", showInsertionMarkers.toString()); } if (basemodFilter != null) { - element.setAttribute("basemodfilter", basemodFilter); + element.setAttribute("basemodfilter", basemodFilter.toString()); } } @@ -1741,11 +1741,14 @@ public void unmarshalXML(Element element, Integer version) { final String attributeValue = element.getAttribute("colorOption"); if("BASE_MODIFICATION_6MA".equals(attributeValue)) { colorOption = ColorOption.BASE_MODIFICATION; - basemodFilter = "a"; + basemodFilter = new BaseModficationFilter("a"); } else if("BASE_MODIFICATION_5MC".equals(attributeValue)) { - colorOption = ColorOption.BASE_MODIFICATION_C; - basemodFilter = "m"; - } else { + colorOption = ColorOption.BASE_MODIFICATION_2COLOR; + basemodFilter = new BaseModficationFilter(null, 'C'); + } else if("BASE_MODIFICATION_C".equals(attributeValue)) { + colorOption = ColorOption.BASE_MODIFICATION; + basemodFilter = new BaseModficationFilter(null, 'C'); + }else { colorOption = ColorOption.valueOf(attributeValue); } } @@ -1822,7 +1825,7 @@ public void unmarshalXML(Element element, Integer version) { showInsertionMarkers = Boolean.parseBoolean(element.getAttribute("showInsertionMarkers")); } if (element.hasAttribute("basemodfilter")) { - basemodFilter = element.getAttribute("basemodfilter"); + basemodFilter = BaseModficationFilter.fromString(element.getAttribute("basemodfilter")); } } diff --git a/src/main/java/org/broad/igv/sam/AlignmentTrackMenu.java b/src/main/java/org/broad/igv/sam/AlignmentTrackMenu.java index 9b5712ef57..e07574cc05 100644 --- a/src/main/java/org/broad/igv/sam/AlignmentTrackMenu.java +++ b/src/main/java/org/broad/igv/sam/AlignmentTrackMenu.java @@ -1,6 +1,5 @@ package org.broad.igv.sam; -import com.google.gson.annotations.Since; import htsjdk.samtools.SAMTag; import htsjdk.samtools.util.Locatable; import org.broad.igv.Globals; @@ -12,6 +11,7 @@ import org.broad.igv.logging.LogManager; import org.broad.igv.logging.Logger; import org.broad.igv.prefs.PreferencesManager; +import org.broad.igv.sam.mods.BaseModficationFilter; import org.broad.igv.sam.mods.BaseModificationUtils; import org.broad.igv.sashimi.SashimiPlot; import org.broad.igv.tools.PFMExporter; @@ -40,8 +40,6 @@ import java.awt.event.MouseEvent; import java.util.*; import java.util.List; -import java.util.concurrent.CompletableFuture; -import java.util.concurrent.Future; import java.util.function.BiConsumer; import java.util.stream.Collectors; @@ -573,9 +571,7 @@ private JRadioButtonMenuItem getColorMenuItem(String label, final AlignmentTrack mi.setSelected(renderOptions.getColorOption() == option); mi.addActionListener(aEvt -> { alignmentTrack.setColorOption(option); - if (option == AlignmentTrack.ColorOption.BASE_MODIFICATION) { - renderOptions.setBasemodFilter(extra); - } + renderOptions.setBasemodFilter(extra == null ? null : new BaseModficationFilter(extra)); alignmentTrack.repaint(); }); @@ -643,42 +639,45 @@ void addColorByMenuItem() { // Base modifications JRadioButtonMenuItem bmMenuItem; - colorMenu.addSeparator(); - Set allModifications = dataManager.getAllBaseModifications(); + Set allModifications = dataManager.getAllBaseModificationKeys().stream().map(bmKey -> bmKey.getModification()).collect(Collectors.toSet()); if (allModifications.size() > 0) { + BaseModficationFilter filter = renderOptions.getBasemodFilter(); - Set cModifications = allModifications.stream().filter(m -> BaseModificationUtils.cModifications.contains(m)).collect(Collectors.toSet()); + colorMenu.addSeparator(); + if (allModifications.size() > 1) { + bmMenuItem = getColorMenuItem("base modification 2-color (all)", AlignmentTrack.ColorOption.BASE_MODIFICATION_2COLOR); + bmMenuItem.setSelected(renderOptions.getColorOption() == AlignmentTrack.ColorOption.BASE_MODIFICATION_2COLOR && filter == null); + colorMenu.add(bmMenuItem); + group.add(bmMenuItem); + } - for (String m : cModifications) { + for (String m : allModifications) { String name = BaseModificationUtils.modificationName(m); - bmMenuItem = getColorMenuItem("base modification (" + name + ")", AlignmentTrack.ColorOption.BASE_MODIFICATION_C, m); - bmMenuItem.setSelected(renderOptions.getColorOption() == AlignmentTrack.ColorOption.BASE_MODIFICATION_C && m.equals(renderOptions.getBasemodFilter())); + bmMenuItem = getColorMenuItem("base modification 2-color (" + name + ")", AlignmentTrack.ColorOption.BASE_MODIFICATION_2COLOR, m); + bmMenuItem.setSelected(renderOptions.getColorOption() == + AlignmentTrack.ColorOption.BASE_MODIFICATION_2COLOR && + (filter != null && filter.pass(m))); colorMenu.add(bmMenuItem); group.add(bmMenuItem); } - if (cModifications.size() > 1) { - bmMenuItem = getColorMenuItem("CpG modification (all C)", AlignmentTrack.ColorOption.BASE_MODIFICATION_C); - bmMenuItem.setSelected(renderOptions.getColorOption() == AlignmentTrack.ColorOption.BASE_MODIFICATION_C); + colorMenu.addSeparator(); + if (allModifications.size() > 1) { + bmMenuItem = getColorMenuItem("base modification (all)", AlignmentTrack.ColorOption.BASE_MODIFICATION); + bmMenuItem.setSelected(renderOptions.getColorOption() == AlignmentTrack.ColorOption.BASE_MODIFICATION && filter == null); colorMenu.add(bmMenuItem); group.add(bmMenuItem); } - colorMenu.addSeparator(); for (String m : allModifications) { String name = BaseModificationUtils.modificationName(m); bmMenuItem = getColorMenuItem("base modification (" + name + ")", AlignmentTrack.ColorOption.BASE_MODIFICATION, m); - bmMenuItem.setSelected(renderOptions.getColorOption() == AlignmentTrack.ColorOption.BASE_MODIFICATION && m.equals(renderOptions.getBasemodFilter())); - colorMenu.add(bmMenuItem); - group.add(bmMenuItem); - } - if (allModifications.size() > 1) { - bmMenuItem = getColorMenuItem("base modification (all)", AlignmentTrack.ColorOption.BASE_MODIFICATION); - bmMenuItem.setSelected(renderOptions.getColorOption() == AlignmentTrack.ColorOption.BASE_MODIFICATION && renderOptions.getBasemodFilter() == null); + bmMenuItem.setSelected(renderOptions.getColorOption() == AlignmentTrack.ColorOption.BASE_MODIFICATION && (filter != null && filter.pass(m))); colorMenu.add(bmMenuItem); group.add(bmMenuItem); } + } diff --git a/src/main/java/org/broad/igv/sam/CoverageTrack.java b/src/main/java/org/broad/igv/sam/CoverageTrack.java index d386f6a2e0..4561f00bf8 100644 --- a/src/main/java/org/broad/igv/sam/CoverageTrack.java +++ b/src/main/java/org/broad/igv/sam/CoverageTrack.java @@ -39,6 +39,7 @@ import org.broad.igv.prefs.IGVPreferences; import org.broad.igv.prefs.PreferencesManager; import org.broad.igv.renderer.*; +import org.broad.igv.sam.mods.BaseModficationFilter; import org.broad.igv.sam.mods.BaseModificationCoverageRenderer; import org.broad.igv.tdf.TDFDataSource; import org.broad.igv.tdf.TDFReader; @@ -505,6 +506,7 @@ private void paint(final RenderContext context, final Rectangle rect, final Alig } // Second pass -- potentially overlay mismatches + Set simplexModifications = dataManager.getSimplexBaseModifications(); for (int idx = 0; idx < nPoints; idx++) { int pos = isSparse ? ((SparseAlignmentCounts) alignmentCounts).getPosition(idx) : start + idx * step; @@ -531,8 +533,8 @@ private void paint(final RenderContext context, final Rectangle rect, final Alig drawBarBisulfite(context, pX, bottomY, dX, barHeight, totalCount, bc); } } else if (colorOption.isBaseMod()) { - String basemodFilter = alignmentTrack != null ? alignmentTrack.getRenderOptions().getBasemodFilter() : null; - BaseModificationCoverageRenderer.drawModifications(context, pX, bottomY, dX, barHeight, pos, alignmentCounts, colorOption, basemodFilter); + BaseModficationFilter basemodFilter = alignmentTrack != null ? alignmentTrack.getRenderOptions().getBasemodFilter() : null; + BaseModificationCoverageRenderer.drawModifications(context, pX, bottomY, dX, barHeight, pos, alignmentCounts, colorOption, basemodFilter, simplexModifications); } else { if (refBases != null) { int refIdx = pos - intervalStart; diff --git a/src/main/java/org/broad/igv/sam/DenseAlignmentCounts.java b/src/main/java/org/broad/igv/sam/DenseAlignmentCounts.java index fefd7a6a37..4aeb31a23d 100644 --- a/src/main/java/org/broad/igv/sam/DenseAlignmentCounts.java +++ b/src/main/java/org/broad/igv/sam/DenseAlignmentCounts.java @@ -134,7 +134,30 @@ public int getTotalCount(int pos) { return 0; } else { return posTotal[offset] + negTotal[offset]; + } + } + public int getTotalPositiveCount(int pos) { + int offset = pos - start; + if (offset < 0 || offset >= posA.length) { + if (log.isDebugEnabled()) { + log.debug("Position out of range: " + pos + " (valid range - " + start + "-" + end); + } + return 0; + } else { + return posTotal[offset]; + } + } + + public int getTotalNegativeCount(int pos) { + int offset = pos - start; + if (offset < 0 || offset >= posA.length) { + if (log.isDebugEnabled()) { + log.debug("Position out of range: " + pos + " (valid range - " + start + "-" + end); + } + return 0; + } else { + return negTotal[offset]; } } diff --git a/src/main/java/org/broad/igv/sam/ReducedMemoryAlignment.java b/src/main/java/org/broad/igv/sam/ReducedMemoryAlignment.java index 90fe24cf65..24421b73d9 100644 --- a/src/main/java/org/broad/igv/sam/ReducedMemoryAlignment.java +++ b/src/main/java/org/broad/igv/sam/ReducedMemoryAlignment.java @@ -468,10 +468,19 @@ public int getTotalCount(int pos) { return 0; } else { return (int) Math.round(total[offset]); - } } + @Override + public int getTotalPositiveCount(int pos) { + throw new RuntimeException(" Method getTotalPositiveCount not implemented"); + } + + @Override + public int getTotalNegativeCount(int pos) { + throw new RuntimeException(" Method getTotalNegativeCount not implemented"); + } + @Override public void incCounts(Alignment alignment) { diff --git a/src/main/java/org/broad/igv/sam/SparseAlignmentCounts.java b/src/main/java/org/broad/igv/sam/SparseAlignmentCounts.java index d2852f0eaf..02de430601 100644 --- a/src/main/java/org/broad/igv/sam/SparseAlignmentCounts.java +++ b/src/main/java/org/broad/igv/sam/SparseAlignmentCounts.java @@ -144,6 +144,32 @@ public int getTotalCount(int pos) { } } + @Override + public int getTotalPositiveCount(int pos) { + if (!indexMap.containsKey(pos)) { + if (log.isDebugEnabled()) { + log.debug("Position out of range: " + pos + " (valid range - " + start + "-" + end); + } + return 0; + } else { + int idx = getIndex(pos); + return getCountFromList(posTotal, idx); + } + } + + @Override + public int getTotalNegativeCount(int pos) { + if (!indexMap.containsKey(pos)) { + if (log.isDebugEnabled()) { + log.debug("Position out of range: " + pos + " (valid range - " + start + "-" + end); + } + return 0; + } else { + int idx = getIndex(pos); + return getCountFromList(negTotal, idx); + + } + } public int getTotalQuality(int pos) { if (!indexMap.containsKey(pos)) { diff --git a/src/main/java/org/broad/igv/sam/mods/BaseModficationFilter.java b/src/main/java/org/broad/igv/sam/mods/BaseModficationFilter.java new file mode 100644 index 0000000000..999501a8ca --- /dev/null +++ b/src/main/java/org/broad/igv/sam/mods/BaseModficationFilter.java @@ -0,0 +1,87 @@ +package org.broad.igv.sam.mods; + +import java.util.Map; +import java.util.stream.Collectors; +import java.util.stream.Stream; + +public class BaseModficationFilter { + + String modification; + char base = 0; + + public BaseModficationFilter(String modification, char base) { + this.modification = modification; + this.base = base; + } + + public BaseModficationFilter(String modification) { + this.modification = modification; + this.base = 0; + } + + public boolean pass(String mod) { + if(mod.startsWith("NONE_")) { + char b = mod.charAt(5); + return (this.base > 0 && this.base == b) || + (modToBaseMap.containsKey(this.modification.substring(0, 1)) && + modToBaseMap.get(this.modification.substring(0, 1)).charAt(0) == b); + } else { + return (modification == null || modification.contains(mod)); + } + } + + public boolean pass(String mod, char b) { + + if(this.base != 0 && this.base == b) { + return true; + } else { + return pass(mod); + } + } + + @Override + public String toString() { + return (modification == null ? "" : modification) + "," + (base == 0 ? "" : base); + } + + public static BaseModficationFilter fromString(String str) { + int idx = str.indexOf(","); + String mod; + char base; + if (idx < 0) { + // Backward compatibility + mod = str; + base = 0; + } else if (idx == 0) { + mod = null; + base = str.charAt(1); + } else { + mod = str.substring(0, idx); + base = idx == str.length() - 1 ? 0 : str.charAt(idx + 1); + } + return new BaseModficationFilter(mod, base); + + } + + public static Map modToBaseMap; + static { + modToBaseMap = Stream.of(new String[][]{ + {"m", "C"}, + {"h", "C"}, + {"f", "C"}, + {"c", "C"}, + {"C", "C"}, + {"g", "T"}, + {"e", "T"}, + {"b", "T"}, + {"T", "T"}, + {"a", "A"}, + {"A", "A"}, + {"o", "G"}, + {"G", "G"}, + {"n", "N"}, + {"N", "N"} + }).collect(Collectors.toMap(data -> data[0], data -> data[1])); + + } +} diff --git a/src/main/java/org/broad/igv/sam/mods/BaseModificationColors.java b/src/main/java/org/broad/igv/sam/mods/BaseModificationColors.java index 44e1343b88..6a1a2ae398 100644 --- a/src/main/java/org/broad/igv/sam/mods/BaseModificationColors.java +++ b/src/main/java/org/broad/igv/sam/mods/BaseModificationColors.java @@ -2,7 +2,6 @@ import org.broad.igv.logging.LogManager; import org.broad.igv.logging.Logger; -import org.broad.igv.prefs.PreferencesManager; import org.broad.igv.sam.AlignmentTrack; import java.awt.*; @@ -29,7 +28,7 @@ public class BaseModificationColors { static HashMap colors5MC = new HashMap<>(); static Color genericColor = new Color(132, 178, 158); - public static Color noModColor5MC = Color.blue; + public static Color noModColor = Color.blue; static Color hColor = new Color(11, 132, 165); static Color oColor = new Color(111, 78, 129); @@ -50,6 +49,10 @@ public class BaseModificationColors { colors.put("e", eColor); colors.put("b", bColor); colors.put("a", aColor); + colors.put("NONE_A", noModColor); + colors.put("NONE_C", noModColor); + colors.put("NONE_T", noModColor); + colors.put("NONE_G", noModColor); colors5MC.put("h", new Color(255, 0, 255)); // Modify h for 5mC to distinguish from blue } @@ -57,73 +60,29 @@ public class BaseModificationColors { * Cache for modified colors */ static Map modColorMap = new HashMap<>(); - static Map modColorMap5MC = new HashMap<>(); - public static Color getModColor(String modification, byte likelihood, AlignmentTrack.ColorOption colorOption) { + public static Color getModColor(String modification, int l, AlignmentTrack.ColorOption colorOption) { // Note the pallete will always return a color, either an initially seeded one if supplied or a random color. Color baseColor = getBaseColor(modification, colorOption); - int l = Byte.toUnsignedInt(likelihood); - if (l > 255) { + if (l > 210) { return baseColor; } String key = modification + "--" + l; - if (colorOption == AlignmentTrack.ColorOption.BASE_MODIFICATION_C) { - - if (!modColorMap5MC.containsKey(key)) { - int alpha = Math.min(255, (int) (l * l / 64f - 4 * l + 256)); // quadratic - if (l >= 128) { - modColorMap5MC.put(key, new Color(baseColor.getRed(), baseColor.getGreen(), baseColor.getBlue(), alpha)); - } else { - modColorMap5MC.put(key, new Color(noModColor5MC.getRed(), noModColor5MC.getGreen(), noModColor5MC.getBlue(), alpha)); - } - } - - return modColorMap5MC.get(key); - - } else { - if (l > 250) { - return baseColor; - } - double threshold = 256 * PreferencesManager.getPreferences().getAsFloat("SAM.BASEMOD_THRESHOLD"); - if (l < threshold) { - l = 0; - } - if (!modColorMap.containsKey(key)) { - modColorMap.put(key, new Color(baseColor.getRed(), baseColor.getGreen(), baseColor.getBlue(), l)); - } - return modColorMap.get(key); - } - } - - public static Color getNoModColor(byte likelihood) { - - // Note the pallete will always return a color, either an initially seeded one if supplied or a random color. - Color baseColor = noModColor5MC; - - int l = Byte.toUnsignedInt(likelihood); - if (l > 255) { - return baseColor; - } - - String key = "NOMOD--" + l; - - if (!modColorMap5MC.containsKey(key)) { - int alpha = Math.min(255, (int) (l * l / 64f - 4 * l + 256)); // quadratic - modColorMap5MC.put(key, new Color(baseColor.getRed(), baseColor.getGreen(), baseColor.getBlue(), alpha)); + if (!modColorMap.containsKey(key)) { + int alpha = Math.max(20, Math.min(255, (int) (l * l / 64f - 4 * l + 256))); + modColorMap.put(key, new Color(baseColor.getRed(), baseColor.getGreen(), baseColor.getBlue(), alpha)); } - - return modColorMap5MC.get(key); + return modColorMap.get(key); } + private static Color getBaseColor(String modification, AlignmentTrack.ColorOption colorOption) { - if ((colorOption == AlignmentTrack.ColorOption.BASE_MODIFICATION_C) && colors5MC.containsKey(modification)) { - return colors5MC.get(modification); - } else if (colors.containsKey(modification)) { + if (colors.containsKey(modification)) { return colors.get(modification); } else { return genericColor; diff --git a/src/main/java/org/broad/igv/sam/mods/BaseModificationCounts.java b/src/main/java/org/broad/igv/sam/mods/BaseModificationCounts.java index e885e636a1..e61846a74e 100644 --- a/src/main/java/org/broad/igv/sam/mods/BaseModificationCounts.java +++ b/src/main/java/org/broad/igv/sam/mods/BaseModificationCounts.java @@ -3,6 +3,7 @@ import org.broad.igv.sam.Alignment; import org.broad.igv.sam.AlignmentBlock; import org.broad.igv.sam.AlignmentTrack; +import org.broad.igv.util.collections.ByteArrayList; import java.util.*; @@ -17,21 +18,13 @@ public class BaseModificationCounts { */ LinkedHashSet allModifications; - /** - * Map for counts for each modification (e.g. m, h, etc). Key is base+modification identifier, value is map of position -> count - */ - Map> counts; - /** - * Map for capturing modification likelihood "pileup", key is modification identifier, - * value is map of position -> sum of likelihoods for modifications at that position - */ - Map> likelihoodSums; + Map> likelihoods; + public BaseModificationCounts() { allModifications = new LinkedHashSet<>(); - counts = new HashMap<>(); - likelihoodSums = new HashMap<>(); + likelihoods = new HashMap<>(); } /** @@ -53,97 +46,109 @@ public void incrementCounts(Alignment alignment) { // Loop through read sequence index ("i") for (int i = block.getBases().startOffset; i < block.getBases().startOffset + block.getBases().length; i++) { - for (BaseModificationSet bmset : baseModificationSets) { - - //String modification = bmset.getModification(); - BaseModificationKey key = BaseModificationKey.getKey(bmset.getBase(), bmset.getStrand(), bmset.getModification()); - Map likelihoods = bmset.getLikelihoods(); - - if (bmset.containsPosition(i)) { - - int lh = Byte.toUnsignedInt(likelihoods.get(i)); - - Map modCounts = counts.get(key); - if (modCounts == null) { - modCounts = new HashMap<>(); - counts.put(key, modCounts); - } - - Map modLikelihoods = likelihoodSums.get(key); - if (modLikelihoods == null) { - modLikelihoods = new HashMap<>(); - likelihoodSums.put(key, modLikelihoods); - } - - int blockIdx = i - block.getBases().startOffset; - int position = block.getStart() + blockIdx; // genomic position - - int c = modCounts.containsKey(position) ? modCounts.get(position) + 1 : 1; - - - int l = modLikelihoods.containsKey(position) ? modLikelihoods.get(position) + lh : lh; - modCounts.put(position, c); - modLikelihoods.put(position, l); - + int blockIdx = i - block.getBases().startOffset; + int position = block.getStart() + blockIdx; +if(position == 119094723) { + System.out.println(); +} + // Loop through base modification sets + + char canonicalBase = 0; + int noModLH = 255; + for (BaseModificationSet bmSet : baseModificationSets) { + final Map bmSetLikelihoods = bmSet.getLikelihoods(); + if (bmSetLikelihoods != null && bmSet.containsPosition(i)) { // TODO or flag == '.' and readbase = bmSet.canonicalbase + byte byteLikelihood = bmSetLikelihoods.get(i); + BaseModificationKey modKey = BaseModificationKey.getKey(bmSet.getBase(), bmSet.getStrand(), bmSet.getModification()); + allModifications.add(modKey); + pushLikelihood(position, byteLikelihood, modKey, likelihoods); + + canonicalBase = bmSet.getCanonicalBase(); // Assumed same for all modifications at this position, modificatons on both bases at a position not supported + noModLH -= Byte.toUnsignedInt(byteLikelihood); } - allModifications.add(key); } + + if (canonicalBase != 0) { + BaseModificationKey noModKey = BaseModificationKey.getKey(canonicalBase, '+', "NONE_" + canonicalBase); + allModifications.add(noModKey); + pushLikelihood(position, (byte) noModLH, noModKey, likelihoods); + } + } } } } - public int getCount(int position, BaseModificationKey key) { - Map modCounts = counts.get(key); - if (modCounts != null && modCounts.containsKey(position)) { - return modCounts.get(position); - } else { - return 0; + private void pushLikelihood(int position, byte byteLikelihood, BaseModificationKey modKey, Map> likelihoods) { + Map t = likelihoods.get(modKey); + if (t == null) { + t = new HashMap<>(); + likelihoods.put(modKey, t); + } + ByteArrayList byteArrayList = t.get(position); + if (byteArrayList == null) { + byteArrayList = new ByteArrayList(100); + t.put(position, byteArrayList); } + byteArrayList.add(byteLikelihood); } - public int getLikelhoodSum(int position, BaseModificationKey key) { - Map modLikelihoods = likelihoodSums.get(key); - if (modLikelihoods != null && modLikelihoods.containsKey(position)) { - return modLikelihoods.get(position); + public int getCount(int position, BaseModificationKey key, int threshold) { + Map t = likelihoods.get(key); + return _getCount(position, threshold, t); + } + + + private int _getCount(int position, int threshold, Map t) { + ByteArrayList byteArrayList = t.get(position); + if (byteArrayList == null) { + return 0; } else { - return getCount(position, key) * 255; + int count = 0; + for (int i = 0; i < byteArrayList.size(); i++) { + int lh = Byte.toUnsignedInt(byteArrayList.get(i)); + if (lh > threshold) { + count++; + } + } + return count; } } + public Set getAllModificationKeys() { return allModifications; } public String getValueString(int position, AlignmentTrack.ColorOption colorOption) { StringBuffer buffer = new StringBuffer(); - for (Map.Entry> entry : counts.entrySet()) { + + // /** + // * Map for capturing modification likelihood "pileup", key is modification identifier, + // * value is map of position -> sum of likelihoods for modifications at that position + // */ + // Map> likelihoodSums; + + + for (Map.Entry> entry : likelihoods.entrySet()) { BaseModificationKey key = entry.getKey(); - Map modCounts = entry.getValue(); - if (modCounts.containsKey(position)) { - final Integer count = modCounts.get(position); - int lh = (int) (((100.0f / 255) * getLikelhoodSum(position, entry.getKey())) / count); + Map t = entry.getValue(); + if (t.containsKey(position)) { + int count = this.getCount(position, key, 255 / 2); String modName = BaseModificationUtils.modificationName(key.modification); - buffer.append("Modification: " + modName + " (" + key.base + key.strand + ", " + count + " @ " + lh + "%)
"); + buffer.append("Modification: " + modName + " (" + key.base + key.strand + "), count > threshold: " + count + "
"); } } +// for (Map.Entry> entry : nomodLikelihoods.entrySet()) { +// BaseModificationKey key = entry.getKey(); +// Map t = entry.getValue(); +// if (t.containsKey(position)) { +// int count = this.getCount(position, key, 255 / 2); +// String modName = BaseModificationUtils.modificationName(key.modification); +// buffer.append("Modification: " + modName + " (" + key.base + key.strand + "), count > threshold: " + count + "
"); +// } +// +// } return buffer.toString(); } - - /** - * For debugging - */ - public void dump() { - for (Map.Entry> entry : counts.entrySet()) { - - String modification = entry.getKey().toString(); - Map modCounts = entry.getValue(); - System.out.println("Modification: " + modification); - for (Map.Entry modKey : modCounts.entrySet()) { - System.out.println(modKey.getKey() + " " + modKey.getValue()); - } - } - } - - } diff --git a/src/main/java/org/broad/igv/sam/mods/BaseModificationCoverageRenderer.java b/src/main/java/org/broad/igv/sam/mods/BaseModificationCoverageRenderer.java index a5ab891961..7df2f184c3 100644 --- a/src/main/java/org/broad/igv/sam/mods/BaseModificationCoverageRenderer.java +++ b/src/main/java/org/broad/igv/sam/mods/BaseModificationCoverageRenderer.java @@ -1,16 +1,16 @@ package org.broad.igv.sam.mods; +import htsjdk.samtools.util.SequenceUtil; import org.broad.igv.sam.AlignmentCounts; import org.broad.igv.sam.AlignmentTrack.ColorOption; import org.broad.igv.track.RenderContext; import java.awt.*; import java.util.*; +import java.util.List; public class BaseModificationCoverageRenderer { - - public static void drawModifications(RenderContext context, int pX, int pBottom, @@ -19,189 +19,53 @@ public static void drawModifications(RenderContext context, int pos, AlignmentCounts alignmentCounts, ColorOption colorOption, - String basemodFilter) { - - switch (colorOption) { - case BASE_MODIFICATION_C: - draw5MC(context, pX, pBottom, dX, barHeight, pos, alignmentCounts, basemodFilter); - break; - default: - draw(context, pX, pBottom, dX, barHeight, pos, alignmentCounts, colorOption, basemodFilter); - } - } + BaseModficationFilter filter, + Set simplexModifications) { - /* - chr11:119,094,758 -Total count: 26 -A : 0 -C : 26 (100%, 16+, 10- ) -G : 0 -T : 0 -N : 0Modification: m (C+, 15 @ 47%) - - */ - - private static void draw(RenderContext context, - int pX, - int pBottom, - int dX, - int barHeight, - int pos, - AlignmentCounts alignmentCounts, - ColorOption colorOption, - String filter) { - - boolean debug = pos == 119094762; + int modThreshold = 255 / 2; BaseModificationCounts modificationCounts = alignmentCounts.getModifiedBaseCounts(); if (modificationCounts != null) { - final BaseModificationKey mCKey = BaseModificationKey.getKey('C', '+', "m"); - final BaseModificationKey mGKey = BaseModificationKey.getKey('G', '-', "m"); + Set allModificationKeys = modificationCounts.getAllModificationKeys(); + List sortedKeys = new ArrayList<>(allModificationKeys); + Collections.sort(sortedKeys); + + // Color bar modification counts (# of modifications > threshold) + int total = alignmentCounts.getTotalCount(pos); + + for (BaseModificationKey key : sortedKeys) { - Graphics2D graphics = context.getGraphics(); + if (filter != null && !filter.pass(key.modification, key.getCanonicalBase())) continue; + if(key.modification.startsWith("NONE_") && colorOption != ColorOption.BASE_MODIFICATION_2COLOR) continue; - Set allModificationKeys = modificationCounts.getAllModificationKeys(); - boolean cpgMode = allModificationKeys.contains(mCKey) && !allModificationKeys.contains(mGKey); + byte base = (byte) key.getBase(); + final byte compl = SequenceUtil.complement(base); + int modifiable = alignmentCounts.getCount(pos, base) + alignmentCounts.getCount(pos, compl); + int detectable = simplexModifications.contains(key.modification) ? + alignmentCounts.getPosCount(pos, base) + alignmentCounts.getNegCount(pos, compl) : + modifiable; - // Merge complementary sets (same modification, opposite strands) - Map likelihoodSums = new LinkedHashMap<>(); - for (BaseModificationKey key : allModificationKeys) { - String modification = key.getModification(); - if(filter != null && !filter.equals(modification)) continue; - float currentCount = likelihoodSums.containsKey(modification) ? likelihoodSums.get(modification) : 0; - likelihoodSums.put(modification, currentCount + modificationCounts.getLikelhoodSum(pos, key)); - } + if (detectable == 0) continue; //No informative reads - // Color bar by likelihood weighted count - int total = alignmentCounts.getTotalCount(pos); + int count = modificationCounts.getCount(pos, key, modThreshold); + float modFraction = (((float) modifiable) / total) * (((float) count) / detectable); + int modHeight = Math.round(modFraction * barHeight); + + int baseY = pBottom - modHeight; + Color modColor = BaseModificationColors.getModColor(key.modification, 255, colorOption); + Graphics2D graphics = context.getGraphics(); + graphics.setColor(modColor); + graphics.fillRect(pX, baseY, dX, modHeight); + pBottom = baseY; - int sumModHeight = 0; - for (String modification : likelihoodSums.keySet()) { - - if (likelihoodSums.get(modification) > 0) { - - float modFraction; - if (cpgMode & "m".equals(modification)) { - // Special mode for out-of-spec 5mC CpG convention. Calls are recorded for C+ only. - // We adjust the height of the bar to account for the missing G- calls. - final byte base = (byte) 'C'; - final byte compl = (byte) 'G'; - final int posCount = alignmentCounts.getPosCount(pos, base); - final int negCount = alignmentCounts.getNegCount(pos, compl); - - // Is this a "C" or "G" reference site? We want to determine this from the counts data - // directly, this should work for all but pathological edge cases - final int cCounts = alignmentCounts.getCount(pos, base); - final int gCounts = alignmentCounts.getCount(pos, compl); - final boolean cSite = cCounts > gCounts; - final int referenceCounts = cSite ? cCounts : gCounts; - final int strandCount = cSite ? posCount : negCount; - modFraction = (((float) referenceCounts) / total) * (likelihoodSums.get(modification) / (strandCount * 255f)); - } else { - modFraction = likelihoodSums.get(modification) / (total * 255f); - } - int modHeight = Math.round(modFraction * barHeight); - - int baseY = pBottom - modHeight; - Color modColor = BaseModificationColors.getModColor(modification, (byte) 255, colorOption); - graphics.setColor(modColor); - graphics.fillRect(pX, baseY, dX, modHeight); - pBottom = baseY; - - if(debug) System.out.println("draw " + modification + " " + modHeight); - } } } } - private static void draw5MC(RenderContext context, - int pX, - int pBottom, - int dX, - int barHeight, - int pos, - AlignmentCounts alignmentCounts, - String basemodFilter) { + static Map modRankOrder; - boolean debug = pos == 119094768; - BaseModificationCounts modificationCounts = alignmentCounts.getModifiedBaseCounts(); - - if (modificationCounts != null) { - - final byte base = (byte) 'C'; - final byte compl = (byte) 'G'; - - Map likelihoodSums = new HashMap<>(); - Map modCounts = new HashMap<>(); - for (BaseModificationKey key : modificationCounts.getAllModificationKeys()) { - - // This coloring mode is exclusively for "C" modifications - if (key.getCanonicalBase() != 'C') continue; - if (key.getModification().equals(basemodFilter) || basemodFilter == null) { - String mod = key.getModification(); - final int count = modificationCounts.getCount(pos, key); - if (count > 0) { - modCounts.put(mod, count); - final int likelhoodSum = modificationCounts.getLikelhoodSum(pos, key); - likelihoodSums.put(mod, likelhoodSum); - } - } - } - if (likelihoodSums.size() > 0) { - - // Special mode for out-of-spec 5mC CpG convention. - // Calls are made for the CG dinucleotide and only recorded on 1 strand. We adjust the height - // of the bar to account for the missing G- calls. This is an approximation and assumes the - // distribution of calls is ~ equal on both strands. - - // Is this a "C" or "G" reference site? We want to determine this from the counts data directly, - // this should work for all but pathological edge cases - final int cCounts = alignmentCounts.getCount(pos, base); - final int gCounts = alignmentCounts.getCount(pos, compl); - final boolean cSite = cCounts > gCounts; - final int referenceCounts = cSite ? cCounts : gCounts; - - // Compute "snp factor", ratio of count of base calls that could be modified to - // total count. This is normally close to 1, but can be less due non CG bases at this location (e.g. snps) - double snpFactor = ((double) referenceCounts) / alignmentCounts.getTotalCount(pos); - - double calledBarHeight = snpFactor * barHeight; - double t = referenceCounts * 255; // If all bases are called this is the total sum of all likelihoods, including "no mod" likelihood - - // Likelihood of no modification. It is assumed that the likelihood of no modification == (1 - sum(likelihood)) - //int c = Collections.max(modCounts.values()); - double noModProb = t; - for (String m : modCounts.keySet()) { - if (likelihoodSums.containsKey(m)) { - noModProb -= likelihoodSums.get(m); - } - } - - // Draw "no mod" bar - int noModHeight = (int) Math.round((noModProb / t) * calledBarHeight); - int baseY = pBottom - noModHeight; - Graphics2D graphics = context.getGraphics(); - graphics.setColor(BaseModificationColors.noModColor5MC); - graphics.fillRect(pX, baseY, dX, noModHeight); - - // Loop through modifications drawing bar for each - String[] orderedMods = likelihoodSums.keySet().toArray(new String[0]); - Arrays.sort(orderedMods, (o1, o2) -> -1 * o1.compareTo(o2)); - for (String m : orderedMods) { - Color mColor = BaseModificationColors.getModColor(m, (byte) 255, ColorOption.BASE_MODIFICATION_C); - int mModHeight = (int) Math.round(((likelihoodSums.get(m)) / t) * calledBarHeight); - if (mModHeight > 0) { - baseY -= mModHeight; - graphics.setColor(mColor); - graphics.fillRect(pX, baseY, dX, mModHeight); - } - if(debug) System.out.println("draw5mC " + m + " " + mModHeight); - } - } - } - } } diff --git a/src/main/java/org/broad/igv/sam/mods/BaseModificationKey.java b/src/main/java/org/broad/igv/sam/mods/BaseModificationKey.java index a83a3c0480..5a6cbbf3af 100644 --- a/src/main/java/org/broad/igv/sam/mods/BaseModificationKey.java +++ b/src/main/java/org/broad/igv/sam/mods/BaseModificationKey.java @@ -7,7 +7,7 @@ import java.util.Map; import java.util.Objects; -public class BaseModificationKey { +public class BaseModificationKey implements Comparable { char base; char strand; String modification; @@ -24,7 +24,10 @@ public static BaseModificationKey getKey(char base, char strand, String modifica keyCache.put(name, key); return key; } + } + public static BaseModificationKey getNomodKey(char base) { + return getKey(base, '+', "NONE"); } private BaseModificationKey(char base, char strand, String modification) { @@ -66,4 +69,39 @@ public int hashCode() { public String toString() { return "" + base + strand + modification; } + + + static Map modificationRankOrder ; + + @Override + public int compareTo(Object o) { + BaseModificationKey otherKey = (BaseModificationKey) o; + String mod1 = this.modification; + String mod2 = otherKey.modification; + + if(mod1.equals(mod2)) { + return (byte) this.strand - (byte) otherKey.strand; + } + + if(modificationRankOrder == null) { + modificationRankOrder = new HashMap<>(); + String [] tmp = {"NONE_C", "NONE_T", "NONE_G", "NONE_A", "m", "h", "f", "c", "C", "g", "e", "b", "T", "U", "a", "A", "o", "G", "n", "N"}; + for(int idx = 0; idx < tmp.length; idx++) { + modificationRankOrder.put(tmp[idx], idx); + } + } + if (modificationRankOrder.containsKey(mod1) & modificationRankOrder.containsKey(mod2)) { + return modificationRankOrder.get(mod1) - modificationRankOrder.get(mod2); + } else if (modificationRankOrder.containsKey(mod1)) { + return 1; + } else if (modificationRankOrder.containsKey(mod2)) { + return -1; + } else { + return mod1.compareTo(mod2); + } + } + + + + } diff --git a/src/main/java/org/broad/igv/sam/mods/BaseModificationRenderer.java b/src/main/java/org/broad/igv/sam/mods/BaseModificationRenderer.java index 36052073ff..da44d8f818 100644 --- a/src/main/java/org/broad/igv/sam/mods/BaseModificationRenderer.java +++ b/src/main/java/org/broad/igv/sam/mods/BaseModificationRenderer.java @@ -16,34 +16,7 @@ public static void drawModifications( Rectangle rowRect, Graphics g, AlignmentTrack.ColorOption colorOption, - String basemodFilter) { - - switch (colorOption) { - case BASE_MODIFICATION_C: - draw5mC(alignment, bpStart, locScale, rowRect, g, basemodFilter); - break; - default: - draw(alignment, bpStart, locScale, rowRect, g, basemodFilter); - } - - } - - /** - * Helper function for AlignmentRenderer. Draw base modifications over alignment. - * - * @param alignment - * @param bpStart - * @param locScale - * @param rowRect - * @param g - */ - private static void draw( - Alignment alignment, - double bpStart, - double locScale, - Rectangle rowRect, - Graphics g, - String filter) { + BaseModficationFilter filter) { List baseModificationSets = alignment.getBaseModificationSets(); @@ -66,123 +39,49 @@ private static void draw( continue; } - // Search all sets for modifications of this base. For now keeps mod with > probability - // TODO -- merge mods in some way - byte lh = 0; + // Search all sets for modifications of this base. Pick modification with largest likelhiood + int maxLh = 0; + int noModLh = 255; String modification = null; + char canonicalBase = 0; + boolean modificationFound = false; for (BaseModificationSet bmSet : baseModificationSets) { - if (filter == null || filter.equals(bmSet.getModification())) { - if (bmSet.containsPosition(i)) { - if (modification == null || Byte.toUnsignedInt(bmSet.getLikelihoods().get(i)) > Byte.toUnsignedInt(lh)) { - modification = bmSet.getModification(); - lh = bmSet.getLikelihoods().get(i); - } + if (bmSet.containsPosition(i)) { + int lh = Byte.toUnsignedInt(bmSet.getLikelihoods().get(i)); + // TODO ELSE if flag == '.' lh = 0 + noModLh -= lh; + if (modification == null || lh > maxLh) { + modification = bmSet.getModification(); + canonicalBase = bmSet.getCanonicalBase(); + maxLh = lh; + modificationFound = modificationFound || filter == null || filter.pass(modification, canonicalBase); } } } - if (modification != null) { - - Color c = BaseModificationColors.getModColor(modification, lh, AlignmentTrack.ColorOption.BASE_MODIFICATION); - g.setColor(c); - - // Expand narrow width to make more visible - if (dX < 3) { - dX = 3; - pX--; + if(modificationFound) { + Color c = null; + if (noModLh > 127 && colorOption == AlignmentTrack.ColorOption.BASE_MODIFICATION_2COLOR) { + c = BaseModificationColors.getModColor("NONE_" + canonicalBase, noModLh, colorOption); + } else if (modification != null && maxLh > 127 && + (filter == null || filter.pass(modification, canonicalBase))) /* || flag == . */ { + c = BaseModificationColors.getModColor(modification, maxLh, colorOption); } - g.fillRect(pX, pY, dX, Math.max(1, dY - 2)); - } - } - } - } - } - - /** - * Helper function for AlignmentRenderer. Draw base modifications over alignment for "5mC" mode. - *

- * Notes: - * Designed primarily for visualization of 5mC modifications compatible with existing bisulfite seq viz - * - 5mC methylated bases colored red - * - Non modified bases colored blue - * - Other modificationc colored as defined in BaseModificationColors - *

- * If multiple modifications are specified for a base the modification with the highest probability is - * drawn. - * - * @param alignment - * @param bpStart - * @param locScale - * @param rowRect - * @param g - */ - private static void draw5mC( - Alignment alignment, - double bpStart, - double locScale, - Rectangle rowRect, - Graphics g, - String basemodFilter) { + if (c != null) { + g.setColor(c); - List baseModificationSets = alignment.getBaseModificationSets(); - if (baseModificationSets != null) { - - for (AlignmentBlock block : alignment.getAlignmentBlocks()) { - // Compute bounds - int pY = (int) rowRect.getY(); - int dY = (int) rowRect.getHeight(); - int dX = (int) Math.max(1, (1.0 / locScale)); - - for (int i = block.getBases().startOffset; i < block.getBases().startOffset + block.getBases().length; i++) { - - int blockIdx = i - block.getBases().startOffset; - int pX = (int) ((block.getStart() + blockIdx - bpStart) / locScale); - - // Don't draw out of clipping rect - if (pX > rowRect.getMaxX()) { - break; - } else if (pX + dX < rowRect.getX()) { - continue; - } - - // Search all sets for modifications of this base, select modification with largest likelihood - int lh = -1; - String modification = null; - - // Compare likelihoods, including likelihood of no modification - int noModificationLikelihood = 255; - for (BaseModificationSet bmSet : baseModificationSets) { - - // This coloring mode is exclusively for "C" modifications, either 5mC or all C mods - if (bmSet.getCanonicalBase() != 'C') continue; - if (bmSet.getModification().equals(basemodFilter) || basemodFilter == null) { - if (bmSet.containsPosition(i)) { - int l = Byte.toUnsignedInt(bmSet.getLikelihoods().get(i)); - noModificationLikelihood -= l; - if (modification == null || l > lh) { - modification = bmSet.getModification(); - lh = l; - } + // Expand narrow width to make more visible + if (dX < 3) { + dX = 3; + pX--; } + g.fillRect(pX, pY, dX, Math.max(1, dY - 2)); } } - - if (modification != null) { - - Color c = noModificationLikelihood > lh ? - BaseModificationColors.getNoModColor((byte) noModificationLikelihood) : - BaseModificationColors.getModColor(modification, (byte) lh, AlignmentTrack.ColorOption.BASE_MODIFICATION_C); - g.setColor(c); - - // Expand narrow width to make more visible - if (dX < 3) { - dX = 3; - pX--; - } - g.fillRect(pX, pY, dX, Math.max(1, dY - 2)); - } } } } } } + + diff --git a/src/main/java/org/broad/igv/sam/mods/BaseModificationSet.java b/src/main/java/org/broad/igv/sam/mods/BaseModificationSet.java index 8be28f8028..6f5c37ee74 100644 --- a/src/main/java/org/broad/igv/sam/mods/BaseModificationSet.java +++ b/src/main/java/org/broad/igv/sam/mods/BaseModificationSet.java @@ -9,15 +9,16 @@ public class BaseModificationSet { char base; char strand; - String modification; Map likelihoods; + char canonicalBase; public BaseModificationSet(char base, char strand, String modification, Map likelihoods) { this.base = base; this.modification = modification; this.strand = strand; this.likelihoods = likelihoods; + this.canonicalBase = strand == '+' ? base : (char) SequenceUtil.complement((byte) base); } public char getBase() { @@ -25,7 +26,7 @@ public char getBase() { } public char getCanonicalBase() { - return strand == '+' ? base : (char) SequenceUtil.complement((byte) base); + return canonicalBase; } public String getModification() { diff --git a/src/main/java/org/broad/igv/util/collections/ByteArrayList.java b/src/main/java/org/broad/igv/util/collections/ByteArrayList.java new file mode 100644 index 0000000000..2a0fbb63eb --- /dev/null +++ b/src/main/java/org/broad/igv/util/collections/ByteArrayList.java @@ -0,0 +1,134 @@ +/* + * The MIT License (MIT) + * + * Copyright (c) 2007-2015 Broad Institute + * + * Permission is hereby granted, free of charge, to any person obtaining a copy + * of this software and associated documentation files (the "Software"), to deal + * in the Software without restriction, including without limitation the rights + * to use, copy, modify, merge, publish, distribute, sublicense, and/or sell + * copies of the Software, and to permit persons to whom the Software is + * furnished to do so, subject to the following conditions: + * + * The above copyright notice and this permission notice shall be included in + * all copies or substantial portions of the Software. + * + * + * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR + * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, + * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE + * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER + * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, + * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN + * THE SOFTWARE. + */ + +package org.broad.igv.util.collections; + + +/** + * Author: jrobinso + * Date: Jul 22, 2010 + *

+ * ArrayList type collection for int types. Purpose is to avoid the need to create an object for each entry + * in the standard java collections. + */ +public class ByteArrayList { + + + private transient byte[] elements; + + private int size; + + + public ByteArrayList() { + this(100); + } + + public ByteArrayList(int initialCapacity) { + if (initialCapacity < 0) + throw new IllegalArgumentException("Illegal Capacity: " + initialCapacity); + this.elements = new byte[initialCapacity]; + } + + public ByteArrayList(byte[] elements) { + this.elements = elements; + size = elements.length; + } + + public void add(byte e) { + if (size + 1 >= elements.length) { + grow(); + } + elements[size++] = e; + } + + public void addAll(byte[] args) { + byte[] newElements = new byte[size + args.length]; + System.arraycopy(elements, 0, newElements, 0, size); + System.arraycopy(args, 0, newElements, size, args.length); + elements = newElements; + size += args.length; + } + + public void addAll(ByteArrayList aList) { + addAll(aList.toArray()); + } + + public byte get(int idx) { + return elements[idx]; + } + + public int size() { + return size; + } + + public boolean isEmpty() { + return size == 0; + } + + /** + * Empty all elements. This logically clears the collection but does not free up any space. + */ + public void clear() { + size = 0; + } + + private void grow() { + int oldCapacity = elements.length; + int newCapacity; + if (oldCapacity < 10000000) { + newCapacity = oldCapacity * 2; + } else { + newCapacity = (oldCapacity * 3) / 2 + 1; + } + byte[] tmp = new byte[newCapacity]; + System.arraycopy(elements, 0, tmp, 0, elements.length); + elements = tmp; + } + + + public byte[] toArray() { + trimToSize(); + return elements; + } + + + private void trimToSize() { + int oldCapacity = elements.length; + if (size < oldCapacity) { + byte[] tmp = new byte[size]; + System.arraycopy(elements, 0, tmp, 0, size); + elements = tmp; + } + } + + public void set(int idx, byte i) { + while(idx >= elements.length) { + grow(); + } + elements[idx] = i; + idx++; + if(idx > size) size = idx; // Tried Math.max here, it showed up in cpu profiles! + } +} diff --git a/src/main/java/org/broad/igv/util/collections/IntArrayList.java b/src/main/java/org/broad/igv/util/collections/IntArrayList.java index 0235c4d4a5..a63a1452b9 100644 --- a/src/main/java/org/broad/igv/util/collections/IntArrayList.java +++ b/src/main/java/org/broad/igv/util/collections/IntArrayList.java @@ -25,7 +25,6 @@ package org.broad.igv.util.collections; -import org.broad.igv.util.IntComparator; /** * Author: jrobinso diff --git a/src/test/java/org/broad/igv/sam/mods/BaseModficationFilterTest.java b/src/test/java/org/broad/igv/sam/mods/BaseModficationFilterTest.java new file mode 100644 index 0000000000..aa59f13558 --- /dev/null +++ b/src/test/java/org/broad/igv/sam/mods/BaseModficationFilterTest.java @@ -0,0 +1,32 @@ +package org.broad.igv.sam.mods; + +import org.junit.Test; + +import static org.junit.Assert.*; + +public class BaseModficationFilterTest { + + @Test + public void pass() { + BaseModficationFilter filter = new BaseModficationFilter(null, 'C'); + assertTrue(filter.pass(null, 'C')); + + String nomod = "NONE_C"; + assertTrue(filter.pass(nomod)); + } + + @Test + public void fromString() { + + String str = "hm,C"; + BaseModficationFilter filter = BaseModficationFilter.fromString(str); + + assertEquals("hm", filter.modification); + assertEquals('C', filter.base); + + assertTrue(filter.pass("m")); + assertTrue(filter.pass("h")); + assertTrue(filter.pass("NONE_C")); + assertFalse(filter.pass("a")); + } +} \ No newline at end of file diff --git a/src/test/java/org/broad/igv/sam/mods/BaseModificationCountsTest.java b/src/test/java/org/broad/igv/sam/mods/BaseModificationCountsTest.java index 54495ae2f9..e26b22c0f0 100644 --- a/src/test/java/org/broad/igv/sam/mods/BaseModificationCountsTest.java +++ b/src/test/java/org/broad/igv/sam/mods/BaseModificationCountsTest.java @@ -41,11 +41,40 @@ public void incrementCounts() throws IOException { BaseModificationKey key = BaseModificationKey.getKey('C', '+', "m"); for(int i=0; i bamiter = bamreader.query(chr, start, end, false); + int readCount = 0; + + BaseModificationCounts counts = new BaseModificationCounts(); + while (bamiter.hasNext()) { + Alignment alignment = bamiter.next(); + counts.incrementCounts(alignment); + readCount++; + } + assertTrue("No data retrieved: " + readCount, readCount > 0); - // counts.dump(); + BaseModificationKey cmKey = BaseModificationKey.getKey('C', '+', "m"); + int aboveThreshold = counts.getCount(119094723, cmKey, 127); + assertEquals("Counts above threshold", 3, aboveThreshold); + BaseModificationKey noModKey = BaseModificationKey.getKey('C', '+', "NONE"); + int belowThreshold = counts.getCount( 119094723, noModKey, 127); + assertEquals("Counts below threshold", 12, belowThreshold); } } \ No newline at end of file diff --git a/src/test/java/org/broad/igv/sam/mods/BaseModificationsTest.java b/src/test/java/org/broad/igv/sam/mods/BaseModificationsTest.java index e05ac4125b..84133d2cf7 100644 --- a/src/test/java/org/broad/igv/sam/mods/BaseModificationsTest.java +++ b/src/test/java/org/broad/igv/sam/mods/BaseModificationsTest.java @@ -70,6 +70,7 @@ public void testOrientTopRev() { // ....45..........6.......4...8.0..... //top-rev 16 * 0 0 * * 0 0 ATATGGCATATCCCCCGCCGATCCGCTAGAGATCCT * Mm:Z:C+m,1,3,0; Ml:B:C,128,153,179 + //sequence actually read TATACCGTATAGGGGGCGGCTAGGCGATCTCTAGGT <= read direction byte[] sequence = "ATATGGCATATCCCCCGCCGATCCGCTAGAGATCCT".getBytes(); String mm = "C+m,1,3,0"; byte[] ml = {(byte) 128, (byte) 153, (byte) 179};