Changeset 4288

Timestamp:

Feb 13, 2012, 9:16:52 PM (12 years ago)

Author:

marianne

Message:

Refs #774. Bug fix: Added two-way alignment. Skipping alignment if there are no matching sequences.
Added: Recall and precision evaluation for features with sequences. Feature number coverage increase computation.

File:

• trunk/plugin/src/org/proteios/plugins/FeatureSequencePropagator.java (modified) (56 diffs)

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

@@ -86 +86 @@
    */
   private static final org.apache.log4j.Logger log = org.apache.log4j.LogManager
-  .getLogger("org.proteios.plugins.featuresequencepropagator");
+      .getLogger("org.proteios.plugins.featuresequencepropagator");
 
 
@@ -98 +98 @@
   {
     return new AboutImpl("Propagate identifications",
-      "Match features and propagate identifications", "0.2",
-      "2011, Marianne Sandin, Fredrik Levander", null, null,
-    "http://www.proteios.org");
+        "Match features and propagate identifications", "0.2",
+        "2011, Marianne Sandin, Fredrik Levander", null, null,
+        "http://www.proteios.org");
   }
 
@@ -107 +107 @@
   public void init(SessionControl sc, ParameterValues configuration,
       ParameterValues job)
-  throws BaseException
-  {
+          throws BaseException
+          {
     super.init(sc, configuration, job);
-  }
+          }
 
 
@@ -123 +123 @@
       Timestamp tStmp = new Timestamp(new Date().getTime());
       long initialClusterID = tStmp.getTime();
-      
+
       double tolerance = 0.01;
       double timeTol = 0.5;
@@ -129 +129 @@
       // dbControl = object used for communicating with database
       DbControl dc = sc.newDbControl();
-      
+
       //used for alignment file
       DbControl dcFile = sc.newDbControl();
@@ -136 +136 @@
       timeTol = ((Float) job.getValue("timeTol")).doubleValue();
       boolean matchBetweenFractions = (Boolean) job
-      .getValue("matchBetweenFractions");
+          .getValue("matchBetweenFractions");
       Project project = (Project) job.getValue("project");
       // qf is used to query the database
@@ -145 +145 @@
         ItemFactory factory = new ItemFactory(dcFile);
 
-        Annotator annotator = new Annotator(factory);
-        AnnotationType fractionAnnotationType = annotator
-        .getAnnotationType("fraction", Type.STRING, File.TYPE);
+        //Annotator annotator = new Annotator(factory);
+        //AnnotationType fractionAnnotationType = annotator
+          //  .getAnnotationType("fraction", Type.STRING, File.TYPE);
         //creates an instance of the template class
         File outCoreFile = factory.create(File.class);
@@ -157 +157 @@
         outCoreFile.setDirectory(project.getProjectDirectory());
         outCoreFile.setName("AlignLog-" + job.getId() + ".txt");
-        
+
         //datFile.setDirectory(project.getProjectDirectory());
         //datFile.setName("AlignDat" + job.getId() + ".txt");
-        
+
         dcFile.saveItem(outCoreFile);
         //dc.saveItem(datFile);
@@ -166 +166 @@
         OutputStream outstream = outCoreFile.getUploadStream(false);
         //OutputStream outstream2 = datFile.getUploadStream(false);
-        
+
         PrintWriter writer = new PrintWriter(outstream);
         //PrintWriter writerDat = new PrintWriter(outstream2);
-        
+
         writer.println("Propagation report for project:" + project
-          .getName());
+            .getName());
         writer.println("m/z tolerance:" + tolerance);
         writer.println("RT (minutes) tolerance:" + timeTol);
@@ -186 +186 @@
         // OBS! has hql object as input
         featureQuery.restrictPermanent(Restrictions.eq(Hql
-          .property("project"), Hql.entity(project)));
+            .property("project"), Hql.entity(project)));
 
         int nbrOfFiles = uniqueMsFiles.size();
@@ -192 +192 @@
         //int totMatchEvents = nbrOfFiles - 1;
         //int currMatchEvents = 0;
-        
+
         FilePoint[] simList = new FilePoint[nbrOfFiles*(nbrOfFiles-1)/2];
         Arrays.fill(simList, new FilePoint());
-        
+
         ArrayList<TreeSet<Feature>> featureList = new ArrayList<TreeSet<Feature>>(nbrOfFiles);
         ArrayList<ArrayList<Feature>> totFeatureList = new ArrayList<ArrayList<Feature>>(nbrOfFiles);
-        
-      //  double[][] mzDist = new double[nbrOfFiles][];
-        PearsonsCorrelation pearCorr = new PearsonsCorrelation();
-        
+
+        //  double[][] mzDist = new double[nbrOfFiles][];
+        //PearsonsCorrelation pearCorr = new PearsonsCorrelation();
+
         double[] meanFeatureLength = new double[nbrOfFiles];
         boolean[] hasSequence = new boolean[nbrOfFiles];
         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();
         
-        //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];
 
         for (int fileNbr = 0; fileNbr<nbrOfFiles; fileNbr++)
         {
           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())
           {
@@ -224 +228 @@
           // Only retrieve features with peptide sequences assigned
           featureQuery.restrict(Restrictions.eq(Hql
-            .property("msFile"), Hql.entity(currentMsFile)));
+              .property("msFile"), Hql.entity(currentMsFile)));
           featureQuery.restrict(Restrictions.gteq(Hql
-            .property("chargeState"),Expressions.integer(2)));
+              .property("chargeState"),Expressions.integer(2)));
           featureQuery.order(Orders.asc(Hql.property("apexRetentionTimeInMinutes")));
-          
-          dc.reconnect();
+
           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);
+
+            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());
-          
-          
+
+
           //setMzDistribution(allFeatures,mzDist,meanFeatureLength,fileNbr);
           featureQuery.reset();
