Correct MS1 isotope error based on the Mercury method.

fcyu · Feb 4, 2018 · 6539f32 · 6539f32
1 parent fd9494a
commit 6539f32
Show file tree

Hide file tree

Showing 6 changed files with 1,108 additions and 52 deletions.
diff --git a/src/main/java/proteomics/ECL2.java b/src/main/java/proteomics/ECL2.java
@@ -133,7 +133,27 @@ private ECL2(String parameter_path, String spectra_path, String dbName) throws I
             System.exit(1);
         }
 
-        new PreSpectra(spectra_parser, build_index_obj, parameter_map, ext, sqlPath);
+        double ms1Tolerance = Double.valueOf(parameter_map.get("ms1_tolerance"));
+        double leftInverseMs1Tolerance = 1 / (1 + ms1Tolerance * 1e-6);
+        double rightInverseMs1Tolerance = 1 / (1 - ms1Tolerance * 1e-6);
+        int ms1ToleranceUnit = Integer.valueOf(parameter_map.get("ms1_tolerance_unit"));
+
+        PreSpectra preSpectra = new PreSpectra(spectra_parser, ms1Tolerance, leftInverseMs1Tolerance, rightInverseMs1Tolerance, ms1ToleranceUnit, build_index_obj, parameter_map, ext, sqlPath);
+
+        if (dev && parameter_map.get("C13_correction").contentEquals("1")) {
+            BufferedWriter writer = new BufferedWriter(new FileWriter("spectrum.dev.csv"));
+            writer.write("scanNum,charge,finalIsotopeCorrectionNum,isotopeCorrectionNum,pearsonCorrelationCoefficient,expMz1,expMz2,expMz3,expInt1,expInt2,expInt3,theoMz1,theoMz2,theoMz3,theoInt1,theoInt2,theoInt3\n");
+            Map<Integer, TreeMap<Integer, TreeSet<PreSpectra.DevEntry>>> scanDevEntryMap = preSpectra.getScanDevEntryMap();
+            for (int scanNum : scanDevEntryMap.keySet()) {
+                TreeMap<Integer, TreeSet<PreSpectra.DevEntry>> chargeDevEntryMap = scanDevEntryMap.get(scanNum);
+                for (int charge : chargeDevEntryMap.keySet()) {
+                    for (PreSpectra.DevEntry devEntry : chargeDevEntryMap.get(charge)) {
+                        writer.write(String.format(Locale.US, "%d,%d,%d,%d,%f,%f,%f,%f,%f,%f,%f,%f,%f,%f,%f,%f,%f\n", scanNum, charge, devEntry.isotopeCorrectionNum, devEntry.isotopeCorrectionNum, devEntry.pearsonCorrelationCoefficient, devEntry.expMatrix[0][0], devEntry.expMatrix[1][0], devEntry.expMatrix[2][0], devEntry.expMatrix[0][1], devEntry.expMatrix[1][1], devEntry.expMatrix[2][1], devEntry.theoMatrix[0][0], devEntry.theoMatrix[1][0], devEntry.theoMatrix[2][0], devEntry.theoMatrix[0][1], devEntry.theoMatrix[1][1], devEntry.theoMatrix[2][1]));
+                    }
+                }
+            }
+            writer.close();
+        }
 
         logger.info("Start searching...");
 

diff --git a/src/main/java/proteomics/Search/Search.java b/src/main/java/proteomics/Search/Search.java
@@ -20,7 +20,6 @@ public class Search {
     private final MassTool mass_tool_obj;
     private final TreeMap<Integer, Set<String>> bin_seq_map;
     private final BuildIndex build_index_obj;
-    private final int[] C13_correction_range;
     final double single_chain_t;
     private final boolean cal_evalue;
 
@@ -40,17 +39,10 @@ public Search(BuildIndex build_index_obj, Map<String, String> parameter_map, dou
         } else {
             cal_evalue = true;
         }
-
-        String[] temp = parameter_map.get("C13_correction_range").split(",");
-        C13_correction_range = new int[temp.length];
-        for (int i = 0; i < temp.length; ++i) {
-            C13_correction_range[i] = Integer.valueOf(temp[i].trim());
-        }
-        Arrays.sort(C13_correction_range);
     }
 
     ResultEntry doSearch(SparseVector xcorrPL, TreeMap<Integer, List<Double>> binScoresMap, double massWithoutLinker, int precursorCharge, double precursorMass, String scanId) throws IOException {
-        int max_chain_bin_idx = Math.min(build_index_obj.massToBin(massWithoutLinker + C13_correction_range[C13_correction_range.length - 1] * MassTool.C13_DIFF) + 1 - bin_seq_map.firstKey(), bin_seq_map.lastKey());
+        int max_chain_bin_idx = Math.min(build_index_obj.massToBin(massWithoutLinker) + 1 - bin_seq_map.firstKey(), bin_seq_map.lastKey());
         int min_chain_bin_idx = bin_seq_map.firstKey();
 
         if (max_chain_bin_idx < min_chain_bin_idx) {
@@ -80,20 +72,11 @@ ResultEntry doSearch(SparseVector xcorrPL, TreeMap<Integer, List<Double>> binSco
                 break;
             }
 
-            double left_mass_2;
-            double right_mass_2;
-            int left_idx_2;
-            int right_idx_2;
-            NavigableMap<Integer, Set<String>> sub_map = new TreeMap<>();
-
-            // consider C13 correction
-            for (int i : C13_correction_range) {
-                left_mass_2 = massWithoutLinker + i * MassTool.C13_DIFF - build_index_obj.binToRightMass(idx_1) - leftMs1Tol;
-                right_mass_2 = massWithoutLinker + i * MassTool.C13_DIFF - build_index_obj.binToLeftMass(idx_1) + rightMs1Tol;
-                left_idx_2 = build_index_obj.massToBin(left_mass_2);
-                right_idx_2 = build_index_obj.massToBin(right_mass_2) + 1;
-                sub_map.putAll(bin_seq_map.subMap(left_idx_2, true, right_idx_2, true));
-            }
+            double left_mass_2 = massWithoutLinker - build_index_obj.binToRightMass(idx_1) - leftMs1Tol;
+            double right_mass_2 = massWithoutLinker - build_index_obj.binToLeftMass(idx_1) + rightMs1Tol;
+            int left_idx_2 = build_index_obj.massToBin(left_mass_2);
+            int right_idx_2 = build_index_obj.massToBin(right_mass_2) + 1;
+            NavigableMap<Integer, Set<String>> sub_map  = bin_seq_map.subMap(left_idx_2, true, right_idx_2, true);
 
             if (!sub_map.isEmpty()) {
                 if (!checkedBinSet.contains(idx_1)) {