Update Readme

readme.md

@@ -12,7 +12,7 @@ Unlike pairSEQ, which calculates p-values for every TCR alpha/beta overlap and c
 against a null distribution, BiGpairSEQ does not do any statistical calculations
 directly.
 
-BiGpairSEQ creates a [weightd bipartite graph](https://en.wikipedia.org/wiki/Bipartite_graph) representing the sample plate.
+BiGpairSEQ creates a [weighted bipartite graph](https://en.wikipedia.org/wiki/Bipartite_graph) representing the sample plate.
 The distinct TCRA and TCRB sequences form the two sets of vertices. Every TCRA/TCRB pair that share a well
 are connected by an edge, with the edge weight set to the number of wells in which both sequences appear.
 (Sequences present in *all* wells are filtered out prior to creating the graph, as there is no signal in their occupancy pattern.)
@@ -29,17 +29,13 @@ Unfortunately, it's a fairly new algorithm, and not yet implemented by the graph
 So this program instead uses the Fibonacci heap-based algorithm of Fredman and Tarjan (1987), which has a worst-case
 runtime of **O(n (n log(n) + m))**. The algorithm is implemented as described in Melhorn and Näher (1999).
 
-The current version of the program uses a pairing heap instead of a Fibonacci heap for its priority queue,
-which has lower theoretical efficiency but also lower complexity overhead, and is often equivalently performant
-in practice.
-
 ## USAGE
 
 ### RUNNING THE PROGRAM
 
 [Download the current version of BiGpairSEQ_Sim.](https://gitea.ejsf.synology.me/efischer/BiGpairSEQ/releases)
 
-BiGpairSEQ_Sim is an executable .jar file. Requires Java 11 or higher. [OpenJDK 17](https://jdk.java.net/17/)
+BiGpairSEQ_Sim is an executable .jar file. Requires Java 14 or higher. [OpenJDK 17](https://jdk.java.net/17/)
 recommended.
 
 Run with the command:
@@ -70,6 +66,18 @@ Please select an option:
 0) Exit
 ```
 
+By default, the Options menu looks like this:
+```
+--------------OPTIONS---------------
+1) Turn on cell sample file caching
+2) Turn on plate file caching
+3) Turn on graph/data file caching
+4) Turn off serialized binary graph output
+5) Turn on GraphML graph output
+6) Maximum weight matching algorithm options
+0) Return to main menu
+```
+
 ### INPUT/OUTPUT
 
 To run the simulation, the program reads and writes 4 kinds of files:
@@ -79,20 +87,25 @@ To run the simulation, the program reads and writes 4 kinds of files:
 * Matching Results files in CSV format
 
 These files are often generated in sequence. When entering filenames, it is not necessary to include the file extension
-(.csv or .ser). When reading or writing files, the program will automatically add the correct extension to any filename without one.
+(.csv or .ser). When reading or writing files, the program will automatically add the correct extension to any filename 
+without one.
 
 To save file I/O time, the most recent instance of each of these four
-files either generated or read from disk can be cached in program memory. This is especially important for Graph/Data files,
+files either generated or read from disk can be cached in program memory. When caching is active, subsequent uses of the 
-which can be several gigabytes in size. Since some simulations may require running multiple, 
+same data file won't need to be read in again until another file of that type is used or generated,
-differently-configured BiGpairSEQ matchings on the same graph, keeping the most recent graph cached can reduce execution time
-
-Subsequent uses of the same data file won't need to be read in again until another file of that type is used or generated,
 or caching is turned off for that file type. The program checks whether it needs to update its cached data by comparing
 filenames as entered by the user. On encountering a new filename, the program flushes its cache and reads in the new file.
 
-The program's caching behavior can be controlled in the Options menu. By default, caching for cell sample and 
+(Note that cached Graph/Data files must be transformed back into their original state after a matching experiment, which
-sample plate files is OFF, and caching for graph/data files is ON.
+may take some time. Whether file I/O or graph transformation takes longer for graph/data files is likely to be
+device-specific.)
 
+The program's caching behavior can be controlled in the Options menu. By default, all caching is OFF.
+
+The program can optionally output Graph/Data files in .GraphML format (.graphml) for data portability. This can be 
+turned on in the Options menu. By default, GraphML output is OFF.
+
+---
 #### Cell Sample Files
 Cell Sample files consist of any number of distinct "T cells." Every cell contains 
 four sequences: Alpha CDR3, Beta CDR3, Alpha CDR1, Beta CDR1. The sequences are represented by
@@ -110,7 +123,6 @@ Comments are preceded by `#`
 
 Structure:
 
----
     # Sample contains 1 unique CDR1 for every 4 unique CDR3s.
 | Alpha CDR3 | Beta CDR3 | Alpha CDR1 | Beta CDR1 |
 |---|---|---|---|
