Changeset 4289

Timestamp:

Feb 16, 2012, 2:49:23 PM (12 years ago)

Author:

marianne

Message:

Refs #774. Distribution tolerance update. ClusterID set only if features are aligned. Alignment if no common sequences.

Location:

trunk

Files:

• api/core/conf/common-queries.xml (modified) (1 diff)
• api/core/src/org/proteios/core/Feature.java (modified) (2 diffs)
• plugin/src/org/proteios/plugins/FeatureSequencePropagator.java (modified) (38 diffs)

trunk/api/core/conf/common-queries.xml

@@ -264 +264 @@
     </description>
   </query>
-  
-  <query id="GET_UNIQUE_SEQUENCES_IN_MSFILE_FOR_PROJECT" type="HQL">
-    <sql>
-      SELECT DISTINCT f.peptideSequence
-      FROM FeatureData f
-      WHERE f.project = :project
-      AND f.msFile = :msFile
-      AND NOT f.peptideSequence IS NULL
-      AND f.chargeState IS 2
-      ORDER by f.apexRetentionTimeInMinutes
-    </sql>
-    <description>
-      Load all distinct (unique) sequences for a ms file in the Feature table for a project 
-    </description>
-  </query>
-  
 
   <query id="GET_IDS_FOR_SEPARATION_METHODS_WITH_SEPARATION_EVENTS" type="HQL">

trunk/api/core/src/org/proteios/core/Feature.java

@@ -21 +21 @@
   public static final Item TYPE = Item.PROTEIOS_FEATURE;
 
-  
-  /*@SuppressWarnings("unchecked")
-  public static List<String> getAllClusterIDs(DbControl dc)
-  {
-    org.hibernate.Query query = HibernateUtil.getPredefinedQuery(dc
-      .getHibernateSession(),
-      "GET_ALL_FEATURES_WITH_SPECIFIC_CLUSTERID");
-    List<String> retval = query.list();
-    return retval;
-  }*/
-  
   
   @SuppressWarnings("unchecked")
