update version number

readme.md

@@ -347,15 +347,21 @@ roughly as though it had a constant well population equal to the plate's average
   * ~~Apache Commons CSV library writes entries a row at a time~~ 
     * _Got this working, but at the cost of a profoundly strange bug in graph occupancy filtering. Have reverted the repo until I can figure out what caused that. Given how easily Thingiverse transposes CSV matrices in R, might not even be worth fixing.
 * ~~Enable GraphML output in addition to serialized object binaries, for data portability~~ DONE
-  * ~~Custom vertex type with attribute for sequence occupancy?~~ ABANDONED
+  * ~~Custom vertex type with attribute for sequence occupancy?~~ DONE
     * Advantage: would eliminate the need to use maps to associate vertices with sequences, which would make the code easier to understand.
-    * Have a branch where this is implemented, but there's a bug that broke matching. Don't currently have time to fix.
+    * ~~Have a branch where this is implemented, but there's a bug that broke matching. Don't currently have time to fix.~~
 * ~~Re-implement command line arguments, to enable scripting and statistical simulation studies~~ DONE
 * Re-implement CDR1 matching method
 * Implement Duan and Su's maximum weight matching algorithm
     * Add controllable algorithm-type parameter?
     * This would be fun and valuable, but probably take more time than I have for a hobby project.
+* Implement an algorithm for approximating a maximum weight matching
+  * Some of these run in linear or near-linear time
+  * given that the underlying biological samples have many, many sources of error, this would probably be the most useful option in practice. It seems less mathematically elegant, though, and so less fun for me.
 * Implement Vose's alias method for arbitrary statistical distributions of cells
+  * Should probably refactor to use apache commons rng for this
+* Use commons JCS for caching
+* Enable post-filtering instead of pre-filtering. Pre-filtering of things like singleton sequences or saturating-occupancy sequences reduces graph size, but could conceivably reduce pairing accuracy by throwing away data. While these sequences have very little signal, it would be interesting to compare unfiltered results to filtered results. This would require a much, much faster MWM algorithm, though, to handle the much larger graphs. Possible one of the linear-time approximation algorithms.
   
   
 ## CITATIONS

src/main/java/BiGpairSEQ.java