@@ -134,11 +146,14 @@ Options when making a Sample Plate file:
     * Standard deviation size 
   * Exponential
     * Lambda value
-      * *(Based on the slope of the graph in Figure 4C of the pairSEQ paper, the distribution of the original experiment was exponential with a lambda of approximately 0.6. (Howie, et al. 2015))*
+      * *(Based on the slope of the graph in Figure 4C of the pairSEQ paper, the distribution of the original experiment was approximately exponential with a lambda ~0.6. (Howie, et al. 2015))*
 * Total number of wells on the plate
-* Number of sections on plate
+* Well populations random or fixed
-* Number of T cells per well
+  * If random, minimum and maximum population sizes
-  * per section, if more than one section
+  * If fixed
+    * Number of sections on plate
+    * Number of T cells per well
+      * per section, if more than one section
 * Dropout rate
 
 Files are in CSV format. There are no header labels. Every row represents a well. 
@@ -152,7 +167,6 @@ Dropout sequences are replaced with the value `-1`. Comments are preceded by `#`
 
 Structure:
 
----
 ```
 # Cell source file name:
 # Each row represents one well on the plate
@@ -181,14 +195,19 @@ Options for creating a Graph/Data file:
 * The Cell Sample file to use
 * The Sample Plate file to use. (This must have been generated from the selected Cell Sample file.)
 
-These files do not have a human-readable structure, and are not portable to other programs. (Export of graphs in a 
+These files do not have a human-readable structure, and are not portable to other programs.
-portable data format may be implemented in the future. The tricky part is encoding the necessary metadata.)
+
+(For portability to other software, turn on GraphML output in the Options menu. This will produce a .graphml file
+for the weighted graph, with vertex attributes sequence, type, and occupancy data.)
 
 ---
 
 #### Matching Results Files 
-Matching results files consist of the results of a BiGpairSEQ matching simulation. Making them requires a Graph and 
+Matching results files consist of the results of a BiGpairSEQ matching simulation. Making them requires a serialized
-Data file. Matching results files are in CSV format. Rows are sequence pairings with extra relevant data. Columns are pairing-specific details.
+binary Graph/Data file (.ser). (Because .graphML files are larger than .ser files, BiGpairSEQ_Sim supports .graphML
+output only. Graph/data input must use a serialized binary.)
+
+Matching results files are in CSV format. Rows are sequence pairings with extra relevant data. Columns are pairing-specific details.
 Metadata about the matching simulation is included as comments. Comments are preceded by `#`.
 
 Options when running a BiGpairSEQ simulation of CDR3 alpha/beta matching:
@@ -203,7 +222,6 @@ Options when running a BiGpairSEQ simulation of CDR3 alpha/beta matching:
 
 Example output:
 
