Add read depth attributes to graphml output

# Conflicts:
#	src/main/java/GraphModificationFunctions.java
#	src/main/java/GraphWithMapData.java
#	src/main/java/Simulator.java
#	src/main/java/Vertex.java

readme.md

@@ -20,8 +20,8 @@ The problem of pairing TCRA/TCRB sequences thus reduces to the "assignment probl
 matching on a bipartite graph--the subset of vertex-disjoint edges whose weights sum to the maximum possible value.
 
 This is a well-studied combinatorial optimization problem, with many known solutions.
-The most efficient algorithm known to the author for maximum weight matching of a bipartite graph with strictly integral weights
-is from Duan and Su (2012). For a graph with m edges, n vertices per side, and maximum integer edge weight N, 
+The most efficient algorithm known to the author for maximum weight matching of a bipartite graph with strictly integral
+weights is from Duan and Su (2012). For a graph with m edges, n vertices per side, and maximum integer edge weight N, 
 their algorithm runs in **O(m sqrt(n) log(N))** time. As the graph representation of a pairSEQ experiment is 
 bipartite with integer weights, this algorithm is ideal for BiGpairSEQ.
 
@@ -43,13 +43,13 @@ Run with the command:
 `java -jar BiGpairSEQ_Sim.jar`
 
 Processing sample plates with tens of thousands of sequences may require large amounts 
-of RAM. It is often desirable to increase the JVM maximum heap allocation with the -Xmx flag.
+of RAM. It is often desirable to increase the JVM maximum heap allocation with the `-Xmx` flag.
 For example, to run the program with 32 gigabytes of memory, use the command:
 
 `java -Xmx32G -jar BiGpairSEQ_Sim.jar`
 
 There are a number of command line options, to allow the program to be used in shell scripts. For a full list,
-use the -help flag:
+use the `-help` flag:
 
 `java -jar BiGpairSEQ_Sim.jar -help`
 
@@ -108,7 +108,7 @@ 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 
+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.
 
 ---
@@ -200,11 +200,20 @@ then use it for multiple different BiGpairSEQ simulations.
 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.)
+* Whether to simulate sequence read depth. If simulated:
+  * The read depth (number of times each sequence is read)
+  * The read error rate (probability a sequence is misread)
+  * The error collision rate (probability two misreads produce the same spurious sequence)
 
 These files do not have a human-readable structure, and are not portable to other programs.
 