-          
+
           log.debug("Mean length: " +meanFeatureLength[fileNbr]);
           //log.debug("Length of mzDist: " +mzDist[fileNbr].length);
-        
-          
+
+
           //ordered according to RT for the treeset so that the first one (least RT) will be picked
           featureQuery.restrict(Restrictions.eq(Hql
-            .property("msFile"), Hql.entity(currentMsFile)));
+              .property("msFile"), Hql.entity(currentMsFile)));
           featureQuery.restrict(Restrictions.neq(Hql
-            .property("peptideSequence"), null));
+              .property("peptideSequence"), null));
           featureQuery.restrict(Restrictions.gteq(Hql
-            .property("chargeState"),Expressions.integer(2)));
+              .property("chargeState"),Expressions.integer(2)));
           featureQuery.order(Orders.asc(Hql.property("apexRetentionTimeInMinutes")));
 
@@ -271 +274 @@
           List<Feature> features =  featureQuery.list(dc);
           featureQuery.reset();
-          
-          featureList.add(new TreeSet(new SeqComparator()));
+
+          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());
           writer.println("Number of features with sequences in file " + fileNbr + ": " + features.size());
-          writer.println("Number of unique features with sequences in file " + fileNbr + ": " + featureList.get(fileNbr).size());
-        
+          writer.println("Number of unique features with sequences in file " + fileNbr + ": " +featureList.get(fileNbr).size());
+
           writer.println();
-          
+
           log.debug("File number " +fileNbr + " is " + currentMsFile);
           log.debug("Number of all features in file " + fileNbr + ": " + allFeatures.size());
           log.debug("Number of features with sequences in file " + fileNbr + ": " + features.size());
           log.debug("Number of unique features with sequences in file " + fileNbr + ": " + featureList.get(fileNbr).size());
-          
-          
-          
+
+
+
           //String frac1 = getFirstAnnotationValue(annotator,
           //fractionAnnotationType, currentMsFile);
@@ -300 +306 @@
 
         writer.println("===============================");
-        
-        
-        
+
+
+
         int countNbrMatches=0;
-        
+
         for(int fileNbr=0; fileNbr<nbrOfFiles-1; fileNbr++ )
         { 
-          
-          
+
+
           for (int secFileNbr=fileNbr+1; secFileNbr<nbrOfFiles; secFileNbr++) 
           { 
-            
+
             countNbrMatches++;
-            
+
             if (progress != null)
             {
               progress.display((100 *countNbrMatches)/((nbrOfFiles*(nbrOfFiles-1)/2)),
-                "Calculating similarities. Processing file nr " + (fileNbr+1) + " and " +(secFileNbr+1));
+                  "Calculating similarities. Processing file nr " + (fileNbr+1) + " and " +(secFileNbr+1));
             }
-            
+
             writer.println("Calculating similarities. Processing file nr " + (fileNbr+1) + " and " +(secFileNbr+1) +". Progress is: " +(100 *countNbrMatches)/((nbrOfFiles*(nbrOfFiles-1)/2)));
             //double simValue = pearCorr.correlation(mzDist[fileNbr], mzDist[secFileNbr]);
             //writer.println("Similarity for file number" +fileNbr + " and " + secFileNbr +": " +simValue);
-            
-            ArrayList<Feature> features = new ArrayList(featureList.get(fileNbr));
-            ArrayList<Feature> features2 = new ArrayList(featureList.get(secFileNbr));
-            
-            //ArrayList<Feature> features = new ArrayList(totFeatureList.get(fileNbr));
-            //ArrayList<Feature> features2 = new ArrayList(totFeatureList.get(secFileNbr));
-            
-            double meanWeight1=meanFeatureLength[fileNbr];
-            double meanWeight2=meanFeatureLength[secFileNbr];
+
+            ArrayList<Feature> features = new ArrayList<Feature>(featureList.get(fileNbr));
+            ArrayList<Feature> features2 = new ArrayList<Feature>(featureList.get(secFileNbr));
+
+            //double meanWeight1=meanFeatureLength[fileNbr];
+            //double meanWeight2=meanFeatureLength[secFileNbr];
             /*  
 
@@ -349 +352 @@
 
             int indexKey = (secFileNbr+1)*(secFileNbr-2)/2+fileNbr+1;
-            
-            
-            
+
+
+
             //  writer.println("The index key :" + indexKey);
             //  
@@ -358 +361 @@
 
             ArrayList<Point> points = getPointMatchBySeq(features, features2);
-            
-            double percentageForRegression = 10;
+
+            //double percentageForRegression = 10;
             //ArrayList<Point> points = getPointMatchByTol(features,features2,percentageForRegression);
-            
+
             int totNbrMatches = points.size();
             writer.println("Total number of matches: " +totNbrMatches);
+
+            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;
+            }
+
 
             Collections.sort(points,new DiffComparator());
@@ -372 +383 @@
             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)
@@ -382 +393 @@
               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;
@@ -448 +459 @@
             //}
 
-            
+
             //  writer.println("Similarity:" + simValues[indexKey]);
-            
+
 
             PolynomialSplineFunction func = estimateSplineFunc(points, writer);
             simList[indexKey] = new FilePoint(fileNbr,secFileNbr,simValue,func,evalPoints);
-            
+
             writer.println();
           }
@@ -460 +471 @@
         }
         //writer.println("Best match is :" + (firstMatch+1) +" and " +(secMatch+1)+ "\n");
-
-        //TODO: What if func is null?
+        
+        
         Arrays.sort(simList,Collections.reverseOrder(new SimComparator()));
-        
+
         // Propagate and get count
-        
+      
         countNbrMatches=0;
         for(FilePoint f:simList)
         { 
-          
+          if (f.getSimValue()>0){
           countNbrMatches++;
-          
+
           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 " +countNbrMatches +" of " +simList.length);
           }
-          
-          
-          
+
+
+
           //TODO: här ska allFeatures skickas in för ytterligare utvärdering
           count += matchFeatureLists(f,factory, project,writer,featureQuery,uniqueMsFiles,tolerance,timeTol);
-          
-        }
-
-        /*seqList.add(ListUtils.sum(seqList.get(firstMatch),seqList.get(secMatch))); 
-          writer.println("Size of union list is :" + seqList.get(seqList.size()-1).size() + "\n");
-          seqList.set(firstMatch,null); 
-          seqList.set(secMatch,null);
-
-          featureSet.add(bestSet); 
-          writer.println("Size of union set is :" + bestSet.size() + "\n");
-          featureSet.set(firstMatch,null); 
-          featureSet.set(secMatch,null);
-
-          currMatchEvents++;*/
-
-
+
+          //TODO: lägg in hur många unika sekvenser filen har EFTER matchning
+          }
+        }
+        
         
         //writerDat.close();
         // Commit to database
         dc.commit();
