custom vc

sa501428 · Nov 7, 2024 · 07cb2df · 07cb2df
1 parent e0b210b
commit 07cb2df
Show file tree

Hide file tree

Showing 3 changed files with 28 additions and 5 deletions.
diff --git a/src/cli/Main.java b/src/cli/Main.java
@@ -22,7 +22,7 @@
 
 public class Main {
 
-    public static final String VERSION_NUM = "0.172.0";
+    public static final String VERSION_NUM = "0.173.0";
     public static boolean printVerboseComments = false;
 
     public static void printGeneralUsageAndExit(int exitCode, String cUsage) {

diff --git a/src/cli/clt/CommandLineParser.java b/src/cli/clt/CommandLineParser.java
@@ -51,6 +51,7 @@ public class CommandLineParser extends CmdLineParser {
     private final Option minDisValOption = addIntegerOption("min-dist");
     private final Option maxDistValOption = addIntegerOption("max-dist");
     private final Option interChromosomalOption = addBooleanOption("include-inter");
+    private final Option buildVCNorms = addBooleanOption("local-vc");
     private final Option onlyOneOption = addBooleanOption("only-one");
     private final Option noOrientationOption = addBooleanOption("no-orientation");
     private final Option aggregateNormalization = addBooleanOption("ag-norm");
@@ -154,6 +155,10 @@ public boolean getIncludeInterChromosomal() {
         return optionToBoolean(interChromosomalOption);
     }
 
+    public boolean getBuildCustomVCNorms() {
+        return optionToBoolean(buildVCNorms);
+    }
+
     public int getPercentileOption(int val) {
         return optionToInteger(percentileOption, val);
     }

diff --git a/src/cli/clt/apa/APA3.java b/src/cli/clt/apa/APA3.java
@@ -66,6 +66,7 @@ public class APA3 {
     private final int minPeakDist; // distance between two bins, can be changed in opts
     private final int maxPeakDist;
     private final boolean includeInterChr;
+    private final boolean useCustomBuiltNorms;
     private final boolean useAgNorm;
     private final int numOutputs;
 
@@ -75,7 +76,6 @@ public class APA3 {
     private final List<double[]> globalRowSums;
     private final List<double[]> globalColSums;
 
-    private final AtomicInteger currChromPair = new AtomicInteger(0);
     private final NormalizationType vcNorm = NormalizationHandler.VC;
     private final Object key = new Object();
     private final HiCZoom zoom;
@@ -126,6 +126,7 @@ public APA3(String[] args, CommandLineParser parser) {
         minPeakDist = parser.getMinDistVal(2 * window);
         maxPeakDist = parser.getMaxDistVal(Integer.MAX_VALUE);
         includeInterChr = parser.getIncludeInterChromosomal();
+        useCustomBuiltNorms = parser.getBuildCustomVCNorms();
 
         matrixWidthL = 2 * window + 1;
         globalAPAMatrices = initializeMatrices(numOutputs, matrixWidthL);
@@ -170,6 +171,7 @@ public void run() {
         Map<Integer, RegionConfiguration> chromosomePairs = new HashMap<>();
         int pairCounter = HiCUtils.populateChromosomePairs(chromosomePairs,
                 handler.getChromosomeArrayWithoutAllByAll(), includeInterChr);
+        final AtomicInteger currChromPair = new AtomicInteger(0);
 
         ParallelizationTools.launchParallelizedCode(numThreads, () -> {
 
@@ -192,9 +194,16 @@ public void run() {
                         MatrixZoomData zd = matrix.getZoomData(zoom);
                         if (zd != null) {
                             try {
-                                processLoopsForRegion(zd, loops, outputs, distanceBoundCalculator);
+                                double[] customVCNormVector1 = new double[(int) (chr1.getLength() / resolution) + 1];
+                                double[] customVCNormVector2 = new double[(int) (chr2.getLength() / resolution) + 1];
+                                processLoopsForRegion(zd, loops, outputs, distanceBoundCalculator,
+                                        customVCNormVector1, customVCNormVector2);
                                 if (useAgNorm) {
-                                    doAggregateNormalization(chr1, chr2, zoom, vcNorm, loops, rowSums, colSums);
+                                    if (useCustomBuiltNorms) {
+                                        doAggregateNormalization(loops, rowSums, colSums, customVCNormVector1, customVCNormVector2);
+                                    } else {
+                                        doAggregateNormalization(chr1, chr2, zoom, vcNorm, loops, rowSums, colSums);
+                                    }
                                 }
                                 System.out.println("Completed " + chr1.getName() + "_" + chr2.getName());
                             } catch (Exception e) {
@@ -243,6 +252,12 @@ private void doAggregateNormalization(Chromosome chr1, Chromosome chr2, HiCZoom
             vector2 = ds.getNormalizationVector(chr2.getIndex(), zoom, vcNorm).getData().getValues().get(0);
         }
 
+        doAggregateNormalization(cloops, rowSums, colSums, vector1, vector2);
+    }
+
+    private void doAggregateNormalization(List<List<Feature2D>> cloops, List<double[]> rowSums, List<double[]> colSums,
+                                          double[] vector1, double[] vector2) {
+
         for (int i = 0; i < cloops.size(); i++) {
             List<Feature2D> loops = cloops.get(i);
             double[] rowSum = rowSums.get(i);
@@ -258,14 +273,17 @@ private void doAggregateNormalization(Chromosome chr1, Chromosome chr2, HiCZoom
     }
 
     protected void processLoopsForRegion(MatrixZoomData zd, List<List<Feature2D>> allLoops,
-                                         List<float[][]> outputs, DistanceBoundCalculator distanceBoundCalculator) {
+                                         List<float[][]> outputs, DistanceBoundCalculator distanceBoundCalculator,
+                                         double[] customVCNormVector1, double[] customVCNormVector2) {
         RegionsOfInterest roi = new RegionsOfInterest(resolution, window, matrixWidthL, allLoops);
         int counter = 0;
 
         Iterator<ContactRecord> it = ExpectedUtils.getIterator(zd, norm);
         while (it.hasNext()) {
             ContactRecord cr = it.next();
             if (cr.getCounts() > 0) {
+                customVCNormVector1[cr.getBinX()] += cr.getCounts();
+                customVCNormVector2[cr.getBinY()] += cr.getCounts();
                 if (distanceBoundCalculator.inDistanceRange(cr)) {
                     if (roi.probablyContainsRecord(cr)) {
                         for (int i = 0; i < allLoops.size(); i++) {