----
 ```
 # Source Sample Plate file: 4MilCellsPlate.csv
 # Source Graph and Data file: 4MilCellsPlateGraph.ser
@@ -254,27 +272,30 @@ slightly less time than the simulation itself. Real elapsed time from start to f
 ## TODO
 
 * ~~Try invoking GC at end of workloads to reduce paging to disk~~ DONE
-* Hold graph data in memory until another graph is read-in? ~~ABANDONED~~ ~~UNABANDONED~~ DONE
+* ~~Hold graph data in memory until another graph is read-in? ABANDONED UNABANDONED~~ DONE
   * ~~*No, this won't work, because BiGpairSEQ simulations alter the underlying graph based on filtering constraints. Changes would cascade with multiple experiments.*~~
-  * Might have figured out a way to do it, by taking edges out and then putting them back into the graph. This may actually be possible. If so, awesome.
+  * Might have figured out a way to do it, by taking edges out and then putting them back into the graph. This may actually be possible.
+  * It is possible, though the modifications to the graph incur their own performance penalties. Need testing to see which option is best.
+* ~~Test whether pairing heap (currently used) or Fibonacci heap is more efficient for priority queue in current matching algorithm~~ DONE
+  * ~~in theory Fibonacci heap should be more efficient, but complexity overhead may eliminate theoretical advantage~~
+  * ~~Add controllable heap-type parameter?~~
+    * Parameter implemented. Fibonacci heap the current default.
+* ~~Implement sample plates with random numbers of T cells per well.~~ DONE
+  * Possible BiGpairSEQ advantage over pairSEQ: BiGpairSEQ is resilient to variations in well population sizes on a sample plate; pairSEQ is not.
+    * preliminary data suggests that BiGpairSEQ behaves roughly as though the whole plate had whatever the *average* well concentration is, but that's still speculative.
 * See if there's a reasonable way to reformat Sample Plate files so that wells are columns instead of rows. 
   * ~~Problem is variable number of cells in a well~~
   * ~~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._
+    * _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
+    * 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
-* Implement sample plates with random numbers of T cells per well.
-  * Possible BiGpairSEQ advantage over pairSEQ: BiGpairSEQ is resilient to variations in well population sizes on a sample plate; pairSEQ is not.
-    * preliminary data suggests that BiGpairSEQ behaves roughly as though the whole plate had whatever the *average* well concentration is, but that's still speculative.
-* Enable GraphML output in addition to serialized object binaries, for data portability
-  * Custom vertex type with attribute for sequence occupancy?
 * Re-implement CDR1 matching method
 * Implement Duan and Su's maximum weight matching algorithm
-  * Add controllable algorithm-type parameter?
+    * Add controllable algorithm-type parameter?
-* Test whether pairing heap (currently used) or Fibonacci heap is more efficient for priority queue in current matching algorithm
+    * This would be fun and valuable, but probably take more time than I have for a hobby project.
-  * in theory Fibonacci heap should be more efficient, but complexity overhead may eliminate theoretical advantage
+  
-  * Add controllable heap-type parameter?
-
-
   
 ## CITATIONS
 * Howie, B., Sherwood, A. M., et al. ["High-throughput pairing of T cell receptor alpha and beta sequences."](https://pubmed.ncbi.nlm.nih.gov/26290413/) Sci. Transl. Med. 7, 301ra131 (2015)

src/main/java/BiGpairSEQ.java

@@ -1,6 +1,6 @@
 import java.util.Random;
 
-//main class. For choosing interface type and caching file data
+//main class. For choosing interface type and holding settings
 public class BiGpairSEQ {
 
     private static final Random rand = new Random();
@@ -12,7 +12,10 @@ public class BiGpairSEQ {
     private static String graphFilename = null;
     private static boolean cacheCells = false;
     private static boolean cachePlate = false;
-    private static boolean cacheGraph = true;
+    private static boolean cacheGraph = false;
+    private static String priorityQueueHeapType = "FIBONACCI";
+    private static boolean outputBinary = true;
+    private static boolean outputGraphML = false;
 
     public static void main(String[] args) {
         if (args.length == 0) {
@@ -151,4 +154,23 @@ public class BiGpairSEQ {
         }
         BiGpairSEQ.cacheGraph = cacheGraph;
     }
+
+    public static String getPriorityQueueHeapType() {
+        return priorityQueueHeapType;
+    }
+
+    public static void setPairingHeap() {
+        priorityQueueHeapType = "PAIRING";
+    }
+
+    public static void setFibonacciHeap() {
+        priorityQueueHeapType = "FIBONACCI";
+    }
+
+    public static boolean outputBinary() {return outputBinary;}
+    public static void setOutputBinary(boolean b) {outputBinary = b;}
+
+    public static boolean outputGraphML() {return outputGraphML;}
+    public static void setOutputGraphML(boolean b) {outputGraphML = b;}
+
 }

src/main/java/CellSample.java

@@ -1,10 +1,37 @@
+import java.util.ArrayList;
+import java.util.Collections;
 import java.util.List;
+import java.util.stream.IntStream;
 
 public class CellSample {
 
     private List<Integer[]> cells;
     private Integer cdr1Freq;
 
+    public CellSample(Integer numDistinctCells, Integer cdr1Freq){
+        this.cdr1Freq = cdr1Freq;
+        List<Integer> numbersCDR3 = new ArrayList<>();
+        List<Integer> numbersCDR1 = new ArrayList<>();
+        Integer numDistCDR3s = 2 * numDistinctCells + 1;
+        IntStream.range(1, numDistCDR3s + 1).forEach(i -> numbersCDR3.add(i));
+        IntStream.range(numDistCDR3s + 1, numDistCDR3s + 1 + (numDistCDR3s / cdr1Freq) + 1).forEach(i -> numbersCDR1.add(i));
+        Collections.shuffle(numbersCDR3);
+        Collections.shuffle(numbersCDR1);
+
+        //Each cell represented by 4 values
+        //two CDR3s, and two CDR1s. First two values are CDR3s (alpha, beta), second two are CDR1s (alpha, beta)
+        List<Integer[]> distinctCells = new ArrayList<>();
+        for(int i = 0; i < numbersCDR3.size() - 1; i = i + 2){
+            Integer tmpCDR3a = numbersCDR3.get(i);
+            Integer tmpCDR3b = numbersCDR3.get(i+1);
+            Integer tmpCDR1a = numbersCDR1.get(i % numbersCDR1.size());
+            Integer tmpCDR1b = numbersCDR1.get((i+1) % numbersCDR1.size());
+            Integer[] tmp = {tmpCDR3a, tmpCDR3b, tmpCDR1a, tmpCDR1b};
+            distinctCells.add(tmp);
+        }
+        this.cells = distinctCells;
+    }
+
     public CellSample(List<Integer[]> cells, Integer cdr1Freq){
         this.cells = cells;
         this.cdr1Freq = cdr1Freq;

src/main/java/CommandLineInterface.java

@@ -288,7 +288,7 @@ public class CommandLineInterface {
 
     //for calling from command line
     public static void makeCells(String filename, Integer numCells, Integer cdr1Freq){
-        CellSample sample = Simulator.generateCellSample(numCells, cdr1Freq);
+        CellSample sample = new CellSample(numCells, cdr1Freq);
         CellFileWriter writer = new CellFileWriter(filename, sample);
         writer.writeCellsToFile();
     }

src/main/java/GraphMLFileReader.java

@@ -1,35 +0,0 @@
-import org.jgrapht.graph.SimpleWeightedGraph;
-import org.jgrapht.nio.graphml.GraphMLImporter;
-
-import java.io.BufferedReader;
-import java.io.IOException;
-import java.nio.file.Files;
-import java.nio.file.Path;
-
-public class GraphMLFileReader {
-
-    private String filename;
-    private SimpleWeightedGraph graph;
-
-    public GraphMLFileReader(String filename, SimpleWeightedGraph graph) {
-        if(!filename.matches(".*\\.graphml")){
-            filename = filename + ".graphml";
-        }
-        this.filename = filename;
-        this.graph = graph;
-
-        try(//don't need to close reader bc of try-with-resources auto-closing
-            BufferedReader reader = Files.newBufferedReader(Path.of(filename));
-        ){
-            GraphMLImporter<SimpleWeightedGraph, BufferedReader> importer = new GraphMLImporter<>();
-            importer.importGraph(graph, reader);
-        }
-        catch (IOException ex) {
-            System.out.println("Graph file " + filename + " not found.");
-            System.err.println(ex);
-        }
-    }
-
-    public SimpleWeightedGraph getGraph() { return graph; }
-
-}

src/main/java/GraphMLFileWriter.java

@@ -1,4 +1,8 @@
+import org.jgrapht.graph.DefaultWeightedEdge;
 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;
 
@@ -7,25 +11,69 @@ import java.io.IOException;
 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;
 
 
-    public GraphMLFileWriter(String filename, SimpleWeightedGraph graph) {
+    public GraphMLFileWriter(String filename, GraphWithMapData data) {
         if(!filename.matches(".*\\.graphml")){
             filename = filename + ".graphml";
         }
         this.filename = filename;
-        this.graph = graph;
+        this.data = data;
     }
 
+//    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 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);
         ){
-            GraphMLExporter<SimpleWeightedGraph, BufferedWriter> exporter = new GraphMLExporter<>();
+            //create exporter. Let the vertex labels be the unique ids for the vertices
+            GraphMLExporter<Integer, SimpleWeightedGraph<Vertex, DefaultWeightedEdge>> exporter = new GraphMLExporter<>(v -> v.toString());
+            //set to export weights
+            exporter.setExportEdgeWeights(true);
+            //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))));
+                }
+                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);
+            //export the graph
             exporter.exportGraph(graph, writer);
         } catch(IOException ex){
             System.out.println("Could not make new file named "+filename);
@@ -33,3 +81,4 @@ public class GraphMLFileWriter {
         }
     }
 }
+

src/main/java/GraphModificationFunctions.java

@@ -4,61 +4,75 @@ import org.jgrapht.graph.SimpleWeightedGraph;
 import java.util.ArrayList;
 import java.util.List;
 import java.util.Map;
-import java.util.Set;
 
-public abstract class GraphModificationFunctions {
+public interface GraphModificationFunctions {
 
     //remove over- and under-weight edges
-    public static List<Integer[]> filterByOverlapThresholds(SimpleWeightedGraph<Integer, DefaultWeightedEdge> graph,
+    static List<Integer[]> filterByOverlapThresholds(SimpleWeightedGraph<Integer, DefaultWeightedEdge> graph,
-                                                            int low, int high) {
+                                                            int low, int high, boolean saveEdges) {
         List<Integer[]> removedEdges = new ArrayList<>();
-        for(DefaultWeightedEdge e: graph.edgeSet()){
+        for (DefaultWeightedEdge e : graph.edgeSet()) {
-            if ((graph.getEdgeWeight(e) > high) || (graph.getEdgeWeight(e) < low)){
+            if ((graph.getEdgeWeight(e) > high) || (graph.getEdgeWeight(e) < low)) {
-                Integer source = graph.getEdgeSource(e);
+                if(saveEdges) {
-                Integer target = graph.getEdgeTarget(e);
+                    Integer source = graph.getEdgeSource(e);
-                Integer weight = (int) graph.getEdgeWeight(e);
+                    Integer target = graph.getEdgeTarget(e);
-                Integer[] edge = {source, target, weight};
+                    Integer weight = (int) graph.getEdgeWeight(e);
-                removedEdges.add(edge);
+                    Integer[] edge = {source, target, weight};
+                    removedEdges.add(edge);
+                }
+                else {
+                    graph.setEdgeWeight(e, 0.0);
+                }
             }
         }
-        for (Integer[] edge : removedEdges) {
+        if(saveEdges) {
-            graph.removeEdge(edge[0], edge[1]);
+            for (Integer[] edge : removedEdges) {
+                graph.removeEdge(edge[0], edge[1]);
+            }
         }
         return removedEdges;
     }
 
     //Remove edges for pairs with large occupancy discrepancy
-    public static List<Integer[]> filterByRelativeOccupancy(SimpleWeightedGraph<Integer, DefaultWeightedEdge> graph,
+    static List<Integer[]> filterByRelativeOccupancy(SimpleWeightedGraph<Integer, DefaultWeightedEdge> graph,
                                                   Map<Integer, Integer> alphaWellCounts,
                                                   Map<Integer, Integer> betaWellCounts,
                                                   Map<Integer, Integer> plateVtoAMap,
                                                   Map<Integer, Integer> plateVtoBMap,
-                                                  Integer maxOccupancyDifference) {
+                                                  Integer maxOccupancyDifference, boolean saveEdges) {
         List<Integer[]> removedEdges = new ArrayList<>();
         for (DefaultWeightedEdge e : graph.edgeSet()) {
             Integer alphaOcc = alphaWellCounts.get(plateVtoAMap.get(graph.getEdgeSource(e)));
             Integer betaOcc = betaWellCounts.get(plateVtoBMap.get(graph.getEdgeTarget(e)));
             if (Math.abs(alphaOcc - betaOcc) >= maxOccupancyDifference) {
-                Integer source = graph.getEdgeSource(e);
+                if (saveEdges) {
-                Integer target = graph.getEdgeTarget(e);
+                    Integer source = graph.getEdgeSource(e);
-                Integer weight = (int) graph.getEdgeWeight(e);
+                    Integer target = graph.getEdgeTarget(e);
-                Integer[] edge = {source, target, weight};
+                    Integer weight = (int) graph.getEdgeWeight(e);
-                removedEdges.add(edge);
+                    Integer[] edge = {source, target, weight};
+                    removedEdges.add(edge);
+                }
+                else {
+                    graph.setEdgeWeight(e, 0.0);
+                }
             }
         }
-        for (Integer[] edge : removedEdges) {
+        if(saveEdges) {
-            graph.removeEdge(edge[0], edge[1]);
+            for (Integer[] edge : removedEdges) {
+                graph.removeEdge(edge[0], edge[1]);
+            }
         }
         return removedEdges;
     }
 
     //Remove edges for pairs where overlap size is significantly lower than the well occupancy
-    public static List<Integer[]> filterByOverlapPercent(SimpleWeightedGraph<Integer, DefaultWeightedEdge> graph,
+    static List<Integer[]> filterByOverlapPercent(SimpleWeightedGraph<Integer, DefaultWeightedEdge> graph,
                                                          Map<Integer, Integer> alphaWellCounts,
                                                          Map<Integer, Integer> betaWellCounts,
                                                          Map<Integer, Integer> plateVtoAMap,
                                                          Map<Integer, Integer> plateVtoBMap,
-                                                         Integer minOverlapPercent) {
+                                                         Integer minOverlapPercent,
+                                                         boolean saveEdges) {
         List<Integer[]> removedEdges = new ArrayList<>();
         for (DefaultWeightedEdge e : graph.edgeSet()) {
             Integer alphaOcc = alphaWellCounts.get(plateVtoAMap.get(graph.getEdgeSource(e)));
@@ -66,20 +80,27 @@ public abstract class GraphModificationFunctions {
             double weight = graph.getEdgeWeight(e);
             double min = minOverlapPercent / 100.0;
             if ((weight / alphaOcc < min) || (weight / betaOcc < min)) {
-                Integer source = graph.getEdgeSource(e);
+                if(saveEdges) {
-                Integer target = graph.getEdgeTarget(e);
+                    Integer source = graph.getEdgeSource(e);
-                Integer intWeight = (int) graph.getEdgeWeight(e);
+                    Integer target = graph.getEdgeTarget(e);
-                Integer[] edge = {source, target, intWeight};
+                    Integer intWeight = (int) graph.getEdgeWeight(e);
-                removedEdges.add(edge);
+                    Integer[] edge = {source, target, intWeight};
+                    removedEdges.add(edge);
+                }
+                else {
+                    graph.setEdgeWeight(e, 0.0);
+                }
             }
         }
-        for (Integer[] edge : removedEdges) {
+        if(saveEdges) {
-            graph.removeEdge(edge[0], edge[1]);
+            for (Integer[] edge : removedEdges) {
+                graph.removeEdge(edge[0], edge[1]);
+            }
         }
         return removedEdges;
     }
 
-    public static void addRemovedEdges(SimpleWeightedGraph<Integer, DefaultWeightedEdge> graph,
+    static void addRemovedEdges(SimpleWeightedGraph<Integer, DefaultWeightedEdge> graph,
                                        List<Integer[]> removedEdges) {
         for (Integer[] edge : removedEdges) {
             DefaultWeightedEdge e = graph.addEdge(edge[0], edge[1]);

src/main/java/InteractiveInterface.java

@@ -38,10 +38,10 @@ public class InteractiveInterface {
                         case 3 -> makeCDR3Graph();
                         case 4 -> matchCDR3s();
                         //case 6 -> matchCellsCDR1();
-                        case 8 -> options();
+                        case 8 -> mainOptions();
                         case 9 -> acknowledge();
                         case 0 -> quit = true;
-                        default -> throw new InputMismatchException("Invalid input.");
+                        default -> System.out.println("Invalid input.");
                     }
                 } catch (InputMismatchException | IOException ex) {
                     System.out.println(ex);
@@ -74,7 +74,7 @@ public class InteractiveInterface {
             System.out.println(ex);
             sc.next();
         }
-        CellSample sample = Simulator.generateCellSample(numCells, cdr1Freq);
+        CellSample sample = new CellSample(numCells, cdr1Freq);
         assert filename != null;
         System.out.println("Writing cells to file");
         CellFileWriter writer = new CellFileWriter(filename, sample);
@@ -252,7 +252,6 @@ public class InteractiveInterface {
         String filename = null;
         String cellFile = null;
         String plateFile = null;
-
         try {
             String str = "\nGenerating bipartite weighted graph encoding occupancy overlap data ";
             str = str.concat("\nrequires a cell sample file and a sample plate file.");
@@ -310,9 +309,16 @@ public class InteractiveInterface {
             List<Integer[]> cells = cellSample.getCells();
             GraphWithMapData data = Simulator.makeGraph(cells, plate, true);
             assert filename != null;
-            GraphDataObjectWriter dataWriter = new GraphDataObjectWriter(filename, data);
+            if(BiGpairSEQ.outputBinary()) {
-            dataWriter.writeDataToFile();
+                GraphDataObjectWriter dataWriter = new GraphDataObjectWriter(filename, data);
-            System.out.println("Graph and Data file written to: " + filename);
+                dataWriter.writeDataToFile();
+                System.out.println("Serialized binary graph/data file written to: " + filename);
+            }
+            if(BiGpairSEQ.outputGraphML()) {
+                GraphMLFileWriter graphMLWriter = new GraphMLFileWriter(filename, data);
+                graphMLWriter.writeGraphToFile();
+                System.out.println("GraphML file written to: " + filename);
+            }
             if(BiGpairSEQ.cacheGraph()) {
                 BiGpairSEQ.setGraphInMemory(data, filename);
 
@@ -493,13 +499,16 @@ public class InteractiveInterface {
 //        }
 //    }
 
-    private static void options(){
+    private static void mainOptions(){
         boolean backToMain = false;
         while(!backToMain) {
             System.out.println("\n--------------OPTIONS---------------");
             System.out.println("1) Turn " + getOnOff(!BiGpairSEQ.cacheCells()) + " cell sample file caching");
             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("6) Maximum weight matching algorithm options");
             System.out.println("0) Return to main menu");
             try {
                 input = sc.nextInt();
@@ -507,8 +516,11 @@ public class InteractiveInterface {
                     case 1 -> BiGpairSEQ.setCacheCells(!BiGpairSEQ.cacheCells());
                     case 2 -> BiGpairSEQ.setCachePlate(!BiGpairSEQ.cachePlate());
                     case 3 -> BiGpairSEQ.setCacheGraph(!BiGpairSEQ.cacheGraph());
+                    case 4 -> BiGpairSEQ.setOutputBinary(!BiGpairSEQ.outputBinary());
+                    case 5 -> BiGpairSEQ.setOutputGraphML(!BiGpairSEQ.outputGraphML());
+                    case 6 -> algorithmOptions();
                     case 0 -> backToMain = true;
-                    default -> throw new InputMismatchException("Invalid input.");
+                    default -> System.out.println("Invalid input");
                 }
             } catch (InputMismatchException ex) {
                 System.out.println(ex);
@@ -517,11 +529,49 @@ public class InteractiveInterface {
         }
     }
 
+    /**
+     * Helper function for printing menu items in mainOptions(). Returns a string based on the value of parameter.
+     *
+     * @param b - a boolean value
+     * @return String "on" if b is true, "off" if b is false
+     */
     private static String getOnOff(boolean b) {
         if (b) { return "on";}
         else { return "off"; }
     }
 
+    private static void algorithmOptions(){
+        boolean backToOptions = false;
+        while(!backToOptions) {
+            System.out.println("\n---------ALGORITHM OPTIONS----------");
+            System.out.println("1) Use scaling algorithm by Duan and Su.");
+            System.out.println("2) Use LEDA book algorithm with Fibonacci heap priority queue");
+            System.out.println("3) Use LEDA book algorithm with pairing heap priority queue");
+            System.out.println("0) Return to Options menu");
+            try {
+                input = sc.nextInt();
+                switch (input) {
+                    case 1 -> System.out.println("This option is not yet implemented. Choose another.");
+                    case 2 -> {
+                        BiGpairSEQ.setFibonacciHeap();
+                        System.out.println("MWM algorithm set to LEDA with Fibonacci heap");
+                        backToOptions = true;
+                    }
+                    case 3 -> {
+                        BiGpairSEQ.setPairingHeap();
+                        System.out.println("MWM algorithm set to LEDA with pairing heap");
+                        backToOptions = true;
+                    }
+                    case 0 -> backToOptions = true;
+                    default -> System.out.println("Invalid input");
+                }
+            } catch (InputMismatchException ex) {
+                System.out.println(ex);
+                sc.next();
+            }
+        }
+    }
+
     private static void acknowledge(){
         System.out.println("This program simulates BiGpairSEQ, a graph theory based adaptation");
         System.out.println("of the pairSEQ algorithm for pairing T cell receptor sequences.");

src/main/java/Simulator.java

@@ -3,6 +3,7 @@ import org.jgrapht.alg.matching.MaximumWeightBipartiteMatching;
 import org.jgrapht.generate.SimpleWeightedBipartiteGraphMatrixGenerator;
 import org.jgrapht.graph.DefaultWeightedEdge;
 import org.jgrapht.graph.SimpleWeightedGraph;
+import org.jheaps.tree.FibonacciHeap;
 import org.jheaps.tree.PairingHeap;
 
 import java.math.BigDecimal;
@@ -16,36 +17,12 @@ import java.util.stream.IntStream;
 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 {
+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;
 
-    public static CellSample generateCellSample(Integer numDistinctCells, Integer cdr1Freq) {
-        //In real T cells, CDR1s have about one third the diversity of CDR3s
-        List<Integer> numbersCDR3 = new ArrayList<>();
-        List<Integer> numbersCDR1 = new ArrayList<>();
-        Integer numDistCDR3s = 2 * numDistinctCells + 1;
-        IntStream.range(1, numDistCDR3s + 1).forEach(i -> numbersCDR3.add(i));
-        IntStream.range(numDistCDR3s + 1, numDistCDR3s + 1 + (numDistCDR3s / cdr1Freq) + 1).forEach(i -> numbersCDR1.add(i));
-        Collections.shuffle(numbersCDR3);
-        Collections.shuffle(numbersCDR1);
-
-        //Each cell represented by 4 values
-        //two CDR3s, and two CDR1s. First two values are CDR3s (alpha, beta), second two are CDR1s (alpha, beta)
-        List<Integer[]> distinctCells = new ArrayList<>();
-        for(int i = 0; i < numbersCDR3.size() - 1; i = i + 2){
-            Integer tmpCDR3a = numbersCDR3.get(i);
-            Integer tmpCDR3b = numbersCDR3.get(i+1);
-            Integer tmpCDR1a = numbersCDR1.get(i % numbersCDR1.size());
-            Integer tmpCDR1b = numbersCDR1.get((i+1) % numbersCDR1.size());
-            Integer[] tmp = {tmpCDR3a, tmpCDR3b, tmpCDR1a, tmpCDR1b};
-            distinctCells.add(tmp);
-        }
-        return new CellSample(distinctCells, cdr1Freq);
-    }
-
     //Make the graph needed for matching CDR3s
     public static GraphWithMapData makeGraph(List<Integer[]> distinctCells, Plate samplePlate, boolean verbose) {
         Instant start = Instant.now();
@@ -146,8 +123,8 @@ public class Simulator {
                                             Integer highThreshold, Integer maxOccupancyDifference,
                                             Integer minOverlapPercent, boolean verbose) {
         Instant start = Instant.now();
-        //Integer arrays will contain TO VERTEX, FROM VERTEX, and WEIGHT (which I'll need to cast to double)
         List<Integer[]> removedEdges = new ArrayList<>();
+        boolean saveEdges = BiGpairSEQ.cacheGraph();
         int numWells = data.getNumWells();
         Integer alphaCount = data.getAlphaCount();
         Integer betaCount = data.getBetaCount();
@@ -160,33 +137,50 @@ public class Simulator {
 
         //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));
+        removedEdges.addAll(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));
+                plateVtoAMap, plateVtoBMap, 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));
+                plateVtoAMap, plateVtoBMap, 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");}
-        //Attempting to use addressable heap to improve performance
+        MaximumWeightBipartiteMatching maxWeightMatching;
-        MaximumWeightBipartiteMatching maxWeightMatching =
+        //Use correct heap type for priority queue
-                new MaximumWeightBipartiteMatching(graph,
+        String heapType = BiGpairSEQ.getPriorityQueueHeapType();
+        switch (heapType) {
+            case "PAIRING" -> {
+                maxWeightMatching = new MaximumWeightBipartiteMatching(graph,
                         plateVtoAMap.keySet(),
                         plateVtoBMap.keySet(),
                         i -> new PairingHeap(Comparator.naturalOrder()));
+            }
+            case "FIBONACCI" -> {
+                maxWeightMatching = new MaximumWeightBipartiteMatching(graph,
+                        plateVtoAMap.keySet(),
+                        plateVtoBMap.keySet(),
+                        i -> new FibonacciHeap(Comparator.naturalOrder()));
+            }
+            default -> {
+                maxWeightMatching = new MaximumWeightBipartiteMatching(graph,
+                        plateVtoAMap.keySet(),
+                        plateVtoBMap.keySet());
+            }
+        }
+        //get the matching
         MatchingAlgorithm.Matching<String, DefaultWeightedEdge> graphMatching = maxWeightMatching.getMatching();
         if(verbose){System.out.println("Matching completed");}
         Instant stop = Instant.now();
@@ -242,6 +236,7 @@ public class Simulator {
         }
 
         //Metadata comments for CSV file
+        String algoType = "LEDA book with heap: " + heapType;
         int min = Math.min(alphaCount, betaCount);
         //rate of attempted matching
         double attemptRate = (double) (trueCount + falseCount) / min;
@@ -272,6 +267,7 @@ public class Simulator {
         Map<String, String> metadata = new LinkedHashMap<>();
         metadata.put("sample plate filename", data.getSourceFilename());
         metadata.put("graph filename", dataFilename);
+        metadata.put("algorithm type", algoType);
         metadata.put("well populations", wellPopulationsString);
         metadata.put("total alphas found", alphaCount.toString());
         metadata.put("total betas found", betaCount.toString());
@@ -292,10 +288,11 @@ public class Simulator {
             }
         }
 
-        //put the removed edges back on the graph
+        if(saveEdges) {
-        System.out.println("Restoring removed edges to graph.");
+            //put the removed edges back on the graph
-        GraphModificationFunctions.addRemovedEdges(graph, removedEdges);
+            System.out.println("Restoring removed edges to graph.");
-
+            GraphModificationFunctions.addRemovedEdges(graph, removedEdges);
+        }
         //return MatchingResult object
         return output;
     }
@@ -671,7 +668,7 @@ public class Simulator {
 
     private static Map<Integer, Integer> makeVertexToSequenceMap(Map<Integer, Integer> sequences, Integer startValue) {
         Map<Integer, Integer> map = new LinkedHashMap<>(); //LinkedHashMap to preserve order of entry
-        Integer index = startValue;
+        Integer index = startValue; //is this necessary? I don't think I use this.
         for (Integer k: sequences.keySet()) {
             map.put(index, k);
             index++;

src/main/java/Vertex.java

@@ -1,14 +1,20 @@
+
+
 public class Vertex {
-    private final Integer peptide;
+    private final Integer vertexLabel;
+    private final Integer sequence;
     private final Integer occupancy;
 
-    public Vertex(Integer peptide, Integer occupancy) {
+    public Vertex(Integer vertexLabel, Integer sequence, Integer occupancy) {
-        this.peptide = peptide;
+        this.vertexLabel = vertexLabel;
+        this.sequence = sequence;
         this.occupancy = occupancy;
     }
 
-    public Integer getPeptide() {
+    public Integer getVertexLabel() { return vertexLabel; }
-        return peptide;
+
+    public Integer getSequence() {
+        return sequence;
     }
 
     public Integer getOccupancy() {