-        response.setDone("Features matched for " + uniqueMsFiles.size() + " MS files. " + count + " features sequences propagated.");
-        writer.println("Features matched for " + uniqueMsFiles.size() + " MS files. " + count + " features sequences propagated.");
+        
+        
+        writer.println("===============================");
+        
+        featureQuery.reset();
+        dc = sc.newDbControl();
+        for(int fileNbr=0; fileNbr<nbrOfFiles; fileNbr++ )
+        {
+          File currentMsFile = uniqueMsFiles.get(fileNbr);
+          featureQuery.restrict(Restrictions.eq(Hql
+              .property("msFile"), Hql.entity(currentMsFile)));
+          featureQuery.restrict(Restrictions.neq(Hql
+              .property("peptideSequence"), null));
+          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]);
+          featureQuery.reset();
+        }
+        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();
@@ -536 +559 @@
     String retval = null;
     Annotation anno2 = annotator.getAnnotation(currentMsFile,
-      annotationType);
+        annotationType);
     List<String> fracList2 = extracted(anno2);
     if (fracList2 != null && !fracList2.isEmpty())
@@ -545 +568 @@
   }
 
-  
+
   private ArrayList<Point> getPointMatchBySeq(ArrayList<Feature> features, ArrayList<Feature> features2)
   {
-    
+
     Iterator<Feature> itr = features.iterator();
     Iterator<Feature> itr2;
 
     ArrayList<Point> points = new ArrayList<Point>();
-    
+
     while(itr.hasNext())
     {
       Feature f = itr.next();
-      
+
       itr2 = features2.iterator();
-      
+
       int nbrHits = 0;
 
@@ -565 +588 @@
 
         Feature f2 = itr2.next();
-        
-        
+
+
         if(f.getPeptideSequence().equals(f2.getPeptideSequence()) && f.getChargeState().equals(f2.getChargeState()))
         { 
@@ -572 +595 @@
 
           Point tempPoint = new Point(f,f2);
-          
+
           //contains is based on the RT x-value, no duplicate x is allowed for the regression
           if(!points.contains(tempPoint))
@@ -579 +602 @@
               if(!(Math.abs(f.getMassToChargeRatio()-f2.getMassToChargeRatio())>0.01) )
               { 
-                
-                
+
+
                 //double weight = Math.max(Math.abs(tempPoint.getXLength()-meanWeight1)/meanWeight1,Math.abs(tempPoint.getYLength()-meanWeight2)/meanWeight2);
                 //weight = 1/weight;
                 //tempPoint.setWeight(weight);
-                
+
                 points.add(tempPoint);
               }
@@ -607 +630 @@
           }
 
-          
+
         }
 
@@ -621 +644 @@
 
     }
-    
+
     return points;
-    
-  }
-  
+
+  }
+
   private ArrayList<Point> getPointMatchByTol(ArrayList<Feature> featuresIn, ArrayList<Feature> featuresIn2, double percentage )
   {
-    
+
     Comparator<Feature> apexC = new ApexIntensityComparator();
     int cutOff1 = (int) Math.floor(featuresIn.size()*(100-percentage)/100);
     int cutOff2 = (int) Math.floor(featuresIn2.size()*(100-percentage)/100);
-    
+
     Collections.sort(featuresIn,apexC);
     Collections.sort(featuresIn2,apexC);
-    
+
     //Collections.sort(featuresIn,Collections.reverseOrder(apexC));
     //Collections.sort(featuresIn2,Collections.reverseOrder(apexC));
-    
+
     List<Feature> features = featuresIn.subList(cutOff1,featuresIn.size()-1);
     List<Feature> features2 = featuresIn2.subList(cutOff2,featuresIn2.size()-1);
-    
+
     log.debug("Tol matching size: " +features.size() +" and " +features2.size());
-    
-    
+
+
     Iterator<Feature> itr = features.iterator();
     Iterator<Feature> itr2;
 
     ArrayList<Point> points = new ArrayList<Point>();
-    
+
     while(itr.hasNext())
     {
       Feature f = itr.next();
-      
+
       itr2 = features2.iterator();
-      
+
       int nbrHits = 0;
 
@@ -685 +708 @@
               //multiple hit for the tolerance
 
-              
+
               log.debug("Duplicate tol hit. RT diff p1: " +points.get(points.size()-1).getDiff() + " p2: " +tempPoint.getDiff());
               //use the one with the "BEST" match, ie the one with the smallest difference
@@ -700 +723 @@
           }
 
-          
-        }
-
-      }
-
-    }
-    
+
+        }
+
+      }
+
+    }
+
     Comparator<Feature> rtC = new RTComparatorF();
     Collections.sort(featuresIn,rtC);
     Collections.sort(featuresIn2,rtC);