@@ -50 +39 @@
   
   
-
-
-  @SuppressWarnings("unchecked")
-  public static List<String> getUniqueFeatureSequences(Project project, DbControl dc, File msFile)
-  {
-    org.hibernate.Query query = HibernateUtil.getPredefinedQuery(dc
-      .getHibernateSession(),
-      "GET_UNIQUE_SEQUENCES_IN_MSFILE_FOR_PROJECT");
-    query.setEntity("project", project.getData());
-    query.setEntity("msFile", msFile.getData());
-    List<String> retval = query.list();
-    return retval;
-  }
-
-  @SuppressWarnings("unchecked")
-  public static List<String> getUniqueFeatureSequnces(Project project, DbControl dc, File featureFile)
-  {
-    org.hibernate.Query query = HibernateUtil.getPredefinedQuery(dc
-      .getHibernateSession(),
-      "GET_UNIQUE_SEQUNCES_IN_FEATURE_FILE_FOR_PROJECT");
-    query.setEntity("project", project.getData());
-    query.setEntity("featureFile", featureFile.getData());
-    List<String> retval = query.list();
-    return retval;
-  }
-
-
   /*
    * From the Identifiable interface

trunk/plugin/src/org/proteios/plugins/FeatureSequencePropagator.java

@@ -123 +123 @@
       Timestamp tStmp = new Timestamp(new Date().getTime());
       long initialClusterID = tStmp.getTime();
-
-      double tolerance = 0.01;
-      double timeTol = 0.5;
+      
+      //DEFAULT VALUES! Used if no sequences in the project at all!
+      double defaultTolerance = 0.01;
+      double defaultTimeTol = 0.5;
+      
+      double tolerance = 0;
+      double timeTol = 0;
       int count = 0;
       // dbControl = object used for communicating with database
@@ -133 +137 @@
       DbControl dcFile = sc.newDbControl();
       // Retrieve the job parameters
-      tolerance = ((Float) job.getValue("tolerance")).doubleValue();
-      timeTol = ((Float) job.getValue("timeTol")).doubleValue();
+      //tolerance = ((Float) job.getValue("tolerance")).doubleValue();
+      //timeTol = ((Float) job.getValue("timeTol")).doubleValue();
       boolean matchBetweenFractions = (Boolean) job
           .getValue("matchBetweenFractions");
@@ -147 +151 @@
         //Annotator annotator = new Annotator(factory);
         //AnnotationType fractionAnnotationType = annotator
-          //  .getAnnotationType("fraction", Type.STRING, File.TYPE);
+        //  .getAnnotationType("fraction", Type.STRING, File.TYPE);
         //creates an instance of the template class
         File outCoreFile = factory.create(File.class);
@@ -172 +176 @@
         writer.println("Propagation report for project:" + project
             .getName());
-        writer.println("m/z tolerance:" + tolerance);
-        writer.println("RT (minutes) tolerance:" + timeTol);
+        writer.println("Default m/z tolerance if no common pair-wise sequences in project: " + defaultTolerance);
+        writer.println("Default RT (minutes) tolerance if no common pair-wise sequences in project: " + defaultTimeTol);
         writer
         .println("Matching between fraction:" + matchBetweenFractions);
+
         //get unique ms files for the project
         List<File> uniqueMsFiles = Feature.getUniqueMsFiles(project, dc);
-        writer.println("Number of files with MS features:" + uniqueMsFiles.size() + "\n");
-        log.debug("Number of files with MS features:" + uniqueMsFiles.size() + "\n");
+        int nbrOfFiles = uniqueMsFiles.size();
+        writer.println("Number of files with MS features:" + nbrOfFiles + "\n");
+        log.debug("Number of files with MS features:" + nbrOfFiles + "\n");
 
         //query for the class
         ItemQuery<Feature> featureQuery = qf.select(Feature.class);
+
         // Only features in project should be retrieved
         // OBS! has hql object as input
@@ -188 +195 @@
             .property("project"), Hql.entity(project)));
 
-        int nbrOfFiles = uniqueMsFiles.size();
-        //int totNbrFiles = 2*nbrOfFiles-2;
-        //int totMatchEvents = nbrOfFiles - 1;
-        //int currMatchEvents = 0;
+
 
         FilePoint[] simList = new FilePoint[nbrOfFiles*(nbrOfFiles-1)/2];
@@ -197 +201 @@
 
         ArrayList<TreeSet<Feature>> featureList = new ArrayList<TreeSet<Feature>>(nbrOfFiles);
-        ArrayList<ArrayList<Feature>> totFeatureList = new ArrayList<ArrayList<Feature>>(nbrOfFiles);
 
         //  double[][] mzDist = new double[nbrOfFiles][];
@@ -206 +209 @@
         Arrays.fill(hasSequence,true);
 
-        //TODO:What to do if hasSequence is false? Filterfunction needed for features when file has no sequences. 
-
         //Random randGen = new Random();
-        
+
         double[] originalFeatWithSeq = new double[nbrOfFiles];
 
@@ -216 +217 @@
           File currentMsFile = uniqueMsFiles.get(fileNbr);
 
-          List<String> currSeq = Feature.getUniqueFeatureSequences(project, dc, currentMsFile);
-          writer.println("Number of unique sequences in file:" + fileNbr + " " + currSeq.size());
-          
-          if (currSeq==null || currSeq.isEmpty())
-          {
-            hasSequence[fileNbr]=false;
+
+          //Clearing clusterID for current project
+          featureQuery.restrict(Restrictions.eq(Hql
+              .property("msFile"), Hql.entity(currentMsFile)));
+          featureQuery.restrict(Restrictions.neq(Hql
+              .property("clusterId"), null));
+
+
+          List<Feature> allFeaturesClusterID = featureQuery.list(dc);
+          featureQuery.reset();
+
+          //Clearing the clusterID:s in current project
+          for (Feature clusterF: allFeaturesClusterID){
+            if (progress != null)
+            {
+              progress.display((fileNbr+1)/nbrOfFiles,
+                  "Clearing clusterID in project. Processing file nr " + (fileNbr+1));
+            }
+
+
+            clusterF.setClusterId(null);  
+            dc.reattachItem(clusterF);
+
           }
-          //seqList.add(currSeq);
+          dc.commit();
+          dc = sc.newDbControl();
 
           // Retrieve features for current ms file
@@ -234 +253 @@
 
           List<Feature> allFeatures = featureQuery.list(dc);
-
-          //TODO: Setting clusterID, should be done together with distribution calculation
-          for (Feature clusterF: allFeatures){
-            //if (clusterF.getClusterId()==null){
-
-            clusterF.setClusterId(initialClusterID++);  
-            dc.reattachItem(clusterF);
-            //}
-          }
-          dc.commit();
-
-          dc = sc.newDbControl();
           allFeatures = featureQuery.list(dc);
-
-          totFeatureList.add(new ArrayList<Feature>());
-          totFeatureList.get(fileNbr).addAll(allFeatures);
-
-          writer.println("Number of all features in file:" + fileNbr + " " + allFeatures.size());
 
 
@@ -277 +279 @@
           featureList.add(new TreeSet<Feature>(new SeqComparator()));
           featureList.get(fileNbr).addAll(features);
-          
-
-          /*for (Feature f:featureList.get(fileNbr))
-          {
-            log.debug(f.getPeptideSequence()+ " "+f.getApexRetentionTimeInMinutes()+ " " +f.getChargeState());
-          }*/
-
-          
+
+
           originalFeatWithSeq[fileNbr]=features.size();
-          
+
           writer.println("File number " +fileNbr + " is " + currentMsFile);
           writer.println("Number of all features in file " + fileNbr + ": " + allFeatures.size());