-(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.)
+*Optional GraphML output*
+
+For portability of graph data to other software, turn on [GraphML](http://graphml.graphdrawing.org/index.html) output 
+in the Options menu in interactive mode, or use the `-graphml`command line argument. This will produce a .graphml file 
+for the weighted graph, with vertex attributes for sequence, type, and occupancy data. This graph contains all the data 
+necessary for the BiGpairSEQ matching algorithm. It does not include the data to measure pairing accuracy; for that, 
+compare the matching results to the original Cell Sample .csv file.
 
 ---
 
@@ -259,49 +268,116 @@ Example output:
 P-values are calculated *after* BiGpairSEQ matching is completed, for purposes of comparison only, 
 using the (2021 corrected) formula from the original pairSEQ paper. (Howie, et al. 2015)
 
-### PERFORMANCE
-Performance details of the example excerpted above:
+
+## PERFORMANCE (old results; need updating to reflect current, improved simulator performance)
 
 On a home computer with a Ryzen 5600X CPU, 64GB of 3200MHz DDR4 RAM (half of which was allocated to the Java Virtual Machine), and a PCIe 3.0 SSD, running Linux Mint 20.3 Edge (5.13 kernel), 
 the author ran a BiGpairSEQ simulation of a 96-well sample plate with 30,000 T cells/well comprising ~11,800 alphas and betas,
-taken from a sample of 4,000,000 distinct cells with an exponential frequency distribution.
+taken from a sample of 4,000,000 distinct cells with an exponential frequency distribution (lambda 0.6).
 
 With min/max occupancy threshold of 3 and 94 wells for matching, and no other pre-filtering, BiGpairSEQ identified 5,151 
 correct pairings and 18 incorrect pairings, for an accuracy of 99.652%.
 
-The simulation time was 14'22". If intermediate results were held in memory, this would be equivalent to the total elapsed time.
+The total simulation time was 14'22". If intermediate results were held in memory, this would be equivalent to the total elapsed time.
 
 Since this implementation of BiGpairSEQ writes intermediate results to disk (to improve the efficiency of *repeated* simulations
 with different filtering options), the actual elapsed time was greater. File I/O time was not measured, but took 
 slightly less time than the simulation itself. Real elapsed time from start to finish was under 30 minutes.
 
+As mentioned in the theory section, performance could be improved by implementing a more efficient algorithm for finding
+the maximum weight matching.
+
+## BEHAVIOR WITH RANDOMIZED WELL POPULATIONS
+
+A series of BiGpairSEQ simulations were conducted using a cell sample file of 3.5 million unique T cells. From these cells,
+10 sample plate files were created. All of these sample plates had 96 wells, used an exponential distribution with a lambda of 0.6, and
+had a sequence dropout rate of 10%.
+
+The well populations of the plates were:
+* One sample plate with 1000 T cells/well
+* One sample plate with 2000 T cells/well
+* One sample plate with 3000 T cells/well
+* One sample plate with 4000 T cells/well
+* One sample plate with 5000 T cells/well
+* Five sample plates with each individual well's population randomized, from 1000 to 5000 T cells. (Average population ~3000 T cells/well.)
+
+All BiGpairSEQ simulations were run with a low overlap threshold of 3 and a high overlap threshold of 94.
+No optional filters were used, so pairing was attempted for all sequences with overlaps within the threshold values.
+
+Constant well population plate results:
+
+| |1000 Cell/Well Plate|2000 Cell/Well Plate|3000 Cell/Well Plate|4000 Cell/Well Plate|5000 Cell/Well Plate
+|---|---|---|---|---|---|
+|Total Alphas Found|6407|7330|7936|8278|8553|
+|Total Betas Found|6405|7333|7968|8269|8582|
+|Pairing Attempt Rate|0.661|0.653|0.600|0.579|0.559|
+|Correct Pairing Count|4231|4749|4723|4761|4750|
+|Incorrect Pairing Count|3|34|40|26|29|
+|Pairing Error Rate|0.000709|0.00711|0.00840|0.00543|0.00607|
+|Simulation Time (Seconds)|500|643|700|589|598|
+
+Randomized well population plate results:
+
+| |Random Plate 1 | Random Plate 2|Random Plate 3|Random Plate 4|Random Plate 5|Average|
+|---|---|---|---|---|---|---|
+Total Alphas Found|7853|7904|7964|7898|7917|7907|
+Total Betas Found|7851|7891|7920|7910|7894|7893|
+Pairing Attempt Rate|0.607|0.610|0.601|0.605|0.603|0.605|
+Correct Pairing Count|4718|4782|4721|4755|4731|4741|
+Incorrect Pairing Count|51|35|42|27|29|37|
+Pairing Error Rate|0.0107|0.00727|0.00882|0.00565|0.00609|0.00771|
+Simulation Time (Seconds)|590|677|730|618|615|646|
+
+The average results for the randomized plates are closest to the constant plate with 3000 T cells/well.
+This and several other tests indicate that BiGpairSEQ treats a sample plate with a highly variable number of T cells/well
+roughly as though it had a constant well population equal to the plate's average well population.
+
 ## 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
   * ~~*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.
-  * It is possible, though the modifications to the graph incur their own performance penalties. Need testing to see which option is best.
+  * It is possible, though the modifications to the graph incur their own performance penalties. Need testing to see which option is best. It may be computer-specific.
 * ~~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.
+  * Possible BiGpairSEQ advantage over pairSEQ: BiGpairSEQ is resilient to variations in well population sizes on a sample plate; pairSEQ is not due to nature of probability calculations.
     * 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. 
+* ~~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.
+  * ~~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
+* ~~Implement custom Vertex class to simplify code and make it easier to implement different MWM algorithms~~ DONE
+  * Advantage: would eliminate the need to use maps to associate vertices with sequences, which would make the code easier to understand.
+  * This also seems to be faster when using the same algorithm than the version with lots of maps, which is a nice bonus!
+* ~~Implement simulation of read depth, and of read errors. Pre-filter graph for difference in read count to eliminate spurious sequences.~~ DONE
+  * Pre-filtering based on comparing (read depth) * (occupancy) to (read count) for each sequence works extremely well
+* ~~Add read depth simulation options to CLI~~ DONE
+* ~~Update graphml output to reflect current Vertex class attributes~~ DONE
+  * Individual well data from the SequenceRecords could be included, if there's ever a reason for it
+* Update matching metadata output options in CLI
+* Update performance data in this readme
 * Re-implement CDR1 matching method
+* Refactor simulator code to collect all needed data in a single scan of the plate
+  * Currently it scans once for the vertices and then again for the edge weights. This made simulating read depth awkward, and incompatible with caching of plate files.
+  * This would be a fairly major rewrite of the simulator code, but could make things faster, and would definitely make them cleaner.
 * 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 auction algorithm for maximum weight matching
+* 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
+* Parameterize pre-filtering. Currently, sequences present in all wells are filtered out before constructing the graph, which massively reduces graph size. But, ideally, no pre-filtering would be necessary.
   
   
 ## CITATIONS
@@ -314,7 +390,7 @@ slightly less time than the simulation itself. Real elapsed time from start to f
 * [JGraphT](https://jgrapht.org) -- Graph theory data structures and algorithms
 * [JHeaps](https://www.jheaps.org) -- For pairing heap priority queue used in maximum weight matching algorithm
 * [Apache Commons CSV](https://commons.apache.org/proper/commons-csv/) -- For CSV file output
-* [Apache Commons CLI](https://commons.apache.org/proper/commons-cli/) -- To enable command line arguments for scripting. (**Awaiting re-implementation**.)
+* [Apache Commons CLI](https://commons.apache.org/proper/commons-cli/) -- To enable command line arguments for scripting.
 
 ## ACKNOWLEDGEMENTS
 BiGpairSEQ was conceived in collaboration with Dr. Alice MacQueen, who brought the original 

src/main/java/BiGpairSEQ.java

@@ -13,9 +13,10 @@ public class BiGpairSEQ {
     private static boolean cacheCells = false;
     private static boolean cachePlate = false;
     private static boolean cacheGraph = false;
-    private static String priorityQueueHeapType = "FIBONACCI";
+    private static HeapType priorityQueueHeapType = HeapType.FIBONACCI;
     private static boolean outputBinary = true;
     private static boolean outputGraphML = false;
+    private static final String version = "version 3.0";
 
     public static void main(String[] args) {
         if (args.length == 0) {
@@ -156,15 +157,15 @@ public class BiGpairSEQ {
     }
 
     public static String getPriorityQueueHeapType() {
-        return priorityQueueHeapType;
+        return priorityQueueHeapType.name();
     }
 
     public static void setPairingHeap() {
-        priorityQueueHeapType = "PAIRING";
+        priorityQueueHeapType = HeapType.PAIRING;
     }
 
     public static void setFibonacciHeap() {
-        priorityQueueHeapType = "FIBONACCI";
+        priorityQueueHeapType = HeapType.FIBONACCI;
     }
 
     public static boolean outputBinary() {return outputBinary;}
@@ -172,5 +173,5 @@ public class BiGpairSEQ {
 
     public static boolean outputGraphML() {return outputGraphML;}
     public static void setOutputGraphML(boolean b) {outputGraphML = b;}
-
+    public static String getVersion() { return version; }
 }

src/main/java/CellFileReader.java

@@ -12,7 +12,7 @@ import java.util.List;
 public class CellFileReader {
 
     private String filename;
-    private List<Integer[]> distinctCells = new ArrayList<>();
+    private List<String[]> distinctCells = new ArrayList<>();
     private Integer cdr1Freq;
 
     public CellFileReader(String filename) {
@@ -32,11 +32,11 @@ public class CellFileReader {
             CSVParser parser = new CSVParser(reader, cellFileFormat);
         ){
             for(CSVRecord record: parser.getRecords()) {
-                Integer[] cell = new Integer[4];
-                cell[0] = Integer.valueOf(record.get("Alpha CDR3"));
-                cell[1] = Integer.valueOf(record.get("Beta CDR3"));
-                cell[2] = Integer.valueOf(record.get("Alpha CDR1"));
-                cell[3] = Integer.valueOf(record.get("Beta CDR1"));
+                String[] cell = new String[4];
+                cell[0] = record.get("Alpha CDR3");
+                cell[1] = record.get("Beta CDR3");
+                cell[2] = record.get("Alpha CDR1");
+                cell[3] = record.get("Beta CDR1");
                 distinctCells.add(cell);
             }
 
@@ -47,8 +47,8 @@ public class CellFileReader {
         }
 
         //get CDR1 frequency
-        ArrayList<Integer> cdr1Alphas = new ArrayList<>();
-        for (Integer[] cell : distinctCells) {
+        ArrayList<String> cdr1Alphas = new ArrayList<>();
+        for (String[] cell : distinctCells) {
             cdr1Alphas.add(cell[3]);
         }
         double count = cdr1Alphas.stream().distinct().count();
@@ -62,14 +62,4 @@ public class CellFileReader {
     }
 
     public String getFilename() { return filename;}
-
-    //Refactor everything that uses this to have access to a Cell Sample and get the cells there instead.
-    public List<Integer[]> getListOfDistinctCellsDEPRECATED(){
-        return distinctCells;
-    }
-
-    public Integer getCellCountDEPRECATED() {
-        //Refactor everything that uses this to have access to a Cell Sample and get the count there instead.
-        return distinctCells.size();
-    }
 }

src/main/java/CellFileWriter.java

@@ -11,7 +11,7 @@ import java.util.List;
 public class CellFileWriter {
 
     private String[] headers = {"Alpha CDR3", "Beta CDR3", "Alpha CDR1", "Beta CDR1"};
-    List<Integer[]> cells;
+    List<String[]> cells;
     String filename;
     Integer cdr1Freq;
 
@@ -35,7 +35,7 @@ public class CellFileWriter {
             printer.printComment("Sample contains 1 unique CDR1 for every " + cdr1Freq + "unique CDR3s.");
             printer.printRecords(cells);
         } catch(IOException ex){
-            System.out.println("Could not make new file named "+filename);
+            System.out.println("Could not make new file named " + filename);
             System.err.println(ex);
         }
     }

src/main/java/CellSample.java

@@ -5,7 +5,7 @@ import java.util.stream.IntStream;
 
 public class CellSample {
 
-    private List<Integer[]> cells;
+    private List<String[]> cells;
     private Integer cdr1Freq;
 
     public CellSample(Integer numDistinctCells, Integer cdr1Freq){
@@ -13,31 +13,39 @@ public class CellSample {
         List<Integer> numbersCDR3 = new ArrayList<>();
         List<Integer> numbersCDR1 = new ArrayList<>();
         Integer numDistCDR3s = 2 * numDistinctCells + 1;
+        //Assign consecutive integers for each CDR3. This ensures they are all unique.
         IntStream.range(1, numDistCDR3s + 1).forEach(i -> numbersCDR3.add(i));
+        //After all CDR3s are assigned, start assigning consecutive integers to CDR1s
+        //There will usually be fewer integers in the CDR1 list, which will allow repeats below
         IntStream.range(numDistCDR3s + 1, numDistCDR3s + 1 + (numDistCDR3s / cdr1Freq) + 1).forEach(i -> numbersCDR1.add(i));
+        //randomize the order of the numbers in the lists
         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<>();
+        List<String[]> 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};
+            //Go through entire CDR3 list once, make pairs of alphas and betas
+            String tmpCDR3a = numbersCDR3.get(i).toString();
+            String tmpCDR3b = numbersCDR3.get(i+1).toString();
+            //Go through the (likely shorter) CDR1 list as many times as necessary, make pairs of alphas and betas
+            String tmpCDR1a = numbersCDR1.get(i % numbersCDR1.size()).toString();
+            String tmpCDR1b = numbersCDR1.get((i+1) % numbersCDR1.size()).toString();
+            //Make the array representing the cell
+            String[] tmp = {tmpCDR3a, tmpCDR3b, tmpCDR1a, tmpCDR1b};
+            //Add the cell to the list of distinct cells
             distinctCells.add(tmp);
         }
         this.cells = distinctCells;
     }
 
-    public CellSample(List<Integer[]> cells, Integer cdr1Freq){
+    public CellSample(List<String[]> cells, Integer cdr1Freq){
         this.cells = cells;
         this.cdr1Freq = cdr1Freq;
     }
 
-    public List<Integer[]> getCells(){
+    public List<String[]> getCells(){
         return cells;
     }
 

src/main/java/CommandLineInterface.java

@@ -35,6 +35,10 @@ import java.util.stream.Stream;
  * output : name of the output file
  * graphml : output a graphml file
  * binary : output a serialized binary object file
+ * IF SIMULATING READ DEPTH, ALL THESE ARE REQUIRED. Absence indicates not simulating read depth
+ * readdepth: number of reads per sequence
+ * readerrorprob: probability of reading a sequence incorrectly
+ * errcollisionprob: probability of two read errors being identical
  *
  * Match flags:
  * graphFile : name of graph and data file to use as input
@@ -62,15 +66,18 @@ public class CommandLineInterface {
 
             if (line.hasOption("help")) {
                 HelpFormatter formatter = new HelpFormatter();
-                formatter.printHelp("BiGpairSEQ_Sim", mainOptions);
+                formatter.printHelp("BiGpairSEQ_Sim.jar", mainOptions);
                 System.out.println();
-                formatter.printHelp("BiGpairSEQ_SIM -cells", cellOptions);
+                formatter.printHelp("BiGpairSEQ_Sim.jar -cells", cellOptions);
                 System.out.println();
-                formatter.printHelp("BiGpairSEQ_Sim -plate", plateOptions);
+                formatter.printHelp("BiGpairSEQ_Sim.jar -plate", plateOptions);
                 System.out.println();
-                formatter.printHelp("BiGpairSEQ_Sim -graph", graphOptions);
+                formatter.printHelp("BiGpairSEQ_Sim.jar -graph", graphOptions);
                 System.out.println();
-                formatter.printHelp("BiGpairSEQ_Sim -match", matchOptions);
+                formatter.printHelp("BiGpairSEQ_Sim.jar -match", matchOptions);
+            }
+            else if (line.hasOption("version")) {
+                System.out.println("BiGpairSEQ_Sim " + BiGpairSEQ.getVersion());
             }
             else if (line.hasOption("cells")) {
                 line = parser.parse(cellOptions, Arrays.copyOfRange(args, 1, args.length));
@@ -139,7 +146,20 @@ public class CommandLineInterface {
                 CellSample cells = getCells(cellFilename);
                 //get plate
                 Plate plate = getPlate(plateFilename);
-                GraphWithMapData graph = Simulator.makeGraph(cells, plate, false);
+                GraphWithMapData graph;
+                Integer readDepth = 1;
+                Double readErrorRate = 0.0;
+                Double errorCollisionRate = 0.0;
+                if (line.hasOption("rd")) {
+                    readDepth = Integer.parseInt(line.getOptionValue("rd"));
+                }
+                if (line.hasOption("err")) {
+                    readErrorRate = Double.parseDouble(line.getOptionValue("err"));
+                }
+                if (line.hasOption("coll")) {
+                    errorCollisionRate = Double.parseDouble(line.getOptionValue("coll"));
+                }
+                graph = Simulator.makeCDR3Graph(cells, plate, readDepth, readErrorRate, errorCollisionRate, false);
                 if (!line.hasOption("no-binary")) { //output binary file unless told not to
                     GraphDataObjectWriter writer = new GraphDataObjectWriter(outputFilename, graph, false);
                     writer.writeDataToFile();
@@ -153,17 +173,24 @@ public class CommandLineInterface {
             else if (line.hasOption("match")) { //can add a flag for which match type in future, spit this in two
                 line = parser.parse(matchOptions, Arrays.copyOfRange(args, 1, args.length));
                 String graphFilename = line.getOptionValue("g");
-                String outputFilename = line.getOptionValue("o");
+
+                String outputFilename;
+                if(line.hasOption("o")) {
+                    outputFilename = line.getOptionValue("o");
+                }
+                else {
+                    outputFilename = null;
+                }
                 Integer minThreshold = Integer.parseInt(line.getOptionValue("min"));
                 Integer maxThreshold = Integer.parseInt(line.getOptionValue("max"));
-                Integer minOverlapPct;
+                int minOverlapPct;
                 if (line.hasOption("minpct")) { //see if this filter is being used
                     minOverlapPct = Integer.parseInt(line.getOptionValue("minpct"));
                 }
                 else {
                     minOverlapPct = 0;
                 }
-                Integer maxOccupancyDiff;
+                int maxOccupancyDiff;
                 if (line.hasOption("maxdiff")) { //see if this filter is being used
                     maxOccupancyDiff = Integer.parseInt(line.getOptionValue("maxdiff"));
                 }
@@ -173,10 +200,38 @@ public class CommandLineInterface {
                 GraphWithMapData graph = getGraph(graphFilename);
                 MatchingResult result = Simulator.matchCDR3s(graph, graphFilename, minThreshold, maxThreshold,
                         maxOccupancyDiff, minOverlapPct, false);
-                MatchingFileWriter writer = new MatchingFileWriter(outputFilename, result);
-                writer.writeResultsToFile();
+                if(outputFilename != null){
+                    MatchingFileWriter writer = new MatchingFileWriter(outputFilename, result);
+                    writer.writeResultsToFile();
+                }
                 //can put a bunch of ifs for outputting various things from the MatchingResult to System.out here
                 //after I put those flags in the matchOptions
+                if(line.hasOption("print-metadata")) {
+                    for (String k : result.getMetadata().keySet()) {
+                        System.out.println(k + ": " + result.getMetadata().get(k));
+                    }
+                }
+                if(line.hasOption("print-error")) {
+                    System.out.println("pairing error rate: " + result.getPairingErrorRate());
+                }
+                if(line.hasOption("print-attempt")) {
+                    System.out.println("pairing attempt rate: " +result.getPairingAttemptRate());
+                }
+                if(line.hasOption("print-correct")) {
+                    System.out.println("correct pairings: " + result.getCorrectPairingCount());
+                }
+                if(line.hasOption("print-incorrect")) {
+                    System.out.println("incorrect pairings: " + result.getIncorrectPairingCount());
+                }
+                if(line.hasOption("print-alphas")) {
+                    System.out.println("total alphas found: " + result.getAlphaCount());
+                }
+                if(line.hasOption("print-betas")) {
+                    System.out.println("total betas found: " + result.getBetaCount());
+                }
+                if(line.hasOption("print-time")) {
+                    System.out.println("simulation time (seconds): " + result.getSimulationTime());
+                }
             }
         }
         catch (ParseException exp) {
@@ -216,8 +271,11 @@ public class CommandLineInterface {
                 .longOpt("match-cdr3")
                 .desc("Matches CDR3s. Requires a graph/data file.")
                 .build();
+        Option printVersion = Option.builder("version")
+                .desc("Prints the program version number to stdout").build();
         OptionGroup mainGroup = new OptionGroup();
         mainGroup.addOption(help);
+        mainGroup.addOption(printVersion);
         mainGroup.addOption(makeCells);
         mainGroup.addOption(makePlate);
         mainGroup.addOption(makeGraph);
@@ -297,7 +355,7 @@ public class CommandLineInterface {
                 .desc("Randomize well populations on sample plate. Takes two arguments: the minimum possible population and the maximum possible population.")
                 .hasArgs()
                 .numberOfArgs(2)
-                .argName("minimum maximum")
+                .argName("min> <max")
                 .build();
         Option specificWellPopulations = Option.builder("pop")
                 .desc("The well populations for each section of the sample plate. There will be as many sections as there are populations given.")
@@ -326,28 +384,49 @@ public class CommandLineInterface {
         Options graphOptions = new Options();
         Option cellFilename = Option.builder("c")
                 .longOpt("cell-file")
-                .desc("Cell sample file to use for checking accuracy")
+                .desc("Cell sample file to use for checking pairing accuracy")
                 .hasArg()
                 .argName("filename")
                 .required().build();
         Option plateFilename = Option.builder("p")
                 .longOpt("plate-filename")
-                .desc("Sample plate file (made from given cell sample file) to construct graph from")
+                .desc("Sample plate file from which to construct graph")
                 .hasArg()
                 .argName("filename")
                 .required().build();
         Option outputGraphML = Option.builder("graphml")
-                .desc("Output GraphML file")
+                .desc("(Optional) Output GraphML file")
                 .build();
         Option outputSerializedBinary = Option.builder("nb")
                 .longOpt("no-binary")
-                .desc("Don't output serialized binary file")
+                .desc("(Optional) Don't output serialized binary file")
+                .build();
+        Option readDepth = Option.builder("rd")
+                .longOpt("read-depth")
+                .desc("(Optional) The number of times to read each sequence.")
+                .hasArg()
+                .argName("depth")
+                .build();
+        Option readErrorProb = Option.builder("err")
+                .longOpt("read-error-prob")
+                .desc("(Optional) The probability that a sequence will be misread. (0.0 - 1.0)")
+                .hasArg()
+                .argName("prob")
+                .build();
+        Option errorCollisionProb = Option.builder("coll")
+                .longOpt("error-collision-prob")
+                .desc("(Optional) The probability that two misreads will produce the same spurious sequence. (0.0 - 1.0)")
+                .hasArg()
+                .argName("prob")
                 .build();
         graphOptions.addOption(cellFilename);
         graphOptions.addOption(plateFilename);
         graphOptions.addOption(outputFileOption());
         graphOptions.addOption(outputGraphML);
         graphOptions.addOption(outputSerializedBinary);
+        graphOptions.addOption(readDepth);
+        graphOptions.addOption(readErrorProb);
+        graphOptions.addOption(errorCollisionProb);
         return graphOptions;
     }
 
@@ -379,15 +458,46 @@ public class CommandLineInterface {
                 .hasArg()
                 .argName("number")
                 .build();
-        matchCDR3options.addOption(graphFilename);
-        matchCDR3options.addOption(minOccupancyOverlap);
-        matchCDR3options.addOption(maxOccupancyOverlap);
-        matchCDR3options.addOption(minOverlapPercent);
-        matchCDR3options.addOption(maxOccupancyDifference);
-        matchCDR3options.addOption(outputFileOption());
-        //options for output to System.out
-       //Option printPairingErrorRate = Option.builder()
+        Option outputFile = Option.builder("o") //can't call the method this time, because this one's optional
+                .longOpt("output-file")
+                .hasArg()
+                .argName("filename")
+                .desc("(Optional) Name of output the output file. If not present, no file will be written.")
+                .build();
+        matchCDR3options.addOption(graphFilename)
+                .addOption(minOccupancyOverlap)
+                .addOption(maxOccupancyOverlap)
+                .addOption(minOverlapPercent)
+                .addOption(maxOccupancyDifference)
+                .addOption(outputFile);
 
+        //options for output to System.out
+        Option printAlphaCount = Option.builder().longOpt("print-alphas")
+                .desc("(Optional) Print the number of distinct alpha sequences to stdout.").build();
+        Option printBetaCount = Option.builder().longOpt("print-betas")
+                .desc("(Optional) Print the number of distinct beta sequences to stdout.").build();
+        Option printTime = Option.builder().longOpt("print-time")
+                .desc("(Optional) Print the total simulation time to stdout.").build();
+        Option printErrorRate = Option.builder().longOpt("print-error")
+               .desc("(Optional) Print the pairing error rate to stdout").build();
+        Option printAttempt = Option.builder().longOpt("print-attempt")
+               .desc("(Optional) Print the pairing attempt rate to stdout").build();
+        Option printCorrect = Option.builder().longOpt("print-correct")
+               .desc("(Optional) Print the number of correct pairs to stdout").build();
+        Option printIncorrect = Option.builder().longOpt("print-incorrect")
+               .desc("(Optional) Print the number of incorrect pairs to stdout").build();
+        Option printMetadata = Option.builder().longOpt("print-metadata")
+                       .desc("(Optional) Print a full summary of the matching results to stdout.").build();
+
+       matchCDR3options
+               .addOption(printErrorRate)
+               .addOption(printAttempt)
+               .addOption(printCorrect)
+               .addOption(printIncorrect)
+               .addOption(printMetadata)
+               .addOption(printAlphaCount)
+               .addOption(printBetaCount)
+               .addOption(printTime);
         return matchCDR3options;
     }
 

src/main/java/GraphDataObjectReader.java

@@ -4,7 +4,7 @@ public class GraphDataObjectReader {
 
     private GraphWithMapData data;
     private String filename;
-    private boolean verbose = true;
+
 
     public GraphDataObjectReader(String filename, boolean verbose) throws IOException {
         if(!filename.matches(".*\\.ser")){
@@ -15,10 +15,13 @@ public class GraphDataObjectReader {
             BufferedInputStream fileIn = new BufferedInputStream(new FileInputStream(filename));
                 ObjectInputStream in = new ObjectInputStream(fileIn))
         {
-            System.out.println("Reading graph data from file. This may take some time");
-            System.out.println("File I/O time is not included in results");
+            if (verbose) {
+                System.out.println("Reading graph data from file. This may take some time");
+                System.out.println("File I/O time is not included in results");
+            }
             data = (GraphWithMapData) in.readObject();
         } catch (FileNotFoundException | ClassNotFoundException ex) {
+            System.out.println("Graph/data file " + filename + " not found.");
             ex.printStackTrace();
         }
     }

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,67 @@ 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));
+        ga.put("read depth", DefaultAttribute.createAttribute(data.getReadDepth().toString()));
+        ga.put("read error rate", DefaultAttribute.createAttribute(data.getReadErrorRate().toString()));
+        ga.put("error collision rate", DefaultAttribute.createAttribute(data.getErrorCollisionRate().toString()));
+        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
+            //NEED TO ADD NEW FIELD FOR READ COUNT
             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()));
+                attributes.put("read count", DefaultAttribute.createAttribute(v.getReadCount()));
                 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);
+            exporter.registerAttribute("read count", AttributeCategory.NODE, AttributeType.STRING);
             //export the graph
             exporter.exportGraph(graph, writer);
         } catch(IOException ex){
@@ -81,4 +97,3 @@ public class GraphMLFileWriter {
         }
     }
 }
-

src/main/java/GraphModificationFunctions.java

@@ -2,23 +2,24 @@ 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,
+    //remove over- and under-weight edges, return removed edges
+    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<>();
         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,31 +27,27 @@ public interface GraphModificationFunctions {
             }
         }
         if(saveEdges) {
-            for (Integer[] edge : removedEdges) {
+            for (Vertex[] edge : removedEdges.keySet()) {
                 graph.removeEdge(edge[0], edge[1]);
             }
         }
         return removedEdges;
     }
 
-    //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,
+    //Remove edges for pairs with large occupancy discrepancy, return removed edges
+    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,34 +55,30 @@ public interface GraphModificationFunctions {
             }
         }
         if(saveEdges) {
-            for (Integer[] edge : removedEdges) {
+            for (Vertex[] edge : removedEdges.keySet()) {
                 graph.removeEdge(edge[0], edge[1]);
             }
         }
         return removedEdges;
     }
 
-    //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,
+    //Remove edges for pairs where overlap size is significantly lower than the well occupancy, return removed edges
+    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);
+                if (saveEdges) {
+                    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,19 +86,53 @@ 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 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()) {
             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;
@@ -14,33 +15,41 @@ public class GraphWithMapData implements java.io.Serializable {
     private Integer[] wellPopulations;
     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 int readDepth;
+    private double readErrorRate;
+    private double errorCollisionRate;
+    private final Map<String, String> 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 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<String, String> distCellsMapAlphaKey, Integer alphaCount, Integer betaCount,
+                            Integer readDepth, Double readErrorRate, Double errorCollisionRate, Duration time){
+
+//                            Map<Integer, Integer> plateVtoAMap,
+//                            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.readDepth = readDepth;
+        this.readErrorRate = readErrorRate;
+        this.errorCollisionRate = errorCollisionRate;
         this.time = time;
     }
 
@@ -64,33 +73,35 @@ public class GraphWithMapData implements java.io.Serializable {
         return betaCount;
     }
 
-    public Map<Integer, Integer> getDistCellsMapAlphaKey() {
+    public Map<String, String> getDistCellsMapAlphaKey() {
         return distCellsMapAlphaKey;
     }
 
-    public Map<Integer, Integer> getPlateVtoAMap() {
-        return plateVtoAMap;
-    }
+//    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> 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 Integer getReadDepth() { return readDepth; }
 
     public Duration getTime() {
         return time;
@@ -103,4 +114,12 @@ public class GraphWithMapData implements java.io.Serializable {
     public String getSourceFilename() {
         return sourceFilename;
     }
+
+    public Double getReadErrorRate() {
+        return readErrorRate;
+    }
+
+    public Double getErrorCollisionRate() {
+        return errorCollisionRate;
+    }
 }

src/main/java/HeapType.java

@@ -0,0 +1,4 @@
+public enum HeapType {
+    FIBONACCI,
+    PAIRING
+}

src/main/java/InteractiveInterface.java

@@ -250,6 +250,11 @@ public class InteractiveInterface {
         String filename = null;
         String cellFile = null;
         String plateFile = null;
+        Boolean simulateReadDepth = false;
+        //number of times to read each sequence in a well
+        int readDepth = 1;
+        double readErrorRate = 0.0;
+        double errorCollisionRate = 0.0;
         try {
             String str = "\nGenerating bipartite weighted graph encoding occupancy overlap data ";
             str = str.concat("\nrequires a cell sample file and a sample plate file.");
@@ -258,7 +263,36 @@ 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("\nEnable simulation of sequence read depth and sequence read errors? (y/n)");
+            System.out.println("NOTE: sample plate data cannot be cached when simulating read errors");
+            String ans = sc.next();
+            Pattern pattern = Pattern.compile("(?:yes|y)", Pattern.CASE_INSENSITIVE);
+            Matcher matcher = pattern.matcher(ans);
+            if(matcher.matches()){
+                simulateReadDepth = true;
+            }
+            if (simulateReadDepth) {
+                BiGpairSEQ.setCachePlate(false);
+                BiGpairSEQ.clearPlateInMemory();
+                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) {
+                    throw new InputMismatchException("The read error rate must be in the range [0.0, 1.0]");
+                }
+                System.out.println("\nPlease enter the probability of read error collision");
+                System.out.println("(the likelihood that two read errors produce the same spurious sequence)");
+                System.out.print("(0.0 to 1.0): ");
+                errorCollisionRate = sc.nextDouble();
+                if(errorCollisionRate < 0.0 || errorCollisionRate > 1.0) {
+                    throw new InputMismatchException("The error collision probability must be an in the range [0.0, 1.0]");
+                }
+            }
+            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) {
@@ -304,7 +338,7 @@ public class InteractiveInterface {
             System.out.println("Returning to main menu.");
         }
         else{
-            GraphWithMapData data = Simulator.makeGraph(cellSample, plate, true);
+            GraphWithMapData data = Simulator.makeCDR3Graph(cellSample, plate, readDepth, readErrorRate, errorCollisionRate, true);
             assert filename != null;
             if(BiGpairSEQ.outputBinary()) {
                 GraphDataObjectWriter dataWriter = new GraphDataObjectWriter(filename, data);
@@ -504,7 +538,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 {
@@ -570,6 +604,8 @@ public class InteractiveInterface {
     }
 
     private static void acknowledge(){
+        System.out.println("BiGpairSEQ_Sim " + BiGpairSEQ.getVersion());
+        System.out.println();
         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.");
         System.out.println();

src/main/java/MatchingResult.java

@@ -9,27 +9,34 @@ public class MatchingResult {
     private final List<String> comments;
     private final List<String> headers;
     private final List<List<String>> allResults;
-    private final Map<Integer, Integer> matchMap;
-    private final Duration time;
+    private final Map<String, String> matchMap;
 
     public MatchingResult(Map<String, String> metadata, List<String> headers,
-                          List<List<String>> allResults, Map<Integer, Integer>matchMap, Duration time){
+                          List<List<String>> allResults, Map<String, String>matchMap){
         /*
          * POSSIBLE KEYS FOR METADATA MAP ARE:
          * sample plate filename *
          * graph filename *
+         * matching weight *
          * well populations *
-         * total alphas found *
-         * total betas found *
-         * high overlap threshold *
-         * low overlap threshold *
-         * maximum occupancy difference *
-         * minimum overlap percent *
+         * sequence read depth *
+         * sequence read error rate *
+         * read error collision rate *
+         * total alphas read from plate *
+         * total betas read from plate *
+         * alphas in graph (after pre-filtering) *
+         * betas in graph (after pre-filtering) *
+         * high overlap threshold for pairing *
+         * low overlap threshold for pairing *
+         * maximum occupancy difference for pairing *
+         * minimum overlap percent for pairing *
          * pairing attempt rate *
          * correct pairing count *
          * incorrect pairing count *
          * pairing error rate *
-         * simulation time (seconds)
+         * time to generate graph (seconds) *
+         * time to pair sequences (seconds) *
+         * total simulation time (seconds) *
          */
         this.metadata = metadata;
         this.comments = new ArrayList<>();
@@ -39,8 +46,6 @@ public class MatchingResult {
         this.headers = headers;
         this.allResults = allResults;
         this.matchMap = matchMap;
-        this.time = time;
-
     }
 
     public Map<String, String> getMetadata() {return metadata;}
@@ -57,13 +62,13 @@ public class MatchingResult {
         return headers;
     }
 
-    public Map<Integer, Integer> getMatchMap() {
+    public Map<String, String> getMatchMap() {
         return matchMap;
     }
 
-    public Duration getTime() {
-        return time;
-    }
+//    public Duration getTime() {
+//        return time;
+//    }
 
     public String getPlateFilename() {
         return metadata.get("sample plate filename");
@@ -84,20 +89,20 @@ public class MatchingResult {
     }
 
     public Integer getAlphaCount() {
-        return Integer.parseInt(metadata.get("total alpha count"));
+        return Integer.parseInt(metadata.get("total alphas read from plate"));
     }
 
     public Integer getBetaCount() {
-        return Integer.parseInt(metadata.get("total beta count"));
+        return Integer.parseInt(metadata.get("total betas read from plate"));
     }
 
-    public Integer getHighOverlapThreshold() { return Integer.parseInt(metadata.get("high overlap threshold"));}
+    public Integer getHighOverlapThreshold() { return Integer.parseInt(metadata.get("high overlap threshold for pairing"));}
 
-    public Integer getLowOverlapThreshold() { return Integer.parseInt(metadata.get("low overlap threshold"));}
+    public Integer getLowOverlapThreshold() { return Integer.parseInt(metadata.get("low overlap threshold for pairing"));}
 
-    public Integer getMaxOccupancyDifference() { return Integer.parseInt(metadata.get("maximum occupancy difference"));}
+    public Integer getMaxOccupancyDifference() { return Integer.parseInt(metadata.get("maximum occupancy difference for pairing"));}
 
-    public Integer getMinOverlapPercent() { return Integer.parseInt(metadata.get("minimum overlap percent"));}
+    public Integer getMinOverlapPercent() { return Integer.parseInt(metadata.get("minimum overlap percent for pairing"));}
 
     public Double getPairingAttemptRate() { return Double.parseDouble(metadata.get("pairing attempt rate"));}
 
@@ -107,6 +112,6 @@ public class MatchingResult {
 
     public Double getPairingErrorRate() { return Double.parseDouble(metadata.get("pairing error rate"));}
 
-    public String getSimulationTime() { return metadata.get("simulation time (seconds)"); }
+    public String getSimulationTime() { return metadata.get("total simulation time (seconds)"); }
 
 }

src/main/java/Plate.java

@@ -5,13 +5,14 @@ TODO: Implement exponential distribution using inversion method - DONE
 TODO: Implement discrete frequency distributions using Vose's Alias Method
  */
 
+
 import java.util.*;
 
 public class Plate {
     private CellSample cells;
     private String sourceFile;
     private String filename;
-    private List<List<Integer[]>> wells;
+    private List<List<String[]>> wells;
     private final Random rand = BiGpairSEQ.getRand();
     private int size;
     private double error;
@@ -48,13 +49,13 @@ public class Plate {
     }
 
     //constructor for returning a Plate from a PlateFileReader
-    public Plate(String filename, List<List<Integer[]>> wells) {
+    public Plate(String filename, List<List<String[]>> wells) {
         this.filename = filename;
         this.wells = wells;
         this.size = wells.size();
 
         List<Integer> concentrations = new ArrayList<>();
-        for (List<Integer[]> w: wells) {
+        for (List<String[]> w: wells) {
             if(!concentrations.contains(w.size())){
                 concentrations.add(w.size());
             }
@@ -65,7 +66,7 @@ public class Plate {
         }
     }
 
-    private void fillWellsExponential(List<Integer[]> cells, double lambda){
+    private void fillWellsExponential(List<String[]> cells, double lambda){
         this.lambda = lambda;
         exponential = true;
         int numSections = populations.length;
@@ -74,17 +75,17 @@ public class Plate {
         int n;
         while (section < numSections){
             for (int i = 0; i < (size / numSections); i++) {
-                List<Integer[]> well = new ArrayList<>();
+                List<String[]> well = new ArrayList<>();
                 for (int j = 0; j < populations[section]; j++) {
                     do {
                         //inverse transform sampling: for random number u in [0,1), x = log(1-u) / (-lambda)
                         m = (Math.log10((1 - rand.nextDouble()))/(-lambda)) * Math.sqrt(cells.size());
                     } while (m >= cells.size() || m < 0);
                     n = (int) Math.floor(m);
-                    Integer[] cellToAdd = cells.get(n).clone();
+                    String[] cellToAdd = cells.get(n).clone();
                     for(int k = 0; k < cellToAdd.length; k++){
-                        if(Math.abs(rand.nextDouble()) < error){//error applied to each seqeunce
-                            cellToAdd[k] = -1;
+                        if(Math.abs(rand.nextDouble()) <= error){//error applied to each sequence
+                            cellToAdd[k] = "-1";
                         }
                     }
                     well.add(cellToAdd);
@@ -95,7 +96,7 @@ public class Plate {
         }
     }
 
-    private void fillWells( List<Integer[]> cells, double stdDev) {
+    private void fillWells( List<String[]> cells, double stdDev) {
         this.stdDev = stdDev;
         int numSections = populations.length;
         int section = 0;
@@ -103,16 +104,16 @@ public class Plate {
         int n;
         while (section < numSections){
             for (int i = 0; i < (size / numSections); i++) {
-                List<Integer[]> well = new ArrayList<>();
+                List<String[]> well = new ArrayList<>();
                 for (int j = 0; j < populations[section]; j++) {
                     do {
                         m = (rand.nextGaussian() * stdDev) + (cells.size() / 2);
                     } while (m >= cells.size() || m < 0);
                     n = (int) Math.floor(m);
-                    Integer[] cellToAdd = cells.get(n).clone();
+                    String[] cellToAdd = cells.get(n).clone();
                     for(int k = 0; k < cellToAdd.length; k++){
                         if(Math.abs(rand.nextDouble()) < error){//error applied to each sequence
-                            cellToAdd[k] = -1;
+                            cellToAdd[k] = "-1";
                         }
                     }
                     well.add(cellToAdd);
@@ -143,39 +144,185 @@ public class Plate {
         return error;
     }
 
-    public List<List<Integer[]>> getWells() {
+    public List<List<String[]>> getWells() {
         return wells;
     }
 
-    //returns a map of the counts of the sequence at cell index sIndex, in all wells
-    public Map<Integer, Integer> assayWellsSequenceS(int... sIndices){
-        return this.assayWellsSequenceS(0, size, sIndices);
-    }
+//    //returns a map of the counts of the sequence at cell index sIndex, in all wells
+//    public void assayWellsSequenceS(Map<String, Integer> sequences, int... sIndices){
+//        this.assayWellsSequenceS(sequences, 0, size, sIndices);
+//    }
+//
+//    //returns a map of the counts of the sequence at cell index sIndex, in a specific well
+//    public void assayWellsSequenceS(Map<String, Integer> sequences, int n, int... sIndices) {
+//        this.assayWellsSequenceS(sequences, n, n+1, sIndices);
+//    }
+//
+//    //returns a map of the counts of the sequence at cell index sIndex, in a range of wells
+//    public void assayWellsSequenceS(Map<String, Integer> sequences, int start, int end, int... sIndices) {
+//        for(int sIndex: sIndices){
+//            for(int i = start; i < end; i++){
+//                countSequences(sequences, wells.get(i), sIndex);
+//            }
+//        }
+//    }
+//    //For the sequences at cell indices sIndices, counts number of unique sequences in the given well into the given map
+//    private void countSequences(Map<String, Integer> wellMap, List<String[]> well, int... sIndices) {
+//        for(String[] cell : well) {
+//            for(int sIndex: sIndices){
+//                //skip dropout sequences, which have value -1
+//                if(!"-1".equals(cell[sIndex])){
+//                    wellMap.merge(cell[sIndex], 1, (oldValue, newValue) -> oldValue + newValue);
+//                }
+//            }
+//        }
+//    }
 
-    //returns a map of the counts of the sequence at cell index sIndex, in a specific well
-    public Map<Integer, Integer> assayWellsSequenceS(int n, int... sIndices) { return this.assayWellsSequenceS(n, n+1, sIndices);}
-
-    //returns a map of the counts of the sequence at cell index sIndex, in a range of wells
-    public Map<Integer, Integer> assayWellsSequenceS(int start, int end, int... sIndices) {
-        Map<Integer,Integer> assay = new HashMap<>();
-        for(int pIndex: sIndices){
-            for(int i = start; i < end; i++){
-                countSequences(assay, wells.get(i), pIndex);
-            }
-        }
-        return assay;
-    }
-    //For the sequences at cell indices sIndices, counts number of unique sequences in the given well into the given map
-    private void countSequences(Map<Integer, Integer> wellMap, List<Integer[]> well, int... sIndices) {
-        for(Integer[] cell : well) {
-            for(int sIndex: sIndices){
-                if(cell[sIndex] != -1){
-                    wellMap.merge(cell[sIndex], 1, (oldValue, newValue) -> oldValue + newValue);
+    //For the sequences at cell indices sIndices, counts number of unique sequences in all well into the given map
+    public Map<String, SequenceRecord> countSequences(Integer readDepth, Double readErrorRate,
+                                                       Double errorCollisionRate, int... sIndices) {
+        SequenceType[] sequenceTypes = EnumSet.allOf(SequenceType.class).toArray(new SequenceType[0]);
+        Map<String, Integer> distinctMisreadCounts = new HashMap<>();
+        Map<String, SequenceRecord> sequenceMap = new LinkedHashMap<>();
+        for (int well = 0; well < size; well++) {
+            for (String[] cell : wells.get(well)) {
+                for (int sIndex : sIndices) {
+                    //skip dropout sequences, which have value -1
+                    if (!"-1".equals(cell[sIndex])) {
+                        for (int j = 0; j < readDepth; j++) {
+                            //Misread sequence
+                            if (rand.nextDouble() < readErrorRate) {
+                                StringBuilder spurious = new StringBuilder(cell[sIndex]);
+                                //if this sequence hasn't been misread before, or the read error is unique,
+                                //append one more "*" than has been appended before
+                                if (rand.nextDouble() > errorCollisionRate || !distinctMisreadCounts.containsKey(cell[sIndex])) {
+                                    distinctMisreadCounts.merge(cell[sIndex], 1, (oldValue, newValue) -> oldValue + newValue);
+                                    for (int k = 0; k < distinctMisreadCounts.get(cell[sIndex]); k++) {
+                                        spurious.append("*");
+                                    }
+                                    SequenceRecord tmp = new SequenceRecord(spurious.toString(), sequenceTypes[sIndex]);
+                                    tmp.addRead(well);
+                                    sequenceMap.put(spurious.toString(), tmp);
+                                }
+                                //if this is a read error collision, randomly choose a number of "*"s that has been appended before
+                                else {
+                                    int starCount = rand.nextInt(distinctMisreadCounts.get(cell[sIndex]));
+                                    for (int k = 0; k < starCount; k++) {
+                                        spurious.append("*");
+                                    }
+                                    sequenceMap.get(spurious.toString()).addRead(well);
+                                }
+                            }
+                            //sequence is read correctly
+                            else {
+                                if (!sequenceMap.containsKey(cell[sIndex])) {
+                                    SequenceRecord tmp = new SequenceRecord(cell[sIndex], sequenceTypes[sIndex]);
+                                    tmp.addRead(well);
+                                    sequenceMap.put(cell[sIndex], tmp);
+                                } else {
+                                    sequenceMap.get(cell[sIndex]).addRead(well);
+                                }
+                            }
+                        }
+                    }
                 }
             }
         }
+        return sequenceMap;
     }
 
+
+//    //returns a map of the counts of the sequence at cell index sIndex, in all 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> misreadCounts, Map<String, Integer> occupancyMap, Map<String, Integer> readCountMap,
+//                                                 int readDepth, double readErrorProb, double errorCollisionProb, int... sIndices) {
+//        this.assayWellsSequenceSWithReadDepth(misreadCounts, 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> misreadCounts, Map<String, Integer> occupancyMap, Map<String, Integer> readCountMap,
+//                                                 int readDepth, double readErrorProb, double errorCollisionProb,
+//                                                 int n, int... sIndices) {
+//       this.assayWellsSequenceSWithReadDepth(misreadCounts, 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 void assayWellsSequenceSWithReadDepth(Map<String, Integer> misreadCounts, Map<String, Integer> occupancyMap, Map<String, Integer> readCountMap,
+//                                                 int readDepth, double readErrorProb, double errorCollisionProb,
+//                                                 int start, int end, int... sIndices) {
+//        for(int sIndex: sIndices){
+//            for(int i = start; i < end; i++){
+//                countSequencesWithReadDepth(misreadCounts, occupancyMap, readCountMap, readDepth, readErrorProb, errorCollisionProb, wells.get(i), sIndex);
+//            }
+//        }
+//    }
+//
+//    //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.
+//    //NOTE: this function changes the content of the well, adding spurious cells to contain the misread sequences
+//    //(this is necessary because, in the simulation, the plate is read multiple times, but random misreads can only
+//    //be simulated once).
+//    //(Possibly I should refactor all of this to only require a single plate assay, to speed things up. Or at least
+//    //to see if it would speed things up.)
+//    private void countSequencesWithReadDepth(Map<String, Integer> distinctMisreadCounts, Map<String, Integer> occupancyMap, Map<String, Integer> readCountMap,
+//                                             int readDepth, double readErrorProb, double errorCollisionProb,
+//                                             List<String[]> well, int... sIndices) {
+//        //list of spurious cells to add to well after counting
+//        List<String[]> spuriousCells = new ArrayList<>();
+//        for(String[] cell : well) {
+//            //new potential spurious cell for each cell that gets read
+//            String[] spuriousCell = new String[SequenceType.values().length];
+//            //initialize spurious cell with all dropout sequences
+//            Arrays.fill(spuriousCell, "-1");
+//            //has a read error occurred?
+//            boolean readError = false;
+//            for(int sIndex: sIndices){
+//                //skip dropout sequences, which have value "-1"
+//                if(!"-1".equals(cell[sIndex])){
+//                    Map<String, Integer> sequencesWithReadCounts = new LinkedHashMap<>();
+//                    for(int i = 0; i < readDepth; i++) {
+//                        if (rand.nextDouble() <= readErrorProb) {
+//                            readError = true;
+//                            //Read errors are represented by appending "*"s to the end of the sequence some number of times
+//                            StringBuilder spurious = new StringBuilder(cell[sIndex]);
+//                            //if this sequence hasn't been misread before, or the read error is unique,
+//                            //append one more "*" than has been appended before
+//                            if (!distinctMisreadCounts.containsKey(cell[sIndex]) || rand.nextDouble() > errorCollisionProb) {
+//                                distinctMisreadCounts.merge(cell[sIndex], 1, (oldValue, newValue) -> oldValue + newValue);
+//                                for (int j = 0; j < distinctMisreadCounts.get(cell[sIndex]); j++) {
+//                                    spurious.append("*");
+//                                }
+//                            }
+//                            //if this is a read error collision, randomly choose a number of "*"s that has been appended before
+//                            else {
+//                                int starCount = rand.nextInt(distinctMisreadCounts.get(cell[sIndex]));
+//                                for (int j = 0; j < starCount; j++) {
+//                                    spurious.append("*");
+//                                }
+//                            }
+//                            sequencesWithReadCounts.merge(spurious.toString(), 1, (oldValue, newValue) -> oldValue + newValue);
+//                            //add spurious sequence to spurious cell
+//                            spuriousCell[sIndex] = spurious.toString();
+//                        }
+//                        else {
+//                            sequencesWithReadCounts.merge(cell[sIndex], 1, (oldValue, newValue) -> oldValue + newValue);
+//                        }
+//                    }
+//                    for(String seq : sequencesWithReadCounts.keySet()) {
+//                        occupancyMap.merge(seq, 1, (oldValue, newValue) -> oldValue + newValue);
+//                        readCountMap.merge(seq, sequencesWithReadCounts.get(seq), (oldValue, newValue) -> oldValue + newValue);
+//                    }
+//                }
+//            }
+//            if (readError) { //only add a new spurious cell if there was a read error
+//                spuriousCells.add(spuriousCell);
+//            }
+//        }
+//        //add all spurious cells to the well
+//        well.addAll(spuriousCells);
+//    }
+
     public String getSourceFileName() {
         return sourceFile;
     }

src/main/java/PlateFileReader.java

@@ -13,7 +13,7 @@ import java.util.regex.Pattern;
 
 public class PlateFileReader {
 
-    private List<List<Integer[]>> wells = new ArrayList<>();
+    private List<List<String[]>> wells = new ArrayList<>();
     private String filename;
 
     public PlateFileReader(String filename){
@@ -32,17 +32,17 @@ public class PlateFileReader {
             CSVParser parser = new CSVParser(reader, plateFileFormat);
         ){
             for(CSVRecord record: parser.getRecords()) {
-                List<Integer[]> well = new ArrayList<>();
+                List<String[]> well = new ArrayList<>();
                 for(String s: record) {
                     if(!"".equals(s)) {
-                        String[] intString = s.replaceAll("\\[", "")
+                        String[] sequences = s.replaceAll("\\[", "")
                                 .replaceAll("]", "")
                                 .replaceAll(" ", "")
                                 .split(",");
-                        //System.out.println(intString);
-                        Integer[] arr = new Integer[intString.length];
-                        for (int i = 0; i < intString.length; i++) {
-                            arr[i] = Integer.valueOf(intString[i]);
+                        //System.out.println(sequences);
+                        String[] arr = new String[sequences.length];
+                        for (int i = 0; i < sequences.length; i++) {
+                            arr[i] = sequences[i];
                         }
                         well.add(arr);
                     }

src/main/java/PlateFileWriter.java

@@ -10,7 +10,7 @@ import java.util.*;
 
 public class PlateFileWriter {
     private int size;
-    private List<List<Integer[]>> wells;
+    private List<List<String[]>> wells;
     private double stdDev;
     private double lambda;
     private Double error;
@@ -40,13 +40,13 @@ public class PlateFileWriter {
     }
 
     public void writePlateFile(){
-        Comparator<List<Integer[]>> listLengthDescending = Comparator.comparingInt(List::size);
+        Comparator<List<String[]>> listLengthDescending = Comparator.comparingInt(List::size);
         wells.sort(listLengthDescending.reversed());
         int maxLength = wells.get(0).size();
         List<List<String>> wellsAsStrings = new ArrayList<>();
-        for (List<Integer[]> w: wells){
+        for (List<String[]> w: wells){
             List<String> tmp = new ArrayList<>();
-            for(Integer[] c: w) {
+            for(String[] c: w) {
                 tmp.add(Arrays.toString(c));
             }
             wellsAsStrings.add(tmp);

src/main/java/SequenceRecord.java

@@ -0,0 +1,65 @@
+/*
+Class to represent individual sequences, holding their well occupancy and read count information.
+Will make a map of these keyed to the sequences themselves.
+Ideally, I'll be able to construct both the Vertices and the weights matrix from this map.
+
+ */
+
+import java.io.Serializable;
+import java.util.*;
+
+public class SequenceRecord implements Serializable {
+    private final String sequence;
+    private final SequenceType type;
+    //keys are well numbers, values are read count in that well
+    private final Map<Integer, Integer> wells;
+
+    public SequenceRecord (String sequence, SequenceType type) {
+        this.sequence = sequence;
+        this.type = type;
+        this.wells = new LinkedHashMap<>();
+    }
+
+    //this shouldn't be necessary, since the sequence will be the map key, but
+    public String getSequence() {
+        return sequence;
+    }
+
+    public SequenceType getSequenceType(){
+        return type;
+    }
+
+    //use this to update the record for each new read
+    public void addRead(Integer wellNumber) {
+        wells.merge(wellNumber,1, Integer::sum);
+    }
+
+    //don't know if I'll ever need this
+    public void addWellData(Integer wellNumber, Integer readCount) {
+        wells.put(wellNumber, readCount);
+    }
+
+    public Set<Integer> getWells() {
+        return wells.keySet();
+    }
+
+    public Map<Integer, Integer> getWellOccupancies() { return wells;}
+
+    public boolean isInWell(Integer wellNumber) {
+        return wells.containsKey(wellNumber);
+    }
+
+    public Integer getOccupancy() {
+        return wells.size();
+    }
+
+    //read count for whole plate
+    public Integer getReadCount(){
+        return wells.values().stream().mapToInt(Integer::valueOf).sum();
+    }
+
+    //read count in a specific well
+    public Integer getReadCount(Integer wellNumber) {
+        return wells.get(wellNumber);
+    }
+ }

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

@@ -12,107 +12,126 @@ import java.text.NumberFormat;
 import java.time.Instant;
 import java.time.Duration;
 import java.util.*;
-import java.util.stream.IntStream;
+/*
+ Refactor notes
+ What would be necessary to do everything with only one scan through the sample plate?
+ I would need to keep a list of sequences (real and spurious), and metadata about each sequence.
+ I would need the data:
+ * # of each well the sequence appears in
+ * Read count in that well
+ */
 
-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
-    public static GraphWithMapData makeGraph(CellSample cellSample, Plate samplePlate, boolean verbose) {
+
+    public static GraphWithMapData makeCDR3Graph(CellSample cellSample, Plate samplePlate, int readDepth,
+                                                 double readErrorRate, double errorCollisionRate, boolean verbose) {
+        //start timing
         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()};
+        List<String[]> distinctCells = cellSample.getCells();
         int numWells = samplePlate.getSize();
 
+        //Make a hashmap keyed to alphas, values are associated betas.
         if(verbose){System.out.println("Making cell maps");}
-        //HashMap keyed to Alphas, values Betas
-        Map<Integer, Integer> distCellsMapAlphaKey = makeSequenceToSequenceMap(distinctCells, 0, 1);
+        Map<String, String> distCellsMapAlphaKey = makeSequenceToSequenceMap(distinctCells,
+                SequenceType.CDR3_ALPHA.ordinal(), SequenceType.CDR3_BETA.ordinal());
         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);
-        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");}
+        //Make linkedHashMap keyed to sequences, values are SequenceRecords reflecting plate statistics
+        if(verbose){System.out.println("Making sample plate sequence maps");}
+        Map<String, SequenceRecord> alphaSequences = samplePlate.countSequences(readDepth, readErrorRate,
+                errorCollisionRate, alphaIndices);
+        int alphaCount = alphaSequences.size();
+        if(verbose){System.out.println("Alphas sequences read: " + alphaCount);}
+        Map<String, SequenceRecord> betaSequences = samplePlate.countSequences(readDepth, readErrorRate,
+                errorCollisionRate, betaIndices);
+        int betaCount = betaSequences.size();
+        if(verbose){System.out.println("Betas sequences read: " + betaCount);}
+        if(verbose){System.out.println("Sample plate sequence maps made");}
 
+        //pre-filter saturating sequences and sequences likely to be misreads
         if(verbose){System.out.println("Removing sequences present in all wells.");}
-        filterByOccupancyThresholds(allAlphas, 1, numWells - 1);
-        filterByOccupancyThresholds(allBetas, 1, numWells - 1);
+        filterByOccupancyThresholds(alphaSequences, 1, numWells - 1);
+        filterByOccupancyThresholds(betaSequences, 1, 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();
-        if(verbose){System.out.println("Remaining betas count: " + pairableBetaCount);}
+        if(verbose){System.out.println("Remaining alpha sequence count: " + alphaSequences.size());}
+        if(verbose){System.out.println("Remaining beta sequence count: " + betaSequences.size());}
+        if (readDepth > 1) {
+            if(verbose){System.out.println("Removing sequences with disparate occupancies and read counts");}
+            filterByOccupancyAndReadCount(alphaSequences, readDepth);
+            filterByOccupancyAndReadCount(betaSequences, readDepth);
+            if(verbose){System.out.println("Sequences removed");}
+            if(verbose){System.out.println("Remaining alpha sequence count: " + alphaSequences.size());}
+            if(verbose){System.out.println("Remaining beta sequence count: " + betaSequences.size());}
+        }
+        int pairableAlphaCount = alphaSequences.size();
+        if(verbose){System.out.println("Remaining alpha sequence count: " + pairableAlphaCount);}
+        int pairableBetaCount = betaSequences.size();
+        if(verbose){System.out.println("Remaining beta sequence count: " + pairableBetaCount);}
 
+        //construct the graph. For simplicity, going to make
         if(verbose){System.out.println("Making vertex maps");}
         //For the SimpleWeightedBipartiteGraphMatrixGenerator, all vertices must have
         //distinct numbers associated with them. Since I'm using a 2D array, that means
         //distinct indices between the rows and columns. vertexStartValue lets me track where I switch
         //from numbering rows to columns, so I can assign unique numbers to every vertex, and then
         //subtract the vertexStartValue from betas to use their vertex labels as array indices
-        Integer vertexStartValue = 0;
+        int vertexStartValue = 0;
         //keys are sequential integer vertices, values are alphas
-        Map<Integer, Integer> plateVtoAMap = makeVertexToSequenceMap(allAlphas, vertexStartValue);
+        Map<String, Integer> plateAtoVMap = makeSequenceToVertexMap(alphaSequences, vertexStartValue);
         //new start value for vertex to beta map should be one more than final vertex value in alpha map
-        vertexStartValue += plateVtoAMap.size();
-        //keys are sequential integers vertices, values are betas
-        Map<Integer, Integer> plateVtoBMap = makeVertexToSequenceMap(allBetas, vertexStartValue);
-        //keys are alphas, values are sequential integer vertices from previous map
-        Map<Integer, Integer> plateAtoVMap = invertVertexMap(plateVtoAMap);
-        //keys are betas, values are sequential integer vertices from previous map
-        Map<Integer, Integer> plateBtoVMap = invertVertexMap(plateVtoBMap);
+        vertexStartValue += plateAtoVMap.size();
+        //keys are betas, values are sequential integers
+        Map<String, Integer> plateBtoVMap = makeSequenceToVertexMap(betaSequences, vertexStartValue);
         if(verbose){System.out.println("Vertex maps made");}
-
         //make adjacency matrix for bipartite graph generator
         //(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
-        Map<Integer, Integer> alphaWellCounts = new HashMap<>();
-        //count how many wells each beta 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);
-        if(verbose){System.out.println("Matrix created");}
-
-        //create bipartite graph
-        if(verbose){System.out.println("Creating graph");}
+        if(verbose){System.out.println("Making adjacency matrix");}
+        double[][] weights = new double[plateAtoVMap.size()][plateBtoVMap.size()];
+        fillAdjacencyMatrix(weights, vertexStartValue, alphaSequences, betaSequences, plateAtoVMap, plateBtoVMap);
+        if(verbose){System.out.println("Adjacency matrix made");}
+        //make bipartite graph
+        if(verbose){System.out.println("Making bipartite weighted 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<Vertex> alphaVertices = new ArrayList<>();
+        for (String seq : plateAtoVMap.keySet()) {
+            Vertex alphaVertex = new Vertex(alphaSequences.get(seq), plateAtoVMap.get(seq));
+            alphaVertices.add(alphaVertex);
+        }
+        //Sort to make sure the order of vertices in list matches the order of the adjacency matrix
+        Collections.sort(alphaVertices);
+        //Add ordered list of vertices to the graph
         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<Vertex> betaVertices = new ArrayList<>();
+        for (String seq : plateBtoVMap.keySet()) {
+            Vertex betaVertex = new Vertex(betaSequences.get(seq), plateBtoVMap.get(seq));
+            betaVertices.add(betaVertex);
+        }
+        //Sort to make sure the order of vertices in list matches the order of the adjacency matrix
+        Collections.sort(betaVertices);
+        //Add ordered list of vertices to the graph
+        graphGenerator.second(betaVertices);
         //use adjacency matrix of weight created previously
         graphGenerator.weights(weights);
         graphGenerator.generateGraph(graph);
         if(verbose){System.out.println("Graph created");}
-
+        //stop timing
         Instant stop = Instant.now();
         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,
+                alphaCount, betaCount, readDepth, readErrorRate, errorCollisionRate, time);
         //Set source file name in graph to name of sample plate
         output.setSourceFilename(samplePlate.getFilename());
         //return GraphWithMapData object
@@ -124,61 +143,67 @@ 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();
-        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();
+        //Integer alphaCount = data.getAlphaCount();
+        //Integer betaCount = data.getBetaCount();
+        Map<String, String> distCellsMapAlphaKey = data.getDistCellsMapAlphaKey();
+        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 graphAlphaCount = alphas.size();
+        Integer graphBetaCount = 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
+        //Find Maximum Weight Matching
         //using jheaps library class PairingHeap for improved efficiency
-        if(verbose){System.out.println("Finding maximum weighted matching");}
+        if(verbose){System.out.println("Finding maximum weight matching");}
         MaximumWeightBipartiteMatching maxWeightMatching;
         //Use correct heap type for priority queue
         String heapType = BiGpairSEQ.getPriorityQueueHeapType();
         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
@@ -205,14 +230,14 @@ public class Simulator implements GraphModificationFunctions {
         int trueCount = 0;
         int falseCount = 0;
         boolean check;
-        Map<Integer, Integer> matchMap = new HashMap<>();
+        Map<String, String> 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);
             //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()));
             if(check) {
                 trueCount++;
             }
@@ -220,17 +245,19 @@ public class Simulator implements GraphModificationFunctions {
                 falseCount++;
             }
             List<String> result = new ArrayList<>();
-            result.add(plateVtoAMap.get(source).toString());
+            //alpha sequence
+            result.add(source.getSequence());
             //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());
             //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,18 +265,20 @@ public class Simulator implements GraphModificationFunctions {
 
         //Metadata comments for CSV file
         String algoType = "LEDA book with heap: " + heapType;
-        int min = Math.min(alphaCount, betaCount);
+        int min = Math.min(graphAlphaCount, graphBetaCount);
+        //matching weight
+        BigDecimal totalMatchingWeight = maxWeightMatching.getMatchingWeight();
         //rate of attempted matching
         double attemptRate = (double) (trueCount + falseCount) / min;
         BigDecimal attemptRateTrunc = new BigDecimal(attemptRate, mc);
         //rate of pairing error
         double pairingErrorRate = (double) falseCount / (trueCount + falseCount);
         BigDecimal pairingErrorRateTrunc;
-        if(pairingErrorRate == NaN || pairingErrorRate == POSITIVE_INFINITY || pairingErrorRate == NEGATIVE_INFINITY) {
-            pairingErrorRateTrunc = new BigDecimal(-1, mc);
+        if(Double.isFinite(pairingErrorRate)) {
+            pairingErrorRateTrunc = new BigDecimal(pairingErrorRate, mc);
         }
         else{
-            pairingErrorRateTrunc = new BigDecimal(pairingErrorRate, mc);
+            pairingErrorRateTrunc = new BigDecimal(-1, mc);
         }
         //get list of well populations
         Integer[] wellPopulations = data.getWellPopulations();
@@ -261,28 +290,44 @@ public class Simulator implements GraphModificationFunctions {
             populationsStringBuilder.append(wellPopulations[i].toString());
         }
         String wellPopulationsString = populationsStringBuilder.toString();
+        //graph generation time
+        Duration graphTime = data.getTime();
+        //MWM run time
+        Duration pairingTime = Duration.between(start, stop);
         //total simulation time
-        Duration time = Duration.between(start, stop);
-        time = time.plus(data.getTime());
+        Duration totalTime = graphTime.plus(pairingTime);
+
 
         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("MWM 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());
-        metadata.put("high overlap threshold", highThreshold.toString());
-        metadata.put("low overlap threshold", lowThreshold.toString());
-        metadata.put("minimum overlap percent", minOverlapPercent.toString());
-        metadata.put("maximum occupancy difference", maxOccupancyDifference.toString());
+        metadata.put("sequence read depth", data.getReadDepth().toString());
+        metadata.put("sequence read error rate", data.getReadErrorRate().toString());
+        metadata.put("read error collision rate", data.getErrorCollisionRate().toString());
+        metadata.put("total alphas read from plate", data.getAlphaCount().toString());
+        metadata.put("total betas read from plate", data.getBetaCount().toString());
+        //HARD CODED, PARAMETERIZE LATER
+        metadata.put("pre-filter sequences present in all wells", "true");
+        //HARD CODED, PARAMETERIZE LATER
+        metadata.put("pre-filter sequences based on occupancy/read count discrepancy", "true");
+        metadata.put("alphas in graph (after pre-filtering)", graphAlphaCount.toString());
+        metadata.put("betas in graph (after pre-filtering)", graphBetaCount.toString());
+        metadata.put("high overlap threshold for pairing", highThreshold.toString());
+        metadata.put("low overlap threshold for pairing", lowThreshold.toString());
+        metadata.put("minimum overlap percent for pairing", minOverlapPercent.toString());
+        metadata.put("maximum occupancy difference for pairing", maxOccupancyDifference.toString());
         metadata.put("pairing attempt rate", attemptRateTrunc.toString());
         metadata.put("correct pairing count", Integer.toString(trueCount));
         metadata.put("incorrect pairing count", Integer.toString(falseCount));
         metadata.put("pairing error rate", pairingErrorRateTrunc.toString());
-        metadata.put("simulation time (seconds)", nf.format(time.toSeconds()));
+        metadata.put("time to generate graph (seconds)", nf.format(graphTime.toSeconds()));
+        metadata.put("time to pair sequences (seconds)",nf.format(pairingTime.toSeconds()));
+        metadata.put("total simulation time (seconds)", nf.format(totalTime.toSeconds()));
         //create MatchingResult object
-        MatchingResult output = new MatchingResult(metadata, header, allResults, matchMap, time);
+        MatchingResult output = new MatchingResult(metadata, header, allResults, matchMap);
         if(verbose){
             for(String s: output.getComments()){
                 System.out.println(s);
@@ -604,81 +649,77 @@ public class Simulator implements GraphModificationFunctions {
 //    }
 
     //Remove sequences based on occupancy
-    public static void filterByOccupancyThresholds(Map<Integer, Integer> wellMap, int low, int high){
-        List<Integer> noise = new ArrayList<>();
-        for(Integer k: wellMap.keySet()){
-            if((wellMap.get(k) > high) || (wellMap.get(k) < low)){
+    public static void filterByOccupancyThresholds(Map<String, SequenceRecord> wellMap, int low, int high){
+        List<String> noise = new ArrayList<>();
+        for(String k: wellMap.keySet()){
+            if((wellMap.get(k).getOccupancy() > high) || (wellMap.get(k).getOccupancy() < low)){
                 noise.add(k);
             }
         }
-        for(Integer k: noise) {
+        for(String k: noise) {
             wellMap.remove(k);
         }
     }
 
-    //Counts the well occupancy of the row peptides and column peptides into given maps, and
-    //fills weights in the given 2D array
-    private static void countSequencesAndFillMatrix(Plate samplePlate,
-                                                    Map<Integer,Integer> allRowSequences,
-                                                    Map<Integer,Integer> allColumnSequences,
-                                                    Map<Integer,Integer> rowSequenceToVertexMap,
-                                                    Map<Integer,Integer> columnSequenceToVertexMap,
-                                                    int[] rowSequenceIndices,
-                                                    int[] colSequenceIndices,
-                                                    Map<Integer, Integer> rowSequenceCounts,
-                                                    Map<Integer,Integer> columnSequenceCounts,
-                                                    double[][] weights){
-        Map<Integer, Integer> wellNRowSequences = null;
-        Map<Integer, Integer> wellNColumnSequences = null;
-        int vertexStartValue = rowSequenceToVertexMap.size();
-        int numWells = samplePlate.getSize();
-        for (int n = 0; n < numWells; n++) {
-            wellNRowSequences = samplePlate.assayWellsSequenceS(n, rowSequenceIndices);
-            for (Integer a : wellNRowSequences.keySet()) {
-                if(allRowSequences.containsKey(a)){
-                    rowSequenceCounts.merge(a, 1, (oldValue, newValue) -> oldValue + newValue);
-                }
+    public static void filterByOccupancyAndReadCount(Map<String, SequenceRecord> sequences, int readDepth) {
+        List<String> noise = new ArrayList<>();
+        for(String k : sequences.keySet()){
+            //occupancy times read depth should be more than half the sequence read count if the read error rate is low
+            Integer threshold = (sequences.get(k).getOccupancy() * readDepth) / 2;
+            if(sequences.get(k).getReadCount() < threshold) {
+                noise.add(k);
             }
-            wellNColumnSequences = samplePlate.assayWellsSequenceS(n, colSequenceIndices);
-            for (Integer b : wellNColumnSequences.keySet()) {
-                if(allColumnSequences.containsKey(b)){
-                    columnSequenceCounts.merge(b, 1, (oldValue, newValue) -> oldValue + newValue);
-                }
-            }
-            for (Integer i : wellNRowSequences.keySet()) {
-                if(allRowSequences.containsKey(i)){
-                    for (Integer j : wellNColumnSequences.keySet()) {
-                        if(allColumnSequences.containsKey(j)){
-                            weights[rowSequenceToVertexMap.get(i)][columnSequenceToVertexMap.get(j) - vertexStartValue] += 1.0;
-                        }
-                    }
-                }
-            }
-
+        }
+        for(String k : noise) {
+            sequences.remove(k);
         }
     }
 
-    private static Map<Integer, Integer> makeSequenceToSequenceMap(List<Integer[]> cells, int keySequenceIndex,
-                                                                   int valueSequenceIndex){
-        Map<Integer, Integer> keySequenceToValueSequenceMap = new HashMap<>();
-        for (Integer[] cell : cells) {
+    private static Map<String, String> makeSequenceToSequenceMap(List<String[]> cells, int keySequenceIndex,
+                                                                 int valueSequenceIndex){
+        Map<String, String> keySequenceToValueSequenceMap = new HashMap<>();
+        for (String[] cell : cells) {
            keySequenceToValueSequenceMap.put(cell[keySequenceIndex], cell[valueSequenceIndex]);
         }
         return keySequenceToValueSequenceMap;
     }
 
-    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; //is this necessary? I don't think I use this.
-        for (Integer k: sequences.keySet()) {
+    private static Map<Integer, String> makeVertexToSequenceMap(Map<String, SequenceRecord> sequences, Integer startValue) {
+        Map<Integer, String> map = new LinkedHashMap<>(); //LinkedHashMap to preserve order of entry
+        Integer index = startValue;
+        for (String k: sequences.keySet()) {
             map.put(index, k);
             index++;
         }
         return map;
     }
 
-    private static Map<Integer, Integer> invertVertexMap(Map<Integer, Integer> map) {
-        Map<Integer, Integer> inverse = new HashMap<>();
+    private static Map<String, Integer> makeSequenceToVertexMap(Map<String, SequenceRecord> sequences, Integer startValue) {
+        Map<String, Integer> map = new LinkedHashMap<>(); //LinkedHashMap to preserve order of entry
+        Integer index = startValue;
+        for (String k: sequences.keySet()) {
+            map.put(k, index);
+            index++;
+        }
+        return map;
+    }
+
+    private static void fillAdjacencyMatrix(double[][] weights, Integer vertexOffsetValue, Map<String, SequenceRecord> rowSequences,
+                                            Map<String, SequenceRecord> columnSequences, Map<String, Integer> rowToVertexMap,
+                                            Map<String, Integer> columnToVertexMap) {
+        for (String rowSeq: rowSequences.keySet()) {
+            for (Integer well: rowSequences.get(rowSeq).getWells()) {
+                for (String colSeq: columnSequences.keySet()) {
+                    if (columnSequences.get(colSeq).isInWell(well)) {
+                        weights[rowToVertexMap.get(rowSeq)][columnToVertexMap.get(colSeq) - vertexOffsetValue] += 1.0;
+                    }
+                }
+            }
+        }
+    }
+
+    private static Map<String, Integer> invertVertexMap(Map<Integer, String> map) {
+        Map<String, Integer> inverse = new HashMap<>();
         for (Integer k : map.keySet()) {
             inverse.put(map.get(k), k);
         }

src/main/java/Vertex.java

@@ -1,23 +1,75 @@
+import org.jheaps.AddressableHeap;
 
+import java.io.Serializable;
+import java.util.Map;
 
-public class Vertex {
-    private final Integer vertexLabel;
-    private final Integer sequence;
-    private final Integer occupancy;
+public class Vertex implements Serializable, Comparable<Vertex> {
+    private SequenceRecord record;
+    private Integer vertexLabel;
+    private Double potential;
+    private AddressableHeap queue;
 
-    public Vertex(Integer vertexLabel, Integer sequence, Integer occupancy) {
+    public Vertex(SequenceRecord record, Integer vertexLabel) {
+        this.record = record;
         this.vertexLabel = vertexLabel;
-        this.sequence = sequence;
-        this.occupancy = occupancy;
     }
 
-    public Integer getVertexLabel() { return vertexLabel; }
+    public SequenceRecord getRecord() { return record; }
 
-    public Integer getSequence() {
-        return sequence;
+    public SequenceType getType() { return record.getSequenceType(); }
+
+    public Integer getVertexLabel() {
+        return vertexLabel;
+    }
+
+    public String getSequence() {
+        return record.getSequence();
     }
 
     public Integer getOccupancy() {
-        return occupancy;
+        return record.getOccupancy();
+    }
+
+    public Integer getReadCount() { return record.getReadCount(); }
+
+    public Map<Integer, Integer> getWellOccupancies() { return record.getWellOccupancies(); }
+
+    @Override //adapted from JGraphT example code
+    public int hashCode()
+    {
+        return (this.getSequence() == null) ? 0 : this.getSequence().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 (this.getSequence() == null) {
+            return other.getSequence() == null;
+        } else {
+            return this.getSequence().equals(other.getSequence());
+        }
+    }
+
+    @Override //adapted from JGraphT example code
+    public String toString()
+    {
+        StringBuilder sb = new StringBuilder();
+        sb.append("(").append(vertexLabel)
+                .append(", Type: ").append(this.getType().name())
+                .append(", Sequence: ").append(this.getSequence())
+                .append(", Occupancy: ").append(this.getOccupancy()).append(")");
+        return sb.toString();
+    }
+
+    @Override
+    public int compareTo(Vertex other) {
+        return this.vertexLabel - other.getVertexLabel();
     }
 }