-    
+
     log.debug("Tol matching points size: " +points.size());
 
-    
+
     return points;
-    
-  }
-  
-  
+
+  }
+
+
   /**
    * Match feature lists and propagate sequences from first to second list.
    */
   private int matchFeatureLists(FilePoint fp, ItemFactory factory,
-      Project project, PrintWriter writer, ItemQuery<Feature> featureQuery,List<File> uniqueMsFiles, double tolerance, double timeTolin)
+      Project project, PrintWriter writer, ItemQuery<Feature> featureQuery,List<File> uniqueMsFiles, double tolerance, double timeTol)
   { 
-    
-    
+
+
     DbControl dc = sc.newDbControl();
-    
+
     ArrayList<Point> features = fp.getEvalList();
     PolynomialSplineFunction func = fp.getSplineFunction();
-    
+
     featureQuery.restrict(Restrictions.eq(Hql
-      .property("msFile"), Hql.entity(uniqueMsFiles.get(fp.getFirstFile()))));
+        .property("msFile"), Hql.entity(uniqueMsFiles.get(fp.getFirstFile()))));
     featureQuery.restrict(Restrictions.gteq(Hql
-      .property("chargeState"),Expressions.integer(2)));
+        .property("chargeState"),Expressions.integer(2)));
     featureQuery.order(Orders.asc(Hql.property("apexRetentionTimeInMinutes")));
     List<Feature> tempAllFeatures = featureQuery.list(dc);
     ArrayList<Feature> allFeatures = new ArrayList<Feature>(tempAllFeatures);
     featureQuery.reset();
-    
+
     featureQuery.restrict(Restrictions.eq(Hql
         .property("msFile"), Hql.entity(uniqueMsFiles.get(fp.getFirstFile()))));
-      featureQuery.restrict(Restrictions.gteq(Hql
+    featureQuery.restrict(Restrictions.gteq(Hql
         .property("chargeState"),Expressions.integer(2)));
-      featureQuery.order(Orders.asc(Hql.property("apexRetentionTimeInMinutes")));
+    featureQuery.order(Orders.asc(Hql.property("apexRetentionTimeInMinutes")));
     featureQuery.restrict(Restrictions.neq(Hql.property("peptideSequence"), null));
     List<Feature> tempAllFeaturesSeq = featureQuery.list(dc);
     ArrayList<Feature> allFeaturesSeq = new ArrayList<Feature>(tempAllFeaturesSeq);
     featureQuery.reset();
-    
+
     featureQuery.restrict(Restrictions.eq(Hql
-      .property("msFile"), Hql.entity(uniqueMsFiles.get(fp.getSecondFile()))));
+        .property("msFile"), Hql.entity(uniqueMsFiles.get(fp.getSecondFile()))));
     featureQuery.restrict(Restrictions.gteq(Hql
-      .property("chargeState"),Expressions.integer(2)));
+        .property("chargeState"),Expressions.integer(2)));
     featureQuery.order(Orders.asc(Hql.property("apexRetentionTimeInMinutes")));
     List<Feature> tempAllFeatures2 = featureQuery.list(dc);
@@ -760 +783 @@
     featureQuery.reset();
     
+    featureQuery.restrict(Restrictions.eq(Hql
+        .property("msFile"), Hql.entity(uniqueMsFiles.get(fp.getSecondFile()))));
+    featureQuery.restrict(Restrictions.gteq(Hql
+        .property("chargeState"),Expressions.integer(2)));
+    featureQuery.order(Orders.asc(Hql.property("apexRetentionTimeInMinutes")));
+    featureQuery.restrict(Restrictions.neq(Hql.property("peptideSequence"), null));
+    List<Feature> tempAllFeaturesSeq2 = featureQuery.list(dc);
+    ArrayList<Feature> allFeaturesSeq2 = new ArrayList<Feature>(tempAllFeaturesSeq);
+    featureQuery.reset();
+
     ArrayList<Feature> bestMatchList = new ArrayList<Feature>();
     ArrayList<Feature> fList = new ArrayList<Feature>();
@@ -767 +800 @@
     //double timeTol = 0.5;
     //double[] timeTol = {0.25,0.5,0.75,1,1.5,2,3,4,5};
-    double[] timeTol = {0, timeTolin};
-    
+
+
     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();
+
+
+    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 " +features.size() +" features with sequences for file " +fp.getFirstFile() +" and " +fp.getSecondFile() +" with similarity score: " +fp.getSimValue());
-
+    String localSampleId = null;
+    String fractionId = null;
+
+    for (Feature f : allFeatures)
+    {
+      if (f.getPeptideSequence() != null)
+      {
+        Iterator<Hit> i = f.getHits().iterator();
+        if (i.hasNext())
+        {
+          Hit h = i.next();
+          localSampleId = h.getLocalSampleId();
+          fractionId = h.getFractionId();
+          log
+          .debug("Found local sample id : " + localSampleId + " fraction: " + fractionId);
+          break;
+        }
+      }
+    }
 
     String localSampleId2 = null;
@@ -791 +844 @@
           fractionId2 = h.getFractionId();
           log
-            .debug("Found local sample id2:" + localSampleId2 + " fraction:" + fractionId2);
+          .debug("Found local sample id2: " + localSampleId2 + " fraction for set 2: " + fractionId2);
           break;
         }
@@ -797 +850 @@
     }
 
