update expected model

sa501428 · Aug 27, 2022 · 3857055 · 3857055
1 parent 29e4cdf
commit 3857055
Show file tree

Hide file tree

Showing 2 changed files with 44 additions and 19 deletions.
diff --git a/src/javastraw/expected/ExpectedModel.java b/src/javastraw/expected/ExpectedModel.java
@@ -1,9 +1,25 @@
 package javastraw.expected;
 
 import javastraw.reader.block.ContactRecord;
+import javastraw.reader.mzd.MatrixZoomData;
+import javastraw.reader.type.NormalizationType;
+
+import java.util.Iterator;
 
 public abstract class ExpectedModel {
 
+    public static int getDist(ContactRecord record) {
+        return Math.abs(record.getBinX() - record.getBinY());
+    }
+
+    public static Iterator<ContactRecord> getIterator(MatrixZoomData zd, NormalizationType norm) {
+        if (norm.getLabel().equalsIgnoreCase("none")) {
+            return zd.getDirectIterator();
+        } else {
+            return zd.getNormalizedIterator(norm);
+        }
+    }
+
     public static float getP(double obs, double expected, double superDiagonal) {
         // P = (O - E)/(SD - E)
         return (float) ((obs - expected) / (superDiagonal - expected));

diff --git a/src/javastraw/expected/LogExpectedSpline.java b/src/javastraw/expected/LogExpectedSpline.java
@@ -1,5 +1,6 @@
 package javastraw.expected;
 
+import javastraw.reader.basics.Chromosome;
 import javastraw.reader.block.ContactRecord;
 import javastraw.reader.mzd.MatrixZoomData;
 import javastraw.reader.type.NormalizationType;
@@ -18,15 +19,17 @@ public class LogExpectedSpline extends ExpectedModel {
     private final double nearDiagonalSignal;
     private final float maxX;
 
-    public LogExpectedSpline(MatrixZoomData zd, NormalizationType norm, int maxBin) {
-        maxX = logp1i(maxBin);
-        function = fitDataToFunction(zd, norm, maxBin);
+    public LogExpectedSpline(MatrixZoomData zd, NormalizationType norm, Chromosome chrom, int res) {
+        int maxBin = (int) (chrom.getLength() / res + 1);
+        int[] maxDist = new int[1];
+        function = fitDataToFunction(zd, norm, maxBin, maxDist);
+        maxX = logp1i(Math.min(maxBin, maxDist[0]));
         nearDiagonalSignal = Math.expm1(function.value(0.5));
     }
 
-    private PolynomialSplineFunction fitDataToFunction(MatrixZoomData zd, NormalizationType norm, int maxBin) {
+    private PolynomialSplineFunction fitDataToFunction(MatrixZoomData zd, NormalizationType norm, int maxBin, int[] maxDist) {
 
-        List<double[]> points = getAverageInEachBin(zd, norm, maxBin);
+        List<double[]> points = getAverageInEachBin(zd, norm, maxBin, maxDist);
         double[] x = new double[points.size()];
         double[] y = new double[points.size()];
         for (int i = 0; i < points.size(); i++) {
@@ -39,25 +42,28 @@ private PolynomialSplineFunction fitDataToFunction(MatrixZoomData zd, Normalizat
         return interpolator.interpolate(x, y);
     }
 
-    private List<double[]> getAverageInEachBin(MatrixZoomData zd, NormalizationType norm, int maxBin) {
+    private List<double[]> getAverageInEachBin(MatrixZoomData zd, NormalizationType norm, int maxBin, int[] maxDist) {
 
         double[] initExpected = new double[maxBin];
         long[] countsPerBin = new long[maxBin];
-        populateWithCounts(zd, norm, initExpected, countsPerBin, maxBin);
+        populateWithCounts(zd, norm, initExpected, countsPerBin, maxBin, maxDist);
 
-        List<double[]> points = collapseToSetOfPoints(initExpected, countsPerBin, maxBin);
+        int maxDistToUse = Math.min(maxBin, maxDist[0]);
 
-        List<double[]> finalPoints = new ArrayList<>(points.size());
+        List<double[]> currentPoints = collapseToSetOfPoints(initExpected, countsPerBin, maxDistToUse);
+        List<double[]> finalPoints = new ArrayList<>(currentPoints.size());
         initExpected = null;
         countsPerBin = null;
 
-        for (double[] current : points) {
-            double[] point = new double[2];
-            point[0] = current[2] / current[3];
-            point[1] = current[0] / current[1];
-            finalPoints.add(point);
+        for (double[] current : currentPoints) {
+            if (current[3] > 0 && current[1] > 0) {
+                double[] finalPoint = new double[2];
+                finalPoint[0] = current[2] / current[3];
+                finalPoint[1] = current[0] / current[1];
+                finalPoints.add(finalPoint);
+            }
         }
-        points.clear();
+        currentPoints.clear();
 
         return finalPoints;
     }
@@ -67,7 +73,9 @@ private List<double[]> collapseToSetOfPoints(double[] initExpected, long[] count
         double[] latest = new double[4]; // vals, counts, distances, num_distances
         int k = 0;
         int numToGroup = 1;
-        while (k < maxBin) {
+        //System.out.println(Arrays.toString(initExpected));
+        //System.out.println(Arrays.toString(countsPerBin));
+        while (k < initExpected.length) {
             latest[0] += initExpected[k];
             latest[1] += countsPerBin[k];
             latest[2] += logp1(k);
@@ -97,12 +105,13 @@ private List<double[]> collapseToSetOfPoints(double[] initExpected, long[] count
     }
 
     private void populateWithCounts(MatrixZoomData zd, NormalizationType norm, double[] initExpected,
-                                    long[] countsPerBin, int maxBin) {
-        Iterator<ContactRecord> records = ExpectedUtils.getIterator(zd, norm);
+                                    long[] countsPerBin, int maxBin, int[] maxDist) {
+        Iterator<ContactRecord> records = getIterator(zd, norm);
         while (records.hasNext()) {
             ContactRecord record = records.next();
-            int dist = ExpectedUtils.getDist(record);
+            int dist = getDist(record);
             if (dist < maxBin) {
+                maxDist[0] = Math.max(maxDist[0], dist);
                 initExpected[dist] += logp1(record.getCounts());
                 countsPerBin[dist]++;
             }