Implement read count for vertices

src/main/java/GraphModificationFunctions.java

@@ -12,7 +12,6 @@ public interface GraphModificationFunctions {
     static Map<Vertex[], Integer> filterByOverlapThresholds(SimpleWeightedGraph<Vertex, DefaultWeightedEdge> graph,
                                                             int low, int high, boolean saveEdges) {
         Map<Vertex[], Integer> removedEdges = new HashMap<>();
-        //List<Integer[]> removedEdges = new ArrayList<>();
         for (DefaultWeightedEdge e : graph.edgeSet()) {
             if ((graph.getEdgeWeight(e) > high) || (graph.getEdgeWeight(e) < low)) {
                 if(saveEdges) {
@@ -94,6 +93,37 @@ public interface GraphModificationFunctions {
         return removedEdges;
     }
 
+    static Map<Vertex[], Integer> filterByRelativeReadCount (SimpleWeightedGraph<Vertex, DefaultWeightedEdge> graph, Integer threshold, boolean saveEdges) {
+        Map<Vertex[], Integer> removedEdges = new HashMap<>();
+        Boolean passes;
+        for (DefaultWeightedEdge e : graph.edgeSet()) {
+            Integer alphaReadCount = graph.getEdgeSource(e).getReadCount();
+            Integer betaReadCount = graph.getEdgeTarget(e).getReadCount();
+            passes = RelativeReadCountFilterFunction(threshold, alphaReadCount, betaReadCount);
+            if (!passes) {
+                if (saveEdges) {
+                    Vertex source = graph.getEdgeSource(e);
+                    Vertex target = graph.getEdgeTarget(e);
+                    Integer intWeight = (int) graph.getEdgeWeight(e);
+                    Vertex[] edge = {source, target};
+                    removedEdges.put(edge, intWeight);
+                }
+                else {
+                    graph.setEdgeWeight(e, 0.0);
+                }
+            }
+        }
+        if(saveEdges) {
+            for (Vertex[] edge : removedEdges.keySet()) {
+                graph.removeEdge(edge[0], edge[1]);
+            }
+        }
+        return removedEdges;
+    }
+
+    static Boolean RelativeReadCountFilterFunction(Integer threshold, Integer alphaReadCount, Integer betaReadCount) {
+        return Math.abs(alphaReadCount - betaReadCount) < threshold;
+    }
     static void addRemovedEdges(SimpleWeightedGraph<Vertex, DefaultWeightedEdge> graph,
                                 Map<Vertex[], Integer> removedEdges) {
         for (Vertex[] edge : removedEdges.keySet()) {
@@ -102,4 +132,6 @@ public interface GraphModificationFunctions {
         }
     }
 
+
+
 }

src/main/java/GraphWithMapData.java

@@ -15,6 +15,7 @@ public class GraphWithMapData implements java.io.Serializable {
     private Integer[] wellPopulations;
     private Integer alphaCount;
     private Integer betaCount;
+    private int readDepth;
     private final Map<String, String> distCellsMapAlphaKey;
 //    private final Map<Integer, Integer> plateVtoAMap;
 //    private final Map<Integer, Integer> plateVtoBMap;
@@ -25,7 +26,7 @@ public class GraphWithMapData implements java.io.Serializable {
     private final Duration time;
 
     public GraphWithMapData(SimpleWeightedGraph graph, Integer numWells, Integer[] wellConcentrations,
-                            Map<String, String> distCellsMapAlphaKey, Integer alphaCount, Integer betaCount, Duration time){
+                            Map<String, String> distCellsMapAlphaKey, Integer alphaCount, Integer betaCount, int readDepth, Duration time){
 
 //                            Map<Integer, Integer> plateVtoAMap,
 //                            Map<Integer,Integer> plateVtoBMap, Map<Integer, Integer> plateAtoVMap,
@@ -94,6 +95,8 @@ public class GraphWithMapData implements java.io.Serializable {
 //        return betaWellCounts;
 //    }
 
+    public int getReadDepth() { return readDepth; }
+
     public Duration getTime() {
         return time;
     }

src/main/java/InteractiveInterface.java

@@ -273,6 +273,9 @@ public class InteractiveInterface {
             if (simulateReadDepth) {
                 System.out.print("\nPlease enter read depth (the integer number of reads per sequence): ");
                 readDepth = sc.nextInt();
+                if(readDepth < 1) {
+                    throw new InputMismatchException("The read depth must be an integer >= 1");
+                }
                 System.out.print("\nPlease enter probability of a sequence read error (0.0 to 1.0): ");
                 readErrorRate = sc.nextDouble();
                 if(readErrorRate < 0.0 || readErrorRate > 1.0) {
@@ -332,7 +335,7 @@ public class InteractiveInterface {
             System.out.println("Returning to main menu.");
         }
         else{
-            GraphWithMapData data = Simulator.makeGraph(cellSample, plate, true);
+            GraphWithMapData data = Simulator.makeGraph(cellSample, plate, readDepth, readErrorRate, errorCollisionRate, true);
             assert filename != null;
             if(BiGpairSEQ.outputBinary()) {
                 GraphDataObjectWriter dataWriter = new GraphDataObjectWriter(filename, data);

src/main/java/Plate.java

@@ -178,17 +178,30 @@ public class Plate {
         }
     }
 
+    //returns a map of the counts of the sequence at cell index sIndex, in a specific of wells
+    //Simulates read depth and read errors, counts the number of reads of a unique sequence into the given map.
+    public void assayWellsSequenceSWithReadDepth(Map<String, Integer> occupancyMap, Map<String, Integer> readCountMap,
+                                                 int readDepth, double readErrorProb, double errorCollisionProb, int... sIndices) {
+        this.assayWellsSequenceSWithReadDepth(occupancyMap, readCountMap, readDepth, readErrorProb, errorCollisionProb, 0, size, sIndices);
+    }
+    //returns a map of the counts of the sequence at cell index sIndex, in a specific of wells
+    //Simulates read depth and read errors, counts the number of reads of a unique sequence into the given map.
+    public void assayWellsSequenceSWithReadDepth(Map<String, Integer> occupancyMap, Map<String, Integer> readCountMap,
+                                                 int readDepth, double readErrorProb, double errorCollisionProb,
+                                                 int n, int... sIndices) {
+       this.assayWellsSequenceSWithReadDepth(occupancyMap, readCountMap, readDepth, readErrorProb, errorCollisionProb, n, n+1, sIndices);
+    }
 
     //returns a map of the counts of the sequence at cell index sIndex, in a range of wells
     //Simulates read depth and read errors, counts the number of reads of a unique sequence into the given map.
-    public Map<String, Integer> assayWellsSequenceSWithReadDepth(int start, int end, int... sIndices) {
+    public void assayWellsSequenceSWithReadDepth(Map<String, Integer> occupancyMap, Map<String, Integer> readCountMap,
-        Map<String, Integer> assay = new HashMap<>();
+                                                 int readDepth, double readErrorProb, double errorCollisionProb,
+                                                 int start, int end, int... sIndices) {
         for(int sIndex: sIndices){
             for(int i = start; i < end; i++){
-                countSequences(assay, wells.get(i), sIndex);
+                countSequencesWithReadDepth(occupancyMap, readCountMap, readDepth, readErrorProb, errorCollisionProb, wells.get(i), sIndex);
             }
         }
-        return assay;
     }
     //For the sequences at cell indices sIndices, counts number of unique sequences in the given well into the given map
     //Simulates read depth and read errors, counts the number of reads of a unique sequence into the given map.

src/main/java/Simulator.java

@@ -21,8 +21,12 @@ public class Simulator implements GraphModificationFunctions {
     //sourceVertexIndices and targetVertexIndices are indices within the cell to use as for the two sets of vertices
     //in the bipartite graph. "Source" and "target" are JGraphT terms for the two vertices an edge touches,
     //even if not directed.
-    public static GraphWithMapData makeGraph(CellSample cellSample, Plate samplePlate, boolean verbose) {
+    public static GraphWithMapData makeGraph(CellSample cellSample, Plate samplePlate, int readDepth, double readErrorRate, double errorCollisionRate, boolean verbose) {
         Instant start = Instant.now();
+        Map<String, Integer> allAlphas;
+        Map<String, Integer> allBetas;
+        Map<String, Integer> alphaReadCounts = null;
+        Map<String, Integer> betaReadCounts = null;
         List<String[]> distinctCells = cellSample.getCells();
         int[] alphaIndices = {SequenceType.CDR3_ALPHA.ordinal()};
         int[] betaIndices = {SequenceType.CDR3_BETA.ordinal()};
@@ -35,11 +39,21 @@ public class Simulator implements GraphModificationFunctions {
         if(verbose){System.out.println("Cell maps made");}
 
         if(verbose){System.out.println("Making well maps");}
+        if(readDepth == 1) {
+            allAlphas = new HashMap<>();
+            samplePlate.assayWellsSequenceS(allAlphas, alphaIndices);
+            allBetas = new HashMap<>();
+            samplePlate.assayWellsSequenceS(allBetas, betaIndices);
+        }
+        else {
+            allAlphas = new HashMap<>();
+            alphaReadCounts = new HashMap<>();
+            samplePlate.assayWellsSequenceSWithReadDepth(allAlphas, alphaReadCounts, readDepth, readErrorRate, errorCollisionRate, alphaIndices);
+            allBetas = new HashMap<>();
+            betaReadCounts = new HashMap<>();
+            samplePlate.assayWellsSequenceSWithReadDepth(allBetas, betaReadCounts, readDepth, readErrorRate, errorCollisionRate, betaIndices);
+        }
 
-        Map<String, Integer> allAlphas = new HashMap<>();
-        samplePlate.assayWellsSequenceS(allAlphas, alphaIndices);
-        Map<String, Integer> allBetas = new HashMap<>();
-        samplePlate.assayWellsSequenceS(allBetas, betaIndices);
         int alphaCount = allAlphas.size();
         if(verbose){System.out.println("All alphas count: " + alphaCount);}
         int betaCount = allBetas.size();
@@ -101,7 +115,13 @@ public class Simulator implements GraphModificationFunctions {
         List<Vertex> alphaVertices = new ArrayList<>();
         //start with map of all alphas mapped to vertex values, get occupancy from the alphaWellCounts map
         for (String seq : plateAtoVMap.keySet()) {
-            Vertex alphaVertex = new Vertex(SequenceType.CDR3_ALPHA, seq, alphaWellCounts.get(seq), plateAtoVMap.get(seq));
+            Vertex alphaVertex;
+            if (readDepth == 1) {
+                alphaVertex = new Vertex(SequenceType.CDR3_ALPHA, seq, alphaWellCounts.get(seq), plateAtoVMap.get(seq));
+            }
+            else {
+                alphaVertex = new Vertex(SequenceType.CDR3_ALPHA, seq, alphaWellCounts.get(seq), plateAtoVMap.get(seq), alphaReadCounts.get(seq));
+            }
             alphaVertices.add(alphaVertex);
         }
         //Sort to make sure the order of vertices in list matches the order of the adjacency matrix
@@ -112,7 +132,13 @@ public class Simulator implements GraphModificationFunctions {
         //List<Integer> betaVertices = new ArrayList<>(plateVtoBMap.keySet());//This will work because LinkedHashMap preserves order of entry
         List<Vertex> betaVertices = new ArrayList<>();
         for (String seq : plateBtoVMap.keySet()) {
-            Vertex betaVertex = new Vertex(SequenceType.CDR3_BETA, seq, betaWellCounts.get(seq), plateBtoVMap.get(seq));
+            Vertex betaVertex;
+            if (readDepth == 1) {
+                betaVertex = new Vertex(SequenceType.CDR3_BETA, seq, betaWellCounts.get(seq), plateBtoVMap.get(seq));
+            }
+            else {
+                betaVertex = new Vertex(SequenceType.CDR3_BETA, seq, betaWellCounts.get(seq), plateBtoVMap.get(seq), betaReadCounts.get(seq));
+            }
             betaVertices.add(betaVertex);
         }
         //Sort to make sure the order of vertices in list matches the order of the adjacency matrix
@@ -128,7 +154,7 @@ public class Simulator implements GraphModificationFunctions {
         Duration time = Duration.between(start, stop);
 
         //create GraphWithMapData object
-        GraphWithMapData output = new GraphWithMapData(graph, numWells, samplePlate.getPopulations(), distCellsMapAlphaKey, alphaCount, betaCount, time);
+        GraphWithMapData output = new GraphWithMapData(graph, numWells, samplePlate.getPopulations(), distCellsMapAlphaKey, alphaCount, betaCount, readDepth, time);
         //Set source file name in graph to name of sample plate
         output.setSourceFilename(samplePlate.getFilename());
         //return GraphWithMapData object

src/main/java/Vertex.java

@@ -23,6 +23,7 @@ public class Vertex implements Serializable, Comparable<Vertex> {
         this.vertexLabel = vertexLabel;
         this.sequence = sequence;
         this.occupancy = occupancy;
+        this.readCount = 1;
     }
 
     public Vertex(SequenceType type, String sequence, Integer occupancy, Integer vertexLabel, Integer readCount) {