-    
-    double[] fineSeq = new double[timeTol.length];
-    double[] notFineSeq = new double[timeTol.length];
-    double[] fineTolSeq = new double[timeTol.length];
-    /*
+
+    double fineSeq = 0;
+    double fineTolSeq = 0;
+
     for (Point p : features)
     {     
@@ -814 +866 @@
         {
           fRT = (float) func.value(f
-            .getApexRetentionTimeInMinutes().doubleValue());
+              .getApexRetentionTimeInMinutes().doubleValue());
         }
         catch (ArgumentOutsideDomainException e)
@@ -824 +876 @@
         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
-        //TODO: skriv detta som funktion så att både allFeatures och Points kan skickas in för utvärdering
-        for (int k=0;k<timeTol.length;k++){
-
-
-
-          if ((Math.abs(f.getMassToChargeRatio() - f2.getMassToChargeRatio()) < tolerance) && f2.getChargeState().equals(f.getChargeState()) && (Math.abs(f2.getApexRetentionTimeInMinutes() - fRT) < timeTol[k]))
+        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))
           {
-
-            if (f2seq!=null && fseq!=null && f2seq.equals(fseq))
-            {
-              fineSeq[k]++;
-
-            }else{
-
-              fineTolSeq[k]++;
-
-              writer.println("Current timetol: " +timeTol[k]);
-              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()))));
-
-            }
-
-          }else if (f2seq!=null && fseq!=null && f2seq.equals(fseq) && f2.getChargeState().equals(f.getChargeState()))
-          {
-            notFineSeq[k]++;
-            log.debug("Current timetol: " +timeTol[k]);
-            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()))));
+            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;
-    double[] fine = new double[timeTol.length];
-    double[] fineTol = new double[timeTol.length];
-    double[] fineTest = new double[timeTol.length];
-    
     int minNdx = 0;
-    
-    
-  
+
+    //double timeTolTemp = 0;
+
     for(Feature f:allFeatures){
-      
-      
+
+
       float fRT = f.getApexRetentionTimeInMinutes();
 
@@ -884 +922 @@
         {
           fRT = (float) func.value(f
-            .getApexRetentionTimeInMinutes().doubleValue());
+              .getApexRetentionTimeInMinutes().doubleValue());
         }
         catch (ArgumentOutsideDomainException e)
@@ -890 +928 @@
       }
 
-      while( minNdx<(allFeatures2.size()-1) && fRT-allFeatures2.get(minNdx).getApexRetentionTimeInMinutes()>timeTol[timeTol.length-1]){
+      while( minNdx<(allFeatures2.size()-1) && fRT-allFeatures2.get(minNdx).getApexRetentionTimeInMinutes()>timeTol){
         minNdx++;
       }
-      
+
       Feature bestMatch = null;
-      
-      
-      boolean[] hit = new boolean[timeTol.length];
-      boolean[] moreHits = new boolean[timeTol.length];
-      
-    //  log.debug("For RT: " +fRT +" Starting at: " +allFeatures2.get(minNdx).getApexRetentionTimeInMinutes());
-      
+
+      boolean hit = false;
+    
+
+      //  log.debug("For RT: " +fRT +" Starting at: " +allFeatures2.get(minNdx).getApexRetentionTimeInMinutes());
+
       for(Feature f2: allFeatures2.subList(minNdx,allFeatures2.size()-1))
       { 
 
-
-        /*if(fRT-f2.getApexRetentionTimeInMinutes()>timeTol[timeTol.length-1]){
-          continue;
-        }*/
-
         //slut på sökning
