Skip to content

Commit

Permalink
Peptide.calculatePeptideFragmentsBetween method
Browse files Browse the repository at this point in the history
  • Loading branch information
@danielgeiszler
danielgeiszler committed May 11, 2024
1 parent 2c17a37 commit 8125eae
Show file tree
Hide file tree
Showing 4 changed files with 133 additions and 19 deletions.
61 changes: 42 additions & 19 deletions src/edu/umich/andykong/ptmshepherd/iterativelocalization/IterativeLocalizer.java
View file
Original file line number Diff line number Diff line change
Expand Up @@ -162,6 +162,9 @@ private void fitMatchedIonDistribution() throws Exception {
float massShift = 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[] reducedMzs = reducedIons[0];
Expand Down Expand Up @@ -360,14 +363,14 @@ private void calculateLocalizationProbabilities() throws Exception {
}
}

//if (specName.equals("02330a_GH3_3991_13_PTM_TrainKit_Pmod_Hydroxyproline_200fmol_3xHCD_R1.7708.7708")) {
//if (specName.equals("02330a_GF2_3991_03_PTM_TrainKit_Kmod_Succinyl_200fmol_3xHCD_R1.20972.20972")) {
// this.debugFlag = true;
//}

// Calculate site-specific localization probabilities
float[] mods = psm.getModsAsArray();
boolean[] allowedPoses = parseAllowedPositions(pep, this.allowedAAs, mods);
double[] siteProbs = localizePsm(psm, spec, pep, mods, dMassApex, cBin, allowedPoses, false); // TODO check whether raw, theoretical, or peakapex is better
double[] siteProbs = localizePsm(psm, spec, pep, mods, dMassApex, cBin, allowedPoses, false);

