From 26ff7f2368a7798adfb26b1a69460afb381488c2 Mon Sep 17 00:00:00 2001 From: Daniel Geiszler Date: Mon, 27 May 2024 17:17:40 +0300 Subject: [PATCH] Peptide.calculatePeptideFragmentsBetween method fix Bump to RC7 Implement intensity-dependent NPV --- build.gradle | 2 +- .../andykong/ptmshepherd/PTMShepherd.java | 2 +- .../IterativeLocalizer.java | 50 +++++++++----- .../MatchedIonDistribution.java | 39 ++++++++--- .../andykong/ptmshepherd/utils/Peptide.java | 66 ++++++++++++++++++- test/utils/PeptideTest.java | 57 ++++++++++++++++ 6 files changed, 185 insertions(+), 31 deletions(-) diff --git a/build.gradle b/build.gradle index d9dfabd..1bef502 100644 --- a/build.gradle +++ b/build.gradle @@ -30,7 +30,7 @@ java { targetCompatibility = JavaVersion.VERSION_1_9 } -version = '3.0.0-rc6' +version = '3.0.0-rc7' application { // Define the main class for the application diff --git a/src/edu/umich/andykong/ptmshepherd/PTMShepherd.java b/src/edu/umich/andykong/ptmshepherd/PTMShepherd.java index 6968237..d8959a4 100644 --- a/src/edu/umich/andykong/ptmshepherd/PTMShepherd.java +++ b/src/edu/umich/andykong/ptmshepherd/PTMShepherd.java @@ -70,7 +70,7 @@ public class PTMShepherd { public static final String name = "PTM-Shepherd"; - public static final String version = "3.0.0-rc6"; + public static final String version = "3.0.0-rc7"; static HashMap params; static TreeMap> datasets; diff --git a/src/edu/umich/andykong/ptmshepherd/iterativelocalization/IterativeLocalizer.java b/src/edu/umich/andykong/ptmshepherd/iterativelocalization/IterativeLocalizer.java index c5359c5..fd058ae 100644 --- a/src/edu/umich/andykong/ptmshepherd/iterativelocalization/IterativeLocalizer.java +++ b/src/edu/umich/andykong/ptmshepherd/iterativelocalization/IterativeLocalizer.java @@ -159,44 +159,57 @@ private void fitMatchedIonDistribution() throws Exception { int mutationSite = decoyPep.mutatedResidue; // Calculate fake PTM mass shift - float massShift = AAMasses.monoisotopic_masses[targetPep.pepSeq.charAt(mutationSite) - 'A'] - + float mutatedMassShift = AAMasses.monoisotopic_masses[targetPep.pepSeq.charAt(mutationSite) - 'A'] - AAMasses.monoisotopic_masses[decoyPep.pepSeq.charAt(mutationSite) - 'A']; // Generate site-determining ions //XXXQQQ // Get set of reduced ions that match at least one configuration - float[][] reducedIons = getReducedIons(spec, decoyPep, massShift); + float[][] reducedIons = getReducedIons(spec, decoyPep, mutatedMassShift); float[] reducedMzs = reducedIons[0]; float[] reducedInts = reducedIons[1]; // Get sites that are/would be shifted - //ArrayList sitePepFrags = DownstreamPepFragGenerator.calculatePeptideFragments( - // targetPep, this.ionTypes, mutationSite, 1); - //ArrayList decoySitePepFrags = DownstreamPepFragGenerator.calculatePeptideFragments( - // decoyPep, this.ionTypes, mutationSite, 1); + ArrayList sitePepFrags = DownstreamPepFragGenerator.calculatePeptideFragments( + targetPep, this.ionTypes, mutationSite, 1); + ArrayList decoySitePepFrags = DownstreamPepFragGenerator.calculatePeptideFragments( + decoyPep, this.ionTypes, mutationSite, 1); // Generate target and decoy fragments - ArrayList sitePepFrags = targetPep.calculatePeptideFragments(this.ionTypes, 1); - ArrayList decoySitePepFrags = decoyPep.calculatePeptideFragments(this.ionTypes, 1); + //ArrayList sitePepFrags = targetPep.calculatePeptideFragments(this.ionTypes, 1); + //ArrayList decoySitePepFrags = decoyPep.calculateComplementaryFragments(this.ionTypes, mutatedMassShift, mutationSite, 1); // Score over all possible sites /** + System.out.println("Decoy\t" + decoyPep.pepSeq + "\t" + Arrays.toString(decoyPep.mods)); + System.out.println("Target\t" + targetPep.pepSeq + "\t" + Arrays.toString(targetPep.mods)); + System.out.println(mutatedMassShift); for (int k = 0; k < decoyPep.pepSeq.length(); k++) { - decoyPep.addMod(dMass, k); - ArrayList sitePepFrags = decoyPep.calculatePeptideFragments(this.ionTypes, 2); - float[][] matchedIons = findMatchedIons(sitePepFrags, reducedMzs, reducedInts); + if (k == mutationSite) + continue; + // Add decoy ions + decoyPep.addMod(mutatedMassShift, k); + ArrayList decoySitePepFrags = decoyPep.calculatePeptideFragmentsBetween(this.ionTypes, k, mutationSite, 1); + System.out.println("Decoy\t" + k + "\t" + mutationSite + "\t" + decoySitePepFrags); + float[][] matchedIons = findMatchedIons(decoySitePepFrags, reducedMzs, reducedInts); float[] matchedIonIntensities = matchedIons[0]; + System.out.println("Decoy\t" + k + "\t" + mutationSite + "\t" + Arrays.toString(matchedIonIntensities)); float[] matchedIonMassErrors = matchedIons[1]; - if (k == mutationSite) { - this.matchedIonDist.addIons(matchedIonIntensities, matchedIonMassErrors, false); - } else { - this.matchedIonDist.addIons(matchedIonIntensities, matchedIonMassErrors, true); - } - decoyPep.addMod(-1*dMass, k); + this.matchedIonDist.addIons(matchedIonIntensities, matchedIonMassErrors, true); + decoyPep.addMod(-1*mutatedMassShift, k); + + // Add target ions + ArrayList sitePepFrags = targetPep.calculatePeptideFragmentsBetween(this.ionTypes, k, mutationSite, 1); + System.out.println("Target\t" + k + "\t" + mutationSite + "\t" + sitePepFrags); + matchedIons = findMatchedIons(sitePepFrags, reducedMzs, reducedInts); + matchedIonIntensities = matchedIons[0]; + System.out.println("Target\t" + k + "\t" + mutationSite + "\t" + Arrays.toString(matchedIonIntensities)); + matchedIonMassErrors = matchedIons[1]; + this.matchedIonDist.addIons(matchedIonIntensities, matchedIonMassErrors, false); } **/ - // Add target peptide values to matched ion histograms + // Add target peptide values to matched ion histograms float[][] matchedIons = findMatchedIons(sitePepFrags, reducedMzs, reducedInts); float[] matchedIonIntensities = matchedIons[0]; float[] matchedIonMassErrors = matchedIons[1]; @@ -207,6 +220,7 @@ private void fitMatchedIonDistribution() throws Exception { float[] decoyMatchedIonIntensities = decoyMatchedIons[0]; float[] decoyMatchedMassErrors = decoyMatchedIons[1]; this.matchedIonDist.addIons(decoyMatchedIonIntensities, decoyMatchedMassErrors, true); + } } } diff --git a/src/edu/umich/andykong/ptmshepherd/iterativelocalization/MatchedIonDistribution.java b/src/edu/umich/andykong/ptmshepherd/iterativelocalization/MatchedIonDistribution.java index aba1905..64250be 100644 --- a/src/edu/umich/andykong/ptmshepherd/iterativelocalization/MatchedIonDistribution.java +++ b/src/edu/umich/andykong/ptmshepherd/iterativelocalization/MatchedIonDistribution.java @@ -133,12 +133,16 @@ public void addIon(float intensity, float massError, boolean isDecoy) { if (intensity > 0.0f) { // Only include matched ions, unmatched ions have negative intensity this.datapoints.addRow(massError, intensity, isDecoy); } else { + int intensityIndex = (int) (-1 * intensity); + this.pdf[intensityIndex]++; this.tNegCount++; } } else { if (intensity > 0.0f) { // Only include matched ions, unmatched ions have negative intensity this.datapoints.addRow(massError, intensity, isDecoy); } else { + int intensityIndex = (int) (-1 * intensity); + this.pdfDecoy[intensityIndex]++; this.dNegCount++; } } @@ -190,7 +194,7 @@ public void calculateIonPosterior() { // Sort table and find min and max this.datapoints.sortRowsByProjVal(); // Set up arrays to hold ion counts for targets and decoys and q-values - int arraySize = (int) (100.0 * (Math.abs(this.datapoints.getMaxProjVal() - + int arraySize = (int) (10.0 * (Math.abs(this.datapoints.getMaxProjVal() - this.datapoints.getMinProjVal()) + 1)); this.projTargetCounts = new int[arraySize]; this.projDecoyCounts = new int[arraySize]; @@ -208,17 +212,21 @@ public void calculateIonPosterior() { boolean leftToRight = checkLdaProjectionLeftToRight(); // Fill q-value array with P(ion == true_positive | LDA score) - int nTargets = 0; + int nTargets = 1; int nDecoys = 1; double cMin = 10000.0; + double cmax = -10000.0; if (leftToRight) { // Fill array for (int i = 0; i < arraySize; i++) { - nTargets += this.projTargetCounts[i]; - nDecoys += this.projDecoyCounts[i]; - this.qVals[i] = (double) nDecoys / (double) (nTargets); + nTargets = this.projTargetCounts[i]; + nDecoys = this.projDecoyCounts[i]; + double q = (double) (nDecoys + 1) / (double) (nDecoys + nTargets + 2); + this.qVals[i] = q; + System.out.println(this.qVals[i]); } // Make array monotonic + /** for (int i = arraySize-1; i >= 0; i--) { if (this.qVals[i] < cMin) { cMin = this.qVals[i]; @@ -230,13 +238,17 @@ public void calculateIonPosterior() { i++; } } + **/ } else { for (int i = arraySize-1; i >= 0; i--) { - nTargets += this.projTargetCounts[i]; - nDecoys += this.projDecoyCounts[i]; - this.qVals[i] = (double) nDecoys / (double) (nTargets); + nTargets = this.projTargetCounts[i]; + nDecoys = this.projDecoyCounts[i]; + double q = (double) (nDecoys + 1) / (double) (nDecoys + nTargets + 2); + this.qVals[i] = q; + System.out.println(this.qVals[i]); } // Make array monotonic + /** for (int i = 0; i < arraySize; i++) { if (this.qVals[i] < cMin) { cMin = this.qVals[i]; @@ -248,6 +260,7 @@ public void calculateIonPosterior() { i--; } } + **/ } } @@ -281,6 +294,12 @@ public void calculateLdaWeights() throws Exception { public void calculateNegativePredictiveValue() { // todo see if adding 1 here helps this.negPredictiveValue = (double) (this.tNegCount + 1) / (double) (this.dNegCount); + + this.ionPosterior = new double[((int) (100.0 * this.binMult + 1))]; + for (int i = 0; i < this.ionPosterior.length; i++) { + this.ionPosterior[i] = (double) ((this.pdf[i] + 1) / (double) (this.pdf[i] + this.pdfDecoy[i] + 2)); + System.out.println(this.ionPosterior[i]); + } //this.negPredictiveValue = Math.max(this.negPredictiveValue, 0.01); } @@ -388,7 +407,7 @@ public double calcIonProbability(float intensity) { } private int translateLdaValToIndex(double projVal) { - int valIndex = (int) ((projVal - this.datapoints.getMinProjVal()) * 100.0); + int valIndex = (int) ((projVal - this.datapoints.getMinProjVal()) * 10.0); if (valIndex < 0) return 0; else if (valIndex >= this.qVals.length) @@ -424,7 +443,7 @@ public double calculateIonProbabilityLda(float intensity, float massError) { double projVal; double prob; if (intensity < 0) - prob = this.negPredictiveValue; + prob = this.ionPosterior[(int) (intensity * -1 * this.binMult)]; else { projVal = this.ldaProcessor.projectData(massError, intensity).getEntry(0,0); int projValIndex = translateLdaValToIndex(projVal); diff --git a/src/edu/umich/andykong/ptmshepherd/utils/Peptide.java b/src/edu/umich/andykong/ptmshepherd/utils/Peptide.java index d40354b..6015787 100644 --- a/src/edu/umich/andykong/ptmshepherd/utils/Peptide.java +++ b/src/edu/umich/andykong/ptmshepherd/utils/Peptide.java @@ -131,7 +131,71 @@ else if (curIonType == 'x' || curIonType == 'y' || curIonType == 'z') float cmass = (cTermMass + ccharge * AAMasses.protMass) / ccharge; for (int i = 0; i < cLen - 1; i++) { cmass += (aaMasses[this.pepSeq.charAt(cLen - 1 - i) - 'A'] + mods[cLen - 1 - i]) / ccharge; - if (((rightI - 1) >= i) && ( i > (leftI - 1))) { + if ((leftI < (cLen - 1 - i)) && ((cLen - 1 - i) <= rightI)) { + knownFrags.add(cmass); + } + } + } + } + + return knownFrags; + } + + /** + * Generate peptide ions that are complementary to a modification site (e.g., unshifted where they should be shifted). + * This is designed to be used on a mono-mutated decoy that does not have the mutated mass differene/pseudo delta-mass + * in its modification list. + * @param ionTypes + * @param dmass + * @param site + * @param maxCharge + * @return + */ + public ArrayList calculateComplementaryFragments(String ionTypes, float dmass, int site, int maxCharge) { + ArrayList knownFrags = new ArrayList<>(this.pepSeq.length() * ionTypes.length()); + + // Set up ion types to be identified + ArrayList nIonTypes = new ArrayList<>(); + ArrayList cIonTypes = new ArrayList<>(); + for (int i = 0; i < ionTypes.length(); i++) { + char curIonType = ionTypes.charAt(i); + if (curIonType == 'a' || curIonType == 'b' || curIonType == 'c') + nIonTypes.add(curIonType); + else if (curIonType == 'x' || curIonType == 'y' || curIonType == 'z') + cIonTypes.add(curIonType); + } + + float [] aaMasses = AAMasses.monoisotopic_masses; + float [] fragTypeShifts = AAMasses.ionTypeShifts; + int cLen = this.pepSeq.length(); + + float nTermMass; + for (Character iType : nIonTypes) { + nTermMass = fragTypeShifts[iType - 'a']; + for (int ccharge = 1; ccharge <= maxCharge; ccharge++) { //loop through charge states + float cmass = AAMasses.protMass + nTermMass; //todo this doesn't account for charge? + for (int i = 0; i < cLen - 1; i++) { //loop through positions on the peptide + cmass += (aaMasses[this.pepSeq.charAt(i) - 'A'] + mods[i]) / ccharge; + if (i != 0) { // todo skip a1/b1 ion, check if this needs to be skipped for c too + if (i < site) { + knownFrags.add(cmass + (dmass / ccharge)); + } else { + knownFrags.add(cmass); + } + } + } + } + } + float cTermMass; + for (Character iType : cIonTypes) { + cTermMass = fragTypeShifts[iType - 'x' + 3]; + for (int ccharge = 1; ccharge <= maxCharge; ccharge++) { + float cmass = (cTermMass + ccharge * AAMasses.protMass) / ccharge; + for (int i = 0; i < cLen - 1; i++) { + cmass += (aaMasses[this.pepSeq.charAt(cLen - 1 - i) - 'A'] + mods[cLen - 1 - i]) / ccharge; + if ((cLen - 1 - i) > site) { + knownFrags.add(cmass + (dmass / ccharge)); + } else { knownFrags.add(cmass); } } diff --git a/test/utils/PeptideTest.java b/test/utils/PeptideTest.java index afa342d..67d9260 100644 --- a/test/utils/PeptideTest.java +++ b/test/utils/PeptideTest.java @@ -76,5 +76,62 @@ void calculatePeptideFragmentsBetween() { for (int i = 0; i < expectedFrags.size(); i++) { assertEquals(expectedFrags.get(i), sitePepFrags.get(i), 0.0001, "Mismatch at index " + i); } + + expectedFrags = new ArrayList<>(Arrays.asList(227.102633f, 324.155397f)); + sitePepFrags = pep.calculatePeptideFragmentsBetween("b", 1, 3, 1); + + assertEquals(expectedFrags.size(), sitePepFrags.size()); + for (int i = 0; i < expectedFrags.size(); i++) { + assertEquals(expectedFrags.get(i), sitePepFrags.get(i), 0.0001, "Mismatch at index " + i); + } + + expectedFrags = new ArrayList<>(Arrays.asList(477.219120f, 574.271884f)); + sitePepFrags = pep.calculatePeptideFragmentsBetween("y", 1, 3, 1); + + assertEquals(expectedFrags.size(), sitePepFrags.size()); + for (int i = 0; i < expectedFrags.size(); i++) { + assertEquals(expectedFrags.get(i), sitePepFrags.get(i), 0.0001, "Mismatch at index " + i); + } + } + + @Test + void calculateComplementaryPeptideFragments() { + String seq = "PEPT"; + float[] mods = new float[]{0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f}; + Peptide pep = new Peptide(seq, mods); + + // Simulate 1.0f dmass on position P3 + ArrayList expectedFrags = new ArrayList<>(Arrays.asList(228.102633f, 324.155397f)); + ArrayList sitePepFrags = pep.calculateComplementaryFragments("b", 1.0f, 2, 1); + + assertEquals(expectedFrags.size(), sitePepFrags.size()); + for (int i = 0; i < expectedFrags.size(); i++) { + assertEquals(expectedFrags.get(i), sitePepFrags.get(i), 0.0001, "Mismatch at index " + i); + } + + // Simulate 1.0f dmass on position P3 + expectedFrags = new ArrayList<>(Arrays.asList(121.065520f, 217.118284f, 346.160877f)); + sitePepFrags = pep.calculateComplementaryFragments("y", 1.0f, 2, 1); + + assertEquals(expectedFrags.size(), sitePepFrags.size()); + for (int i = 0; i < expectedFrags.size(); i++) { + assertEquals(expectedFrags.get(i), sitePepFrags.get(i), 0.0001, "Mismatch at index " + i); + } + + } + + @Test + void AddModTest() { + String seq = "PEP"; + float[] mods = new float[]{1.0f, 0.0f, 0.0f}; + float[] expectedMods = new float[]{0.0f, 0.0f, 0.0f}; + Peptide pep = new Peptide(seq, mods); + + pep.addMod(-1*1.0f, 0); + assertArrayEquals(pep.mods, expectedMods); + + pep.addMod(1.0f, 0); + expectedMods = new float[]{1.0f, 0.0f, 0.0f}; + assertArrayEquals(pep.mods, expectedMods); } }