@@ -299 +295 @@
           log.debug("Number of unique features with sequences in file " + fileNbr + ": " + featureList.get(fileNbr).size());
 
-
-
           //String frac1 = getFirstAnnotationValue(annotator,
           //fractionAnnotationType, currentMsFile);
         }
+
+
 
         writer.println("===============================");
@@ -353 +349 @@
             int indexKey = (secFileNbr+1)*(secFileNbr-2)/2+fileNbr+1;
 
-
-
-            //  writer.println("The index key :" + indexKey);
-            //  
-            //if (simValues[indexKey]<0)
-            //{ 
-
             ArrayList<Point> points = getPointMatchBySeq(features, features2);
+            ArrayList<Point> evalPoints = new ArrayList<Point>();
+            double simValue = 0;
+            PolynomialSplineFunction func = null;
 
             //double percentageForRegression = 10;
@@ -370 +362 @@
             if(totNbrMatches==0){
 
-              writer.println("There are no common sequences for file: " + (fileNbr+1) + " and " +(secFileNbr+1));
-              writer.println("These files will not be aligned.");
-              continue;
+              writer.println("There are no common sequences for file: " + (fileNbr+1) + " and " +(secFileNbr+1) +": " +totNbrMatches);
+              writer.println("Alignment will be performed without retention time correction.");
+
+            }else{
+
+
+              Collections.sort(points,new DiffComparator());
+
+              double lowerQuartile = points.get((int) Math.round(0.25*points.size())).getDiff();
+              double upperQuartile = points.get((int) Math.round(0.75*points.size())).getDiff();
+              double medianQuartile = points.get((int) Math.round(0.5*points.size())).getDiff();
+              double interQuartileRange = upperQuartile-lowerQuartile;
+              double upperFence = upperQuartile + 3*interQuartileRange;
+
+              writer.println("Median RT difference before alignment: " +medianQuartile);
+              writer.println("Largest RT difference before alignment: " +points.get(points.size()-1).getDiff());
+
+              int maxCutOff = points.size()-1;
+              while(points.get(maxCutOff).getDiff()>upperFence)
+              { 
+                points.remove(maxCutOff);
+                maxCutOff--;
+              }
+
+
+              writer.println("Removing " +(totNbrMatches-points.size()) + " outliers. Total number of sequence matches is now: " +points.size());
+
+              medianQuartile = points.get((int) Math.round(0.5*points.size())).getDiff();
+
+              writer.println("Median RT difference is now: " +medianQuartile);
+              writer.println("Largest RT difference is now: " +points.get(points.size()-1).getDiff());
+
+              Collections.sort(points,new MzDiffComparator());
+              double medianQuartileMz = points.get((int) Math.round(0.5*points.size())).getMzDiff();
+              writer.println("Median mz difference before alignment: " +medianQuartileMz);
+              writer.println("Largest mz difference before alignment: " +points.get(points.size()-1).getMzDiff());
+
+              tolerance=points.get(points.size()-1).getMzDiff();
+              timeTol=medianQuartile;
+
+
+
+
+
+              //no duplicates in x-values
+              Collections.sort(points, new RTComparator()); 
+
+              log.debug("Points size: " +points.size());
+
+
+              //TODO: bortkommentera när mzDist ska användas som likhetsvärde
+              simValue = points.size()*2.0/(features2.size()+features.size());
+              writer.println("Similarity for file number" +fileNbr + " and " + secFileNbr +": " +simValue);
+
+              /*
+               * Weight order
+              Collections.sort(points, new WeightComparator()); 
+
+              int halfLength = (int) Math.ceil(points.size()/2);
+
+              evalPoints.addAll(points.subList(0,halfLength-1));
+
+              points.subList(0,halfLength).clear();
+
+              Collections.sort(points, new RTComparator()); 
+              Collections.sort(evalPoints, new RTComparator());*/
+
+
+              //Partitioning matching sequences in alignment set and evaluation set
+              int ndx = 0;
+              while(ndx<points.size())
+              { 
+                //TODO: avkommentera när random ska användas
+                //features.remove(randGen.nextInt(features.size()));
+
+                evalPoints.add(points.remove(ndx));
+                ndx++;
+
+              }
+
+              log.debug("Points size: " +points.size());
+              log.debug("Eval size: " +evalPoints.size());
+
+              ArrayList<Point> tempPoints = new ArrayList<Point>();
+              tempPoints = evalPoints;
+              evalPoints = points;
+              points = tempPoints;
+
+              func = estimateSplineFunc(points, writer);
+              if (func==null){
+                writer.println("Adding evaluation points to spline function estimation. No precision and recall will be computed for this pair-wise alignment.");
+                points.addAll(evalPoints);
+                Collections.sort(points, new RTComparator()); 
+                evalPoints.clear();
+                func = estimateSplineFunc(points, writer);
+                if (func==null){
+                  writer.println("Too few points for alignment.");
+                  continue;
+                }
+
+              }
             }
