From ec6713a1c09d121c05a4ad7dd29774f49b78df5e Mon Sep 17 00:00:00 2001
From: eugenefischer <66030419+eugenefischer@users.noreply.github.com>
Date: Thu, 29 Sep 2022 15:50:57 -0500
Subject: [PATCH] Implement filtering for wells with anomalous read counts

---
 src/main/java/Simulator.java | 23 +++++++++++++++++++++++
 1 file changed, 23 insertions(+)

diff --git a/src/main/java/Simulator.java b/src/main/java/Simulator.java
index dcea030..286e009 100644
--- a/src/main/java/Simulator.java
+++ b/src/main/java/Simulator.java
@@ -69,6 +69,12 @@ public class Simulator implements GraphModificationFunctions {
             if(verbose){System.out.println("Remaining alpha sequence count: " + alphaSequences.size());}
             if(verbose){System.out.println("Remaining beta sequence count: " + betaSequences.size());}
         }
+        if (realSequenceCollisionRate > 0.0) {
+            if(verbose){System.out.println("Removing wells with anomalous read counts from sequence records");}
+            filterWellsByReadCount(alphaSequences);
+            filterWellsByReadCount(betaSequences);
+            if(verbose){System.out.println("Wells with anomalous read counts removed from sequence records");}
+        }
 
         //construct the graph. For simplicity, going to make
         if(verbose){System.out.println("Making vertex maps");}
@@ -673,6 +679,23 @@ public class Simulator implements GraphModificationFunctions {
         }
     }
 
+    public static void filterWellsByReadCount(Map<String, SequenceRecord> sequences) {
+        for (String k: sequences.keySet()) {
+            //If a sequence has read count R and appears in W wells, then on average its read count in each
+            //well should be R/W. Delete any wells where the read count is less than R/2W.
+            Integer threshold = sequences.get(k).getReadCount() / (2 * sequences.get(k).getOccupancy());
+            List<Integer> noise = new ArrayList<>();
+            for (Integer well: sequences.get(k).getWells()) {
+                if (sequences.get(k).getReadCount(well) < threshold) {
+                    noise.add(well);
+                }
+            }
+            for (Integer well: noise) {
+                sequences.get(k).deleteWell(well);
+            }
+        }
+    }
+
     private static Map<String, String> makeSequenceToSequenceMap(List<String[]> cells, int keySequenceIndex,
                                                                  int valueSequenceIndex){
         Map<String, String> keySequenceToValueSequenceMap = new HashMap<>();