// Update prior probabilities
if (!finalPass)
Expand All @@ -383,6 +386,7 @@ private void calculateLocalizationProbabilities() throws Exception {
strEntropies.add(entropyToString(locEntropy));

// TODO check decoy usefulness
/**
Peptide decoyPep = Peptide.generateDecoy(pep, mods, maxProbI, this.rng, "mono-swapped");
boolean[] decoyAllowedPoses = parseAllowedPositions(decoyPep.pepSeq,
this.allowedAAs, decoyPep.mods);
Expand All @@ -398,6 +402,7 @@ private void calculateLocalizationProbabilities() throws Exception {
maxProbToString(decoyMaxProb, decoyMaxProbAA));
}
}
**/

}
}
Expand All @@ -409,8 +414,8 @@ private void calculateLocalizationProbabilities() throws Exception {
"PTM-Shepherd Localization", specNames, strOutputProbs);
psmf.addColumn(psmf.getColumn("PTM-Shepherd Localization") + 1,
"PTM-Shepherd Best Localization", specNames, strMaxProbs);
psmf.addColumn(psmf.getColumn("PTM-Shepherd Best Localization") + 1,
"PTM-Shepherd Best Decoy Localization", specNames, strMaxProbsDecoy);
//psmf.addColumn(psmf.getColumn("PTM-Shepherd Best Localization") + 1,
// "PTM-Shepherd Best Decoy Localization", specNames, strMaxProbsDecoy);
//psmf.addColumn(psmf.getColumn("delta_mass_maxloc") + 1,
// "delta_mass_entropy", specNames, strEntropies);
//psmf.addColumn(psmf.getColumn("PTM-Shepherd Best Localization") + 1,
Expand Down Expand Up @@ -807,7 +812,7 @@ public int compareTo(SpecProbQ o) {

// Assign each q-Val
ArrayList<String> probModelQVals = new ArrayList<>(specNames.size());
ArrayList<String> decoyModelQVals = new ArrayList<>(specNames.size());
//ArrayList<String> decoyModelQVals = new ArrayList<>(specNames.size());
//ArrayList<String> probDecoyModelVals = new ArrayList<>(specNames.size());
//ArrayList<String> entropyDecoyModelVals = new ArrayList<>(specNames.size());
for (int j = 0; j < specNames.size(); j++) {
Expand All @@ -816,15 +821,15 @@ public int compareTo(SpecProbQ o) {
if (unmodFlag) {
// prob model
probModelQVals.add("");
decoyModelQVals.add("");
//decoyModelQVals.add("");
// prob model with decoys
//probDecoyModelVals.add("");
// entropy model
//entropyDecoyModelVals.add("");
} else {
// prob model
probModelQVals.add(new DecimalFormat("0.0000").format(qValsMap.get(specNames.get(j))));
decoyModelQVals.add(new DecimalFormat("0.0000").format(qValsDecoyMap.get(specNames.get(j))));
//decoyModelQVals.add(new DecimalFormat("0.0000").format(qValsDecoyMap.get(specNames.get(j))));
// prob model with decoys
//probDecoyModelVals.add(new DecimalFormat("0.0000").format(qProbDecoyModel[(int) (maxProb * 1000)]));
// entropy model
Expand All @@ -836,8 +841,8 @@ public int compareTo(SpecProbQ o) {
// Send to PSM file
psmf.addColumn(psmf.getColumn("PTM-Shepherd Best Localization") + 1, "PTM-Shepherd q-val",
specNames, probModelQVals);
psmf.addColumn(psmf.getColumn("PTM-Shepherd q-val") + 1, "PTM-Shepherd decoy q-val",
specNames, decoyModelQVals);
//psmf.addColumn(psmf.getColumn("PTM-Shepherd q-val") + 1, "PTM-Shepherd decoy q-val",
// specNames, decoyModelQVals);
/** //TODO figure out what's going on with these before implementing them, assuming they're even worth doing
psmf.addColumn(psmf.getColumn("delta_mass_BH_loc_q") + 1, "delta_mass_prob_decoyAA_q",
specNames, probDecoyModelVals);
Expand Down Expand Up @@ -950,11 +955,13 @@ private double[] computePoissonBinomialLikelihood(String pep, float[] mods, floa
ArrayList<Float> shiftedPepFrags = new ArrayList<Float>(pepFrags.size());
for (Float frag : pepFrags)
shiftedPepFrags.add(frag + dMass);
pepFrags.addAll(shiftedPepFrags);
pepFrags.addAll(shiftedPepFrags); // todo this wont work with charge state 2

if (debugFlag)
if (debugFlag) {
System.out.println("All pep frags");
System.out.println(pepFrags.stream().map(Object::toString)
.collect(Collectors.joining(", ")));
}

// Filter peakMzs and peakInts to only those that match at least one ion
float[] peakMzs = spec.getPeakMZ();
Expand Down Expand Up @@ -997,14 +1004,16 @@ private double[] computePoissonBinomialLikelihood(String pep, float[] mods, floa
float[] matchedIonMassErrors = matchedIons[1];
// Map matched ion intensities to MatchedIonDistribution, negative intensities will be returned as -1
double[] matchedIonProbabilities = this.matchedIonDist.calcIonProbabilities(matchedIonIntensities, matchedIonMassErrors);
if (debugFlag) {
if (debugFlag && (i == 8 || i == 10)) {
System.out.println("\n");
System.out.println("Start");
System.out.println(pep);
System.out.println("Position " + (i + 1));;
System.out.println("Matched ions");
System.out.println(Arrays.toString(matchedIonIntensities));
System.out.println(Arrays.toString(matchedIonProbabilities));
System.out.println("Matched ions done");
System.out.println("End\n");
}
// Format these for Poisson Binomial input
ionProbs[cAllowedPos] = matchedIonProbabilities;
Expand Down Expand Up @@ -1045,9 +1054,12 @@ private double[] computeLikelihoods(String pep, float[] mods, float dMass, boole
shiftedPepFrags.add(frag + dMass);
pepFrags.addAll(shiftedPepFrags);

if (debugFlag)
if (debugFlag) {
System.out.println("All pep frags");
System.out.println(pepFrags.stream().map(Object::toString)
.collect(Collectors.joining(", ")));
System.out.println();
}

// Filter peakMzs and peakInts to only those that match at least one ion
float[] peakMzs = spec.getPeakMZ();
Expand Down Expand Up @@ -1241,9 +1253,19 @@ public float[][] findMatchedIons(ArrayList<Float> pepFrags, float[] peakMzs, flo
float[] matchedIonErrors = new float[peakMzs.length];
Collections.sort(pepFrags);

if (debugFlag)
if (debugFlag) {
System.out.println("PepFrags:");
System.out.println(pepFrags.stream().map(Object::toString)
.collect(Collectors.joining(", ")));
System.out.println();

System.out.println("Spectrum or reduced spectrum");
for (int i = 0; i < peakMzs.length; i++) {
System.out.println(peakMzs[i] + "\t\t" + peakInts[i]);
}
System.out.println();
}


int specIdx = 0;
int maxSpecIdx = peakMzs.length;
Expand Down Expand Up @@ -1277,15 +1299,11 @@ public float[][] findMatchedIons(ArrayList<Float> pepFrags, float[] peakMzs, flo
}
wIndex--;
if (matchedIons == 1) { // If only matched one ion in minVal->maxVal window
if (debugFlag)
System.out.println(peakMzs[wIndex]);
matchedIonIntensities[startIdx] = maxInt;
matchedIonErrors[startIdx] = Math.abs(frag - peakMzs[startIdx]) / peakMzs[startIdx] * 1000000.0f;
} else { // If more than one ion in minVal->maxVal window, max is matched and rest are unmatched
for (int i = startIdx; i <= wIndex; i++) {
if (Math.abs(peakInts[i] - maxInt) < 0.001) {
if (debugFlag)
System.out.println(peakMzs[wIndex] + "\t" + peakInts[wIndex]);
matchedIonIntensities[i] = maxInt;
matchedIonErrors[i] = Math.abs(frag - peakMzs[i]) / peakMzs[i] * 1000000.0f;
break;
Expand All @@ -1299,8 +1317,13 @@ public float[][] findMatchedIons(ArrayList<Float> pepFrags, float[] peakMzs, flo
specIdx++;
}

if (debugFlag)
if (debugFlag) {
System.out.println();
System.out.println("Matched ions");
System.out.println(Arrays.toString(matchedIonIntensities));
System.out.println(Arrays.toString(matchedIonErrors));
System.out.println();
}

float[][] matchedIons = new float[][]{
matchedIonIntensities,
Expand Down
10 changes: 10 additions & 0 deletions src/edu/umich/andykong/ptmshepherd/iterativelocalization/LocalizationLikelihood.java
View file
Original file line number Diff line number Diff line change
Expand Up @@ -31,4 +31,14 @@ public double[] getSiteLikelihoods() {
return this.siteLikelihoods;
}
}

public String toString() {
StringBuffer sb = new StringBuffer("(");
for (Mod m : this.mods) {
for (int i = 0; i < m.siteLikelihoods.length; i++) {
sb.append(Double.toString(m.siteLikelihoods[i]) + ")");
}
}
return sb.toString();
}
}
66 changes: 66 additions & 0 deletions src/edu/umich/andykong/ptmshepherd/utils/Peptide.java
View file
Original file line number Diff line number Diff line change
Expand Up @@ -75,6 +75,72 @@ else if (curIonType == 'x' || curIonType == 'y' || curIonType == 'z')
return knownFrags;
}

/**
* Generate peptide ions between two sites, exclusive of the final site (e.g., left for y-ions, right for b-ions).
* @param ionTypes
* @param siteA
* @param siteB
* @param maxCharge
* @return
*/
public ArrayList<Float> calculatePeptideFragmentsBetween(String ionTypes, int siteA, int siteB, int maxCharge) {
ArrayList<Float> knownFrags = new ArrayList<>(this.pepSeq.length() * ionTypes.length());

// If the sites are the same, there are no fragments between
if (siteA == siteB)
return knownFrags;

// Make sure sites are in correct order
int leftI = Math.min(siteA, siteB);
int rightI = Math.max(siteA, siteB);

// Set up ion types to be identified
ArrayList<Character> nIonTypes = new ArrayList<>();
ArrayList<Character> 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;
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 ((leftI <= i) && (i < rightI)) {
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 (((rightI - 1) >= i) && ( i > (leftI - 1))) {
knownFrags.add(cmass);
}
}
}
}

return knownFrags;
}

public Peptide generateDecoy(Random rng, String method) {
if (method.equals("shuffled"))
return generateShuffledDecoy(this.pepSeq, this.mods, rng);
Expand Down
15 changes: 15 additions & 0 deletions test/utils/PeptideTest.java
View file
Original file line number Diff line number Diff line change
Expand Up @@ -62,4 +62,19 @@ void calculatePeptideFragments() {
System.out.println(sitePepFrags.toString());

}

@Test
void calculatePeptideFragmentsBetween() {
String seq = "PEPTIDE";
float[] mods = new float[]{0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f};
Peptide pep = new Peptide(seq, mods);

ArrayList<Float> expectedFrags = new ArrayList<>(Arrays.asList(324.1554f, 425.20306f, 376.17145f, 477.21912f));
ArrayList<Float> sitePepFrags = pep.calculatePeptideFragmentsBetween("by", 2, 4, 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);
}
}
}

0 comments on commit 8125eae

Please sign in to comment.