-
-
-            Collections.sort(points,new DiffComparator());
-
-            double lowerQuartile = points.get((int) Math.round(0.25*points.size())).getDiff();
-            double upperQuartile = points.get((int) Math.round(0.75*points.size())).getDiff();
-            double medianQuartile = points.get((int) Math.round(0.5*points.size())).getDiff();
-            double interQuartileRange = upperQuartile-lowerQuartile;
-            double upperFence = upperQuartile + 3*interQuartileRange;
-
-            writer.println("Median RT difference before alignment: " +medianQuartile);
-            writer.println("Largest RT difference before alignment: " +points.get(points.size()-1).getDiff());
-
-            int maxCutOff = points.size()-1;
-            while(points.get(maxCutOff).getDiff()>upperFence && points.get(maxCutOff).getDiff()>timeTol)
-            { 
-              points.remove(maxCutOff);
-              maxCutOff--;
-            }
-
-
-            writer.println("Removing " +(totNbrMatches-points.size()) + " outliers. Total number of sequence matches is now: " +points.size());
-
-            medianQuartile = points.get((int) Math.round(0.5*points.size())).getDiff();
-
-            writer.println("Median RT difference is now: " +medianQuartile);
-            writer.println("Largest RT difference is now: " +points.get(points.size()-1).getDiff());
-
-            Collections.sort(points,new MzDiffComparator());
-            double medianQuartileMz = points.get((int) Math.round(0.5*points.size())).getMzDiff();
-            writer.println("Median mz difference before alignment: " +medianQuartileMz);
-            writer.println("Largest mz difference before alignment: " +points.get(points.size()-1).getMzDiff());
-
-            tolerance=points.get(points.size()-1).getMzDiff();
-            timeTol=medianQuartile;
-
-
-            ArrayList<Point> evalPoints = new ArrayList<Point>();
-
-
-            //no duplicates in x-values
-            Collections.sort(points, new RTComparator()); 
-
-            log.debug("Points size: " +points.size());
-
-
-            //TODO: bortkommentera när mzDist ska användas som likhetsvärde
-            double simValue = points.size()*2.0/(features2.size()+features.size());
-            writer.println("Similarity for file number" +fileNbr + " and " + secFileNbr +": " +simValue);
-
-            /*
-            Collections.sort(points, new WeightComparator()); 
-
-            int halfLength = (int) Math.ceil(points.size()/2);
-
-            evalPoints.addAll(points.subList(0,halfLength-1));
-
-            points.subList(0,halfLength).clear();
-
-            Collections.sort(points, new RTComparator()); 
-            Collections.sort(evalPoints, new RTComparator());*/
-
-            int ndx = 0;
-
-            while(ndx<points.size())
-            { 
-              //TODO: avkommentera när random ska användas
-              //features.remove(randGen.nextInt(features.size()));
-
-              evalPoints.add(points.remove(ndx));
-              ndx++;
-
-            }
-
-            log.debug("Points size: " +points.size());
-            log.debug("Eval size: " +evalPoints.size());
-
-            ArrayList<Point> tempPoints = new ArrayList<Point>();
-            tempPoints = evalPoints;
-            evalPoints = points;
-            points = tempPoints;
-
-
-            //}
-
-
-            //  writer.println("Similarity:" + simValues[indexKey]);
-
-
-            PolynomialSplineFunction func = estimateSplineFunc(points, writer);
-            simList[indexKey] = new FilePoint(fileNbr,secFileNbr,simValue,func,evalPoints);
+            simList[indexKey] = new FilePoint(fileNbr,secFileNbr,simValue, tolerance, timeTol, func,evalPoints);
 
             writer.println();
@@ -471 +471 @@
         }
         //writer.println("Best match is :" + (firstMatch+1) +" and " +(secMatch+1)+ "\n");
+
+
+        Arrays.sort(simList,Collections.reverseOrder(new SimComparator()));
+
+        // Propagate and get count
         
+        double noSeqMzTol = 0;
+        double noSeqRtTol = 0;
         