-        if(f2.getApexRetentionTimeInMinutes()-fRT>timeTol[timeTol.length-1]){
-      //    log.debug("For RT: " +fRT +"Ending at: " +allFeatures2.get(minNdx).getApexRetentionTimeInMinutes());
+        if(f2.getApexRetentionTimeInMinutes()-fRT>timeTol){
           break;
         }
 
-
-
         if ((Math.abs(f.getMassToChargeRatio() - f2.getMassToChargeRatio()) < tolerance) && f2.getChargeState().equals(f.getChargeState()))
 
         { 
-
-          for (int k=0;k<timeTol.length;k++){
-
-            
-
-            if(Math.abs(f2.getApexRetentionTimeInMinutes() - fRT) < timeTol[k])
-            { 
-              
-              fineTest[k]++;
-              
-              if(!hit[k]){
-                //there is a single hit
-                fine[k]++;
-                hit[k] = true;
-                bestMatch = f2;
-                continue;
-              }else if(!moreHits[k]){
-
-                //there is more than one hit for this peptide sequence
-                
-                fineTol[k]++;
-              //  log.debug("There is a second hit: " +fineTol[k]);
-                //moreHits[k] = true;
-              }
-
+          if(Math.abs(f2.getApexRetentionTimeInMinutes() - fRT) < timeTol)
+          { 
+            if(!hit){
+              //there is a single hit             
+              hit = true;
+              bestMatch = f2; 
+            }else{
               //check if the second and further hits are better than the current bestMatch, will always contain the best match of all features, even split peaks
               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) + "new is: " +Math.abs(f2.getApexRetentionTimeInMinutes() - fRT));
+                log.debug("Old rt diff is: " +Math.abs(bestMatch.getApexRetentionTimeInMinutes() - fRT));
                 bestMatch = f2;
+                log.debug("New rt diff is: " +Math.abs(bestMatch.getApexRetentionTimeInMinutes() - fRT));
+                
               }
 
             }
-            
 
           }
-
-        }
-        
-      }
-
-        //TODO: Same with sequence. Write back to feature table as well as hit table. May not be restricted at all to project in future.
-      /*QueryFactory qf2 = new QueryFactory();
-      ItemQuery<Feature> featureQuery2 = qf2.select(Feature.class);
-      featureQuery2.restrict(Restrictions.eq(Hql.property("project"), Hql.entity(project)));*/
-      
-      /*ItemQuery query = ...
-      List<Integer> selected = ...
-
-      query.restrict(
-        Restrictions.in(
-          Hql.property("id"),
-          Expressions.parameter("selectedItems", selected, Type.INT)
-      ));
-*/
-      
-      
+        }
+
+      }
+
+
       if(bestMatch!=null)
       { 
+
+
+        //log.debug("Matched features with RT " + f.getApexRetentionTimeInMinutes() + " and RT:" + bestMatch.getApexRetentionTimeInMinutes());
         
-        
-        log.debug("Matched features with RT " + f.getApexRetentionTimeInMinutes() + " and RT:" + bestMatch.getApexRetentionTimeInMinutes());
-        count++;
-        
-        log.debug("First feature: " +f.getClusterId() + " bestMatch feature: " +bestMatch.getClusterId());
+
+        //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());
-          
-          //All features that have the same clusterID as bestMatch will have their clusterID changed to that of f
-          //TODO: Restrict to project to start with - LATER REMOVE SO ALL FEATURES IN DATABASE WILL BE CHANGED
-          
-          //long tempID = bestMatch.getClusterId();
-          
-          /*dc.reattachItem(bestMatch);
-          bestMatch.setClusterId(f.getClusterId());
-          
-          dc.commit();  
-          dc = sc.newDbControl();
-          dc.reattachItem(f);*/
-          
-          /*
-          log.debug("First feature: " +f.getClusterId() + " bestMatch feature: " +bestMatch.getClusterId());
-          
-          
-          //feature should find itself in the database
-          
-          featureQuery2.restrict(Restrictions.eq(Hql.property("clusterId"), Expressions.aLong(f.getClusterId())));
-          List<Feature> clusterIDFeatures =  featureQuery2.list(dc);
-          
-          //size of clusterF is larger than 0 only if there exists features that have the same clusterID as f
-          log.debug("Size of clusterIdFeatures: " +clusterIDFeatures.size());
-          featureQuery2.reset();
-          
-          for (Feature clusterF: clusterIDFeatures){
-            dc.reattachItem(clusterF);
-            clusterF.setClusterId(f.getClusterId());
-            
-          }*/
-        }
-        
-        
-        //TODO: Propagate the sequence (within the project for now)-> if both have a sequence or both are null, do nothing (CHANGE LATER)
-        //Update hits table
-        
-        /*if (f.getPeptideSequence()!=null && bestMatch.getPeptideSequence()==null){  
-          updateHitsList(f,bestMatch,localSampleId2, fractionId2,allFeatures,writer,dc,factory,project);
-        }else if(f.getPeptideSequence()==null && bestMatch.getPeptideSequence()!=null){ 
-          updateHitsList(bestMatch,f,localSampleId2,fractionId2,allFeatures2,writer,dc,factory,project);
-        }*/
-        
-        
-      }
-      
-
-    }
+
+        }
+
+      }
+
+
+    }
+    
+    int hitsCounter = 0;
     
     if (!bestMatchClusterIDList.isEmpty()){
@@ -1055 +1012 @@
       log.debug("Size of fList: " +fList.size());
 
-
+      count=bestMatchList.size();
 
 
@@ -1064 +1021 @@
 
         int featureNbr = bestMatchClusterIDList.indexOf(bestMatchClusterID);
-        log.debug("Featurenr: " +featureNbr);
+        //log.debug("Featurenr: " +featureNbr);
 
         Feature f = fList.get(featureNbr);
@@ -1071 +1028 @@
         dc.reattachItem(bestMatch);
 
-        log.debug("First feature: " +f.getClusterId() + " bestMatch feature: " +bestMatch.getClusterId());
+        //log.debug("First feature: " +f.getClusterId() + " bestMatch feature: " +bestMatch.getClusterId());
 
         bestMatch.setClusterId(f.getClusterId());
-        
+
         //updating hits
         if (f.getPeptideSequence()!= null && bestMatch.getPeptideSequence()==null){
-          updateHitsList(f,bestMatch,localSampleId2,fractionId2,allFeaturesSeq,writer,dc,factory,project);
-        }
-
-        log.debug("First feature after set: " +f.getClusterId() + " bestMatch feature: " +bestMatch.getClusterId());
+          //log.debug("Updating hit for sequence: " +f.getPeptideSequence());
+          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);
+        }
+        //TODO: what to do if both have sequences?
+
+        //log.debug("First feature after set: " +f.getClusterId() + " bestMatch feature: " +bestMatch.getClusterId());
 
       }
@@ -1087 +1048 @@
     }
     dc.commit();