@@ -16,7 +16,7 @@ public class BiGpairSEQ {
     private static String priorityQueueHeapType = "FIBONACCI";
     private static boolean outputBinary = true;
     private static boolean outputGraphML = false;
-    private static final String version = "version 2.0";
+    private static final String version = "version 3.0";
 
     public static void main(String[] args) {
         if (args.length == 0) {

src/main/java/GraphMLFileWriter.java

@@ -3,8 +3,9 @@ import org.jgrapht.graph.SimpleWeightedGraph;
 import org.jgrapht.nio.Attribute;
 import org.jgrapht.nio.AttributeType;
 import org.jgrapht.nio.DefaultAttribute;
-import org.jgrapht.nio.dot.DOTExporter;
 import org.jgrapht.nio.graphml.GraphMLExporter;
+import org.jgrapht.nio.graphml.GraphMLExporter.AttributeCategory;
+import org.w3c.dom.Attr;
 
 import java.io.BufferedWriter;
 import java.io.IOException;
@@ -12,14 +13,14 @@ import java.nio.file.Files;
 import java.nio.file.Path;
 import java.nio.file.StandardOpenOption;
 import java.util.HashMap;
-import java.util.LinkedHashMap;
 import java.util.Map;
 
 public class GraphMLFileWriter {
 
     String filename;
+    SimpleWeightedGraph graph;
     GraphWithMapData data;
-
+    Map<String, Attribute> graphAttributes;
 
     public GraphMLFileWriter(String filename, GraphWithMapData data) {
         if(!filename.matches(".*\\.graphml")){
@@ -27,52 +28,61 @@ public class GraphMLFileWriter {
         }
         this.filename = filename;
         this.data = data;
+        this.graph = data.getGraph();
+        graphAttributes = createGraphAttributes();
     }
 
-//    public void writeGraphToFile() {
-//        try(BufferedWriter writer = Files.newBufferedWriter(Path.of(filename), StandardOpenOption.CREATE_NEW);
-//        ){
-//            GraphMLExporter<SimpleWeightedGraph, BufferedWriter> exporter = new GraphMLExporter<>();
-//            exporter.exportGraph(graph, writer);
-//        } catch(IOException ex){
-//            System.out.println("Could not make new file named "+filename);
-//            System.err.println(ex);
-//        }
-//    }
+    public GraphMLFileWriter(String filename, SimpleWeightedGraph<Vertex, DefaultWeightedEdge> graph) {
+        if(!filename.matches(".*\\.graphml")){
+            filename = filename + ".graphml";
+        }
+        this.filename = filename;
+        this.graph = graph;
+    }
+
+    private Map<String, Attribute> createGraphAttributes(){
+        Map<String, Attribute> ga = new HashMap<>();
+        //Sample plate filename
+        ga.put("sample plate filename", DefaultAttribute.createAttribute(data.getSourceFilename()));
+        // Number of wells
+        ga.put("well count", DefaultAttribute.createAttribute(data.getNumWells().toString()));
+        //Well populations
+        Integer[] wellPopulations = data.getWellPopulations();
+        StringBuilder populationsStringBuilder = new StringBuilder();
+        populationsStringBuilder.append(wellPopulations[0].toString());
+        for(int i = 1; i < wellPopulations.length; i++){
+            populationsStringBuilder.append(", ");
+            populationsStringBuilder.append(wellPopulations[i].toString());
+        }
+        String wellPopulationsString = populationsStringBuilder.toString();
+        ga.put("well populations", DefaultAttribute.createAttribute(wellPopulationsString));
+        return ga;
+    }
 
     public void writeGraphToFile() {
-        SimpleWeightedGraph graph = data.getGraph();
-        Map<Integer, Integer> vertexToAlphaMap = data.getPlateVtoAMap();
-        Map<Integer, Integer> vertexToBetaMap = data.getPlateVtoBMap();
-        Map<Integer, Integer> alphaOccs = data.getAlphaWellCounts();
-        Map<Integer, Integer> betaOccs = data.getBetaWellCounts();
         try(BufferedWriter writer = Files.newBufferedWriter(Path.of(filename), StandardOpenOption.CREATE_NEW);
         ){
             //create exporter. Let the vertex labels be the unique ids for the vertices
-            GraphMLExporter<Integer, SimpleWeightedGraph<Vertex, DefaultWeightedEdge>> exporter = new GraphMLExporter<>(v -> v.toString());
+            GraphMLExporter<Vertex, SimpleWeightedGraph<Vertex, DefaultWeightedEdge>> exporter = new GraphMLExporter<>(v -> v.getVertexLabel().toString());
             //set to export weights
             exporter.setExportEdgeWeights(true);
+            //Set graph attributes
+            exporter.setGraphAttributeProvider( () -> graphAttributes);
             //set type, sequence, and occupancy attributes for each vertex
             exporter.setVertexAttributeProvider( v -> {
                 Map<String, Attribute> attributes = new HashMap<>();
-                if(vertexToAlphaMap.containsKey(v)) {
-                    attributes.put("type", DefaultAttribute.createAttribute("CDR3 Alpha"));
-                    attributes.put("sequence", DefaultAttribute.createAttribute(vertexToAlphaMap.get(v)));
-                    attributes.put("occupancy", DefaultAttribute.createAttribute(
-                            alphaOccs.get(vertexToAlphaMap.get(v))));
-                }
-                else if(vertexToBetaMap.containsKey(v)) {
-                    attributes.put("type", DefaultAttribute.createAttribute("CDR3 Beta"));
-                    attributes.put("sequence", DefaultAttribute.createAttribute(vertexToBetaMap.get(v)));
-                    attributes.put("occupancy", DefaultAttribute.createAttribute(
-                            betaOccs.get(vertexToBetaMap.get(v))));
-                }
+                attributes.put("type", DefaultAttribute.createAttribute(v.getType().name()));
+                attributes.put("sequence", DefaultAttribute.createAttribute(v.getSequence()));
+                attributes.put("occupancy", DefaultAttribute.createAttribute(v.getOccupancy()));
                 return attributes;
             });
             //register the attributes
-            exporter.registerAttribute("type", GraphMLExporter.AttributeCategory.NODE, AttributeType.STRING);
-            exporter.registerAttribute("sequence", GraphMLExporter.AttributeCategory.NODE, AttributeType.STRING);
-            exporter.registerAttribute("occupancy", GraphMLExporter.AttributeCategory.NODE, AttributeType.STRING);
+            for(String s : graphAttributes.keySet()) {
+                exporter.registerAttribute(s, AttributeCategory.GRAPH, AttributeType.STRING);
+            }
+            exporter.registerAttribute("type", AttributeCategory.NODE, AttributeType.STRING);
+            exporter.registerAttribute("sequence", AttributeCategory.NODE, AttributeType.STRING);
+            exporter.registerAttribute("occupancy", AttributeCategory.NODE, AttributeType.STRING);
             //export the graph
             exporter.exportGraph(graph, writer);
         } catch(IOException ex){
@@ -81,4 +91,3 @@ public class GraphMLFileWriter {
         }
     }
 }
-

src/main/java/GraphModificationFunctions.java

@@ -2,23 +2,25 @@ import org.jgrapht.graph.DefaultWeightedEdge;
 import org.jgrapht.graph.SimpleWeightedGraph;
 
 import java.util.ArrayList;
+import java.util.HashMap;
 import java.util.List;
 import java.util.Map;
 
 public interface GraphModificationFunctions {
 
     //remove over- and under-weight edges
-    static List<Integer[]> filterByOverlapThresholds(SimpleWeightedGraph<Integer, DefaultWeightedEdge> graph,
+    static Map<Vertex[], Integer> filterByOverlapThresholds(SimpleWeightedGraph<Vertex, DefaultWeightedEdge> graph,
                                                             int low, int high, boolean saveEdges) {
-        List<Integer[]> removedEdges = new ArrayList<>();
+        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) {
-                    Integer source = graph.getEdgeSource(e);
-                    Integer target = graph.getEdgeTarget(e);
+                    Vertex source = graph.getEdgeSource(e);
+                    Vertex target = graph.getEdgeTarget(e);
                     Integer weight = (int) graph.getEdgeWeight(e);
-                    Integer[] edge = {source, target, weight};
-                    removedEdges.add(edge);
+                    Vertex[] edge = {source, target};
+                    removedEdges.put(edge, weight);
                 }
                 else {
                     graph.setEdgeWeight(e, 0.0);
@@ -26,7 +28,7 @@ public interface GraphModificationFunctions {
             }
         }
         if(saveEdges) {
-            for (Integer[] edge : removedEdges) {
+            for (Vertex[] edge : removedEdges.keySet()) {
                 graph.removeEdge(edge[0], edge[1]);
             }
         }
@@ -34,23 +36,19 @@ public interface GraphModificationFunctions {
     }
 
     //Remove edges for pairs with large occupancy discrepancy
-    static List<Integer[]> filterByRelativeOccupancy(SimpleWeightedGraph<Integer, DefaultWeightedEdge> graph,
-                                                  Map<Integer, Integer> alphaWellCounts,
-                                                  Map<Integer, Integer> betaWellCounts,
-                                                  Map<Integer, Integer> plateVtoAMap,
-                                                  Map<Integer, Integer> plateVtoBMap,
+    static Map<Vertex[], Integer> filterByRelativeOccupancy(SimpleWeightedGraph<Vertex, DefaultWeightedEdge> graph,
                                                   Integer maxOccupancyDifference, boolean saveEdges) {
-        List<Integer[]> removedEdges = new ArrayList<>();
+        Map<Vertex[], Integer> removedEdges = new HashMap<>();
         for (DefaultWeightedEdge e : graph.edgeSet()) {
-            Integer alphaOcc = alphaWellCounts.get(plateVtoAMap.get(graph.getEdgeSource(e)));
-            Integer betaOcc = betaWellCounts.get(plateVtoBMap.get(graph.getEdgeTarget(e)));
+            Integer alphaOcc = graph.getEdgeSource(e).getOccupancy();
+            Integer betaOcc = graph.getEdgeTarget(e).getOccupancy();
             if (Math.abs(alphaOcc - betaOcc) >= maxOccupancyDifference) {
                 if (saveEdges) {
-                    Integer source = graph.getEdgeSource(e);
-                    Integer target = graph.getEdgeTarget(e);
+                    Vertex source = graph.getEdgeSource(e);
+                    Vertex target = graph.getEdgeTarget(e);
                     Integer weight = (int) graph.getEdgeWeight(e);
-                    Integer[] edge = {source, target, weight};
-                    removedEdges.add(edge);
+                    Vertex[] edge = {source, target};
+                    removedEdges.put(edge, weight);
                 }
                 else {
                     graph.setEdgeWeight(e, 0.0);
@@ -58,7 +56,7 @@ public interface GraphModificationFunctions {
             }
         }
         if(saveEdges) {
-            for (Integer[] edge : removedEdges) {
+            for (Vertex[] edge : removedEdges.keySet()) {
                 graph.removeEdge(edge[0], edge[1]);
             }
         }
@@ -66,26 +64,22 @@ public interface GraphModificationFunctions {
     }
 
     //Remove edges for pairs where overlap size is significantly lower than the well occupancy
-    static List<Integer[]> filterByOverlapPercent(SimpleWeightedGraph<Integer, DefaultWeightedEdge> graph,
-                                                         Map<Integer, Integer> alphaWellCounts,
-                                                         Map<Integer, Integer> betaWellCounts,
-                                                         Map<Integer, Integer> plateVtoAMap,
-                                                         Map<Integer, Integer> plateVtoBMap,
+    static Map<Vertex[], Integer> filterByOverlapPercent(SimpleWeightedGraph<Vertex, DefaultWeightedEdge> graph,
                                                          Integer minOverlapPercent,
                                                          boolean saveEdges) {
-        List<Integer[]> removedEdges = new ArrayList<>();
+        Map<Vertex[], Integer> removedEdges = new HashMap<>();
         for (DefaultWeightedEdge e : graph.edgeSet()) {
-            Integer alphaOcc = alphaWellCounts.get(plateVtoAMap.get(graph.getEdgeSource(e)));
-            Integer betaOcc = betaWellCounts.get(plateVtoBMap.get(graph.getEdgeTarget(e)));
+            Integer alphaOcc = graph.getEdgeSource(e).getOccupancy();
+            Integer betaOcc = graph.getEdgeTarget(e).getOccupancy();
             double weight = graph.getEdgeWeight(e);
             double min = minOverlapPercent / 100.0;
             if ((weight / alphaOcc < min) || (weight / betaOcc < min)) {
                 if (saveEdges) {
-                    Integer source = graph.getEdgeSource(e);
-                    Integer target = graph.getEdgeTarget(e);
+                    Vertex source = graph.getEdgeSource(e);
+                    Vertex target = graph.getEdgeTarget(e);
                     Integer intWeight = (int) graph.getEdgeWeight(e);
-                    Integer[] edge = {source, target, intWeight};
-                    removedEdges.add(edge);
+                    Vertex[] edge = {source, target};
+                    removedEdges.put(edge, intWeight);
                 }
                 else {
                     graph.setEdgeWeight(e, 0.0);
@@ -93,18 +87,18 @@ public interface GraphModificationFunctions {
             }
         }
         if(saveEdges) {
-            for (Integer[] edge : removedEdges) {
+            for (Vertex[] edge : removedEdges.keySet()) {
                 graph.removeEdge(edge[0], edge[1]);
             }
         }
         return removedEdges;
     }
 
-    static void addRemovedEdges(SimpleWeightedGraph<Integer, DefaultWeightedEdge> graph,
-                                       List<Integer[]> removedEdges) {
-        for (Integer[] edge : removedEdges) {
+    static void addRemovedEdges(SimpleWeightedGraph<Vertex, DefaultWeightedEdge> graph,
+                                Map<Vertex[], Integer> removedEdges) {
+        for (Vertex[] edge : removedEdges.keySet()) {
             DefaultWeightedEdge e = graph.addEdge(edge[0], edge[1]);
-            graph.setEdgeWeight(e, (double) edge[2]);
+            graph.setEdgeWeight(e, removedEdges.get(edge));
         }
     }
 

src/main/java/GraphWithMapData.java

@@ -6,6 +6,7 @@ import java.util.Map;
 //Can't just write the graph, because I need the occupancy data too.
 //Makes most sense to serialize object and write that to a file.
 //Which means there's no reason to split map data and graph data up.
+//Custom vertex class means a lot of the map data can now be encoded in the graph itself
 public class GraphWithMapData implements java.io.Serializable {
 
     private String sourceFilename;
@@ -15,32 +16,33 @@ public class GraphWithMapData implements java.io.Serializable {
     private Integer alphaCount;
     private Integer betaCount;
     private final Map<Integer, Integer> distCellsMapAlphaKey;
-    private final Map<Integer, Integer> plateVtoAMap;
-    private final Map<Integer, Integer> plateVtoBMap;
-    private final Map<Integer, Integer> plateAtoVMap;
-    private final Map<Integer, Integer> plateBtoVMap;
-    private final Map<Integer, Integer> alphaWellCounts;
-    private final Map<Integer, Integer> betaWellCounts;
+//    private final Map<Integer, Integer> plateVtoAMap;
+//    private final Map<Integer, Integer> plateVtoBMap;
+//    private final Map<Integer, Integer> plateAtoVMap;
+//    private final Map<Integer, Integer> plateBtoVMap;
+//    private final Map<Integer, Integer> alphaWellCounts;
+//    private final Map<Integer, Integer> betaWellCounts;
     private final Duration time;
 
     public GraphWithMapData(SimpleWeightedGraph graph, Integer numWells, Integer[] wellConcentrations,
-                            Integer alphaCount, Integer betaCount,
-                            Map<Integer, Integer> distCellsMapAlphaKey, Map<Integer, Integer> plateVtoAMap,
-                            Map<Integer,Integer> plateVtoBMap, Map<Integer, Integer> plateAtoVMap,
-                            Map<Integer, Integer> plateBtoVMap, Map<Integer, Integer> alphaWellCounts,
-                            Map<Integer, Integer> betaWellCounts, Duration time) {
+                            Map<Integer, Integer> distCellsMapAlphaKey, Duration time){
+
+//                            Map<Integer, Integer> plateVtoAMap, Integer alphaCount, Integer betaCount,
+//                            Map<Integer,Integer> plateVtoBMap, Map<Integer, Integer> plateAtoVMap,
+//                            Map<Integer, Integer> plateBtoVMap, Map<Integer, Integer> alphaWellCounts,
+//                            Map<Integer, Integer> betaWellCounts,) {
         this.graph = graph;
         this.numWells = numWells;
         this.wellPopulations = wellConcentrations;
         this.alphaCount = alphaCount;
         this.betaCount = betaCount;
         this.distCellsMapAlphaKey = distCellsMapAlphaKey;
-        this.plateVtoAMap = plateVtoAMap;
-        this.plateVtoBMap = plateVtoBMap;
-        this.plateAtoVMap = plateAtoVMap;
-        this.plateBtoVMap = plateBtoVMap;
-        this.alphaWellCounts = alphaWellCounts;
-        this.betaWellCounts = betaWellCounts;
+//        this.plateVtoAMap = plateVtoAMap;
+//        this.plateVtoBMap = plateVtoBMap;
+//        this.plateAtoVMap = plateAtoVMap;
+//        this.plateBtoVMap = plateBtoVMap;
+//        this.alphaWellCounts = alphaWellCounts;
+//        this.betaWellCounts = betaWellCounts;
         this.time = time;
     }
 
@@ -56,41 +58,41 @@ public class GraphWithMapData implements java.io.Serializable {
         return wellPopulations;
     }
 
-    public Integer getAlphaCount() {
-        return alphaCount;
-    }
-
-    public Integer getBetaCount() {
-        return betaCount;
-    }
+//    public Integer getAlphaCount() {
+//        return alphaCount;
+//    }
+//
+//    public Integer getBetaCount() {
+//        return betaCount;
+//    }
 
     public Map<Integer, Integer> getDistCellsMapAlphaKey() {
         return distCellsMapAlphaKey;
     }
 
-    public Map<Integer, Integer> getPlateVtoAMap() {
-        return plateVtoAMap;
-    }
-
-    public Map<Integer, Integer> getPlateVtoBMap() {
-        return plateVtoBMap;
-    }
-
-    public Map<Integer, Integer> getPlateAtoVMap() {
-        return plateAtoVMap;
-    }
-
-    public Map<Integer, Integer> getPlateBtoVMap() {
-        return plateBtoVMap;
-    }
-
-    public Map<Integer, Integer> getAlphaWellCounts() {
-        return alphaWellCounts;
-    }
-
-    public Map<Integer, Integer> getBetaWellCounts() {
-        return betaWellCounts;
-    }
+//    public Map<Integer, Integer> getPlateVtoAMap() {
+//        return plateVtoAMap;
+//    }
+//
+//    public Map<Integer, Integer> getPlateVtoBMap() {
+//        return plateVtoBMap;
+//    }
+//
+//    public Map<Integer, Integer> getPlateAtoVMap() {
+//        return plateAtoVMap;
+//    }
+//
+//    public Map<Integer, Integer> getPlateBtoVMap() {
+//        return plateBtoVMap;
+//    }
+//
+//    public Map<Integer, Integer> getAlphaWellCounts() {
+//        return alphaWellCounts;
+//    }
+//
+//    public Map<Integer, Integer> getBetaWellCounts() {
+//        return betaWellCounts;
+//    }
 
     public Duration getTime() {
         return time;

src/main/java/InteractiveInterface.java

@@ -258,7 +258,7 @@ public class InteractiveInterface {
             cellFile = sc.next();
             System.out.print("\nPlease enter name of an existing sample plate file: ");
             plateFile = sc.next();
-            System.out.println("\nThe graph and occupancy data will be written to a serialized binary file.");
+            System.out.println("\nThe graph and occupancy data will be written to a file.");
             System.out.print("Please enter a name for the output file: ");
             filename = sc.next();
         } catch (InputMismatchException ex) {
@@ -504,7 +504,7 @@ public class InteractiveInterface {
             System.out.println("2) Turn " + getOnOff(!BiGpairSEQ.cachePlate()) + " plate file caching");
             System.out.println("3) Turn " + getOnOff(!BiGpairSEQ.cacheGraph()) + " graph/data file caching");
             System.out.println("4) Turn " + getOnOff(!BiGpairSEQ.outputBinary()) + " serialized binary graph output");
-            System.out.println("5) Turn " + getOnOff(!BiGpairSEQ.outputGraphML()) + " GraphML graph output");
+            System.out.println("5) Turn " + getOnOff(!BiGpairSEQ.outputGraphML()) + " GraphML graph output (for data portability to other programs)");
             System.out.println("6) Maximum weight matching algorithm options");
             System.out.println("0) Return to main menu");
             try {

src/main/java/SequenceType.java

@@ -0,0 +1,8 @@
+//enum for tagging types of sequences
+//Listed in order that they appear in a cell array, so ordinal() method will return correct index
+public enum SequenceType {
+    CDR3_ALPHA,
+    CDR3_BETA,
+    CDR1_ALPHA,
+    CDR1_BETA
+}

src/main/java/Simulator.java

@@ -18,17 +18,17 @@ import static java.lang.Float.*;
 
 //NOTE: "sequence" in method and variable names refers to a peptide sequence from a simulated T cell
 public class Simulator implements GraphModificationFunctions {
-    private static final int cdr3AlphaIndex = 0;
-    private static final int cdr3BetaIndex = 1;
-    private static final int cdr1AlphaIndex = 2;
-    private static final int cdr1BetaIndex = 3;
 
-    //Make the graph needed for matching CDR3s
+
+    //Make the graph needed for matching sequences.
+    //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) {
         Instant start = Instant.now();
         List<Integer[]> distinctCells = cellSample.getCells();
-        int[] alphaIndex = {cdr3AlphaIndex};
-        int[] betaIndex = {cdr3BetaIndex};
+        int[] alphaIndices = {SequenceType.CDR3_ALPHA.ordinal()};
+        int[] betaIndices = {SequenceType.CDR3_BETA.ordinal()};
 
         int numWells = samplePlate.getSize();
 
@@ -38,18 +38,20 @@ public class Simulator implements GraphModificationFunctions {
         if(verbose){System.out.println("Cell maps made");}
 
         if(verbose){System.out.println("Making well maps");}
-        Map<Integer, Integer> allAlphas = samplePlate.assayWellsSequenceS(alphaIndex);
-        Map<Integer, Integer> allBetas = samplePlate.assayWellsSequenceS(betaIndex);
+
+        Map<Integer, Integer> allAlphas = samplePlate.assayWellsSequenceS(alphaIndices);
+        Map<Integer, Integer> allBetas = samplePlate.assayWellsSequenceS(betaIndices);
         int alphaCount = allAlphas.size();
         if(verbose){System.out.println("All alphas count: " + alphaCount);}
         int betaCount = allBetas.size();
         if(verbose){System.out.println("All betas count: " + betaCount);}
         if(verbose){System.out.println("Well maps made");}
 
-        if(verbose){System.out.println("Removing singleton sequences and sequences present in all wells.");}
-        filterByOccupancyThresholds(allAlphas, 2, numWells - 1);
-        filterByOccupancyThresholds(allBetas, 2, numWells - 1);
-        if(verbose){System.out.println("Sequences removed");}
+
+//        if(verbose){System.out.println("Removing singleton sequences and sequences present in all wells.");}
+//        filterByOccupancyThresholds(allAlphas, 2, numWells - 1);
+//        filterByOccupancyThresholds(allBetas, 2, numWells - 1);
+//        if(verbose){System.out.println("Sequences removed");}
         int pairableAlphaCount = allAlphas.size();
         if(verbose){System.out.println("Remaining alphas count: " + pairableAlphaCount);}
         int pairableBetaCount = allBetas.size();
@@ -78,29 +80,40 @@ public class Simulator implements GraphModificationFunctions {
         //(technically this is only 1/4 of an adjacency matrix, but that's all you need
         //for a bipartite graph, and all the SimpleWeightedBipartiteGraphMatrixGenerator class expects.)
         if(verbose){System.out.println("Creating adjacency matrix");}
-        //Count how many wells each alpha appears in
+        //Count how many wells each alpha sequence appears in
         Map<Integer, Integer> alphaWellCounts = new HashMap<>();
-        //count how many wells each beta appears in
+        //count how many wells each beta sequence appears in
         Map<Integer, Integer> betaWellCounts = new HashMap<>();
         //the adjacency matrix to be used by the graph generator
         double[][] weights = new double[plateVtoAMap.size()][plateVtoBMap.size()];
         countSequencesAndFillMatrix(samplePlate, allAlphas, allBetas, plateAtoVMap,
-                plateBtoVMap, alphaIndex, betaIndex, alphaWellCounts, betaWellCounts, weights);
+                plateBtoVMap, alphaIndices, betaIndices, alphaWellCounts, betaWellCounts, weights);
         if(verbose){System.out.println("Matrix created");}
 
         //create bipartite graph
         if(verbose){System.out.println("Creating graph");}
         //the graph object
-        SimpleWeightedGraph<Integer, DefaultWeightedEdge> graph =
+        SimpleWeightedGraph<Vertex, DefaultWeightedEdge> graph =
                 new SimpleWeightedGraph<>(DefaultWeightedEdge.class);
         //the graph generator
         SimpleWeightedBipartiteGraphMatrixGenerator graphGenerator = new SimpleWeightedBipartiteGraphMatrixGenerator();
         //the list of alpha vertices
-        List<Integer> alphaVertices = new ArrayList<>(plateVtoAMap.keySet()); //This will work because LinkedHashMap preserves order of entry
+        //List<Integer> alphaVertices = new ArrayList<>(plateVtoAMap.keySet()); //This will work because LinkedHashMap preserves order of entry
+        List<Vertex> alphaVertices = new ArrayList<>();
+        //start with map of all alphas mapped to vertex values, get occupancy from the alphaWellCounts map
+        for (Integer seq : plateAtoVMap.keySet()) {
+            Vertex alphaVertex = new Vertex(SequenceType.CDR3_ALPHA, seq, alphaWellCounts.get(seq), plateAtoVMap.get(seq));
+            alphaVertices.add(alphaVertex);
+        }
         graphGenerator.first(alphaVertices);
         //the list of beta vertices
-        List<Integer> betaVertices = new ArrayList<>(plateVtoBMap.keySet());
-        graphGenerator.second(betaVertices); //This will work because LinkedHashMap preserves order of entry
+        //List<Integer> betaVertices = new ArrayList<>(plateVtoBMap.keySet());//This will work because LinkedHashMap preserves order of entry
+        List<Vertex> betaVertices = new ArrayList<>();
+        for (Integer seq : plateBtoVMap.keySet()) {
+            Vertex betaVertex = new Vertex(SequenceType.CDR3_BETA, seq, betaWellCounts.get(seq), plateBtoVMap.get(seq));
+            betaVertices.add(betaVertex);
+        }
+        graphGenerator.second(betaVertices);
         //use adjacency matrix of weight created previously
         graphGenerator.weights(weights);
         graphGenerator.generateGraph(graph);
@@ -110,9 +123,7 @@ public class Simulator implements GraphModificationFunctions {
         Duration time = Duration.between(start, stop);
 
         //create GraphWithMapData object
-        GraphWithMapData output = new GraphWithMapData(graph, numWells, samplePlate.getPopulations(), alphaCount, betaCount,
-                distCellsMapAlphaKey, plateVtoAMap, plateVtoBMap, plateAtoVMap,
-                plateBtoVMap, alphaWellCounts, betaWellCounts, time);
+        GraphWithMapData output = new GraphWithMapData(graph, numWells, samplePlate.getPopulations(), distCellsMapAlphaKey, time);
         //Set source file name in graph to name of sample plate
         output.setSourceFilename(samplePlate.getFilename());
         //return GraphWithMapData object
@@ -124,39 +135,46 @@ public class Simulator implements GraphModificationFunctions {
                                             Integer highThreshold, Integer maxOccupancyDifference,
                                             Integer minOverlapPercent, boolean verbose) {
         Instant start = Instant.now();
-        List<Integer[]> removedEdges = new ArrayList<>();
+        SimpleWeightedGraph<Vertex, DefaultWeightedEdge> graph = data.getGraph();
+        Map<Vertex[], Integer> removedEdges = new HashMap<>();
         boolean saveEdges = BiGpairSEQ.cacheGraph();
         int numWells = data.getNumWells();
-        Integer alphaCount = data.getAlphaCount();
-        Integer betaCount = data.getBetaCount();
+        //Integer alphaCount = data.getAlphaCount();
+        //Integer betaCount = data.getBetaCount();
         Map<Integer, Integer> distCellsMapAlphaKey = data.getDistCellsMapAlphaKey();
-        Map<Integer, Integer> plateVtoAMap = data.getPlateVtoAMap();
-        Map<Integer, Integer> plateVtoBMap = data.getPlateVtoBMap();
-        Map<Integer, Integer> alphaWellCounts = data.getAlphaWellCounts();
-        Map<Integer, Integer> betaWellCounts = data.getBetaWellCounts();
-        SimpleWeightedGraph<Integer, DefaultWeightedEdge> graph = data.getGraph();
+        Set<Vertex> alphas = new HashSet<>();
+        Set<Vertex> betas = new HashSet<>();
+        for(Vertex v: graph.vertexSet()) {
+            if (SequenceType.CDR3_ALPHA.equals(v.getType())){
+                alphas.add(v);
+            }
+            else {
+                betas.add(v);
+            }
+        }
+        Integer alphaCount = alphas.size();
+        Integer betaCount = betas.size();
 
         //remove edges with weights outside given overlap thresholds, add those to removed edge list
         if(verbose){System.out.println("Eliminating edges with weights outside overlap threshold values");}
-        removedEdges.addAll(GraphModificationFunctions.filterByOverlapThresholds(graph, lowThreshold, highThreshold, saveEdges));
+        removedEdges.putAll(GraphModificationFunctions.filterByOverlapThresholds(graph, lowThreshold, highThreshold, saveEdges));
         if(verbose){System.out.println("Over- and under-weight edges removed");}
 
         //remove edges between vertices with too small an overlap size, add those to removed edge list
         if(verbose){System.out.println("Eliminating edges with weights less than " + minOverlapPercent.toString() +
                 " percent of vertex occupancy value.");}
-        removedEdges.addAll(GraphModificationFunctions.filterByOverlapPercent(graph, alphaWellCounts, betaWellCounts,
-                plateVtoAMap, plateVtoBMap, minOverlapPercent, saveEdges));
+        removedEdges.putAll(GraphModificationFunctions.filterByOverlapPercent(graph, minOverlapPercent, saveEdges));
         if(verbose){System.out.println("Edges with weights too far below a vertex occupancy value removed");}
 
         //Filter by relative occupancy
         if(verbose){System.out.println("Eliminating edges between vertices with occupancy difference > "
                 + maxOccupancyDifference);}
-        removedEdges.addAll(GraphModificationFunctions.filterByRelativeOccupancy(graph, alphaWellCounts, betaWellCounts,
-                plateVtoAMap, plateVtoBMap, maxOccupancyDifference, saveEdges));
+        removedEdges.putAll(GraphModificationFunctions.filterByRelativeOccupancy(graph, maxOccupancyDifference, saveEdges));
         if(verbose){System.out.println("Edges between vertices of with excessively different occupancy values " +
                 "removed");}
 
         //Find Maximum Weighted Matching
+        //using jheaps library class PairingHeap for improved efficiency
         if(verbose){System.out.println("Finding maximum weighted matching");}
         MaximumWeightBipartiteMatching maxWeightMatching;
         //Use correct heap type for priority queue
@@ -164,20 +182,20 @@ public class Simulator implements GraphModificationFunctions {
         switch (heapType) {
             case "PAIRING" -> {
                 maxWeightMatching = new MaximumWeightBipartiteMatching(graph,
-                        plateVtoAMap.keySet(),
-                        plateVtoBMap.keySet(),
+                        alphas,
+                        betas,
                         i -> new PairingHeap(Comparator.naturalOrder()));
             }
             case "FIBONACCI" -> {
                 maxWeightMatching = new MaximumWeightBipartiteMatching(graph,
-                        plateVtoAMap.keySet(),
-                        plateVtoBMap.keySet(),
+                        alphas,
+                        betas,
                         i -> new FibonacciHeap(Comparator.naturalOrder()));
             }
             default -> {
                 maxWeightMatching = new MaximumWeightBipartiteMatching(graph,
-                        plateVtoAMap.keySet(),
-                        plateVtoBMap.keySet());
+                        alphas,
+                        betas);
             }
         }
         //get the matching
@@ -207,11 +225,14 @@ public class Simulator implements GraphModificationFunctions {
         Map<Integer, Integer> matchMap = new HashMap<>();
         while(weightIter.hasNext()) {
             e = weightIter.next();
-            Integer source = graph.getEdgeSource(e);
-            Integer target = graph.getEdgeTarget(e);
+            Vertex source = graph.getEdgeSource(e);
+            Vertex target = graph.getEdgeTarget(e);
+            //Integer source = graph.getEdgeSource(e);
+            //Integer target = graph.getEdgeTarget(e);
             //The match map is all matches found, not just true matches!
-            matchMap.put(plateVtoAMap.get(source), plateVtoBMap.get(target));
-            check = plateVtoBMap.get(target).equals(distCellsMapAlphaKey.get(plateVtoAMap.get(source)));
+            matchMap.put(source.getSequence(), target.getSequence());
+            check = target.getSequence().equals(distCellsMapAlphaKey.get(source.getSequence()));
+            //check = plateVtoBMap.get(target).equals(distCellsMapAlphaKey.get(plateVtoAMap.get(source)));
             if(check) {
                 trueCount++;
             }
@@ -219,17 +240,19 @@ public class Simulator implements GraphModificationFunctions {
                 falseCount++;
             }
             List<String> result = new ArrayList<>();
-            result.add(plateVtoAMap.get(source).toString());
+            //alpha sequence
+            result.add(source.getSequence().toString());
             //alpha well count
-            result.add(alphaWellCounts.get(plateVtoAMap.get(source)).toString());
-            result.add(plateVtoBMap.get(target).toString());
+            result.add(source.getOccupancy().toString());
+            //beta sequence
+            result.add(target.getSequence().toString());
             //beta well count
-            result.add(betaWellCounts.get(plateVtoBMap.get(target)).toString());
+            result.add(target.getOccupancy().toString());
             //overlap count
             result.add(Double.toString(graph.getEdgeWeight(e)));
             result.add(Boolean.toString(check));
-            double pValue = Equations.pValue(numWells, alphaWellCounts.get(plateVtoAMap.get(source)),
-                    betaWellCounts.get(plateVtoBMap.get(target)), graph.getEdgeWeight(e));
+            double pValue = Equations.pValue(numWells, source.getOccupancy(),
+                    target.getOccupancy(), graph.getEdgeWeight(e));
             BigDecimal pValueTrunc = new BigDecimal(pValue, mc);
             result.add(pValueTrunc.toString());
             allResults.add(result);
@@ -238,6 +261,8 @@ public class Simulator implements GraphModificationFunctions {
         //Metadata comments for CSV file
         String algoType = "LEDA book with heap: " + heapType;
         int min = Math.min(alphaCount, betaCount);
+        //matching weight
+        BigDecimal totalMatchingWeight = maxWeightMatching.getMatchingWeight();
         //rate of attempted matching
         double attemptRate = (double) (trueCount + falseCount) / min;
         BigDecimal attemptRateTrunc = new BigDecimal(attemptRate, mc);
@@ -268,6 +293,7 @@ public class Simulator implements GraphModificationFunctions {
         metadata.put("sample plate filename", data.getSourceFilename());
         metadata.put("graph filename", dataFilename);
         metadata.put("algorithm type", algoType);
+        metadata.put("matching weight", totalMatchingWeight.toString());
         metadata.put("well populations", wellPopulationsString);
         metadata.put("total alphas found", alphaCount.toString());
         metadata.put("total betas found", betaCount.toString());

src/main/java/Vertex.java

@@ -1,23 +1,92 @@
+import java.io.Serializable;
 
+public class Vertex implements Serializable {
+    private  SequenceType type;
+    private  Integer vertexLabel;
+    private  Integer sequence;
+    private  Integer occupancy;
 
-public class Vertex {
-    private final Integer vertexLabel;
-    private final Integer sequence;
-    private final Integer occupancy;
+    public Vertex(Integer vertexLabel) {
+        this.vertexLabel = vertexLabel;
+    }
+    public Vertex(String vertexLabel) {
+        this.vertexLabel = Integer.parseInt((vertexLabel));
+    }
 
-    public Vertex(Integer vertexLabel, Integer sequence, Integer occupancy) {
+    public Vertex(SequenceType type, Integer sequence, Integer occupancy, Integer vertexLabel) {
+        this.type = type;
         this.vertexLabel = vertexLabel;
         this.sequence = sequence;
         this.occupancy = occupancy;
     }
 
-    public Integer getVertexLabel() { return vertexLabel; }
+
+    public SequenceType getType() {
+        return type;
+    }
+
+    public void setType(String type) {
+        this.type = SequenceType.valueOf(type);
+    }
+
+    public Integer getVertexLabel() {
+        return vertexLabel;
+    }
+
+    public void setVertexLabel(String label) {
+        this.vertexLabel = Integer.parseInt(label);
+    }
 
     public Integer getSequence() {
+
         return sequence;
     }
 
+    public void setSequence(String sequence) {
+        this.sequence = Integer.parseInt(sequence);
+    }
+
     public Integer getOccupancy() {
         return occupancy;
     }
+
+    public void setOccupancy(String occupancy) {
+        this.occupancy = Integer.parseInt(occupancy);
+    }
+
+    @Override //adapted from JGraphT example code
+    public int hashCode()
+    {
+        return (sequence == null) ? 0 : sequence.hashCode();
+    }
+
+    @Override //adapted from JGraphT example code
+    public boolean equals(Object obj)
+    {
+        if (this == obj)
+            return true;
+        if (obj == null)
+            return false;
+        if (getClass() != obj.getClass())
+            return false;
+        Vertex other = (Vertex) obj;
+        if (sequence == null) {
+            return other.sequence == null;
+        } else {
+            return sequence.equals(other.sequence);
+        }
+    }
+
+
+    @Override //adapted from JGraphT example code
+    public String toString()
+    {
+        StringBuilder sb = new StringBuilder();
+        sb.append("(").append(vertexLabel)
+                .append(", Type: ").append(type.name())
+                .append(", Sequence: ").append(sequence)
+                .append(", Occupancy: ").append(occupancy).append(")");
+        return sb.toString();
+    }
+
 }