-        Arrays.sort(simList,Collections.reverseOrder(new SimComparator()));
-
-        // Propagate and get count
-      
-        countNbrMatches=0;
-        for(FilePoint f:simList)
+        for(int i=0; i<simList.length; i++)
         { 
-          if (f.getSimValue()>0){
-          countNbrMatches++;
-
+
+          FilePoint f = simList[i];
+          
           if (progress != null)
           {
             progress.display((100 *countNbrMatches/simList.length),
-                "Aligning files " +(f.getFirstFile()+1) + " and " + (f.getSecondFile()+1) +". File-pair " +countNbrMatches +" of " +simList.length);
+                "Aligning files " +(f.getFirstFile()+1) + " and " + (f.getSecondFile()+1) +". File-pair " +(i+1) +" of " +simList.length +". Similarity score: " +f.getSimValue());
           }
 
-
+          //there are no common features at all, setting the mzTol and rtTol to the one with the least similar files
+          if (f.getMzTol()==0 && f.getRtTol()==0){
+            
+          
+            if (noSeqMzTol == 0 && noSeqRtTol==0){
+
+              if (i-1>=0){
+                
+                noSeqMzTol = simList[i-1].getMzTol();
+                noSeqRtTol = simList[i-1].getRtTol();
+                
+                writer.println("===============================");
+                writer.println("There are no common sequences for these files. Alignment will be performed with mz and rt tolerance for least similar files.");
+                
+              }else{
+                
+                writer.println("===============================");
+                writer.println("There are no sequences in this project. Alignment will be performed with default values for mz and rt tolerance.");
+                
+                noSeqMzTol = defaultTolerance;
+                noSeqRtTol = defaultTimeTol;
+                
+              }
+
+            }else{
+              f.setMzTol(noSeqMzTol);
+              f.setRtTol(noSeqRtTol);
+              
+            }
+
+          }
+          
 
           //TODO: här ska allFeatures skickas in för ytterligare utvärdering
-          count += matchFeatureLists(f,factory, project,writer,featureQuery,uniqueMsFiles,tolerance,timeTol);
-
-          //TODO: lägg in hur många unika sekvenser filen har EFTER matchning
-          }
-        }
-        
-        
+          count += matchFeatureLists(f,factory, project,writer,featureQuery,uniqueMsFiles);
+
+
+
+        }
+
+
         //writerDat.close();
         // Commit to database
         dc.commit();
-        
-        
+
+
         writer.println("===============================");
-        
+
         featureQuery.reset();
         dc = sc.newDbControl();
@@ -517 +550 @@
           featureQuery.restrict(Restrictions.gteq(Hql
               .property("chargeState"),Expressions.integer(2)));
-          
+
           writer.println("Number of features with sequences in file " +currentMsFile +" after alignment: " +featureQuery.list(dc).size());
           writer.println("Total increase in sequence coverage for file " +currentMsFile + ": " +(featureQuery.list(dc).size()-originalFeatWithSeq[fileNbr])/originalFeatWithSeq[fileNbr]);
@@ -523 +556 @@
         }
         dc.close();
-        
+
         writer.println("Total number of features matched for " + uniqueMsFiles.size() + " MS files: " + count);
         response.setDone("Total number of features matched for " + uniqueMsFiles.size() + " MS files: " + count);
-        
+
         // Close log file
         writer.close();
@@ -624 +657 @@
                 log.debug("Point replaced.");
               }
-              //TODO: lägg in hits sorteringen..????
+
             }
 