+
+
+    writer.println();
+
+    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();
+
+
+    return count;
+  }
+
+
+  private int updateHitsList(Feature f, Feature f2, String localSampleId2, String fractionId2, ArrayList<Feature> allFeaturesSeq, PrintWriter writer, DbControl dc, ItemFactory factory, Project project){
+    //create a hit in database for f2 (the feature whose sequence has been updated) 
+
+    f2.setPeptideSequence(f.getPeptideSequence());
+    dc.reattachItem(f);
+    dc.reattachItem(f2);
     
+    int newHitCounter = 0;
     
-    writer.println();
-    //writer.println("Peptide sequences propagated:" + count);
-    writer.println("Alignment recall and precision based on sequences with mzTol " +tolerance +": ");
-    for (int k=0;k<timeTol.length;k++){
-      writer.println("For timeTol " +timeTol[k] +": " +fineSeq[k]/features.size() +" " +(1-fineTolSeq[k]/features.size()));
-    }
-    writer.println("Precision for all features with mzTol " +tolerance +": ");
-    for (int k=0;k<timeTol.length;k++){
-      writer.println("For timeTol " +timeTol[k] +": ");
-      writer.println("Number of total matched features: " +fineTest[k]);
-      writer.println("Number of matched features: " +fine[k] +" and recall value:  " +(fine[k]*2.0/(allFeatures.size()+allFeatures2.size())));
-      writer.println("Number of features matching more than one: " +fineTol[k] +" Precision vaule: " +(1-(fineTol[k]/(fine[k]+fineTol[k]))));
-    }
-    writer.println();
-    
-    
-    return count;
-  }
-
-  
-  private void updateHitsList(Feature f, Feature f2, String localSampleId2, String fractionId2, ArrayList<Feature> allFeaturesSeq, PrintWriter writer, DbControl dc, ItemFactory factory, Project project){
-        //create a hit in database for f2 (the feature whose sequence has been updated) 
-        
-          f2.setPeptideSequence(f.getPeptideSequence());
-          dc.reattachItem(f);
-          dc.reattachItem(f2);
-          
-            Hit h = null;
-            for (Hit h1 : f.getHits())
-            {
-              // Take the best hit matching the feature
-              if (h == null || h1.getCombinedFDR() != null && (h
-                  .getCombinedFDR() == null || (h1
-                      .getCombinedFDR() < h.getCombinedFDR())))
-              {
-                h = h1;
-              }
-            }
+    Hit h = null;
+    for (Hit h1 : f.getHits())
+    {
+      // Take the best hit matching the feature
+      if (h == null || h1.getCombinedFDR() != null && (h
+          .getCombinedFDR() == null || (h1
+              .getCombinedFDR() < h.getCombinedFDR())))
+      {
+        h = h1;
+      }
+    }
 
     /* * If there was no hit for the feature, then take
      * from another feature with the same peptide
      * sequence (ignoring delta of modifications).
-*/
-            if (h == null)
+     */
+    if (h == null)
+    {
+      String fseq = f.getPeptideSequence();
+      if (fseq.contains(" "))
+        fseq = fseq.substring(0, fseq
+            .lastIndexOf(' '));
+      for (Feature f3 : allFeaturesSeq)
+      {
+        if ((fseq.equals(f3.getPeptideSequence()) || (fseq
+            .contains(" ") && f3
+            .getPeptideSequence().startsWith(fseq))) && !f3
+            .getHits().isEmpty() && f.getChargeState().equals(f3.getChargeState()))
+        {
+          for (Hit h1 : f3.getHits())
+          {
+            if (h == null || h1
+                .getCombinedFDR() != null && (h
+                    .getCombinedFDR() == null || (h1
+                        .getCombinedFDR() < h
+                        .getCombinedFDR())))
             {
-              String fseq = f.getPeptideSequence();
-              if (fseq.contains(" "))
-                fseq = fseq.substring(0, fseq
-                  .lastIndexOf(' '));
-              for (Feature f3 : allFeaturesSeq)
-              {
-                  if ((fseq.equals(f3.getPeptideSequence()) || (fseq
-                      .contains(" ") && f3
-                      .getPeptideSequence().startsWith(fseq))) && !f3
-                      .getHits().isEmpty() && f.getChargeState().equals(f3.getChargeState()))
-                  {
-                    for (Hit h1 : f3.getHits())
-                    {
-                      if (h == null || h1
-                          .getCombinedFDR() != null && (h
-                              .getCombinedFDR() == null || (h1
-                                  .getCombinedFDR() < h
-                                  .getCombinedFDR())))
-                      {
-                        h = h1;
-                        log
-                        .debug("Assigning hit from feature at RT " + f3
-                          .getApexRetentionTimeInMinutes() + " with sequence:" + f3
-                          .getPeptideSequence());
-                      }
-                    }
-                  }
-                
-              }
+              h = h1;
+              log
+              .debug("Assigning hit from feature at RT " + f3
+                  .getApexRetentionTimeInMinutes() + " with sequence:" + f3
+                  .getPeptideSequence());
             }
-            if (h != null && localSampleId2 != null)
-            {
-              Hit newHit = factory.create(Hit.class);
-              dc.saveItem(newHit);
-              // To allow the user to trace the identification:
-              newHit.setPeakListFile(h.getPeakListFile());
-              newHit.setIdentificationResultFile(h.getIdentificationResultFile());
-              newHit.setSpectrumStringId(h.getSpectrumStringId());
-              newHit.setSpectrumId(h.getSpectrumId());
-              newHit.setCombinedFDR(h.getCombinedFDR());
-              newHit.setDescription(h.getDescription());
-              newHit.setExternalId(h.getExternalId());
-              newHit.setCharge(h.getCharge());
-              // From the sample we're matching with:
-              newHit.setFractionId(fractionId2);
-              newHit.setLocalSampleId(localSampleId2);
-              newHit.setScoreType("Proteios aligned");
-              newHit.setFeature(f2);
-              newHit.setProject(project);
-              newHit.setPrecursorQuantity(f2
-                .getIntegratedIntensity());
-              newHit.setRetentionTimeInMinutes(f2
-                .getApexRetentionTimeInMinutes());
-              newHit.setPrimaryCombined(true);
-              newHit.setProtein(false);
-              newHit.setFirstInCombination(h);
-              writer.println("New hit generated:" + newHit
-                .getDescription());
-            }
-            else
-            {
-              writer
-              .println("Feature: " + f2
-                .getMassToChargeRatio() + " sequence " + f2
-                .getPeptideSequence() + ". Local sample id: " + localSampleId2);
-              f2.setHits(f.getHits());
-            }
-            
-  }
-  
+          }
+        }
+
+      }
+    }
+    if (h != null && localSampleId2 != null)
+    {
+      Hit newHit = factory.create(Hit.class);
+      dc.saveItem(newHit);
+      // To allow the user to trace the identification:
+      newHit.setPeakListFile(h.getPeakListFile());
+      newHit.setIdentificationResultFile(h.getIdentificationResultFile());
+      newHit.setSpectrumStringId(h.getSpectrumStringId());
+      newHit.setSpectrumId(h.getSpectrumId());
+      newHit.setCombinedFDR(h.getCombinedFDR());
+      newHit.setDescription(h.getDescription());
+      newHit.setExternalId(h.getExternalId());
+      newHit.setCharge(h.getCharge());
+      // From the sample we're matching with:
+      newHit.setFractionId(fractionId2);
+      newHit.setLocalSampleId(localSampleId2);
+      newHit.setScoreType("Proteios aligned");
+      newHit.setFeature(f2);
+      newHit.setProject(project);
+      newHit.setPrecursorQuantity(f2
+          .getIntegratedIntensity());
+      newHit.setRetentionTimeInMinutes(f2
+          .getApexRetentionTimeInMinutes());
+      newHit.setPrimaryCombined(true);
+      newHit.setProtein(false);
+      newHit.setFirstInCombination(h);
+      log.debug("New hit generated:" + newHit
+          .getDescription());
+      newHitCounter++;
+    }
+    else
+    {
+      log.debug("Feature: " + f2
+          .getMassToChargeRatio() + " sequence " + f2
+          .getPeptideSequence() + ". Local sample id: " + localSampleId2);
+      f2.setHits(f.getHits());
+    }
+    
+    return newHitCounter;
+
+  }
+
   private PolynomialSplineFunction estimateSplineFunc(ArrayList<Point> points, PrintWriter writer)
   {
@@ -1216 +1174 @@
       double[] y = new double[points.size()];
       double[] weights = new double[points.size()];
-      
+
 
       int index = 0;
@@ -1227 +1185 @@
         y[index]=tempPoint.getY();
         weights[index]=tempPoint.getWeight();
-        
+
         index++;
       }