@@ -746 +779 @@
    */
   private int matchFeatureLists(FilePoint fp, ItemFactory factory,
-      Project project, PrintWriter writer, ItemQuery<Feature> featureQuery,List<File> uniqueMsFiles, double tolerance, double timeTol)
+      Project project, PrintWriter writer, ItemQuery<Feature> featureQuery,List<File> uniqueMsFiles)
   { 
 
+    int hitsCounter = 0;
+    //timestamp for clusterID
+    Timestamp tStmp = new Timestamp(new Date().getTime());
+    long initialClusterID = tStmp.getTime();
 
     DbControl dc = sc.newDbControl();
 
-    ArrayList<Point> features = fp.getEvalList();
+    ArrayList<Point> evalList = fp.getEvalList();
     PolynomialSplineFunction func = fp.getSplineFunction();
+    double tolerance = fp.getMzTol();
+    double timeTol = fp.getRtTol();
 
     featureQuery.restrict(Restrictions.eq(Hql
@@ -782 +821 @@
     ArrayList<Feature> allFeatures2 = new ArrayList<Feature>(tempAllFeatures2);
     featureQuery.reset();
-    
+
     featureQuery.restrict(Restrictions.eq(Hql
         .property("msFile"), Hql.entity(uniqueMsFiles.get(fp.getSecondFile()))));
@@ -790 +829 @@
     featureQuery.restrict(Restrictions.neq(Hql.property("peptideSequence"), null));
     List<Feature> tempAllFeaturesSeq2 = featureQuery.list(dc);
-    ArrayList<Feature> allFeaturesSeq2 = new ArrayList<Feature>(tempAllFeaturesSeq);
+    ArrayList<Feature> allFeaturesSeq2 = new ArrayList<Feature>(tempAllFeaturesSeq2);
     featureQuery.reset();
+
 
     ArrayList<Feature> bestMatchList = new ArrayList<Feature>();
@@ -803 +843 @@
 
     writer.println("===============================");
-    writer.println("Evaluating match with " +features.size() +" features with sequences for file " +fp.getFirstFile() +" and " +fp.getSecondFile() +" with similarity score: " +fp.getSimValue());
+    writer.println("Evaluating match with " +evalList.size() +" features with sequences for file " +fp.getFirstFile() +" and " +fp.getSecondFile() +" with similarity score: " +fp.getSimValue());
     writer.println();
 
 
-    log.debug("Evaluating match with " +features.size() +" features with sequences for file " +fp.getFirstFile() +" and " +fp.getSecondFile() +" with similarity score: " +fp.getSimValue());
-    
-    
+    log.debug("Evaluating match with " +evalList.size() +" features with sequences for file " +fp.getFirstFile() +" and " +fp.getSecondFile() +" with similarity score: " +fp.getSimValue());
+
+
     String localSampleId = null;
     String fractionId = null;
@@ -854 +894 @@
     double fineTolSeq = 0;
 
-    for (Point p : features)
-    {     
-      Feature f = p.getXFeature();
-      String fseq = f.getPeptideSequence();
+    if (!evalList.isEmpty()){
+
+      for (Point p : evalList)
+      {     
+        Feature f = p.getXFeature();
+        String fseq = f.getPeptideSequence();
+
+        float fRT = f.getApexRetentionTimeInMinutes();
+
+        if (func != null)
+        {
+          try
+          {
+            fRT = (float) func.value(f
+                .getApexRetentionTimeInMinutes().doubleValue());
+          }
+          catch (ArgumentOutsideDomainException e)
+          {}
+        }
+
+        for(Point p2: evalList)
+        {
+          Feature f2 = p2.getYFeature();
+          String f2seq = f2.getPeptideSequence();
+
+          if ((Math.abs(f.getMassToChargeRatio() - f2.getMassToChargeRatio()) < tolerance) && f2.getChargeState().equals(f.getChargeState()) && (Math.abs(f2.getApexRetentionTimeInMinutes() - fRT) < timeTol))
+          {
+            if (f2seq.equals(fseq))
+            {
+              fineSeq++;
+
+            }else{
+
+              fineTolSeq++;
+
+              writer.println("Current timetol: " +timeTol);
+              writer.println("Feature with sequence " + f.getPeptideSequence() + " matching " + f2.getPeptideSequence() + " time diff: " +(Math.abs(f2.getApexRetentionTimeInMinutes() - fRT) + " mz diff: " +(Math.abs(f.getMassToChargeRatio() - f2.getMassToChargeRatio()))));
+
+            }
+
+            //sequences that match but are not within tolerance limits => conflict  
+          }else if (f2seq.equals(fseq) && f2.getChargeState().equals(f.getChargeState())){
+
+            log.debug("Current timetol: " +timeTol +" and current mz tol: " +tolerance);
+            log.debug("Feature with sequence " + f.getPeptideSequence() + " not matching " + f2.getPeptideSequence() + " time diff: " +(Math.abs(f2.getApexRetentionTimeInMinutes() - fRT) + " mz diff: " +(Math.abs(f.getMassToChargeRatio() - f2.getMassToChargeRatio()))));
+
+          }
+        }
+      }
+    }
+
+
+
+    int count = 0;
+    int minNdx = 0;
+
+    //double timeTolTemp = 0;
+    Feature bestMatchPrev = null;
+    for(Feature f:allFeatures){
+
 
       float fRT = f.getApexRetentionTimeInMinutes();
@@ -872 +968 @@
       }
 
-      for(Point p2: features)
-      {
-        Feature f2 = p2.getYFeature();
-        String f2seq = f2.getPeptideSequence();
-        //  log.debug("Feature to match RT:" + f.getApexRetentionTimeInMinutes() + "after shift:" + fRT + " diff:" + (fRT - f.getApexRetentionTimeInMinutes()));
-        //sequence evaluation for recall and tol evaluation for precision
-        if ((Math.abs(f.getMassToChargeRatio() - f2.getMassToChargeRatio()) < tolerance) && f2.getChargeState().equals(f.getChargeState()) && (Math.abs(f2.getApexRetentionTimeInMinutes() - fRT) < timeTol))
-        {
-          
-          //fseq och f2seq ska inte kunna vara null
-          if (f2seq.equals(fseq))
-          {
-            fineSeq++;
-
-          }else{
-
-            fineTolSeq++;
-
-            writer.println("Current timetol: " +timeTol);
-            writer.println("Feature with sequence " + f.getPeptideSequence() + " matching " + f2.getPeptideSequence() + " time diff: " +(Math.abs(f2.getApexRetentionTimeInMinutes() - fRT) + " mz diff: " +(Math.abs(f.getMassToChargeRatio() - f2.getMassToChargeRatio()))));
-
-          }
-        
-        //sequences that should be matching 
-        }else if (f2seq.equals(fseq) && f2.getChargeState().equals(f.getChargeState())){
-          
-          log.debug("Current timetol: " +timeTol +" and current mz tol: " +tolerance);
-          log.debug("Feature with sequence " + f.getPeptideSequence() + " not matching " + f2.getPeptideSequence() + " time diff: " +(Math.abs(f2.getApexRetentionTimeInMinutes() - fRT) + " mz diff: " +(Math.abs(f.getMassToChargeRatio() - f2.getMassToChargeRatio()))));
-      
-        }
-      }
-    }
-
-
-
-    int count = 0;
-    int minNdx = 0;
-
-    //double timeTolTemp = 0;
-
-    for(Feature f:allFeatures){
-
-
-      float fRT = f.getApexRetentionTimeInMinutes();
-
-      if (func != null)
-      {
-        try
-        {
-          fRT = (float) func.value(f
-              .getApexRetentionTimeInMinutes().doubleValue());
-        }
-        catch (ArgumentOutsideDomainException e)
-        {}
-      }
-
       while( minNdx<(allFeatures2.size()-1) && fRT-allFeatures2.get(minNdx).getApexRetentionTimeInMinutes()>timeTol){
         minNdx++;
@@ -935 +975 @@
 
       boolean hit = false;
-    
+
 
       //  log.debug("For RT: " +fRT +" Starting at: " +allFeatures2.get(minNdx).getApexRetentionTimeInMinutes());
@@ -942 +982 @@
       { 
 
-        //slut på sökning
+        //end of search
         if(f2.getApexRetentionTimeInMinutes()-fRT>timeTol){
           break;
@@ -960 +1000 @@
               if(Math.abs(bestMatch.getApexRetentionTimeInMinutes() - fRT)>Math.abs(f2.getApexRetentionTimeInMinutes() - fRT))
               { 
-                log.debug("Best match is updated.");
-                log.debug("Old rt diff is: " +Math.abs(bestMatch.getApexRetentionTimeInMinutes() - fRT));
                 bestMatch = f2;
-                log.debug("New rt diff is: " +Math.abs(bestMatch.getApexRetentionTimeInMinutes() - fRT));
-                
+                //log.debug("Best match is updated.");
+                //log.debug("Old rt diff is: " +Math.abs(bestMatch.getApexRetentionTimeInMinutes() - fRT));
+
+                //log.debug("New rt diff is: " +Math.abs(bestMatch.getApexRetentionTimeInMinutes() - fRT));
+
               }
 
@@ -975 +1016 @@
 
 
-      if(bestMatch!=null)
-      { 
-
-
-        //log.debug("Matched features with RT " + f.getApexRetentionTimeInMinutes() + " and RT:" + bestMatch.getApexRetentionTimeInMinutes());
-        
-
-        //log.debug("First feature: " +f.getClusterId() + " bestMatch feature: " +bestMatch.getClusterId());
-        //All features have clusterID:s, they are set at start of alignment -> all features that have a clusterID have undergone alignment
-        if(!f.getClusterId().equals(bestMatch.getClusterId())){
-          //log.debug("Setting new clusterID");
-
-          bestMatchList.add(bestMatch);
-          fList.add(f);
-          bestMatchClusterIDList.add(bestMatch.getClusterId());
-
-        }
-
-      }
-
-
-    }
-    
-    int hitsCounter = 0;
-    
+      if(bestMatch!=null){  
+
+        if(f.getClusterId()!=null && bestMatch.getClusterId()!=null || f.getClusterId()==null && bestMatch.getClusterId()!=null){
+          if (!bestMatch.getClusterId().equals(f.getClusterId())){
+            //Feature pair have undergone matching before and are not in same cluster => add to bestMatchClusterIDList
+            bestMatchList.add(bestMatch);
+            fList.add(f);
+            bestMatchClusterIDList.add(bestMatch.getClusterId());
+          }
+        }else{
+
+          //single match
+
+          count++;
+
+          if(f.getClusterId()==null && bestMatch.getClusterId()==null){
+            //Feature-pair has not undergone alignment before, no need to add to bestMatchClusterIDList
+            long clusterMatchID=initialClusterID++;
+            bestMatch.setClusterId(clusterMatchID); 
+            f.setClusterId(clusterMatchID);
+
+          }else if (f.getClusterId()!=null && bestMatch.getClusterId()==null){
+            //only f has undergone alignment before but not bestMatch, no need to add to bestMatchClusterIDList
+            bestMatch.setClusterId(f.getClusterId()); 
+          }
+
+          dc.reattachItem(f);
+          dc.reattachItem(bestMatch);
+
+          //Updating hits list for single match
+          if (f.getPeptideSequence()!= null && bestMatch.getPeptideSequence()==null){
+            hitsCounter+=updateHitsList(f,bestMatch,localSampleId2,fractionId2,allFeaturesSeq,writer,dc,factory,project);
+          }else if(f.getPeptideSequence()== null && bestMatch.getPeptideSequence()!=null){
+            hitsCounter+=updateHitsList(bestMatch,f,localSampleId,fractionId,allFeaturesSeq2,writer,dc,factory,project);
+          }
+
+
+        }
+
+        if (bestMatchPrev!=null && bestMatchPrev.getMassToChargeRatio()==bestMatch.getMassToChargeRatio() && Math.abs(bestMatchPrev.getApexRetentionTimeInMinutes()-bestMatch.getApexRetentionTimeInMinutes())<0.00001){
+          log.debug("Previous bestMatch equals new bestMatch.");
+
+        }
+        bestMatchPrev=bestMatch;
+
+      }
+
+
+
+
+
+    }
+
+
     if (!bestMatchClusterIDList.isEmpty()){
 
@@ -1012 +1082 @@
       log.debug("Size of fList: " +fList.size());
 
-      count=bestMatchList.size();
+      count+=bestMatchList.size();
 
 
@@ -1025 +1095 @@
         Feature f = fList.get(featureNbr);
 
+        //log.debug("First feature: " +f.getClusterId() + " bestMatch feature: " +bestMatch.getClusterId());
+        if (f.getClusterId()!=null){
+          bestMatch.setClusterId(f.getClusterId());
+        }else{
+          f.setClusterId(bestMatch.getClusterId());
+        }
+
         dc.reattachItem(f);
         dc.reattachItem(bestMatch);
-
-        //log.debug("First feature: " +f.getClusterId() + " bestMatch feature: " +bestMatch.getClusterId());
-
-        bestMatch.setClusterId(f.getClusterId());
 
         //updating hits
@@ -1039 +1112 @@
           hitsCounter+=updateHitsList(bestMatch,f,localSampleId,fractionId,allFeaturesSeq2,writer,dc,factory,project);
         }
-        //TODO: what to do if both have sequences?
-
-        //log.debug("First feature after set: " +f.getClusterId() + " bestMatch feature: " +bestMatch.getClusterId());
+        //TODO: what to do if both have sequences??? CONFILCT resolution needed!
 
       }
@@ -1054 +1125 @@
     writer.println("Number of features matched:" + count);
     writer.println("Number of hits created:" + hitsCounter);
-    writer.println("Pairwsie alignment recall based on sequences with mzTol " +tolerance +" and  timeTol " +timeTol +": " +fineSeq/features.size());
-    writer.println("Pairwsie alignment precision based on sequences with mzTol " +tolerance +" and  timeTol " +timeTol +": " +(1-fineTolSeq/features.size()));
+    writer.println("Pairwsie alignment recall based on sequences with mzTol " +tolerance +" and  timeTol " +timeTol +": " +fineSeq/evalList.size());
+    writer.println("Pairwsie alignment precision based on sequences with mzTol " +tolerance +" and  timeTol " +timeTol +": " +(1-fineTolSeq/evalList.size()));
     writer.println();
 
@@ -1069 +1140 @@
     dc.reattachItem(f);
     dc.reattachItem(f2);
-    
+
     int newHitCounter = 0;
-    
+
     Hit h = null;
     for (Hit h1 : f.getHits())
@@ -1095 +1166 @@
             .lastIndexOf(' '));
       for (Feature f3 : allFeaturesSeq)
-      {
+      { 
+
         if ((fseq.equals(f3.getPeptideSequence()) || (fseq
             .contains(" ") && f3
@@ -1157 +1229 @@
       f2.setHits(f.getHits());
     }
-    
+
     return newHitCounter;
 
@@ -1196 +1268 @@
       catch (MathException e1)
       {
-        throw new BaseException(e1.getMessage());
+        writer.println("Too few points to perform spline interpolation: " +points.size());
       }
     }
@@ -1478 +1550 @@
   private class FilePoint{
 
-    private double simValue = 0.0;
+    private double simValue = 0;
+    private double mzTol = 0;
+    private double rtTol = 0;
     private PolynomialSplineFunction func = null;
     private int file1 = 0;
@@ -1484 +1558 @@
     private ArrayList<Point> evalPoint = new ArrayList<Point>();
 
+
     public FilePoint()
     {}
 
-    public FilePoint(int fileNbr, int fileNbr2, double simValue, PolynomialSplineFunction func, ArrayList<Point> evalPoint)
+    public FilePoint(int fileNbr, int fileNbr2, double simValue, double mzTol, double rtTol, PolynomialSplineFunction func, ArrayList<Point> evalPoint)
     { 
       this.file1 = fileNbr;
       this.file2 = fileNbr2;
       this.simValue = simValue;
+      this.mzTol = mzTol;
+      this.rtTol = rtTol;
       this.func = func;
       this.evalPoint.addAll(evalPoint);
@@ -1512 +1589 @@
     }
 
+    public double getMzTol(){
+
+      return this.mzTol;
+    }
+    
+    public void setMzTol(double mzTol){
+
+      this.mzTol = mzTol;
+    }
+
+    public double getRtTol(){
+
+      return this.rtTol;
+    }
+    
+    public void setRtTol(double rtTol){
+
+      this.rtTol = rtTol;
+    }
+
     public PolynomialSplineFunction getSplineFunction() {