-      
+
       try
       { 
@@ -1245 +1203 @@
   }
 
-  
+
   private void setMzDistribution(List<Feature> features, double[][] mzDist, double[] meanLength, int fileNbr){
-    
+
     //TODO: get the minMz, maxMz, rtMax values automatically
-    
+
     int minMz=400;
     int maxMz=2000;
-    
+
     int mult = 100;
-    
+
     mzDist[fileNbr] = new double[(mult)*(maxMz-minMz)];
-    
+
     double sum = 0;
     double nbrSeq = 0;
-    
+
     //Double min = Double.POSITIVE_INFINITY;
     //Double max = Double.NEGATIVE_INFINITY;
-    
+
     //int minRT = Integer.MAX_VALUE;
     //int maxRT = Integer.MIN_VALUE;
-    
+
     for(Feature f: features)
     {
       int mz = (int) (Math.round((f.getMassToChargeRatio()-minMz)*mult));
-      
-      
+
+
       mzDist[fileNbr][mz]++;
-      
+
       if(f.getPeptideSequence()!=null)
       {
         nbrSeq++;
         sum+=f.getEndRetentionTimeInMinutes().doubleValue()-f.getStartRetentionTimeInMinutes().doubleValue();
-        
-      }
-    
-    }
-    
+
+      }
+
+    }
+
     meanLength[fileNbr] = sum/nbrSeq;
-    
-  }
-  
+
+  }
+
   private class SeqComparator implements Comparator<Feature>{
 
@@ -1295 +1253 @@
         return f.getChargeState().compareTo(f2.getChargeState());
       }else{
-        
+
         return f.getPeptideSequence().compareTo(f2.getPeptideSequence());
       }
-      
-    }
-    
-  }
-  
+
+    }
+
+  }
+
   private class FeatureSorter implements Comparator<Feature>{
 
@@ -1310 +1268 @@
     { 
       if (f.getPeptideSequence().equals(f2.getPeptideSequence()) && f.getChargeState().equals(f2.getChargeState())){
-        
+
         return f.getApexRetentionTimeInMinutes().compareTo(f2.getApexRetentionTimeInMinutes());
-        
+
       }else if(f.getPeptideSequence().equals(f2.getPeptideSequence())){
-        
+
         return f.getChargeState().compareTo(f2.getChargeState());
-        
+
       }else{
-        
+
         return f.getPeptideSequence().compareTo(f2.getPeptideSequence());
-        
-      }
-    
-      
-      
+
+      }
+
+
+
 
     }
@@ -1341 +1299 @@
 
   }
-  
+
   private class MzDiffComparator implements Comparator<Point>{
 
@@ -1353 +1311 @@
 
   }
-  
+
   private class ApexIntensityComparator implements Comparator<Feature>{
 
@@ -1365 +1323 @@
 
   }
-  
+
   private class IntensityComparator implements Comparator<Feature>{
 
@@ -1377 +1335 @@
 
   }
-  
-  
+
+
   private class RTComparator implements Comparator<Point>{
 
@@ -1390 +1348 @@
 
   }
-  
+
   private class RTComparatorF implements Comparator<Feature>{
 
@@ -1429 +1387 @@
 
     }
-    
-    
+
+
     public Feature getXFeature()
     {
       return this.x;
     }
-    
+
     public Feature getYFeature()
     {
       return this.y;
     }
-    
+
     public double getXLength()
     {
@@ -1460 +1418 @@
       return this.y.getApexRetentionTimeInMinutes().doubleValue();
     }
-    
+
     public double getMzX()
     {
@@ -1486 +1444 @@
       return Math.abs(this.getX()-this.getY());
     }
-    
+
     public double getMzDiff()
     {
@@ -1525 +1483 @@
     private int file2 = 0;
     private ArrayList<Point> evalPoint = new ArrayList<Point>();
-    
+
     public FilePoint()
     {}
-    
+
     public FilePoint(int fileNbr, int fileNbr2, double simValue, PolynomialSplineFunction func, ArrayList<Point> evalPoint)
     { 
@@ -1558 +1516 @@
       return this.func;
     }
-    
+
     public ArrayList<Point> getEvalList() {
-      
+
       return this.evalPoint;
     }