From 9fcc20343d7354368740112d4f7bc9e21a7ce7a3 Mon Sep 17 00:00:00 2001 From: efischer Date: Sat, 26 Feb 2022 10:36:00 -0600 Subject: [PATCH 01/59] Fix GraphML writer --- src/main/java/GraphMLFileWriter.java | 57 +++++++++++++++++++++++-- src/main/java/InteractiveInterface.java | 2 +- 2 files changed, 54 insertions(+), 5 deletions(-) diff --git a/src/main/java/GraphMLFileWriter.java b/src/main/java/GraphMLFileWriter.java index ca5e89d..96ad358 100644 --- a/src/main/java/GraphMLFileWriter.java +++ b/src/main/java/GraphMLFileWriter.java @@ -1,4 +1,8 @@ +import org.jgrapht.graph.DefaultWeightedEdge; import org.jgrapht.graph.SimpleWeightedGraph; +import org.jgrapht.nio.Attribute; +import org.jgrapht.nio.AttributeType; +import org.jgrapht.nio.DefaultAttribute; import org.jgrapht.nio.dot.DOTExporter; import org.jgrapht.nio.graphml.GraphMLExporter; @@ -7,25 +11,69 @@ import java.io.IOException; import java.nio.file.Files; import java.nio.file.Path; import java.nio.file.StandardOpenOption; +import java.util.HashMap; +import java.util.LinkedHashMap; +import java.util.Map; public class GraphMLFileWriter { String filename; - SimpleWeightedGraph graph; + GraphWithMapData data; - public GraphMLFileWriter(String filename, SimpleWeightedGraph graph) { + public GraphMLFileWriter(String filename, GraphWithMapData data) { if(!filename.matches(".*\\.graphml")){ filename = filename + ".graphml"; } this.filename = filename; - this.graph = graph; + this.data = data; } +// public void writeGraphToFile() { +// try(BufferedWriter writer = Files.newBufferedWriter(Path.of(filename), StandardOpenOption.CREATE_NEW); +// ){ +// GraphMLExporter exporter = new GraphMLExporter<>(); +// exporter.exportGraph(graph, writer); +// } catch(IOException ex){ +// System.out.println("Could not make new file named "+filename); +// System.err.println(ex); +// } +// } + public void writeGraphToFile() { + SimpleWeightedGraph graph = data.getGraph(); + Map vertexToAlphaMap = data.getPlateVtoAMap(); + Map vertexToBetaMap = data.getPlateVtoBMap(); + Map alphaOccs = data.getAlphaWellCounts(); + Map betaOccs = data.getBetaWellCounts(); try(BufferedWriter writer = Files.newBufferedWriter(Path.of(filename), StandardOpenOption.CREATE_NEW); ){ - GraphMLExporter exporter = new GraphMLExporter<>(); + //create exporter. Let the vertex labels be the unique ids for the vertices + GraphMLExporter> exporter = new GraphMLExporter<>(v -> v.toString()); + //set to export weights + exporter.setExportEdgeWeights(true); + //set type, sequence, and occupancy attributes for each vertex + exporter.setVertexAttributeProvider( v -> { + Map attributes = new HashMap<>(); + if(vertexToAlphaMap.containsKey(v)) { + attributes.put("type", DefaultAttribute.createAttribute("CDR3 Alpha")); + attributes.put("sequence", DefaultAttribute.createAttribute(vertexToAlphaMap.get(v))); + attributes.put("occupancy", DefaultAttribute.createAttribute( + alphaOccs.get(vertexToAlphaMap.get(v)))); + } + else if(vertexToBetaMap.containsKey(v)) { + attributes.put("type", DefaultAttribute.createAttribute("CDR3 Beta")); + attributes.put("sequence", DefaultAttribute.createAttribute(vertexToBetaMap.get(v))); + attributes.put("occupancy", DefaultAttribute.createAttribute( + betaOccs.get(vertexToBetaMap.get(v)))); + } + return attributes; + }); + //register the attributes + exporter.registerAttribute("type", GraphMLExporter.AttributeCategory.NODE, AttributeType.STRING); + exporter.registerAttribute("sequence", GraphMLExporter.AttributeCategory.NODE, AttributeType.STRING); + exporter.registerAttribute("occupancy", GraphMLExporter.AttributeCategory.NODE, AttributeType.STRING); + //export the graph exporter.exportGraph(graph, writer); } catch(IOException ex){ System.out.println("Could not make new file named "+filename); @@ -33,3 +81,4 @@ public class GraphMLFileWriter { } } } + diff --git a/src/main/java/InteractiveInterface.java b/src/main/java/InteractiveInterface.java index 044e03b..9293cfd 100644 --- a/src/main/java/InteractiveInterface.java +++ b/src/main/java/InteractiveInterface.java @@ -315,7 +315,7 @@ public class InteractiveInterface { System.out.println("Serialized binary graph/data file written to: " + filename); } if(BiGpairSEQ.outputGraphML()) { - GraphMLFileWriter graphMLWriter = new GraphMLFileWriter(filename, data.getGraph()); + GraphMLFileWriter graphMLWriter = new GraphMLFileWriter(filename, data); graphMLWriter.writeGraphToFile(); System.out.println("GraphML file written to: " + filename); } From 8935407ade89af8facc83bfb5dec437df16380f4 Mon Sep 17 00:00:00 2001 From: efischer Date: Sat, 26 Feb 2022 10:38:10 -0600 Subject: [PATCH 02/59] Get rid of GraphML reader, those files are larger than serialized files --- src/main/java/GraphMLFileReader.java | 35 ---------------------------- 1 file changed, 35 deletions(-) delete mode 100644 src/main/java/GraphMLFileReader.java diff --git a/src/main/java/GraphMLFileReader.java b/src/main/java/GraphMLFileReader.java deleted file mode 100644 index 151c078..0000000 --- a/src/main/java/GraphMLFileReader.java +++ /dev/null @@ -1,35 +0,0 @@ -import org.jgrapht.graph.SimpleWeightedGraph; -import org.jgrapht.nio.graphml.GraphMLImporter; - -import java.io.BufferedReader; -import java.io.IOException; -import java.nio.file.Files; -import java.nio.file.Path; - -public class GraphMLFileReader { - - private String filename; - private SimpleWeightedGraph graph; - - public GraphMLFileReader(String filename, SimpleWeightedGraph graph) { - if(!filename.matches(".*\\.graphml")){ - filename = filename + ".graphml"; - } - this.filename = filename; - this.graph = graph; - - try(//don't need to close reader bc of try-with-resources auto-closing - BufferedReader reader = Files.newBufferedReader(Path.of(filename)); - ){ - GraphMLImporter importer = new GraphMLImporter<>(); - importer.importGraph(graph, reader); - } - catch (IOException ex) { - System.out.println("Graph file " + filename + " not found."); - System.err.println(ex); - } - } - - public SimpleWeightedGraph getGraph() { return graph; } - -} From c2db4f87c1999c99d9f89c046a5a864624b14dd5 Mon Sep 17 00:00:00 2001 From: efischer Date: Sat, 26 Feb 2022 11:00:18 -0600 Subject: [PATCH 03/59] Update Readme --- readme.md | 77 +++++++++++++++++++++++++++++++++++-------------------- 1 file changed, 49 insertions(+), 28 deletions(-) diff --git a/readme.md b/readme.md index fa8e894..7eee6c3 100644 --- a/readme.md +++ b/readme.md @@ -66,6 +66,18 @@ Please select an option: 0) Exit ``` +By default, the Options menu looks like this: +``` +--------------OPTIONS--------------- +1) Turn on cell sample file caching +2) Turn on plate file caching +3) Turn on graph/data file caching +4) Turn off serialized binary graph output +5) Turn on GraphML graph output +6) Maximum weight matching algorithm options +0) Return to main menu +``` + ### INPUT/OUTPUT To run the simulation, the program reads and writes 4 kinds of files: @@ -75,21 +87,24 @@ To run the simulation, the program reads and writes 4 kinds of files: * Matching Results files in CSV format These files are often generated in sequence. When entering filenames, it is not necessary to include the file extension -(.csv or .ser). When reading or writing files, the program will automatically add the correct extension to any filename without one. +(.csv or .ser). When reading or writing files, the program will automatically add the correct extension to any filename +without one. To save file I/O time, the most recent instance of each of these four -files either generated or read from disk can be cached in program memory. This is could be important for Graph/Data files, -which can be several gigabytes in size. Since some simulations may require running multiple, -differently-configured BiGpairSEQ matchings on the same graph, keeping the most recent graph cached may reduce execution time. -(The manipulation necessary to re-use a graph incurs its own performance overhead, though, which may scale with graph -size faster than file I/O does. If so, caching is best for smaller graphs.) - -When caching is active, subsequent uses of the same data file won't need to be read in again until another file of that type is used or generated, +files either generated or read from disk can be cached in program memory. When caching is active, subsequent uses of the +same data file won't need to be read in again until another file of that type is used or generated, or caching is turned off for that file type. The program checks whether it needs to update its cached data by comparing filenames as entered by the user. On encountering a new filename, the program flushes its cache and reads in the new file. +(Note that cached Graph/Data files must be transformed back into their original state after a matching experiment, which +may take some time. Whether file I/O or graph transformation takes longer for graph/data files is likely to be +device-specific.) + The program's caching behavior can be controlled in the Options menu. By default, all caching is OFF. +The program can optionally output Graph/Data files in .GraphML format (.graphml) for data portability. This can be +turned on in the Options menu. By default, GraphML output is OFF. + #### Cell Sample Files Cell Sample files consist of any number of distinct "T cells." Every cell contains four sequences: Alpha CDR3, Beta CDR3, Alpha CDR1, Beta CDR1. The sequences are represented by @@ -181,14 +196,19 @@ Options for creating a Graph/Data file: * The Cell Sample file to use * The Sample Plate file to use. (This must have been generated from the selected Cell Sample file.) -These files do not have a human-readable structure, and are not portable to other programs. (Export of graphs in a -portable data format may be implemented in the future. The tricky part is encoding the necessary metadata.) +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.) --- #### Matching Results Files -Matching results files consist of the results of a BiGpairSEQ matching simulation. Making them requires a Graph and -Data file. Matching results files are in CSV format. Rows are sequence pairings with extra relevant data. Columns are pairing-specific details. +Matching results files consist of the results of a BiGpairSEQ matching simulation. Making them requires a serialized +binary Graph/Data file (.ser). (Because .graphML files are larger than .ser files, BiGpairSEQ_Sim supports .graphML +output only. Graph/data input must use a serialized binary.) + +Matching results files are in CSV format. Rows are sequence pairings with extra relevant data. Columns are pairing-specific details. Metadata about the matching simulation is included as comments. Comments are preceded by `#`. Options when running a BiGpairSEQ simulation of CDR3 alpha/beta matching: @@ -258,25 +278,26 @@ slightly less time than the simulation itself. Real elapsed time from start to f * ~~*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. -* 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._ -* Re-implement command line arguments, to enable scripting and statistical simulation studies -* ~~Implement sample plates with random numbers of T cells per well.~~ DONE - * Possible BiGpairSEQ advantage over pairSEQ: BiGpairSEQ is resilient to variations in well population sizes on a sample plate; pairSEQ is not. - * preliminary data suggests that BiGpairSEQ behaves roughly as though the whole plate had whatever the *average* well concentration is, but that's still speculative. -* Enable GraphML output in addition to serialized object binaries, for data portability - * Custom vertex type with attribute for sequence occupancy? -* Re-implement CDR1 matching method -* Implement Duan and Su's maximum weight matching algorithm - * Add controllable algorithm-type parameter? * ~~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. For large graphs, Fibonacci heap wins. Now the new default. - - + * Parameter implemented. Fibonacci heap the current default. +* ~~Implement sample plates with random numbers of T cells per well.~~ DONE + * Possible BiGpairSEQ advantage over pairSEQ: BiGpairSEQ is resilient to variations in well population sizes on a sample plate; pairSEQ is not. + * preliminary data suggests that BiGpairSEQ behaves roughly as though the whole plate had whatever the *average* well concentration is, but that's still speculative. +* See if there's a reasonable way to reformat Sample Plate files so that wells are columns instead of rows. + * ~~Problem is variable number of cells in a well~~ + * ~~Apache Commons CSV library writes entries a row at a time~~ + * _Got this working, but at the cost of a profoundly strange bug in graph occupancy filtering. Have reverted the repo until I can figure out what caused that. Given how easily Thingiverse transposes CSV matrices in R, might not even be worth fixing. +* ~~Enable GraphML output in addition to serialized object binaries, for data portability~~ DONE + * ~~Custom vertex type with attribute for sequence occupancy?~~ABANDONED + * Have a branch where this is implemented, but there's a bug that broke matching. Don't currently have time to fix. +* Re-implement command line arguments, to enable scripting and statistical simulation studies +* Re-implement CDR1 matching method +* Implement Duan and Su's maximum weight matching algorithm + * Add controllable algorithm-type parameter? + * This would be fun and valuable, but probably take more time than I have for a hobby project. + ## CITATIONS * Howie, B., Sherwood, A. M., et al. ["High-throughput pairing of T cell receptor alpha and beta sequences."](https://pubmed.ncbi.nlm.nih.gov/26290413/) Sci. Transl. Med. 7, 301ra131 (2015) From 98bf452891dfe22fe6dd093194f238d022005555 Mon Sep 17 00:00:00 2001 From: efischer Date: Sat, 26 Feb 2022 11:01:20 -0600 Subject: [PATCH 04/59] Update Readme --- readme.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/readme.md b/readme.md index 7eee6c3..5b05c33 100644 --- a/readme.md +++ b/readme.md @@ -104,7 +104,7 @@ The program's caching behavior can be controlled in the Options menu. By default The program can optionally output Graph/Data files in .GraphML format (.graphml) for data portability. This can be turned on in the Options menu. By default, GraphML output is OFF. - +--- #### Cell Sample Files Cell Sample files consist of any number of distinct "T cells." Every cell contains four sequences: Alpha CDR3, Beta CDR3, Alpha CDR1, Beta CDR1. The sequences are represented by From 88eb8aca50be7cb8945efe40d0371f22f2f826c5 Mon Sep 17 00:00:00 2001 From: efischer Date: Sat, 26 Feb 2022 11:01:44 -0600 Subject: [PATCH 05/59] Update Readme --- readme.md | 1 + 1 file changed, 1 insertion(+) diff --git a/readme.md b/readme.md index 5b05c33..f33dfb5 100644 --- a/readme.md +++ b/readme.md @@ -104,6 +104,7 @@ The program's caching behavior can be controlled in the Options menu. By default The program can optionally output Graph/Data files in .GraphML format (.graphml) for data portability. This can be turned on in the Options menu. By default, GraphML output is OFF. + --- #### Cell Sample Files Cell Sample files consist of any number of distinct "T cells." Every cell contains From ff72c9b359ddaee371449db86856910853e2c848 Mon Sep 17 00:00:00 2001 From: efischer Date: Sat, 26 Feb 2022 11:02:23 -0600 Subject: [PATCH 06/59] Update Readme --- readme.md | 3 --- 1 file changed, 3 deletions(-) diff --git a/readme.md b/readme.md index f33dfb5..b7726ed 100644 --- a/readme.md +++ b/readme.md @@ -123,7 +123,6 @@ Comments are preceded by `#` Structure: ---- # Sample contains 1 unique CDR1 for every 4 unique CDR3s. | Alpha CDR3 | Beta CDR3 | Alpha CDR1 | Beta CDR1 | |---|---|---|---| @@ -168,7 +167,6 @@ Dropout sequences are replaced with the value `-1`. Comments are preceded by `#` Structure: ---- ``` # Cell source file name: # Each row represents one well on the plate @@ -224,7 +222,6 @@ Options when running a BiGpairSEQ simulation of CDR3 alpha/beta matching: Example output: ---- ``` # Source Sample Plate file: 4MilCellsPlate.csv # Source Graph and Data file: 4MilCellsPlateGraph.ser From b4cc240048c3bacd2abb9ad582c3368755f7ac68 Mon Sep 17 00:00:00 2001 From: efischer Date: Sat, 26 Feb 2022 11:03:31 -0600 Subject: [PATCH 07/59] Update Readme --- readme.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/readme.md b/readme.md index b7726ed..2ba6163 100644 --- a/readme.md +++ b/readme.md @@ -288,7 +288,7 @@ slightly less time than the simulation itself. Real elapsed time from start to f * ~~Apache Commons CSV library writes entries a row at a time~~ * _Got this working, but at the cost of a profoundly strange bug in graph occupancy filtering. Have reverted the repo until I can figure out what caused that. Given how easily Thingiverse transposes CSV matrices in R, might not even be worth fixing. * ~~Enable GraphML output in addition to serialized object binaries, for data portability~~ DONE - * ~~Custom vertex type with attribute for sequence occupancy?~~ABANDONED + * ~~Custom vertex type with attribute for sequence occupancy?~~ ABANDONED * Have a branch where this is implemented, but there's a bug that broke matching. Don't currently have time to fix. * Re-implement command line arguments, to enable scripting and statistical simulation studies * Re-implement CDR1 matching method From 6b5837e6ce4f99f7e8d06c4a3723f433787696c7 Mon Sep 17 00:00:00 2001 From: efischer Date: Sun, 27 Feb 2022 11:46:11 -0600 Subject: [PATCH 08/59] Add Vose's alias method to to-dos --- readme.md | 1 + 1 file changed, 1 insertion(+) diff --git a/readme.md b/readme.md index 2ba6163..6a35f81 100644 --- a/readme.md +++ b/readme.md @@ -295,6 +295,7 @@ slightly less time than the simulation itself. Real elapsed time from start to f * 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 Vose's alias method for arbitrary statistical distributions of cells ## CITATIONS From 974d2d650cf53d2e3bb9909a6d8f08ce736ac463 Mon Sep 17 00:00:00 2001 From: efischer Date: Sun, 27 Feb 2022 14:17:11 -0600 Subject: [PATCH 09/59] Refactor plate to fill its own wells in its constructor --- src/main/java/Plate.java | 25 +++++++++++++++++++++---- 1 file changed, 21 insertions(+), 4 deletions(-) diff --git a/src/main/java/Plate.java b/src/main/java/Plate.java index 5a13eff..af37ac0 100644 --- a/src/main/java/Plate.java +++ b/src/main/java/Plate.java @@ -8,6 +8,7 @@ TODO: Implement discrete frequency distributions using Vose's Alias Method import java.util.*; public class Plate { + private CellSample cells; private String sourceFile; private List> wells; private final Random rand = BiGpairSEQ.getRand(); @@ -18,6 +19,24 @@ public class Plate { private double lambda; boolean exponential = false; + public Plate(CellSample cells, String cellFilename, int size, Integer[] populations, + double dropoutRate, double stdDev_or_lambda, boolean exponential){ + this.cells = cells; + this.sourceFile = cellFilename; + this.size = size; + this.error = dropoutRate; + this.populations = populations; + this.exponential = exponential; + if (this.exponential) { + this.lambda = stdDev_or_lambda; + fillWellsExponential(cells.getCells(), this.lambda); + } + else { + this.stdDev = stdDev_or_lambda; + fillWells(cells.getCells(), this.stdDev); + } + } + public Plate(int size, double error, Integer[] populations) { this.size = size; @@ -43,10 +62,9 @@ public class Plate { } } - public void fillWellsExponential(String sourceFileName, List cells, double lambda){ + private void fillWellsExponential(List cells, double lambda){ this.lambda = lambda; exponential = true; - sourceFile = sourceFileName; int numSections = populations.length; int section = 0; double m; @@ -74,9 +92,8 @@ public class Plate { } } - public void fillWells(String sourceFileName, List cells, double stdDev) { + private void fillWells( List cells, double stdDev) { this.stdDev = stdDev; - sourceFile = sourceFileName; int numSections = populations.length; int section = 0; double m; From b53f5f1cc0fae99823e64062c5dc04f15afc8bd3 Mon Sep 17 00:00:00 2001 From: efischer Date: Sun, 27 Feb 2022 14:17:16 -0600 Subject: [PATCH 10/59] Refactor plate to fill its own wells in its constructor --- src/main/java/InteractiveInterface.java | 6 ++---- 1 file changed, 2 insertions(+), 4 deletions(-) diff --git a/src/main/java/InteractiveInterface.java b/src/main/java/InteractiveInterface.java index 9293cfd..b4bf965 100644 --- a/src/main/java/InteractiveInterface.java +++ b/src/main/java/InteractiveInterface.java @@ -227,16 +227,14 @@ public class InteractiveInterface { Plate samplePlate; PlateFileWriter writer; if(exponential){ - samplePlate = new Plate(numWells, dropOutRate, populations); - samplePlate.fillWellsExponential(cellFile, cells.getCells(), lambda); + samplePlate = new Plate(cells, cellFile, numWells, populations, dropOutRate, lambda, true); writer = new PlateFileWriter(filename, samplePlate); } else { if (poisson) { stdDev = Math.sqrt(cells.getCellCount()); //gaussian with square root of elements approximates poisson } - samplePlate = new Plate(numWells, dropOutRate, populations); - samplePlate.fillWells(cellFile, cells.getCells(), stdDev); + samplePlate = new Plate(cells, cellFile, numWells, populations, dropOutRate, stdDev, false); writer = new PlateFileWriter(filename, samplePlate); } System.out.println("Writing Sample Plate to file"); From 8ebfc1469f05c7c080b0abeab9b06185f97cd656 Mon Sep 17 00:00:00 2001 From: efischer Date: Sun, 27 Feb 2022 14:25:53 -0600 Subject: [PATCH 11/59] Refactor plate to fill its own wells in its constructor --- src/main/java/InteractiveInterface.java | 2 +- src/main/java/Plate.java | 8 ++++++-- src/main/java/PlateFileReader.java | 7 ++----- src/main/java/Simulator.java | 2 +- 4 files changed, 10 insertions(+), 9 deletions(-) diff --git a/src/main/java/InteractiveInterface.java b/src/main/java/InteractiveInterface.java index b4bf965..0de7a8e 100644 --- a/src/main/java/InteractiveInterface.java +++ b/src/main/java/InteractiveInterface.java @@ -290,7 +290,7 @@ public class InteractiveInterface { else { System.out.println("Reading Sample Plate file: " + plateFile); PlateFileReader plateReader = new PlateFileReader(plateFile); - plate = new Plate(plateReader.getFilename(), plateReader.getWells()); + plate = plateReader.getSamplePlate(); if(BiGpairSEQ.cachePlate()) { BiGpairSEQ.setPlateInMemory(plate, plateFile); } diff --git a/src/main/java/Plate.java b/src/main/java/Plate.java index af37ac0..d0bc1aa 100644 --- a/src/main/java/Plate.java +++ b/src/main/java/Plate.java @@ -10,6 +10,7 @@ import java.util.*; public class Plate { private CellSample cells; private String sourceFile; + private String filename; private List> wells; private final Random rand = BiGpairSEQ.getRand(); private int size; @@ -45,8 +46,9 @@ public class Plate { wells = new ArrayList<>(); } - public Plate(String sourceFileName, List> wells) { - this.sourceFile = sourceFileName; + //constructor for returning a Plate from a PlateFileReader + public Plate(String filename, List> wells) { + this.filename = filename; this.wells = wells; this.size = wells.size(); @@ -176,4 +178,6 @@ public class Plate { public String getSourceFileName() { return sourceFile; } + + public String getFilename() { return filename; } } diff --git a/src/main/java/PlateFileReader.java b/src/main/java/PlateFileReader.java index 7b95cee..27e98b0 100644 --- a/src/main/java/PlateFileReader.java +++ b/src/main/java/PlateFileReader.java @@ -56,11 +56,8 @@ public class PlateFileReader { } - public List> getWells() { - return wells; + public Plate getSamplePlate() { + return new Plate(filename, wells); } - public String getFilename() { - return filename; - } } \ No newline at end of file diff --git a/src/main/java/Simulator.java b/src/main/java/Simulator.java index 6d3f6eb..0cc6df1 100644 --- a/src/main/java/Simulator.java +++ b/src/main/java/Simulator.java @@ -113,7 +113,7 @@ public class Simulator implements GraphModificationFunctions { distCellsMapAlphaKey, plateVtoAMap, plateVtoBMap, plateAtoVMap, plateBtoVMap, alphaWellCounts, betaWellCounts, time); //Set source file name in graph to name of sample plate - output.setSourceFilename(samplePlate.getSourceFileName()); + output.setSourceFilename(samplePlate.getFilename()); //return GraphWithMapData object return output; } From a5f7c0641d30a368f6527e13fa7b4ccec91561b9 Mon Sep 17 00:00:00 2001 From: efischer Date: Sun, 27 Feb 2022 14:51:53 -0600 Subject: [PATCH 12/59] Refactor for better encapsulation with CellSamples --- src/main/java/InteractiveInterface.java | 3 +-- src/main/java/Simulator.java | 3 ++- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/src/main/java/InteractiveInterface.java b/src/main/java/InteractiveInterface.java index 0de7a8e..d826b1f 100644 --- a/src/main/java/InteractiveInterface.java +++ b/src/main/java/InteractiveInterface.java @@ -304,8 +304,7 @@ public class InteractiveInterface { System.out.println("Returning to main menu."); } else{ - List cells = cellSample.getCells(); - GraphWithMapData data = Simulator.makeGraph(cells, plate, true); + GraphWithMapData data = Simulator.makeGraph(cellSample, plate, true); assert filename != null; if(BiGpairSEQ.outputBinary()) { GraphDataObjectWriter dataWriter = new GraphDataObjectWriter(filename, data); diff --git a/src/main/java/Simulator.java b/src/main/java/Simulator.java index 0cc6df1..da13878 100644 --- a/src/main/java/Simulator.java +++ b/src/main/java/Simulator.java @@ -24,8 +24,9 @@ public class Simulator implements GraphModificationFunctions { private static final int cdr1BetaIndex = 3; //Make the graph needed for matching CDR3s - public static GraphWithMapData makeGraph(List distinctCells, Plate samplePlate, boolean verbose) { + public static GraphWithMapData makeGraph(CellSample cellSample, Plate samplePlate, boolean verbose) { Instant start = Instant.now(); + List distinctCells = cellSample.getCells(); int[] alphaIndex = {cdr3AlphaIndex}; int[] betaIndex = {cdr3BetaIndex}; From 5246cc4a0ccf9ba41acb6b0bff70edf47b5243e7 Mon Sep 17 00:00:00 2001 From: efischer Date: Sun, 27 Feb 2022 15:35:07 -0600 Subject: [PATCH 13/59] Re-implement command line options --- src/main/java/BiGpairSEQ.java | 4 +- src/main/java/CommandLineInterface.java | 526 ++++++++++++++---------- src/main/java/Plate.java | 4 +- 3 files changed, 307 insertions(+), 227 deletions(-) diff --git a/src/main/java/BiGpairSEQ.java b/src/main/java/BiGpairSEQ.java index 935888c..3e67a68 100644 --- a/src/main/java/BiGpairSEQ.java +++ b/src/main/java/BiGpairSEQ.java @@ -23,8 +23,8 @@ public class BiGpairSEQ { } else { //This will be uncommented when command line arguments are re-implemented. - //CommandLineInterface.startCLI(args); - System.out.println("Command line arguments are still being re-implemented."); + CommandLineInterface.startCLI(args); + //System.out.println("Command line arguments are still being re-implemented."); } } diff --git a/src/main/java/CommandLineInterface.java b/src/main/java/CommandLineInterface.java index 0c527b1..aadfb3a 100644 --- a/src/main/java/CommandLineInterface.java +++ b/src/main/java/CommandLineInterface.java @@ -1,5 +1,8 @@ import org.apache.commons.cli.*; +import java.io.IOException; +import java.util.stream.Stream; + /* * Class for parsing options passed to program from command line * @@ -29,6 +32,8 @@ import org.apache.commons.cli.*; * cellfile : name of the cell sample file to use as input * platefile : name of the sample plate file to use as input * output : name of the output file + * graphml : output a graphml file + * binary : output a serialized binary object file * * Match flags: * graphFile : name of graph and data file to use as input @@ -49,23 +54,156 @@ public class CommandLineInterface { //should be three separate programs. //There might be a way to do it with option parameters? - //main options set + //main options set - for the four different program modes + Options mainOptions = buildMainOptions(); + CommandLineParser parser = new DefaultParser(); + try{ + CommandLine line = parser.parse(mainOptions, args); + + if (line.hasOption("cells")) { + Options cellOptions = buildCellOptions(); + line = parser.parse(cellOptions, args); + Integer number = Integer.valueOf(line.getOptionValue("n")); + Integer diversity = Integer.valueOf(line.getOptionValue("d")); + String filename = line.getOptionValue("o"); + makeCells(filename, number, diversity); + } + + else if (line.hasOption("plate")) { + Options plateOptions = buildPlateOptions(); + line = parser.parse(plateOptions, args); + //get the cells + String cellFilename = line.getOptionValue("c"); + CellSample cells = getCells(cellFilename); + //get the rest of the parameters + Integer[] populations; + String outputFilename = line.getOptionValue("o"); + Integer numWells = Integer.parseInt(line.getOptionValue("w")); + Double dropoutRate = Double.parseDouble(line.getOptionValue("err")); + if (line.hasOption("random")) { + //Array holding values of minimum and maximum populations + Integer[] min_max = Stream.of(line.getOptionValues("random")) + .mapToInt(Integer::parseInt) + .boxed() + .toArray(Integer[]::new); + populations = BiGpairSEQ.getRand().ints(min_max[0], min_max[1] + 1) + .limit(numWells) + .boxed() + .toArray(Integer[]::new); + } + else if (line.hasOption("pop")) { + populations = Stream.of(line.getOptionValues("pop")) + .mapToInt(Integer::parseInt) + .boxed() + .toArray(Integer[]::new); + } + else{ + populations = new Integer[1]; + populations[0] = 1; + } + //make the plate + Plate plate; + if (line.hasOption("poisson")) { + Double stdDev = Math.sqrt(numWells); + plate = new Plate(cells, cellFilename, numWells, populations, dropoutRate, stdDev, false); + } + else if (line.hasOption("gaussian")) { + Double stdDev = Double.parseDouble(line.getOptionValue("stddev")); + plate = new Plate(cells, cellFilename, numWells, populations, dropoutRate, stdDev, false); + } + else { + assert line.hasOption("exponential"); + Double lambda = Double.parseDouble(line.getOptionValue("lambda")); + plate = new Plate(cells, cellFilename, numWells, populations, dropoutRate, lambda, true); + } + PlateFileWriter writer = new PlateFileWriter(outputFilename, plate); + writer.writePlateFile(); + } + + else if (line.hasOption("graph")) { //Making a graph + Options graphOptions = buildGraphOptions(); + line = parser.parse(graphOptions, args); + String cellFilename = line.getOptionValue("c"); + String plateFilename = line.getOptionValue("p"); + String outputFilename = line.getOptionValue("o"); + //get cells + CellSample cells = getCells(cellFilename); + //get plate + Plate plate = getPlate(plateFilename); + GraphWithMapData graph = Simulator.makeGraph(cells, plate, false); + if (!line.hasOption("no-binary")) { //output binary file unless told not to + GraphDataObjectWriter writer = new GraphDataObjectWriter(outputFilename, graph); + writer.writeDataToFile(); + } + if (line.hasOption("graphml")) { //if told to, output graphml file + GraphMLFileWriter gmlwriter = new GraphMLFileWriter(outputFilename, graph); + gmlwriter.writeGraphToFile(); + } + } + + else if (line.hasOption("match")) { //can add a flag for which match type in future, spit this in two + Options matchOptions = buildMatchCDR3options(); + line = parser.parse(matchOptions, args); + String graphFilename = line.getOptionValue("g"); + String outputFilename = line.getOptionValue("o"); + Integer minThreshold = Integer.parseInt(line.getOptionValue("min")); + Integer maxThreshold = Integer.parseInt(line.getOptionValue("max")); + Integer minOverlapPct; + if (line.hasOption("minpct")) { //see if this filter is being used + minOverlapPct = Integer.parseInt(line.getOptionValue("minpct")); + } + else { + minOverlapPct = 0; + } + Integer maxOccupancyDiff; + if (line.hasOption("maxdiff")) { //see if this filter is being used + maxOccupancyDiff = Integer.parseInt(line.getOptionValue("maxdiff")); + } + else { + maxOccupancyDiff = Integer.MAX_VALUE; + } + GraphWithMapData graph = getGraph(graphFilename); + MatchingResult result = Simulator.matchCDR3s(graph, graphFilename, minThreshold, maxThreshold, + maxOccupancyDiff, minOverlapPct, false); + 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 + } + } + catch (ParseException exp) { + System.err.println("Parsing failed. Reason: " + exp.getMessage()); + } + } + + private static Option outputFileOption() { + Option outputFile = Option.builder("o") + .longOpt("output-file") + .hasArg() + .argName("filename") + .desc("Name of output file") + .required() + .build(); + return outputFile; + } + + private static Options buildMainOptions() { Options mainOptions = new Options(); Option makeCells = Option.builder("cells") .longOpt("make-cells") - .desc("Makes a file of distinct cells") + .desc("Makes a cell sample file of distinct T cells") .build(); Option makePlate = Option.builder("plates") .longOpt("make-plates") - .desc("Makes a sample plate file") + .desc("Makes a sample plate file. Requires a cell sample file.") .build(); Option makeGraph = Option.builder("graph") .longOpt("make-graph") - .desc("Makes a graph and data file") + .desc("Makes a graph/data file. Requires a cell sample file and a sample plate file") .build(); Option matchCDR3 = Option.builder("match") .longOpt("match-cdr3") - .desc("Match CDR3s. Requires a cell sample file and any number of plate files.") + .desc("Matches CDR3s. Requires a graph/data file.") .build(); OptionGroup mainGroup = new OptionGroup(); mainGroup.addOption(makeCells); @@ -74,255 +212,197 @@ public class CommandLineInterface { mainGroup.addOption(matchCDR3); mainGroup.setRequired(true); mainOptions.addOptionGroup(mainGroup); + return mainOptions; + } - //Reuse clones of this for other options groups, rather than making it lots of times - Option outputFile = Option.builder("o") - .longOpt("output-file") - .hasArg() - .argName("filename") - .desc("Name of output file") - .build(); - mainOptions.addOption(outputFile); - - //Options cellOptions = new Options(); - Option numCells = Option.builder("nc") + private static Options buildCellOptions() { + Options cellOptions = new Options(); + Option numCells = Option.builder("n") .longOpt("num-cells") .desc("The number of distinct cells to generate") .hasArg() .argName("number") - .build(); - mainOptions.addOption(numCells); - Option cdr1Freq = Option.builder("d") - .longOpt("peptide-diversity-factor") + .required().build(); + Option cdr3Diversity = Option.builder("d") + .longOpt("diversity-factor") + .desc("The factor by which unique CDR3s outnumber unique CDR1s") .hasArg() - .argName("number") - .desc("Number of distinct CDR3s for every CDR1") - .build(); - mainOptions.addOption(cdr1Freq); - //Option cellOutput = (Option) outputFile.clone(); - //cellOutput.setRequired(true); - //mainOptions.addOption(cellOutput); + .argName("factor") + .required().build(); + cellOptions.addOption(numCells); + cellOptions.addOption(cdr3Diversity); + cellOptions.addOption(outputFileOption()); + return cellOptions; + } - //Options plateOptions = new Options(); - Option inputCells = Option.builder("c") + private static Options buildPlateOptions() { + Options plateOptions = new Options(); + Option cellFile = Option.builder("c") // add this to plate options .longOpt("cell-file") + .desc("The cell sample file to use") .hasArg() - .argName("file") - .desc("The cell sample file used for filling wells") - .build(); - mainOptions.addOption(inputCells); - Option numWells = Option.builder("w") - .longOpt("num-wells") + .required().build(); + Option numWells = Option.builder("w")// add this to plate options + .longOpt("wells") + .desc("The number of wells on the sample plate") .hasArg() - .argName("number") - .desc("The number of wells on each plate") + .required().build(); + //options group for choosing with distribution to use + OptionGroup distributions = new OptionGroup();// add this to plate options + distributions.setRequired(true); + Option poisson = Option.builder("poisson") + .desc("Use a Poisson distribution for cell sample") .build(); - mainOptions.addOption(numWells); - Option numPlates = Option.builder("np") - .longOpt("num-plates") + Option gaussian = Option.builder("gaussian") + .desc("Use a Gaussian distribution for cell sample") + .build(); + Option exponential = Option.builder("exponential") + .desc("Use an exponential distribution for cell sample") + .build(); + distributions.addOption(poisson); + distributions.addOption(gaussian); + distributions.addOption(exponential); + //options group for statistical distribution parameters + OptionGroup statParams = new OptionGroup();// add this to plate options + Option stdDev = Option.builder("stddev") + .desc("Standard deviation for Gaussian distribution") .hasArg() - .argName("number") - .desc("The number of plate files to output") .build(); - mainOptions.addOption(numPlates); - //Option plateOutput = (Option) outputFile.clone(); - //plateOutput.setRequired(true); - //plateOutput.setDescription("Prefix for plate output filenames"); - //mainOptions.addOption(plateOutput); - Option plateErr = Option.builder("err") - .longOpt("drop-out-rate") + Option lambda = Option.builder("lambda") + .desc("Lambda value for exponential distribution") .hasArg() - .argName("number") - .desc("Well drop-out rate. (Probability between 0 and 1)") .build(); - mainOptions.addOption(plateErr); - Option plateConcentrations = Option.builder("t") - .longOpt("t-cells-per-well") + statParams.addOption(stdDev); + statParams.addOption(lambda); + //Option group for random plate or set populations + OptionGroup wellPopOptions = new OptionGroup(); // add this to plate options + wellPopOptions.setRequired(true); + Option randomWellPopulations = Option.builder("random") + .desc("Randomize well populations on sample plate.") .hasArgs() - .argName("number 1, number 2, ...") - .desc("Number of T cells per well for each plate section") + .argName("MIN_POP MAX_POP") .build(); - mainOptions.addOption(plateConcentrations); - -//different distributions, mutually exclusive - OptionGroup plateDistributions = new OptionGroup(); - Option plateExp = Option.builder("exponential") - .desc("Sample from distinct cells with exponential frequency distribution") - .build(); - plateDistributions.addOption(plateExp); - Option plateGaussian = Option.builder("gaussian") - .desc("Sample from distinct cells with gaussain frequency distribution") - .build(); - plateDistributions.addOption(plateGaussian); - Option platePoisson = Option.builder("poisson") - .desc("Sample from distinct cells with poisson frequency distribution") - .build(); - plateDistributions.addOption(platePoisson); - mainOptions.addOptionGroup(plateDistributions); - - Option plateStdDev = Option.builder("stddev") - .desc("Standard deviation for gaussian distribution") - .hasArg() - .argName("number") - .build(); - mainOptions.addOption(plateStdDev); - - Option plateLambda = Option.builder("lambda") - .desc("Lambda for exponential distribution") - .hasArg() - .argName("number") - .build(); - mainOptions.addOption(plateLambda); - - - -// -// String cellFile, String filename, Double stdDev, -// Integer numWells, Integer numSections, -// Integer[] concentrations, Double dropOutRate -// - - //Options matchOptions = new Options(); - inputCells.setDescription("The cell sample file to be used for matching."); - mainOptions.addOption(inputCells); - Option lowThresh = Option.builder("low") - .longOpt("low-threshold") - .hasArg() - .argName("number") - .desc("Sets the minimum occupancy overlap to attempt matching") - .build(); - mainOptions.addOption(lowThresh); - Option highThresh = Option.builder("high") - .longOpt("high-threshold") - .hasArg() - .argName("number") - .desc("Sets the maximum occupancy overlap to attempt matching") - .build(); - mainOptions.addOption(highThresh); - Option occDiff = Option.builder("occdiff") - .longOpt("occupancy-difference") - .hasArg() - .argName("Number") - .desc("Maximum difference in alpha/beta occupancy to attempt matching") - .build(); - mainOptions.addOption(occDiff); - Option overlapPer = Option.builder("ovper") - .longOpt("overlap-percent") - .hasArg() - .argName("Percent") - .desc("Minimum overlap percent to attempt matching (0 -100)") - .build(); - mainOptions.addOption(overlapPer); - Option inputPlates = Option.builder("p") - .longOpt("plate-files") + Option specificWellPopulations = Option.builder("pop") + .longOpt("populations") + .desc("The well populations for each section of the sample plate") .hasArgs() - .desc("Plate files to match") + .argName("SECTION_1_POP [SECTION_2_POP] [SECTION_3_POP] ...") .build(); - mainOptions.addOption(inputPlates); + Option dropoutRate = Option.builder("err") //add this to plate options + .longOpt("dropout-rate") + .hasArg() + .desc("The sequence dropout rate due to amplification error") + .argName("DROPOUT_RATE (value between 0.0 and 1.0)") + .required() + .build(); + wellPopOptions.addOption(randomWellPopulations); + wellPopOptions.addOption(specificWellPopulations); + plateOptions.addOption(cellFile); + plateOptions.addOption(numWells); + plateOptions.addOptionGroup(distributions); + plateOptions.addOptionGroup(statParams); + plateOptions.addOptionGroup(wellPopOptions); + plateOptions.addOption(dropoutRate); + plateOptions.addOption(outputFileOption()); + return plateOptions; + } + + private static Options buildGraphOptions() { + Options graphOptions = new Options(); + Option cellFilename = Option.builder("c") + .longOpt("cell-file") + .desc("Cell sample file to use for checking accuracy") + .hasArg() + .argName("CELL_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") + .hasArg() + .argName("PLATE_FILENAME") + .required().build(); + Option outputGraphML = Option.builder("graphml") + .desc("Output GraphML file") + .build(); + Option outputSerializedBinary = Option.builder("no-binary") + .desc("Don't output serialized binary file") + .build(); + graphOptions.addOption(cellFilename); + graphOptions.addOption(plateFilename); + graphOptions.addOption(outputFileOption()); + return graphOptions; + } + + private static Options buildMatchCDR3options() { + Options matchCDR3options = new Options(); + Option graphFilename = Option.builder("g") + .longOpt("graph-file") + .desc("Graph/data file to use") + .hasArg() + .argName("GRAPH/DATA_FILENAME") + .required().build(); + Option minOccupancyOverlap = Option.builder("min") + .longOpt("min-overlap-size") + .desc("The minimum number of shared wells to attempt to match a sequence pair") + .hasArg() + .argName("MIN_OVERLAP") + .required().build(); + Option maxOccupancyOverlap = Option.builder("max") + .longOpt("max_overlap_size") + .desc("The maximum number of shared wells to attempt to match a sequence pair") + .hasArg() + .argName("MAX_OVERLAP") + .required().build(); + Option minOverlapPercent = Option.builder("minpct") + .longOpt("min-overlap-percent") + .desc("The minimum percentage of a sequence's total occupancy shared by another sequence to attempt matching") + .hasArg() + .build(); + Option maxOccupancyDifference = Option.builder("maxdiff") + .longOpt("max-occupancy-difference") + .desc("The maximum difference in total occupancy between two sequences to attempt matching") + .hasArg() + .build(); + matchCDR3options.addOption(graphFilename); + matchCDR3options.addOption(minOccupancyOverlap); + matchCDR3options.addOption(maxOccupancyOverlap); + matchCDR3options.addOption(minOverlapPercent); + matchCDR3options.addOption(maxOccupancyDifference); + matchCDR3options.addOption(outputFileOption()); + return matchCDR3options; + } - CommandLineParser parser = new DefaultParser(); - try { - CommandLine line = parser.parse(mainOptions, args); - if(line.hasOption("match")){ - //line = parser.parse(mainOptions, args); - //String cellFile = line.getOptionValue("c"); - String graphFile = line.getOptionValue("g"); - Integer lowThreshold = Integer.valueOf(line.getOptionValue(lowThresh)); - Integer highThreshold = Integer.valueOf(line.getOptionValue(highThresh)); - Integer occupancyDifference = Integer.valueOf(line.getOptionValue(occDiff)); - Integer overlapPercent = Integer.valueOf(line.getOptionValue(overlapPer)); - for(String plate: line.getOptionValues("p")) { - matchCDR3s(graphFile, lowThreshold, highThreshold, occupancyDifference, overlapPercent); - } - } - else if(line.hasOption("cells")){ - //line = parser.parse(mainOptions, args); - String filename = line.getOptionValue("o"); - Integer numDistCells = Integer.valueOf(line.getOptionValue("nc")); - Integer freq = Integer.valueOf(line.getOptionValue("d")); - makeCells(filename, numDistCells, freq); - } - else if(line.hasOption("plates")){ - //line = parser.parse(mainOptions, args); - String cellFile = line.getOptionValue("c"); - String filenamePrefix = line.getOptionValue("o"); - Integer numWellsOnPlate = Integer.valueOf(line.getOptionValue("w")); - Integer numPlatesToMake = Integer.valueOf(line.getOptionValue("np")); - String[] concentrationsToUseString = line.getOptionValues("t"); - Integer numSections = concentrationsToUseString.length; + private static CellSample getCells(String cellFilename) { + assert cellFilename != null; + CellFileReader reader = new CellFileReader(cellFilename); + return reader.getCellSample(); + } - Integer[] concentrationsToUse = new Integer[numSections]; - for(int i = 0; i Date: Sun, 27 Feb 2022 15:54:47 -0600 Subject: [PATCH 14/59] Initialize wells --- src/main/java/Plate.java | 1 + 1 file changed, 1 insertion(+) diff --git a/src/main/java/Plate.java b/src/main/java/Plate.java index aae0a07..3bdc713 100644 --- a/src/main/java/Plate.java +++ b/src/main/java/Plate.java @@ -25,6 +25,7 @@ public class Plate { this.cells = cells; this.sourceFile = cellFilename; this.size = numWells; + this.wells = new ArrayList<>(); this.error = dropoutRate; this.populations = populations; this.exponential = exponential; From 3d1f8668ee3b9a052af525ab72a4371036242496 Mon Sep 17 00:00:00 2001 From: efischer Date: Sun, 27 Feb 2022 16:03:57 -0600 Subject: [PATCH 15/59] Control verbose output --- src/main/java/GraphDataObjectWriter.java | 18 ++++++++++++++++-- 1 file changed, 16 insertions(+), 2 deletions(-) diff --git a/src/main/java/GraphDataObjectWriter.java b/src/main/java/GraphDataObjectWriter.java index 3cf79e8..6b6a832 100644 --- a/src/main/java/GraphDataObjectWriter.java +++ b/src/main/java/GraphDataObjectWriter.java @@ -1,3 +1,5 @@ +import org.jgrapht.Graph; + import java.io.BufferedOutputStream; import java.io.FileOutputStream; import java.io.IOException; @@ -7,6 +9,7 @@ public class GraphDataObjectWriter { private GraphWithMapData data; private String filename; + private boolean verbose = true; public GraphDataObjectWriter(String filename, GraphWithMapData data) { if(!filename.matches(".*\\.ser")){ @@ -16,13 +19,24 @@ public class GraphDataObjectWriter { this.data = data; } + public GraphDataObjectWriter(String filename, GraphWithMapData data, boolean verbose) { + this.verbose = verbose; + if(!filename.matches(".*\\.ser")){ + filename = filename + ".ser"; + } + this.filename = filename; + this.data = data; + } + public void writeDataToFile() { try (BufferedOutputStream bufferedOut = new BufferedOutputStream(new FileOutputStream(filename)); ObjectOutputStream out = new ObjectOutputStream(bufferedOut); ){ - System.out.println("Writing graph and occupancy data to file. This may take some time."); - System.out.println("File I/O time is not included in results."); + if(verbose) { + System.out.println("Writing graph and occupancy data to file. This may take some time."); + System.out.println("File I/O time is not included in results."); + } out.writeObject(data); } catch (IOException ex) { ex.printStackTrace(); From a822f69ea4650ed19ccad2ceadb0874d2e31398e Mon Sep 17 00:00:00 2001 From: efischer Date: Sun, 27 Feb 2022 16:07:17 -0600 Subject: [PATCH 16/59] Control verbose output --- src/main/java/GraphDataObjectReader.java | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/src/main/java/GraphDataObjectReader.java b/src/main/java/GraphDataObjectReader.java index c522be8..fea0bff 100644 --- a/src/main/java/GraphDataObjectReader.java +++ b/src/main/java/GraphDataObjectReader.java @@ -1,10 +1,12 @@ import java.io.*; public class GraphDataObjectReader { + private GraphWithMapData data; private String filename; + private boolean verbose = true; - public GraphDataObjectReader(String filename) throws IOException { + public GraphDataObjectReader(String filename, boolean verbose) throws IOException { if(!filename.matches(".*\\.ser")){ filename = filename + ".ser"; } From 05556bce0c2c7f29a6811b709f59a187ed05a903 Mon Sep 17 00:00:00 2001 From: efischer Date: Sun, 27 Feb 2022 16:08:59 -0600 Subject: [PATCH 17/59] Add units to metadata --- src/main/java/Simulator.java | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/main/java/Simulator.java b/src/main/java/Simulator.java index da13878..9818aa5 100644 --- a/src/main/java/Simulator.java +++ b/src/main/java/Simulator.java @@ -280,7 +280,7 @@ public class Simulator implements GraphModificationFunctions { 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", nf.format(time.toSeconds())); + metadata.put("simulation time (seconds)", nf.format(time.toSeconds())); //create MatchingResult object MatchingResult output = new MatchingResult(metadata, header, allResults, matchMap, time); if(verbose){ From 2485ac4cf627ce73c7bcd0029581d277e59cf807 Mon Sep 17 00:00:00 2001 From: efischer Date: Sun, 27 Feb 2022 16:15:26 -0600 Subject: [PATCH 18/59] Add getters to MatchingResult --- src/main/java/MatchingResult.java | 36 ++++++++++++++++++++++--------- 1 file changed, 26 insertions(+), 10 deletions(-) diff --git a/src/main/java/MatchingResult.java b/src/main/java/MatchingResult.java index 5942c82..4502f51 100644 --- a/src/main/java/MatchingResult.java +++ b/src/main/java/MatchingResult.java @@ -21,15 +21,15 @@ public class MatchingResult { * well populations * * total alphas found * * total betas found * - * high overlap threshold - * low overlap threshold - * maximum occupancy difference - * minimum overlap percent - * pairing attempt rate - * correct pairing count - * incorrect pairing count - * pairing error rate - * simulation time + * high overlap threshold * + * low overlap threshold * + * maximum occupancy difference * + * minimum overlap percent * + * pairing attempt rate * + * correct pairing count * + * incorrect pairing count * + * pairing error rate * + * simulation time (seconds) */ this.metadata = metadata; this.comments = new ArrayList<>(); @@ -91,6 +91,22 @@ public class MatchingResult { return Integer.parseInt(metadata.get("total beta count")); } - //put in the rest of these methods following the same pattern + public Integer getHighOverlapThreshold() { return Integer.parseInt(metadata.get("high overlap threshold"));} + + public Integer getLowOverlapThreshold() { return Integer.parseInt(metadata.get("low overlap threshold"));} + + public Integer getMaxOccupancyDifference() { return Integer.parseInt(metadata.get("maximum occupancy difference"));} + + public Integer getMinOverlapPercent() { return Integer.parseInt(metadata.get("minimum overlap percent"));} + + public Double getPairingAttemptRate() { return Double.parseDouble(metadata.get("pairing attempt rate"));} + + public Integer getCorrectPairingCount() { return Integer.parseInt(metadata.get("correct pairing count"));} + + public Integer getIncorrectPairingCount() { return Integer.parseInt(metadata.get("incorrect pairing count"));} + + public Double getPairingErrorRate() { return Double.parseDouble(metadata.get("pairing error rate"));} + + public String getSimulationTime() { return metadata.get("simulation time (seconds)"); } } From 32c5bcaaff94aa85bcf4ab4cc15a2c062b1a4199 Mon Sep 17 00:00:00 2001 From: efischer Date: Sun, 27 Feb 2022 16:16:24 -0600 Subject: [PATCH 19/59] Deactivate file I/O announcement for CLI --- src/main/java/CommandLineInterface.java | 30 ++++++++++++------------- src/main/java/InteractiveInterface.java | 2 +- 2 files changed, 15 insertions(+), 17 deletions(-) diff --git a/src/main/java/CommandLineInterface.java b/src/main/java/CommandLineInterface.java index aadfb3a..75943f9 100644 --- a/src/main/java/CommandLineInterface.java +++ b/src/main/java/CommandLineInterface.java @@ -1,6 +1,7 @@ import org.apache.commons.cli.*; import java.io.IOException; +import java.util.Arrays; import java.util.stream.Stream; /* @@ -48,21 +49,15 @@ import java.util.stream.Stream; public class CommandLineInterface { public static void startCLI(String[] args) { - //These command line options are a big mess - //Really, I don't think command line tools are expected to work in this many different modes - //making cells, making plates, and matching are the sort of thing that UNIX philosophy would say - //should be three separate programs. - //There might be a way to do it with option parameters? - //main options set - for the four different program modes Options mainOptions = buildMainOptions(); CommandLineParser parser = new DefaultParser(); try{ - CommandLine line = parser.parse(mainOptions, args); + CommandLine line = parser.parse(mainOptions, Arrays.copyOfRange(args, 0, 1)); if (line.hasOption("cells")) { Options cellOptions = buildCellOptions(); - line = parser.parse(cellOptions, args); + line = parser.parse(cellOptions, Arrays.copyOfRange(args, 1, args.length)); Integer number = Integer.valueOf(line.getOptionValue("n")); Integer diversity = Integer.valueOf(line.getOptionValue("d")); String filename = line.getOptionValue("o"); @@ -71,7 +66,7 @@ public class CommandLineInterface { else if (line.hasOption("plate")) { Options plateOptions = buildPlateOptions(); - line = parser.parse(plateOptions, args); + line = parser.parse(plateOptions, Arrays.copyOfRange(args, 1, args.length)); //get the cells String cellFilename = line.getOptionValue("c"); CellSample cells = getCells(cellFilename); @@ -122,7 +117,7 @@ public class CommandLineInterface { else if (line.hasOption("graph")) { //Making a graph Options graphOptions = buildGraphOptions(); - line = parser.parse(graphOptions, args); + line = parser.parse(graphOptions, Arrays.copyOfRange(args, 1, args.length)); String cellFilename = line.getOptionValue("c"); String plateFilename = line.getOptionValue("p"); String outputFilename = line.getOptionValue("o"); @@ -132,7 +127,7 @@ public class CommandLineInterface { Plate plate = getPlate(plateFilename); GraphWithMapData graph = Simulator.makeGraph(cells, plate, false); if (!line.hasOption("no-binary")) { //output binary file unless told not to - GraphDataObjectWriter writer = new GraphDataObjectWriter(outputFilename, graph); + GraphDataObjectWriter writer = new GraphDataObjectWriter(outputFilename, graph, false); writer.writeDataToFile(); } if (line.hasOption("graphml")) { //if told to, output graphml file @@ -143,7 +138,7 @@ public class CommandLineInterface { else if (line.hasOption("match")) { //can add a flag for which match type in future, spit this in two Options matchOptions = buildMatchCDR3options(); - line = parser.parse(matchOptions, args); + line = parser.parse(matchOptions, Arrays.copyOfRange(args, 1, args.length)); String graphFilename = line.getOptionValue("g"); String outputFilename = line.getOptionValue("o"); Integer minThreshold = Integer.parseInt(line.getOptionValue("min")); @@ -193,8 +188,8 @@ public class CommandLineInterface { .longOpt("make-cells") .desc("Makes a cell sample file of distinct T cells") .build(); - Option makePlate = Option.builder("plates") - .longOpt("make-plates") + Option makePlate = Option.builder("plate") + .longOpt("make-plate") .desc("Makes a sample plate file. Requires a cell sample file.") .build(); Option makeGraph = Option.builder("graph") @@ -324,12 +319,15 @@ public class CommandLineInterface { Option outputGraphML = Option.builder("graphml") .desc("Output GraphML file") .build(); - Option outputSerializedBinary = Option.builder("no-binary") + Option outputSerializedBinary = Option.builder("nb") + .longOpt("no-binary") .desc("Don't output serialized binary file") .build(); graphOptions.addOption(cellFilename); graphOptions.addOption(plateFilename); graphOptions.addOption(outputFileOption()); + graphOptions.addOption(outputGraphML); + graphOptions.addOption(outputSerializedBinary); return graphOptions; } @@ -389,7 +387,7 @@ public class CommandLineInterface { private static GraphWithMapData getGraph(String graphFilename) { assert graphFilename != null; try{ - GraphDataObjectReader reader = new GraphDataObjectReader(graphFilename); + GraphDataObjectReader reader = new GraphDataObjectReader(graphFilename, false); return reader.getData(); } diff --git a/src/main/java/InteractiveInterface.java b/src/main/java/InteractiveInterface.java index d826b1f..bcedcbc 100644 --- a/src/main/java/InteractiveInterface.java +++ b/src/main/java/InteractiveInterface.java @@ -375,7 +375,7 @@ public class InteractiveInterface { data = BiGpairSEQ.getGraphInMemory(); } else { - GraphDataObjectReader dataReader = new GraphDataObjectReader(graphFilename); + GraphDataObjectReader dataReader = new GraphDataObjectReader(graphFilename, true); data = dataReader.getData(); if(BiGpairSEQ.cacheGraph()) { BiGpairSEQ.setGraphInMemory(data, graphFilename); From 12b003a69f7af869680a60c7af4a8c4124d91330 Mon Sep 17 00:00:00 2001 From: efischer Date: Sun, 27 Feb 2022 16:45:30 -0600 Subject: [PATCH 20/59] Add -help CLI option --- src/main/java/CommandLineInterface.java | 63 +++++++++++++++++-------- 1 file changed, 44 insertions(+), 19 deletions(-) diff --git a/src/main/java/CommandLineInterface.java b/src/main/java/CommandLineInterface.java index 75943f9..a6d84b0 100644 --- a/src/main/java/CommandLineInterface.java +++ b/src/main/java/CommandLineInterface.java @@ -49,14 +49,30 @@ import java.util.stream.Stream; public class CommandLineInterface { public static void startCLI(String[] args) { - //main options set - for the four different program modes + //Options sets for the different modes Options mainOptions = buildMainOptions(); + Options cellOptions = buildCellOptions(); + Options plateOptions = buildPlateOptions(); + Options graphOptions = buildGraphOptions(); + Options matchOptions = buildMatchCDR3options(); + CommandLineParser parser = new DefaultParser(); try{ CommandLine line = parser.parse(mainOptions, Arrays.copyOfRange(args, 0, 1)); - if (line.hasOption("cells")) { - Options cellOptions = buildCellOptions(); + if (line.hasOption("help")) { + HelpFormatter formatter = new HelpFormatter(); + formatter.printHelp("BiGpairSEQ_Sim", mainOptions); + System.out.println(); + formatter.printHelp("BiGpairSEQ_SIM -cells", cellOptions); + System.out.println(); + formatter.printHelp("BiGpairSEQ_Sim -plate", plateOptions); + System.out.println(); + formatter.printHelp("BiGpairSEQ_Sim -graph", graphOptions); + System.out.println(); + formatter.printHelp("BiGpairSEQ_Sim -match", matchOptions); + } + else if (line.hasOption("cells")) { line = parser.parse(cellOptions, Arrays.copyOfRange(args, 1, args.length)); Integer number = Integer.valueOf(line.getOptionValue("n")); Integer diversity = Integer.valueOf(line.getOptionValue("d")); @@ -65,7 +81,6 @@ public class CommandLineInterface { } else if (line.hasOption("plate")) { - Options plateOptions = buildPlateOptions(); line = parser.parse(plateOptions, Arrays.copyOfRange(args, 1, args.length)); //get the cells String cellFilename = line.getOptionValue("c"); @@ -116,7 +131,6 @@ public class CommandLineInterface { } else if (line.hasOption("graph")) { //Making a graph - Options graphOptions = buildGraphOptions(); line = parser.parse(graphOptions, Arrays.copyOfRange(args, 1, args.length)); String cellFilename = line.getOptionValue("c"); String plateFilename = line.getOptionValue("p"); @@ -137,7 +151,6 @@ public class CommandLineInterface { } else if (line.hasOption("match")) { //can add a flag for which match type in future, spit this in two - Options matchOptions = buildMatchCDR3options(); line = parser.parse(matchOptions, Arrays.copyOfRange(args, 1, args.length)); String graphFilename = line.getOptionValue("g"); String outputFilename = line.getOptionValue("o"); @@ -184,6 +197,7 @@ public class CommandLineInterface { private static Options buildMainOptions() { Options mainOptions = new Options(); + Option help = Option.builder("help").build(); Option makeCells = Option.builder("cells") .longOpt("make-cells") .desc("Makes a cell sample file of distinct T cells") @@ -201,6 +215,7 @@ public class CommandLineInterface { .desc("Matches CDR3s. Requires a graph/data file.") .build(); OptionGroup mainGroup = new OptionGroup(); + mainGroup.addOption(help); mainGroup.addOption(makeCells); mainGroup.addOption(makePlate); mainGroup.addOption(makeGraph); @@ -236,11 +251,13 @@ public class CommandLineInterface { .longOpt("cell-file") .desc("The cell sample file to use") .hasArg() + .argName("filename") .required().build(); Option numWells = Option.builder("w")// add this to plate options .longOpt("wells") .desc("The number of wells on the sample plate") .hasArg() + .argName("number") .required().build(); //options group for choosing with distribution to use OptionGroup distributions = new OptionGroup();// add this to plate options @@ -260,12 +277,14 @@ public class CommandLineInterface { //options group for statistical distribution parameters OptionGroup statParams = new OptionGroup();// add this to plate options Option stdDev = Option.builder("stddev") - .desc("Standard deviation for Gaussian distribution") + .desc("If using -gaussian flag, standard deviation for distrbution") .hasArg() + .argName("value") .build(); Option lambda = Option.builder("lambda") - .desc("Lambda value for exponential distribution") + .desc("If using -exponential flag, lambda value for distribution") .hasArg() + .argName("value") .build(); statParams.addOption(stdDev); statParams.addOption(lambda); @@ -275,19 +294,20 @@ public class CommandLineInterface { Option randomWellPopulations = Option.builder("random") .desc("Randomize well populations on sample plate.") .hasArgs() - .argName("MIN_POP MAX_POP") + .numberOfArgs(2) + .argName("max_pop min_pop") .build(); Option specificWellPopulations = Option.builder("pop") .longOpt("populations") - .desc("The well populations for each section of the sample plate") + .desc("The well populations for each section of the sample plate. There will be as many sections as populations given.") .hasArgs() - .argName("SECTION_1_POP [SECTION_2_POP] [SECTION_3_POP] ...") + .argName("1st_pop [2nd_pop]...") .build(); Option dropoutRate = Option.builder("err") //add this to plate options .longOpt("dropout-rate") .hasArg() .desc("The sequence dropout rate due to amplification error") - .argName("DROPOUT_RATE (value between 0.0 and 1.0)") + .argName("rate") .required() .build(); wellPopOptions.addOption(randomWellPopulations); @@ -308,13 +328,13 @@ public class CommandLineInterface { .longOpt("cell-file") .desc("Cell sample file to use for checking accuracy") .hasArg() - .argName("CELL_FILENAME") + .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") .hasArg() - .argName("PLATE_FILENAME") + .argName("filename") .required().build(); Option outputGraphML = Option.builder("graphml") .desc("Output GraphML file") @@ -337,29 +357,31 @@ public class CommandLineInterface { .longOpt("graph-file") .desc("Graph/data file to use") .hasArg() - .argName("GRAPH/DATA_FILENAME") + .argName("filename") .required().build(); Option minOccupancyOverlap = Option.builder("min") .longOpt("min-overlap-size") .desc("The minimum number of shared wells to attempt to match a sequence pair") .hasArg() - .argName("MIN_OVERLAP") + .argName("minimum") .required().build(); Option maxOccupancyOverlap = Option.builder("max") .longOpt("max_overlap_size") .desc("The maximum number of shared wells to attempt to match a sequence pair") .hasArg() - .argName("MAX_OVERLAP") + .argName("maximum") .required().build(); Option minOverlapPercent = Option.builder("minpct") .longOpt("min-overlap-percent") - .desc("The minimum percentage of a sequence's total occupancy shared by another sequence to attempt matching") + .desc("(Optional) The minimum percentage of a sequence's total occupancy shared by another sequence to attempt matching. (0 - 100) ") .hasArg() + .argName("minimum") .build(); Option maxOccupancyDifference = Option.builder("maxdiff") .longOpt("max-occupancy-difference") - .desc("The maximum difference in total occupancy between two sequences to attempt matching") + .desc("(Optional) The maximum difference in total occupancy between two sequences to attempt matching.") .hasArg() + .argName("maximum") .build(); matchCDR3options.addOption(graphFilename); matchCDR3options.addOption(minOccupancyOverlap); @@ -367,6 +389,9 @@ public class CommandLineInterface { matchCDR3options.addOption(minOverlapPercent); matchCDR3options.addOption(maxOccupancyDifference); matchCDR3options.addOption(outputFileOption()); + //options for output to System.out + //Option printPairingErrorRate = Option.builder() + return matchCDR3options; } From b1155f8100e17697ef9d840777db55fcabd86e80 Mon Sep 17 00:00:00 2001 From: efischer Date: Sun, 27 Feb 2022 16:53:46 -0600 Subject: [PATCH 21/59] Format -help CLI option --- src/main/java/CommandLineInterface.java | 30 +++++++++++-------------- 1 file changed, 13 insertions(+), 17 deletions(-) diff --git a/src/main/java/CommandLineInterface.java b/src/main/java/CommandLineInterface.java index a6d84b0..ea5ca23 100644 --- a/src/main/java/CommandLineInterface.java +++ b/src/main/java/CommandLineInterface.java @@ -197,7 +197,9 @@ public class CommandLineInterface { private static Options buildMainOptions() { Options mainOptions = new Options(); - Option help = Option.builder("help").build(); + Option help = Option.builder("help") + .desc("Displays this help menu") + .build(); Option makeCells = Option.builder("cells") .longOpt("make-cells") .desc("Makes a cell sample file of distinct T cells") @@ -292,21 +294,19 @@ public class CommandLineInterface { OptionGroup wellPopOptions = new OptionGroup(); // add this to plate options wellPopOptions.setRequired(true); Option randomWellPopulations = Option.builder("random") - .desc("Randomize well populations on sample plate.") + .desc("Randomize well populations on sample plate. Takes two arguments: the minimum possible population and the maximum possible population.") .hasArgs() .numberOfArgs(2) - .argName("max_pop min_pop") + .argName("minimum maximum") .build(); Option specificWellPopulations = Option.builder("pop") - .longOpt("populations") - .desc("The well populations for each section of the sample plate. There will be as many sections as populations given.") + .desc("The well populations for each section of the sample plate. There will be as many sections as there are populations given.") .hasArgs() - .argName("1st_pop [2nd_pop]...") + .argName("number [number]...") .build(); Option dropoutRate = Option.builder("err") //add this to plate options - .longOpt("dropout-rate") .hasArg() - .desc("The sequence dropout rate due to amplification error") + .desc("The sequence dropout rate due to amplification error. (0.0 - 1.0)") .argName("rate") .required() .build(); @@ -355,33 +355,29 @@ public class CommandLineInterface { Options matchCDR3options = new Options(); Option graphFilename = Option.builder("g") .longOpt("graph-file") - .desc("Graph/data file to use") + .desc("The graph/data file to use") .hasArg() .argName("filename") .required().build(); Option minOccupancyOverlap = Option.builder("min") - .longOpt("min-overlap-size") .desc("The minimum number of shared wells to attempt to match a sequence pair") .hasArg() - .argName("minimum") + .argName("number") .required().build(); Option maxOccupancyOverlap = Option.builder("max") - .longOpt("max_overlap_size") .desc("The maximum number of shared wells to attempt to match a sequence pair") .hasArg() - .argName("maximum") + .argName("number") .required().build(); Option minOverlapPercent = Option.builder("minpct") - .longOpt("min-overlap-percent") .desc("(Optional) The minimum percentage of a sequence's total occupancy shared by another sequence to attempt matching. (0 - 100) ") .hasArg() - .argName("minimum") + .argName("percent") .build(); Option maxOccupancyDifference = Option.builder("maxdiff") - .longOpt("max-occupancy-difference") .desc("(Optional) The maximum difference in total occupancy between two sequences to attempt matching.") .hasArg() - .argName("maximum") + .argName("number") .build(); matchCDR3options.addOption(graphFilename); matchCDR3options.addOption(minOccupancyOverlap); From e10350c214182dd05a5becb809b0113e182dfec9 Mon Sep 17 00:00:00 2001 From: efischer Date: Sun, 27 Feb 2022 16:56:58 -0600 Subject: [PATCH 22/59] Update readme with CLI arguments --- readme.md | 12 +++++++++--- 1 file changed, 9 insertions(+), 3 deletions(-) diff --git a/readme.md b/readme.md index 6a35f81..918b3ca 100644 --- a/readme.md +++ b/readme.md @@ -48,8 +48,13 @@ For example, to run the program with 32 gigabytes of memory, use the command: `java -Xmx32G -jar BiGpairSEQ_Sim.jar` -Once running, BiGpairSEQ_Sim has an interactive, menu-driven CLI for generating files and simulating TCR pairing. The -main menu looks like this: +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: + +`java -jar BiGpairSEQ_SIM -help` + +If no command line arguments are given, BiGpairSEQ_Sim will launch with an interactive, menu-driven CLI for +generating files and simulating TCR pairing. The main menu looks like this: ``` --------BiGPairSEQ SIMULATOR-------- @@ -78,6 +83,7 @@ By default, the Options menu looks like this: 0) Return to main menu ``` + ### INPUT/OUTPUT To run the simulation, the program reads and writes 4 kinds of files: @@ -290,7 +296,7 @@ slightly less time than the simulation itself. Real elapsed time from start to f * ~~Enable GraphML output in addition to serialized object binaries, for data portability~~ DONE * ~~Custom vertex type with attribute for sequence occupancy?~~ ABANDONED * Have a branch where this is implemented, but there's a bug that broke matching. Don't currently have time to fix. -* Re-implement command line arguments, to enable scripting and statistical simulation studies +* ~~Re-implement command line arguments, to enable scripting and statistical simulation studies~~ DONE * Re-implement CDR1 matching method * Implement Duan and Su's maximum weight matching algorithm * Add controllable algorithm-type parameter? From fb4d22e7a40d8c19b4f1fd686d33af68202b61ab Mon Sep 17 00:00:00 2001 From: efischer Date: Sun, 27 Feb 2022 17:00:54 -0600 Subject: [PATCH 23/59] Update readme with CLI arguments --- readme.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/readme.md b/readme.md index 918b3ca..48e51e2 100644 --- a/readme.md +++ b/readme.md @@ -51,7 +51,7 @@ For example, to run the program with 32 gigabytes of memory, use the command: 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: -`java -jar BiGpairSEQ_SIM -help` +`java -jar BiGpairSEQ_SIM.jar -help` If no command line arguments are given, BiGpairSEQ_Sim will launch with an interactive, menu-driven CLI for generating files and simulating TCR pairing. The main menu looks like this: From 6f5afbc6eca4440c2ad143347f36091a2d7d97e0 Mon Sep 17 00:00:00 2001 From: efischer Date: Sun, 27 Feb 2022 17:01:12 -0600 Subject: [PATCH 24/59] Update readme with CLI arguments --- readme.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/readme.md b/readme.md index 48e51e2..cc6cb03 100644 --- a/readme.md +++ b/readme.md @@ -51,7 +51,7 @@ For example, to run the program with 32 gigabytes of memory, use the command: 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: -`java -jar BiGpairSEQ_SIM.jar -help` +`java -jar BiGpairSEQ_Sim.jar -help` If no command line arguments are given, BiGpairSEQ_Sim will launch with an interactive, menu-driven CLI for generating files and simulating TCR pairing. The main menu looks like this: From c8364d8a6e9d3a8323827d5ffb8dd620fa035315 Mon Sep 17 00:00:00 2001 From: efischer Date: Sun, 27 Feb 2022 17:34:20 -0600 Subject: [PATCH 25/59] check verbose flag --- src/main/java/GraphDataObjectReader.java | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/src/main/java/GraphDataObjectReader.java b/src/main/java/GraphDataObjectReader.java index fea0bff..c8940c8 100644 --- a/src/main/java/GraphDataObjectReader.java +++ b/src/main/java/GraphDataObjectReader.java @@ -15,8 +15,10 @@ 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) { ex.printStackTrace(); From f2347e8fc22033a2da63c4dea258e47a41784be0 Mon Sep 17 00:00:00 2001 From: efischer Date: Sun, 27 Feb 2022 17:35:50 -0600 Subject: [PATCH 26/59] check verbose flag --- src/main/java/GraphDataObjectReader.java | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/main/java/GraphDataObjectReader.java b/src/main/java/GraphDataObjectReader.java index c8940c8..8e883ee 100644 --- a/src/main/java/GraphDataObjectReader.java +++ b/src/main/java/GraphDataObjectReader.java @@ -15,7 +15,7 @@ public class GraphDataObjectReader { BufferedInputStream fileIn = new BufferedInputStream(new FileInputStream(filename)); ObjectInputStream in = new ObjectInputStream(fileIn)) { - if(verbose) { + if (this.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"); } From bc5d67680dfe776dd5a432b41e881149706a4c40 Mon Sep 17 00:00:00 2001 From: efischer Date: Sun, 27 Feb 2022 17:36:23 -0600 Subject: [PATCH 27/59] Add flag to print metadata to stdout --- src/main/java/CommandLineInterface.java | 71 ++++++++++++++++++------- 1 file changed, 53 insertions(+), 18 deletions(-) diff --git a/src/main/java/CommandLineInterface.java b/src/main/java/CommandLineInterface.java index ea5ca23..d10fb3d 100644 --- a/src/main/java/CommandLineInterface.java +++ b/src/main/java/CommandLineInterface.java @@ -153,17 +153,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 +180,17 @@ 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)); + } + } } } catch (ParseException exp) { @@ -297,7 +311,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> Date: Sun, 27 Feb 2022 19:08:29 -0600 Subject: [PATCH 28/59] bugfix --- src/main/java/GraphDataObjectReader.java | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/main/java/GraphDataObjectReader.java b/src/main/java/GraphDataObjectReader.java index 8e883ee..ff7edb4 100644 --- a/src/main/java/GraphDataObjectReader.java +++ b/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,7 +15,7 @@ public class GraphDataObjectReader { BufferedInputStream fileIn = new BufferedInputStream(new FileInputStream(filename)); ObjectInputStream in = new ObjectInputStream(fileIn)) { - if (this.verbose) { + 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"); } From e04d2d677790b89ee9399f59ef534a778602d5dd Mon Sep 17 00:00:00 2001 From: efischer Date: Tue, 1 Mar 2022 08:16:06 -0600 Subject: [PATCH 29/59] Fix typos in help menu --- src/main/java/CommandLineInterface.java | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/src/main/java/CommandLineInterface.java b/src/main/java/CommandLineInterface.java index d10fb3d..b0d2a64 100644 --- a/src/main/java/CommandLineInterface.java +++ b/src/main/java/CommandLineInterface.java @@ -62,15 +62,15 @@ 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("cells")) { line = parser.parse(cellOptions, Arrays.copyOfRange(args, 1, args.length)); From f301327693ff0277e229dc34052b2b121f8194d7 Mon Sep 17 00:00:00 2001 From: efischer Date: Tue, 1 Mar 2022 08:24:43 -0600 Subject: [PATCH 30/59] Update readme with -graphml flag --- readme.md | 9 +++++---- 1 file changed, 5 insertions(+), 4 deletions(-) diff --git a/readme.md b/readme.md index cc6cb03..88a4240 100644 --- a/readme.md +++ b/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. @@ -203,8 +203,9 @@ Options for creating a Graph/Data file: 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.) +(For portability to other software, turn on GraphML 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.) --- From 25acf920c231d618a60beb82cf13cd19dad14ab1 Mon Sep 17 00:00:00 2001 From: efischer Date: Tue, 1 Mar 2022 08:34:35 -0600 Subject: [PATCH 31/59] Add version information --- src/main/java/BiGpairSEQ.java | 3 ++- src/main/java/CommandLineInterface.java | 6 ++++++ src/main/java/InteractiveInterface.java | 2 ++ 3 files changed, 10 insertions(+), 1 deletion(-) diff --git a/src/main/java/BiGpairSEQ.java b/src/main/java/BiGpairSEQ.java index 3e67a68..7b3cb95 100644 --- a/src/main/java/BiGpairSEQ.java +++ b/src/main/java/BiGpairSEQ.java @@ -16,6 +16,7 @@ public class BiGpairSEQ { private static String priorityQueueHeapType = "FIBONACCI"; private static boolean outputBinary = true; private static boolean outputGraphML = false; + private static final String version = "version 2.0"; public static void main(String[] args) { if (args.length == 0) { @@ -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; } } diff --git a/src/main/java/CommandLineInterface.java b/src/main/java/CommandLineInterface.java index b0d2a64..8549ee7 100644 --- a/src/main/java/CommandLineInterface.java +++ b/src/main/java/CommandLineInterface.java @@ -72,6 +72,9 @@ public class CommandLineInterface { System.out.println(); 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)); Integer number = Integer.valueOf(line.getOptionValue("n")); @@ -230,8 +233,11 @@ public class CommandLineInterface { .longOpt("match-cdr3") .desc("Matches CDR3s. Requires a graph/data file.") .build(); + Option printVersion = Option.builder("version") + .desc("Print the version number").build(); OptionGroup mainGroup = new OptionGroup(); mainGroup.addOption(help); + mainGroup.addOption(printVersion); mainGroup.addOption(makeCells); mainGroup.addOption(makePlate); mainGroup.addOption(makeGraph); diff --git a/src/main/java/InteractiveInterface.java b/src/main/java/InteractiveInterface.java index bcedcbc..7080807 100644 --- a/src/main/java/InteractiveInterface.java +++ b/src/main/java/InteractiveInterface.java @@ -570,6 +570,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(); From 9fb3095f0ff0a7111e652f6f11289f53f5755939 Mon Sep 17 00:00:00 2001 From: efischer Date: Tue, 1 Mar 2022 08:40:34 -0600 Subject: [PATCH 32/59] Clarify help text --- src/main/java/CommandLineInterface.java | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/main/java/CommandLineInterface.java b/src/main/java/CommandLineInterface.java index 8549ee7..43fac13 100644 --- a/src/main/java/CommandLineInterface.java +++ b/src/main/java/CommandLineInterface.java @@ -234,7 +234,7 @@ public class CommandLineInterface { .desc("Matches CDR3s. Requires a graph/data file.") .build(); Option printVersion = Option.builder("version") - .desc("Print the version number").build(); + .desc("Prints the program version number to stdout").build(); OptionGroup mainGroup = new OptionGroup(); mainGroup.addOption(help); mainGroup.addOption(printVersion); From b5a8b7e2d50e072b331d6c35f383ba21d00d7f85 Mon Sep 17 00:00:00 2001 From: efischer Date: Tue, 1 Mar 2022 08:41:57 -0600 Subject: [PATCH 33/59] update readme --- readme.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/readme.md b/readme.md index 88a4240..7bcbac8 100644 --- a/readme.md +++ b/readme.md @@ -203,7 +203,7 @@ Options for creating a Graph/Data file: 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 in interactive mode, or use the -graphml +(For portability to other software, turn on GraphML 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.) From 67ec3f37642c4b418f42857224ce987f519850cc Mon Sep 17 00:00:00 2001 From: efischer Date: Tue, 1 Mar 2022 08:43:18 -0600 Subject: [PATCH 34/59] update readme --- readme.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/readme.md b/readme.md index 7bcbac8..75f6bdf 100644 --- a/readme.md +++ b/readme.md @@ -49,7 +49,7 @@ 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` From 2ac34518421c70db248b57c7ea5b5db130749d5f Mon Sep 17 00:00:00 2001 From: efischer Date: Tue, 1 Mar 2022 08:43:48 -0600 Subject: [PATCH 35/59] update readme --- readme.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/readme.md b/readme.md index 75f6bdf..f171759 100644 --- a/readme.md +++ b/readme.md @@ -43,7 +43,7 @@ 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` From bedf0894bccf4cc5c08730d2ebc1ec51a06096f2 Mon Sep 17 00:00:00 2001 From: efischer Date: Tue, 1 Mar 2022 08:45:40 -0600 Subject: [PATCH 36/59] update readme --- readme.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/readme.md b/readme.md index f171759..6d123af 100644 --- a/readme.md +++ b/readme.md @@ -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. --- From 1886800873ecc8662de28b8accd333f17c3ffeb6 Mon Sep 17 00:00:00 2001 From: efischer Date: Tue, 1 Mar 2022 08:54:32 -0600 Subject: [PATCH 37/59] update readme --- readme.md | 10 +++++++--- 1 file changed, 7 insertions(+), 3 deletions(-) diff --git a/readme.md b/readme.md index 6d123af..2c07a33 100644 --- a/readme.md +++ b/readme.md @@ -203,9 +203,13 @@ Options for creating a Graph/Data file: 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 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.) +*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. --- From 64209691f0fd6cc7a6e5b76e2b060777e91f3b79 Mon Sep 17 00:00:00 2001 From: efischer Date: Tue, 1 Mar 2022 09:00:58 -0600 Subject: [PATCH 38/59] Check for finite pairing error rate --- src/main/java/Simulator.java | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/main/java/Simulator.java b/src/main/java/Simulator.java index 9818aa5..f74312e 100644 --- a/src/main/java/Simulator.java +++ b/src/main/java/Simulator.java @@ -245,7 +245,7 @@ public class Simulator implements GraphModificationFunctions { //rate of pairing error double pairingErrorRate = (double) falseCount / (trueCount + falseCount); BigDecimal pairingErrorRateTrunc; - if(pairingErrorRate == NaN || pairingErrorRate == POSITIVE_INFINITY || pairingErrorRate == NEGATIVE_INFINITY) { + if(!Double.isFinite(pairingErrorRate)) { pairingErrorRateTrunc = new BigDecimal(-1, mc); } else{ From 8275cf7740b9901af83265b3f999705bd614f296 Mon Sep 17 00:00:00 2001 From: efischer Date: Tue, 1 Mar 2022 09:01:53 -0600 Subject: [PATCH 39/59] Check for finite pairing error rate --- src/main/java/Simulator.java | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/main/java/Simulator.java b/src/main/java/Simulator.java index f74312e..a1b2609 100644 --- a/src/main/java/Simulator.java +++ b/src/main/java/Simulator.java @@ -245,11 +245,11 @@ public class Simulator implements GraphModificationFunctions { //rate of pairing error double pairingErrorRate = (double) falseCount / (trueCount + falseCount); BigDecimal pairingErrorRateTrunc; - if(!Double.isFinite(pairingErrorRate)) { - 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(); From 578bdc0fbf6e5cb4d15d368baf885f5d91c85a49 Mon Sep 17 00:00:00 2001 From: efischer Date: Tue, 1 Mar 2022 15:08:43 -0600 Subject: [PATCH 40/59] clarify help menu text --- src/main/java/CommandLineInterface.java | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/main/java/CommandLineInterface.java b/src/main/java/CommandLineInterface.java index 43fac13..7a5a7c0 100644 --- a/src/main/java/CommandLineInterface.java +++ b/src/main/java/CommandLineInterface.java @@ -422,7 +422,7 @@ public class CommandLineInterface { // 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 all metadata to stdout").build(); + .desc("(Optional) Print summary of matching results to stdout.").build(); matchCDR3options // .addOption(printErrorRate) // .addOption(printAttempt) From 3fc39302c78546db1574d51e29215ee649217a3c Mon Sep 17 00:00:00 2001 From: efischer Date: Tue, 1 Mar 2022 15:24:14 -0600 Subject: [PATCH 41/59] Add detail to error message --- src/main/java/GraphDataObjectReader.java | 1 + 1 file changed, 1 insertion(+) diff --git a/src/main/java/GraphDataObjectReader.java b/src/main/java/GraphDataObjectReader.java index ff7edb4..3db45d8 100644 --- a/src/main/java/GraphDataObjectReader.java +++ b/src/main/java/GraphDataObjectReader.java @@ -21,6 +21,7 @@ public class GraphDataObjectReader { } data = (GraphWithMapData) in.readObject(); } catch (FileNotFoundException | ClassNotFoundException ex) { + System.out.println("Graph/data file " + filename + " not found."); ex.printStackTrace(); } } From 4c872ed48e0e9742b6c7c807fc706ca185ff39bf Mon Sep 17 00:00:00 2001 From: efischer Date: Tue, 1 Mar 2022 15:27:04 -0600 Subject: [PATCH 42/59] Add optional stdout print flags --- src/main/java/CommandLineInterface.java | 61 +++++++++++++++++++------ 1 file changed, 46 insertions(+), 15 deletions(-) diff --git a/src/main/java/CommandLineInterface.java b/src/main/java/CommandLineInterface.java index 7a5a7c0..1bf4dbc 100644 --- a/src/main/java/CommandLineInterface.java +++ b/src/main/java/CommandLineInterface.java @@ -194,6 +194,27 @@ public class CommandLineInterface { 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) { @@ -413,22 +434,32 @@ public class CommandLineInterface { .addOption(outputFile); //options for output to System.out -// 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 summary of matching results to stdout.").build(); + 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(printErrorRate) + .addOption(printAttempt) + .addOption(printCorrect) + .addOption(printIncorrect) + .addOption(printMetadata) + .addOption(printAlphaCount) + .addOption(printBetaCount) + .addOption(printTime); return matchCDR3options; } From 582dc3ef4086516723e1372805952a2ac116c8ba Mon Sep 17 00:00:00 2001 From: efischer Date: Wed, 2 Mar 2022 12:39:40 -0600 Subject: [PATCH 43/59] Update readme --- readme.md | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/readme.md b/readme.md index 2c07a33..adc07ae 100644 --- a/readme.md +++ b/readme.md @@ -264,17 +264,16 @@ 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 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 @@ -286,7 +285,7 @@ slightly less time than the simulation itself. Real elapsed time from start to f * ~~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?~~ @@ -300,6 +299,7 @@ slightly less time than the simulation itself. Real elapsed time from start to f * _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 + * Advantage: would eliminate the need to use maps to associate vertices with sequences, which would make the code easier to understand. * Have a branch where this is implemented, but there's a bug that broke matching. Don't currently have time to fix. * ~~Re-implement command line arguments, to enable scripting and statistical simulation studies~~ DONE * Re-implement CDR1 matching method @@ -319,7 +319,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 From 03e8d3121048d00183d6eb5b0bcf4cfc97ab3fd1 Mon Sep 17 00:00:00 2001 From: efischer Date: Wed, 2 Mar 2022 18:55:19 -0600 Subject: [PATCH 44/59] Add data on randomized well population behavior --- readme.md | 43 +++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 43 insertions(+) diff --git a/readme.md b/readme.md index adc07ae..ece0552 100644 --- a/readme.md +++ b/readme.md @@ -279,6 +279,49 @@ Since this implementation of BiGpairSEQ writes intermediate results to disk (to 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. +## 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. + +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| + +From these results, it can be seen 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 average well population. + ## TODO * ~~Try invoking GC at end of workloads to reduce paging to disk~~ DONE From df047267ee8c40c63dcf64857842dd1f5d550e01 Mon Sep 17 00:00:00 2001 From: efischer Date: Wed, 2 Mar 2022 22:54:17 -0600 Subject: [PATCH 45/59] Add data on randomized well population behavior --- readme.md | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/readme.md b/readme.md index ece0552..894ea2d 100644 --- a/readme.md +++ b/readme.md @@ -264,6 +264,7 @@ 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 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), @@ -279,6 +280,9 @@ Since this implementation of BiGpairSEQ writes intermediate results to disk (to 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 weighted 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, @@ -294,6 +298,7 @@ The well populations of the plates were: * 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: From 325e1ebe2bda75bfff2d01db68f09a6825942494 Mon Sep 17 00:00:00 2001 From: efischer Date: Wed, 2 Mar 2022 23:21:56 -0600 Subject: [PATCH 46/59] Add data on randomized well population behavior --- readme.md | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/readme.md b/readme.md index 894ea2d..5c3210b 100644 --- a/readme.md +++ b/readme.md @@ -324,8 +324,9 @@ 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| -From these results, it can be seen 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 average well population. +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 From b602fb02f1fb339caeaabeb2d068b4fa99bdf814 Mon Sep 17 00:00:00 2001 From: efischer Date: Wed, 2 Mar 2022 23:35:24 -0600 Subject: [PATCH 47/59] Remove obsolete comments --- src/main/java/Simulator.java | 1 - 1 file changed, 1 deletion(-) diff --git a/src/main/java/Simulator.java b/src/main/java/Simulator.java index a1b2609..e0e186c 100644 --- a/src/main/java/Simulator.java +++ b/src/main/java/Simulator.java @@ -157,7 +157,6 @@ public class Simulator implements GraphModificationFunctions { "removed");} //Find Maximum Weighted Matching - //using jheaps library class PairingHeap for improved efficiency if(verbose){System.out.println("Finding maximum weighted matching");} MaximumWeightBipartiteMatching maxWeightMatching; //Use correct heap type for priority queue From 96ba57d65358e6c38412be0414e2fb65387c19cf Mon Sep 17 00:00:00 2001 From: efischer Date: Fri, 4 Mar 2022 16:14:17 -0600 Subject: [PATCH 48/59] Remove singleton sequences from wells in initial filtering --- src/main/java/Simulator.java | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/main/java/Simulator.java b/src/main/java/Simulator.java index e0e186c..8bf462b 100644 --- a/src/main/java/Simulator.java +++ b/src/main/java/Simulator.java @@ -46,9 +46,9 @@ public class Simulator implements GraphModificationFunctions { if(verbose){System.out.println("All betas count: " + betaCount);} if(verbose){System.out.println("Well maps made");} - if(verbose){System.out.println("Removing sequences present in all wells.");} - filterByOccupancyThresholds(allAlphas, 1, numWells - 1); - filterByOccupancyThresholds(allBetas, 1, numWells - 1); + if(verbose){System.out.println("Removing singleton sequences and sequences present in all wells.");} + filterByOccupancyThresholds(allAlphas, 2, numWells - 1); + filterByOccupancyThresholds(allBetas, 2, numWells - 1); if(verbose){System.out.println("Sequences removed");} int pairableAlphaCount = allAlphas.size(); if(verbose){System.out.println("Remaining alphas count: " + pairableAlphaCount);} From 1df86f01dfb9f71f472e4e10cd73caa4db56c27c Mon Sep 17 00:00:00 2001 From: efischer Date: Sat, 5 Mar 2022 12:02:58 -0600 Subject: [PATCH 49/59] parameterized sequence indices --- src/main/java/Simulator.java | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/main/java/Simulator.java b/src/main/java/Simulator.java index 8bf462b..1c62872 100644 --- a/src/main/java/Simulator.java +++ b/src/main/java/Simulator.java @@ -23,6 +23,7 @@ public class Simulator implements GraphModificationFunctions { 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) { Instant start = Instant.now(); @@ -34,7 +35,7 @@ public class Simulator implements GraphModificationFunctions { if(verbose){System.out.println("Making cell maps");} //HashMap keyed to Alphas, values Betas - Map distCellsMapAlphaKey = makeSequenceToSequenceMap(distinctCells, 0, 1); + Map distCellsMapAlphaKey = makeSequenceToSequenceMap(distinctCells, cdr3AlphaIndex, cdr3BetaIndex); if(verbose){System.out.println("Cell maps made");} if(verbose){System.out.println("Making well maps");} From f980722b56034f37ed9a888edffa9e457f7530b9 Mon Sep 17 00:00:00 2001 From: efischer Date: Wed, 21 Sep 2022 18:09:37 +0000 Subject: [PATCH 50/59] update TODO --- readme.md | 1 + 1 file changed, 1 insertion(+) diff --git a/readme.md b/readme.md index 5c3210b..d568daa 100644 --- a/readme.md +++ b/readme.md @@ -330,6 +330,7 @@ roughly as though it had a constant well population equal to the plate's average ## TODO +* Enable post-filtering instead of pre-filtering. Pre-filtering of things like singleton sequences or saturating-occupancy sequences reduces graph size, but could conceivably reduce pairing accuracy by throwing away data. While these sequences have very little signal, it would be interesting to compare unfiltered results to filtered results. This would require a much, much faster MWM algorithm, though, to handle the much larger graphs. Possible one of the linear-time approximation algorithms. * ~~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.*~~ From cf771ce57427ac0c81285252daff90db305ca813 Mon Sep 17 00:00:00 2001 From: eugenefischer <66030419+eugenefischer@users.noreply.github.com> Date: Wed, 21 Sep 2022 13:15:49 -0500 Subject: [PATCH 51/59] parameterized sequence indices --- readme.md | 2 ++ src/main/java/Simulator.java | 12 ++++++++---- 2 files changed, 10 insertions(+), 4 deletions(-) diff --git a/readme.md b/readme.md index 5c3210b..99da1a3 100644 --- a/readme.md +++ b/readme.md @@ -356,6 +356,8 @@ roughly as though it had a constant well population equal to the plate's average * Add controllable algorithm-type parameter? * This would be fun and valuable, but probably take more time than I have for a hobby project. * 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 ## CITATIONS diff --git a/src/main/java/Simulator.java b/src/main/java/Simulator.java index 1c62872..933ad0e 100644 --- a/src/main/java/Simulator.java +++ b/src/main/java/Simulator.java @@ -47,10 +47,11 @@ public class Simulator implements GraphModificationFunctions { if(verbose){System.out.println("All betas count: " + betaCount);} if(verbose){System.out.println("Well maps made");} - if(verbose){System.out.println("Removing singleton sequences and sequences present in all wells.");} - filterByOccupancyThresholds(allAlphas, 2, numWells - 1); - filterByOccupancyThresholds(allBetas, 2, numWells - 1); - if(verbose){System.out.println("Sequences removed");} + +// if(verbose){System.out.println("Removing singleton sequences and sequences present in all wells.");} +// filterByOccupancyThresholds(allAlphas, 2, numWells - 1); +// filterByOccupancyThresholds(allBetas, 2, numWells - 1); +// if(verbose){System.out.println("Sequences removed");} int pairableAlphaCount = allAlphas.size(); if(verbose){System.out.println("Remaining alphas count: " + pairableAlphaCount);} int pairableBetaCount = allBetas.size(); @@ -239,6 +240,8 @@ public class Simulator implements GraphModificationFunctions { //Metadata comments for CSV file String algoType = "LEDA book with heap: " + heapType; int min = Math.min(alphaCount, betaCount); + //matching weight + BigDecimal totalMatchingWeight = maxWeightMatching.getMatchingWeight(); //rate of attempted matching double attemptRate = (double) (trueCount + falseCount) / min; BigDecimal attemptRateTrunc = new BigDecimal(attemptRate, mc); @@ -269,6 +272,7 @@ public class Simulator implements GraphModificationFunctions { metadata.put("sample plate filename", data.getSourceFilename()); metadata.put("graph filename", dataFilename); metadata.put("algorithm type", algoType); + metadata.put("matching weight", totalMatchingWeight.toString()); metadata.put("well populations", wellPopulationsString); metadata.put("total alphas found", alphaCount.toString()); metadata.put("total betas found", betaCount.toString()); From 3ba305abdb2324ef42fe4d87e29804c181b04e35 Mon Sep 17 00:00:00 2001 From: eugenefischer <66030419+eugenefischer@users.noreply.github.com> Date: Wed, 21 Sep 2022 13:30:30 -0500 Subject: [PATCH 52/59] Update ToDo --- readme.md | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/readme.md b/readme.md index 6347ba4..4908318 100644 --- a/readme.md +++ b/readme.md @@ -330,7 +330,6 @@ roughly as though it had a constant well population equal to the plate's average ## TODO -* Enable post-filtering instead of pre-filtering. Pre-filtering of things like singleton sequences or saturating-occupancy sequences reduces graph size, but could conceivably reduce pairing accuracy by throwing away data. While these sequences have very little signal, it would be interesting to compare unfiltered results to filtered results. This would require a much, much faster MWM algorithm, though, to handle the much larger graphs. Possible one of the linear-time approximation algorithms. * ~~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.*~~ @@ -356,9 +355,13 @@ roughly as though it had a constant well population equal to the plate's average * Implement Duan and Su's maximum weight matching algorithm * Add controllable algorithm-type parameter? * This would be fun and valuable, but probably take more time than I have for a hobby project. +* Implement an algorithm for approximating a maximum weight matching + * Some of these run in linear or near-linear time + * given that the underlying biological samples have many, many sources of error, this would probably be the most useful option in practice. It seems less mathematically elegant, though, and so less fun for me. * Implement Vose's alias method for arbitrary statistical distributions of cells * Should probably refactor to use apache commons rng for this * Use commons JCS for caching +* Enable post-filtering instead of pre-filtering. Pre-filtering of things like singleton sequences or saturating-occupancy sequences reduces graph size, but could conceivably reduce pairing accuracy by throwing away data. While these sequences have very little signal, it would be interesting to compare unfiltered results to filtered results. This would require a much, much faster MWM algorithm, though, to handle the much larger graphs. Possible one of the linear-time approximation algorithms. ## CITATIONS From 29b844afd2436962accbcbae91f7bcff1d04e046 Mon Sep 17 00:00:00 2001 From: eugenefischer <66030419+eugenefischer@users.noreply.github.com> Date: Wed, 21 Sep 2022 16:48:26 -0500 Subject: [PATCH 53/59] by-hand merge of needed code from custom vertex branch --- src/main/java/GraphMLFileWriter.java | 79 ++++++------ src/main/java/GraphModificationFunctions.java | 70 +++++------ src/main/java/GraphWithMapData.java | 95 +++++++------- src/main/java/SequenceType.java | 8 ++ src/main/java/Simulator.java | 119 ++++++++++-------- src/main/java/Vertex.java | 81 +++++++++++- 6 files changed, 277 insertions(+), 175 deletions(-) create mode 100644 src/main/java/SequenceType.java diff --git a/src/main/java/GraphMLFileWriter.java b/src/main/java/GraphMLFileWriter.java index 96ad358..8411502 100644 --- a/src/main/java/GraphMLFileWriter.java +++ b/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 graphAttributes; public GraphMLFileWriter(String filename, GraphWithMapData data) { if(!filename.matches(".*\\.graphml")){ @@ -27,52 +28,61 @@ public class GraphMLFileWriter { } this.filename = filename; this.data = data; + this.graph = data.getGraph(); + graphAttributes = createGraphAttributes(); } -// public void writeGraphToFile() { -// try(BufferedWriter writer = Files.newBufferedWriter(Path.of(filename), StandardOpenOption.CREATE_NEW); -// ){ -// GraphMLExporter 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 graph) { + if(!filename.matches(".*\\.graphml")){ + filename = filename + ".graphml"; + } + this.filename = filename; + this.graph = graph; + } + + private Map createGraphAttributes(){ + Map ga = new HashMap<>(); + //Sample plate filename + ga.put("sample plate filename", DefaultAttribute.createAttribute(data.getSourceFilename())); + // Number of wells + ga.put("well count", DefaultAttribute.createAttribute(data.getNumWells().toString())); + //Well populations + Integer[] wellPopulations = data.getWellPopulations(); + StringBuilder populationsStringBuilder = new StringBuilder(); + populationsStringBuilder.append(wellPopulations[0].toString()); + for(int i = 1; i < wellPopulations.length; i++){ + populationsStringBuilder.append(", "); + populationsStringBuilder.append(wellPopulations[i].toString()); + } + String wellPopulationsString = populationsStringBuilder.toString(); + ga.put("well populations", DefaultAttribute.createAttribute(wellPopulationsString)); + return ga; + } public void writeGraphToFile() { - SimpleWeightedGraph graph = data.getGraph(); - Map vertexToAlphaMap = data.getPlateVtoAMap(); - Map vertexToBetaMap = data.getPlateVtoBMap(); - Map alphaOccs = data.getAlphaWellCounts(); - Map 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> exporter = new GraphMLExporter<>(v -> v.toString()); + GraphMLExporter> exporter = new GraphMLExporter<>(v -> v.getVertexLabel().toString()); //set to export weights exporter.setExportEdgeWeights(true); + //Set graph attributes + exporter.setGraphAttributeProvider( () -> graphAttributes); //set type, sequence, and occupancy attributes for each vertex exporter.setVertexAttributeProvider( v -> { Map attributes = new HashMap<>(); - if(vertexToAlphaMap.containsKey(v)) { - attributes.put("type", DefaultAttribute.createAttribute("CDR3 Alpha")); - attributes.put("sequence", DefaultAttribute.createAttribute(vertexToAlphaMap.get(v))); - attributes.put("occupancy", DefaultAttribute.createAttribute( - alphaOccs.get(vertexToAlphaMap.get(v)))); - } - else if(vertexToBetaMap.containsKey(v)) { - attributes.put("type", DefaultAttribute.createAttribute("CDR3 Beta")); - attributes.put("sequence", DefaultAttribute.createAttribute(vertexToBetaMap.get(v))); - attributes.put("occupancy", DefaultAttribute.createAttribute( - betaOccs.get(vertexToBetaMap.get(v)))); - } + attributes.put("type", DefaultAttribute.createAttribute(v.getType().name())); + attributes.put("sequence", DefaultAttribute.createAttribute(v.getSequence())); + attributes.put("occupancy", DefaultAttribute.createAttribute(v.getOccupancy())); return attributes; }); //register the attributes - exporter.registerAttribute("type", GraphMLExporter.AttributeCategory.NODE, AttributeType.STRING); - exporter.registerAttribute("sequence", GraphMLExporter.AttributeCategory.NODE, AttributeType.STRING); - exporter.registerAttribute("occupancy", GraphMLExporter.AttributeCategory.NODE, AttributeType.STRING); + for(String s : graphAttributes.keySet()) { + exporter.registerAttribute(s, AttributeCategory.GRAPH, AttributeType.STRING); + } + exporter.registerAttribute("type", AttributeCategory.NODE, AttributeType.STRING); + exporter.registerAttribute("sequence", AttributeCategory.NODE, AttributeType.STRING); + exporter.registerAttribute("occupancy", AttributeCategory.NODE, AttributeType.STRING); //export the graph exporter.exportGraph(graph, writer); } catch(IOException ex){ @@ -81,4 +91,3 @@ public class GraphMLFileWriter { } } } - diff --git a/src/main/java/GraphModificationFunctions.java b/src/main/java/GraphModificationFunctions.java index a4ab1a0..68cc9d2 100644 --- a/src/main/java/GraphModificationFunctions.java +++ b/src/main/java/GraphModificationFunctions.java @@ -2,23 +2,25 @@ import org.jgrapht.graph.DefaultWeightedEdge; import org.jgrapht.graph.SimpleWeightedGraph; import java.util.ArrayList; +import java.util.HashMap; import java.util.List; import java.util.Map; public interface GraphModificationFunctions { //remove over- and under-weight edges - static List filterByOverlapThresholds(SimpleWeightedGraph graph, + static Map filterByOverlapThresholds(SimpleWeightedGraph graph, int low, int high, boolean saveEdges) { - List removedEdges = new ArrayList<>(); + Map removedEdges = new HashMap<>(); + //List removedEdges = new ArrayList<>(); for (DefaultWeightedEdge e : graph.edgeSet()) { if ((graph.getEdgeWeight(e) > high) || (graph.getEdgeWeight(e) < low)) { if(saveEdges) { - Integer source = graph.getEdgeSource(e); - Integer target = graph.getEdgeTarget(e); + Vertex source = graph.getEdgeSource(e); + Vertex target = graph.getEdgeTarget(e); Integer weight = (int) graph.getEdgeWeight(e); - Integer[] edge = {source, target, weight}; - removedEdges.add(edge); + Vertex[] edge = {source, target}; + removedEdges.put(edge, weight); } else { graph.setEdgeWeight(e, 0.0); @@ -26,7 +28,7 @@ public interface GraphModificationFunctions { } } if(saveEdges) { - for (Integer[] edge : removedEdges) { + for (Vertex[] edge : removedEdges.keySet()) { graph.removeEdge(edge[0], edge[1]); } } @@ -34,23 +36,19 @@ public interface GraphModificationFunctions { } //Remove edges for pairs with large occupancy discrepancy - static List filterByRelativeOccupancy(SimpleWeightedGraph graph, - Map alphaWellCounts, - Map betaWellCounts, - Map plateVtoAMap, - Map plateVtoBMap, + static Map filterByRelativeOccupancy(SimpleWeightedGraph graph, Integer maxOccupancyDifference, boolean saveEdges) { - List removedEdges = new ArrayList<>(); + Map removedEdges = new HashMap<>(); for (DefaultWeightedEdge e : graph.edgeSet()) { - Integer alphaOcc = alphaWellCounts.get(plateVtoAMap.get(graph.getEdgeSource(e))); - Integer betaOcc = betaWellCounts.get(plateVtoBMap.get(graph.getEdgeTarget(e))); + Integer alphaOcc = graph.getEdgeSource(e).getOccupancy(); + Integer betaOcc = graph.getEdgeTarget(e).getOccupancy(); if (Math.abs(alphaOcc - betaOcc) >= maxOccupancyDifference) { if (saveEdges) { - Integer source = graph.getEdgeSource(e); - Integer target = graph.getEdgeTarget(e); + Vertex source = graph.getEdgeSource(e); + Vertex target = graph.getEdgeTarget(e); Integer weight = (int) graph.getEdgeWeight(e); - Integer[] edge = {source, target, weight}; - removedEdges.add(edge); + Vertex[] edge = {source, target}; + removedEdges.put(edge, weight); } else { graph.setEdgeWeight(e, 0.0); @@ -58,7 +56,7 @@ public interface GraphModificationFunctions { } } if(saveEdges) { - for (Integer[] edge : removedEdges) { + for (Vertex[] edge : removedEdges.keySet()) { graph.removeEdge(edge[0], edge[1]); } } @@ -66,26 +64,22 @@ public interface GraphModificationFunctions { } //Remove edges for pairs where overlap size is significantly lower than the well occupancy - static List filterByOverlapPercent(SimpleWeightedGraph graph, - Map alphaWellCounts, - Map betaWellCounts, - Map plateVtoAMap, - Map plateVtoBMap, + static Map filterByOverlapPercent(SimpleWeightedGraph graph, Integer minOverlapPercent, boolean saveEdges) { - List removedEdges = new ArrayList<>(); + Map 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,18 +87,18 @@ public interface GraphModificationFunctions { } } if(saveEdges) { - for (Integer[] edge : removedEdges) { + for (Vertex[] edge : removedEdges.keySet()) { graph.removeEdge(edge[0], edge[1]); } } return removedEdges; } - static void addRemovedEdges(SimpleWeightedGraph graph, - List removedEdges) { - for (Integer[] edge : removedEdges) { + static void addRemovedEdges(SimpleWeightedGraph graph, + Map 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)); } } diff --git a/src/main/java/GraphWithMapData.java b/src/main/java/GraphWithMapData.java index 3795190..0e4c09b 100644 --- a/src/main/java/GraphWithMapData.java +++ b/src/main/java/GraphWithMapData.java @@ -15,32 +15,33 @@ public class GraphWithMapData implements java.io.Serializable { private Integer alphaCount; private Integer betaCount; private final Map distCellsMapAlphaKey; - private final Map plateVtoAMap; - private final Map plateVtoBMap; - private final Map plateAtoVMap; - private final Map plateBtoVMap; - private final Map alphaWellCounts; - private final Map betaWellCounts; +// private final Map plateVtoAMap; +// private final Map plateVtoBMap; +// private final Map plateAtoVMap; +// private final Map plateBtoVMap; +// private final Map alphaWellCounts; +// private final Map betaWellCounts; private final Duration time; public GraphWithMapData(SimpleWeightedGraph graph, Integer numWells, Integer[] wellConcentrations, - Integer alphaCount, Integer betaCount, - Map distCellsMapAlphaKey, Map plateVtoAMap, - Map plateVtoBMap, Map plateAtoVMap, - Map plateBtoVMap, Map alphaWellCounts, - Map betaWellCounts, Duration time) { + Map distCellsMapAlphaKey, Duration time){ + +// Map plateVtoAMap, Integer alphaCount, Integer betaCount, +// Map plateVtoBMap, Map plateAtoVMap, +// Map plateBtoVMap, Map alphaWellCounts, +// Map betaWellCounts,) { this.graph = graph; this.numWells = numWells; this.wellPopulations = wellConcentrations; this.alphaCount = alphaCount; this.betaCount = betaCount; this.distCellsMapAlphaKey = distCellsMapAlphaKey; - this.plateVtoAMap = plateVtoAMap; - this.plateVtoBMap = plateVtoBMap; - this.plateAtoVMap = plateAtoVMap; - this.plateBtoVMap = plateBtoVMap; - this.alphaWellCounts = alphaWellCounts; - this.betaWellCounts = betaWellCounts; +// this.plateVtoAMap = plateVtoAMap; +// this.plateVtoBMap = plateVtoBMap; +// this.plateAtoVMap = plateAtoVMap; +// this.plateBtoVMap = plateBtoVMap; +// this.alphaWellCounts = alphaWellCounts; +// this.betaWellCounts = betaWellCounts; this.time = time; } @@ -56,41 +57,41 @@ public class GraphWithMapData implements java.io.Serializable { return wellPopulations; } - public Integer getAlphaCount() { - return alphaCount; - } - - public Integer getBetaCount() { - return betaCount; - } +// public Integer getAlphaCount() { +// return alphaCount; +// } +// +// public Integer getBetaCount() { +// return betaCount; +// } public Map getDistCellsMapAlphaKey() { return distCellsMapAlphaKey; } - public Map getPlateVtoAMap() { - return plateVtoAMap; - } - - public Map getPlateVtoBMap() { - return plateVtoBMap; - } - - public Map getPlateAtoVMap() { - return plateAtoVMap; - } - - public Map getPlateBtoVMap() { - return plateBtoVMap; - } - - public Map getAlphaWellCounts() { - return alphaWellCounts; - } - - public Map getBetaWellCounts() { - return betaWellCounts; - } +// public Map getPlateVtoAMap() { +// return plateVtoAMap; +// } +// +// public Map getPlateVtoBMap() { +// return plateVtoBMap; +// } +// +// public Map getPlateAtoVMap() { +// return plateAtoVMap; +// } +// +// public Map getPlateBtoVMap() { +// return plateBtoVMap; +// } +// +// public Map getAlphaWellCounts() { +// return alphaWellCounts; +// } +// +// public Map getBetaWellCounts() { +// return betaWellCounts; +// } public Duration getTime() { return time; diff --git a/src/main/java/SequenceType.java b/src/main/java/SequenceType.java new file mode 100644 index 0000000..54de083 --- /dev/null +++ b/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 +} diff --git a/src/main/java/Simulator.java b/src/main/java/Simulator.java index 933ad0e..c94db0c 100644 --- a/src/main/java/Simulator.java +++ b/src/main/java/Simulator.java @@ -18,29 +18,29 @@ import static java.lang.Float.*; //NOTE: "sequence" in method and variable names refers to a peptide sequence from a simulated T cell public class Simulator implements GraphModificationFunctions { - private static final int cdr3AlphaIndex = 0; - private static final int cdr3BetaIndex = 1; - private static final int cdr1AlphaIndex = 2; - private static final int cdr1BetaIndex = 3; - //Make the graph needed for matching CDR3s + //Make the graph needed for matching sequences. + //sourceVertexIndices and targetVertexIndices are indices within the cell to use as for the two sets of vertices + //in the bipartite graph. "Source" and "target" are JGraphT terms for the two vertices an edge touches, + //even if not directed. public static GraphWithMapData makeGraph(CellSample cellSample, Plate samplePlate, boolean verbose) { Instant start = Instant.now(); List distinctCells = cellSample.getCells(); - int[] alphaIndex = {cdr3AlphaIndex}; - int[] betaIndex = {cdr3BetaIndex}; + int[] alphaIndices = {SequenceType.CDR3_ALPHA.ordinal()}; + int[] betaIndices = {SequenceType.CDR3_BETA.ordinal()}; int numWells = samplePlate.getSize(); if(verbose){System.out.println("Making cell maps");} //HashMap keyed to Alphas, values Betas - Map distCellsMapAlphaKey = makeSequenceToSequenceMap(distinctCells, cdr3AlphaIndex, cdr3BetaIndex); + Map distCellsMapAlphaKey = makeSequenceToSequenceMap(distinctCells, 0, 1); if(verbose){System.out.println("Cell maps made");} if(verbose){System.out.println("Making well maps");} - Map allAlphas = samplePlate.assayWellsSequenceS(alphaIndex); - Map allBetas = samplePlate.assayWellsSequenceS(betaIndex); + + Map allAlphas = samplePlate.assayWellsSequenceS(alphaIndices); + Map allBetas = samplePlate.assayWellsSequenceS(betaIndices); int alphaCount = allAlphas.size(); if(verbose){System.out.println("All alphas count: " + alphaCount);} int betaCount = allBetas.size(); @@ -80,29 +80,40 @@ public class Simulator implements GraphModificationFunctions { //(technically this is only 1/4 of an adjacency matrix, but that's all you need //for a bipartite graph, and all the SimpleWeightedBipartiteGraphMatrixGenerator class expects.) if(verbose){System.out.println("Creating adjacency matrix");} - //Count how many wells each alpha appears in + //Count how many wells each alpha sequence appears in Map alphaWellCounts = new HashMap<>(); - //count how many wells each beta appears in + //count how many wells each beta sequence appears in Map betaWellCounts = new HashMap<>(); //the adjacency matrix to be used by the graph generator double[][] weights = new double[plateVtoAMap.size()][plateVtoBMap.size()]; countSequencesAndFillMatrix(samplePlate, allAlphas, allBetas, plateAtoVMap, - plateBtoVMap, alphaIndex, betaIndex, alphaWellCounts, betaWellCounts, weights); + plateBtoVMap, alphaIndices, betaIndices, alphaWellCounts, betaWellCounts, weights); if(verbose){System.out.println("Matrix created");} //create bipartite graph if(verbose){System.out.println("Creating graph");} //the graph object - SimpleWeightedGraph graph = + SimpleWeightedGraph graph = new SimpleWeightedGraph<>(DefaultWeightedEdge.class); //the graph generator SimpleWeightedBipartiteGraphMatrixGenerator graphGenerator = new SimpleWeightedBipartiteGraphMatrixGenerator(); //the list of alpha vertices - List alphaVertices = new ArrayList<>(plateVtoAMap.keySet()); //This will work because LinkedHashMap preserves order of entry + //List alphaVertices = new ArrayList<>(plateVtoAMap.keySet()); //This will work because LinkedHashMap preserves order of entry + List alphaVertices = new ArrayList<>(); + //start with map of all alphas mapped to vertex values, get occupancy from the alphaWellCounts map + for (Integer seq : plateAtoVMap.keySet()) { + Vertex alphaVertex = new Vertex(SequenceType.CDR3_ALPHA, seq, alphaWellCounts.get(seq), plateAtoVMap.get(seq)); + alphaVertices.add(alphaVertex); + } graphGenerator.first(alphaVertices); //the list of beta vertices - List betaVertices = new ArrayList<>(plateVtoBMap.keySet()); - graphGenerator.second(betaVertices); //This will work because LinkedHashMap preserves order of entry + //List betaVertices = new ArrayList<>(plateVtoBMap.keySet());//This will work because LinkedHashMap preserves order of entry + List betaVertices = new ArrayList<>(); + for (Integer seq : plateBtoVMap.keySet()) { + Vertex betaVertex = new Vertex(SequenceType.CDR3_BETA, seq, betaWellCounts.get(seq), plateBtoVMap.get(seq)); + betaVertices.add(betaVertex); + } + graphGenerator.second(betaVertices); //use adjacency matrix of weight created previously graphGenerator.weights(weights); graphGenerator.generateGraph(graph); @@ -112,9 +123,7 @@ public class Simulator implements GraphModificationFunctions { Duration time = Duration.between(start, stop); //create GraphWithMapData object - GraphWithMapData output = new GraphWithMapData(graph, numWells, samplePlate.getPopulations(), alphaCount, betaCount, - distCellsMapAlphaKey, plateVtoAMap, plateVtoBMap, plateAtoVMap, - plateBtoVMap, alphaWellCounts, betaWellCounts, time); + GraphWithMapData output = new GraphWithMapData(graph, numWells, samplePlate.getPopulations(), distCellsMapAlphaKey, time); //Set source file name in graph to name of sample plate output.setSourceFilename(samplePlate.getFilename()); //return GraphWithMapData object @@ -126,39 +135,46 @@ public class Simulator implements GraphModificationFunctions { Integer highThreshold, Integer maxOccupancyDifference, Integer minOverlapPercent, boolean verbose) { Instant start = Instant.now(); - List removedEdges = new ArrayList<>(); + SimpleWeightedGraph graph = data.getGraph(); + Map removedEdges = new HashMap<>(); boolean saveEdges = BiGpairSEQ.cacheGraph(); int numWells = data.getNumWells(); - Integer alphaCount = data.getAlphaCount(); - Integer betaCount = data.getBetaCount(); + //Integer alphaCount = data.getAlphaCount(); + //Integer betaCount = data.getBetaCount(); Map distCellsMapAlphaKey = data.getDistCellsMapAlphaKey(); - Map plateVtoAMap = data.getPlateVtoAMap(); - Map plateVtoBMap = data.getPlateVtoBMap(); - Map alphaWellCounts = data.getAlphaWellCounts(); - Map betaWellCounts = data.getBetaWellCounts(); - SimpleWeightedGraph graph = data.getGraph(); + Set alphas = new HashSet<>(); + Set betas = new HashSet<>(); + for(Vertex v: graph.vertexSet()) { + if (SequenceType.CDR3_ALPHA.equals(v.getType())){ + alphas.add(v); + } + else { + betas.add(v); + } + } + Integer alphaCount = alphas.size(); + Integer betaCount = betas.size(); //remove edges with weights outside given overlap thresholds, add those to removed edge list if(verbose){System.out.println("Eliminating edges with weights outside overlap threshold values");} - removedEdges.addAll(GraphModificationFunctions.filterByOverlapThresholds(graph, lowThreshold, highThreshold, saveEdges)); + removedEdges.putAll(GraphModificationFunctions.filterByOverlapThresholds(graph, lowThreshold, highThreshold, saveEdges)); if(verbose){System.out.println("Over- and under-weight edges removed");} //remove edges between vertices with too small an overlap size, add those to removed edge list if(verbose){System.out.println("Eliminating edges with weights less than " + minOverlapPercent.toString() + " percent of vertex occupancy value.");} - removedEdges.addAll(GraphModificationFunctions.filterByOverlapPercent(graph, alphaWellCounts, betaWellCounts, - plateVtoAMap, plateVtoBMap, minOverlapPercent, saveEdges)); + removedEdges.putAll(GraphModificationFunctions.filterByOverlapPercent(graph, minOverlapPercent, saveEdges)); if(verbose){System.out.println("Edges with weights too far below a vertex occupancy value removed");} //Filter by relative occupancy if(verbose){System.out.println("Eliminating edges between vertices with occupancy difference > " + maxOccupancyDifference);} - removedEdges.addAll(GraphModificationFunctions.filterByRelativeOccupancy(graph, alphaWellCounts, betaWellCounts, - plateVtoAMap, plateVtoBMap, maxOccupancyDifference, saveEdges)); + removedEdges.putAll(GraphModificationFunctions.filterByRelativeOccupancy(graph, maxOccupancyDifference, saveEdges)); if(verbose){System.out.println("Edges between vertices of with excessively different occupancy values " + "removed");} //Find Maximum Weighted Matching + //using jheaps library class PairingHeap for improved efficiency if(verbose){System.out.println("Finding maximum weighted matching");} MaximumWeightBipartiteMatching maxWeightMatching; //Use correct heap type for priority queue @@ -166,20 +182,20 @@ public class Simulator implements GraphModificationFunctions { switch (heapType) { case "PAIRING" -> { maxWeightMatching = new MaximumWeightBipartiteMatching(graph, - plateVtoAMap.keySet(), - plateVtoBMap.keySet(), + alphas, + betas, i -> new PairingHeap(Comparator.naturalOrder())); } case "FIBONACCI" -> { maxWeightMatching = new MaximumWeightBipartiteMatching(graph, - plateVtoAMap.keySet(), - plateVtoBMap.keySet(), + alphas, + betas, i -> new FibonacciHeap(Comparator.naturalOrder())); } default -> { maxWeightMatching = new MaximumWeightBipartiteMatching(graph, - plateVtoAMap.keySet(), - plateVtoBMap.keySet()); + alphas, + betas); } } //get the matching @@ -209,11 +225,14 @@ public class Simulator implements GraphModificationFunctions { Map matchMap = new HashMap<>(); while(weightIter.hasNext()) { e = weightIter.next(); - Integer source = graph.getEdgeSource(e); - Integer target = graph.getEdgeTarget(e); + Vertex source = graph.getEdgeSource(e); + Vertex target = graph.getEdgeTarget(e); + //Integer source = graph.getEdgeSource(e); + //Integer target = graph.getEdgeTarget(e); //The match map is all matches found, not just true matches! - matchMap.put(plateVtoAMap.get(source), plateVtoBMap.get(target)); - check = plateVtoBMap.get(target).equals(distCellsMapAlphaKey.get(plateVtoAMap.get(source))); + matchMap.put(source.getSequence(), target.getSequence()); + check = target.getSequence().equals(distCellsMapAlphaKey.get(source.getSequence())); + //check = plateVtoBMap.get(target).equals(distCellsMapAlphaKey.get(plateVtoAMap.get(source))); if(check) { trueCount++; } @@ -221,17 +240,19 @@ public class Simulator implements GraphModificationFunctions { falseCount++; } List result = new ArrayList<>(); - result.add(plateVtoAMap.get(source).toString()); + //alpha sequence + result.add(source.getSequence().toString()); //alpha well count - result.add(alphaWellCounts.get(plateVtoAMap.get(source)).toString()); - result.add(plateVtoBMap.get(target).toString()); + result.add(source.getOccupancy().toString()); + //beta sequence + result.add(target.getSequence().toString()); //beta well count - result.add(betaWellCounts.get(plateVtoBMap.get(target)).toString()); + result.add(target.getOccupancy().toString()); //overlap count result.add(Double.toString(graph.getEdgeWeight(e))); result.add(Boolean.toString(check)); - double pValue = Equations.pValue(numWells, alphaWellCounts.get(plateVtoAMap.get(source)), - betaWellCounts.get(plateVtoBMap.get(target)), graph.getEdgeWeight(e)); + double pValue = Equations.pValue(numWells, source.getOccupancy(), + target.getOccupancy(), graph.getEdgeWeight(e)); BigDecimal pValueTrunc = new BigDecimal(pValue, mc); result.add(pValueTrunc.toString()); allResults.add(result); diff --git a/src/main/java/Vertex.java b/src/main/java/Vertex.java index ef962ae..a7ddda9 100644 --- a/src/main/java/Vertex.java +++ b/src/main/java/Vertex.java @@ -1,23 +1,92 @@ +import java.io.Serializable; +public class Vertex implements Serializable { + private SequenceType type; + private Integer vertexLabel; + private Integer sequence; + private Integer occupancy; -public class Vertex { - private final Integer vertexLabel; - private final Integer sequence; - private final Integer occupancy; + public Vertex(Integer vertexLabel) { + this.vertexLabel = vertexLabel; + } + public Vertex(String vertexLabel) { + this.vertexLabel = Integer.parseInt((vertexLabel)); + } - public Vertex(Integer vertexLabel, Integer sequence, Integer occupancy) { + public Vertex(SequenceType type, Integer sequence, Integer occupancy, Integer vertexLabel) { + this.type = type; this.vertexLabel = vertexLabel; this.sequence = sequence; this.occupancy = occupancy; } - public Integer getVertexLabel() { return vertexLabel; } + + public SequenceType getType() { + return type; + } + + public void setType(String type) { + this.type = SequenceType.valueOf(type); + } + + public Integer getVertexLabel() { + return vertexLabel; + } + + public void setVertexLabel(String label) { + this.vertexLabel = Integer.parseInt(label); + } public Integer getSequence() { + return sequence; } + public void setSequence(String sequence) { + this.sequence = Integer.parseInt(sequence); + } + public Integer getOccupancy() { return occupancy; } + + public void setOccupancy(String occupancy) { + this.occupancy = Integer.parseInt(occupancy); + } + + @Override //adapted from JGraphT example code + public int hashCode() + { + return (sequence == null) ? 0 : sequence.hashCode(); + } + + @Override //adapted from JGraphT example code + public boolean equals(Object obj) + { + if (this == obj) + return true; + if (obj == null) + return false; + if (getClass() != obj.getClass()) + return false; + Vertex other = (Vertex) obj; + if (sequence == null) { + return other.sequence == null; + } else { + return sequence.equals(other.sequence); + } + } + + + @Override //adapted from JGraphT example code + public String toString() + { + StringBuilder sb = new StringBuilder(); + sb.append("(").append(vertexLabel) + .append(", Type: ").append(type.name()) + .append(", Sequence: ").append(sequence) + .append(", Occupancy: ").append(occupancy).append(")"); + return sb.toString(); + } + } From a054c0c20a191e3bc8b4709b5d7419fe3aad0710 Mon Sep 17 00:00:00 2001 From: eugenefischer <66030419+eugenefischer@users.noreply.github.com> Date: Wed, 21 Sep 2022 16:50:00 -0500 Subject: [PATCH 54/59] update TODO --- readme.md | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/readme.md b/readme.md index 4908318..b6f2251 100644 --- a/readme.md +++ b/readme.md @@ -347,9 +347,9 @@ roughly as though it had a constant well population equal to the plate's average * ~~Apache Commons CSV library writes entries a row at a time~~ * _Got this working, but at the cost of a profoundly strange bug in graph occupancy filtering. Have reverted the repo until I can figure out what caused that. Given how easily Thingiverse transposes CSV matrices in R, might not even be worth fixing. * ~~Enable GraphML output in addition to serialized object binaries, for data portability~~ DONE - * ~~Custom vertex type with attribute for sequence occupancy?~~ ABANDONED + * ~~Custom vertex type with attribute for sequence occupancy?~~ DONE * Advantage: would eliminate the need to use maps to associate vertices with sequences, which would make the code easier to understand. - * Have a branch where this is implemented, but there's a bug that broke matching. Don't currently have time to fix. + * ~~Have a branch where this is implemented, but there's a bug that broke matching. Don't currently have time to fix.~~ * ~~Re-implement command line arguments, to enable scripting and statistical simulation studies~~ DONE * Re-implement CDR1 matching method * Implement Duan and Su's maximum weight matching algorithm From 63317f2aa05cc341eb96bc0bf563f2f4f7aabf9e Mon Sep 17 00:00:00 2001 From: eugenefischer <66030419+eugenefischer@users.noreply.github.com> Date: Wed, 21 Sep 2022 18:08:52 -0500 Subject: [PATCH 55/59] Reword output message --- src/main/java/InteractiveInterface.java | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/main/java/InteractiveInterface.java b/src/main/java/InteractiveInterface.java index 7080807..b0bc636 100644 --- a/src/main/java/InteractiveInterface.java +++ b/src/main/java/InteractiveInterface.java @@ -258,7 +258,7 @@ public class InteractiveInterface { cellFile = sc.next(); System.out.print("\nPlease enter name of an existing sample plate file: "); plateFile = sc.next(); - System.out.println("\nThe graph and occupancy data will be written to a serialized binary file."); + System.out.println("\nThe graph and occupancy data will be written to a file."); System.out.print("Please enter a name for the output file: "); filename = sc.next(); } catch (InputMismatchException ex) { From 06e72314b0e2924836dee19a152297e259306920 Mon Sep 17 00:00:00 2001 From: eugenefischer <66030419+eugenefischer@users.noreply.github.com> Date: Wed, 21 Sep 2022 21:43:47 -0500 Subject: [PATCH 56/59] Reword option menu item --- src/main/java/InteractiveInterface.java | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/main/java/InteractiveInterface.java b/src/main/java/InteractiveInterface.java index b0bc636..5055e00 100644 --- a/src/main/java/InteractiveInterface.java +++ b/src/main/java/InteractiveInterface.java @@ -504,7 +504,7 @@ public class InteractiveInterface { System.out.println("2) Turn " + getOnOff(!BiGpairSEQ.cachePlate()) + " plate file caching"); System.out.println("3) Turn " + getOnOff(!BiGpairSEQ.cacheGraph()) + " graph/data file caching"); System.out.println("4) Turn " + getOnOff(!BiGpairSEQ.outputBinary()) + " serialized binary graph output"); - System.out.println("5) Turn " + getOnOff(!BiGpairSEQ.outputGraphML()) + " GraphML graph output"); + System.out.println("5) Turn " + getOnOff(!BiGpairSEQ.outputGraphML()) + " GraphML graph output (for data portability to other programs)"); System.out.println("6) Maximum weight matching algorithm options"); System.out.println("0) Return to main menu"); try { From 73c83bf35d9748c5505f6fa28d2190b050b5228f Mon Sep 17 00:00:00 2001 From: eugenefischer <66030419+eugenefischer@users.noreply.github.com> Date: Wed, 21 Sep 2022 21:46:00 -0500 Subject: [PATCH 57/59] Add comment on map data encodng --- src/main/java/GraphWithMapData.java | 1 + 1 file changed, 1 insertion(+) diff --git a/src/main/java/GraphWithMapData.java b/src/main/java/GraphWithMapData.java index 0e4c09b..8f9ab03 100644 --- a/src/main/java/GraphWithMapData.java +++ b/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; From e4e5a1f97918fdf08468817af309441b03992817 Mon Sep 17 00:00:00 2001 From: eugenefischer <66030419+eugenefischer@users.noreply.github.com> Date: Thu, 22 Sep 2022 00:00:02 -0500 Subject: [PATCH 58/59] update version number --- src/main/java/BiGpairSEQ.java | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/main/java/BiGpairSEQ.java b/src/main/java/BiGpairSEQ.java index 7b3cb95..e3f46b8 100644 --- a/src/main/java/BiGpairSEQ.java +++ b/src/main/java/BiGpairSEQ.java @@ -16,7 +16,7 @@ public class BiGpairSEQ { private static String priorityQueueHeapType = "FIBONACCI"; private static boolean outputBinary = true; private static boolean outputGraphML = false; - private static final String version = "version 2.0"; + private static final String version = "version 3.0"; public static void main(String[] args) { if (args.length == 0) { From 5c03909a11ecf4e476943491f5647e3b9b33d19b Mon Sep 17 00:00:00 2001 From: eugenefischer <66030419+eugenefischer@users.noreply.github.com> Date: Thu, 22 Sep 2022 01:39:13 -0500 Subject: [PATCH 59/59] Restore pre-filtering of singleton and saturating sequences --- src/main/java/Simulator.java | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/main/java/Simulator.java b/src/main/java/Simulator.java index c94db0c..6ba2cbf 100644 --- a/src/main/java/Simulator.java +++ b/src/main/java/Simulator.java @@ -48,10 +48,10 @@ public class Simulator implements GraphModificationFunctions { if(verbose){System.out.println("Well maps made");} -// if(verbose){System.out.println("Removing singleton sequences and sequences present in all wells.");} -// filterByOccupancyThresholds(allAlphas, 2, numWells - 1); -// filterByOccupancyThresholds(allBetas, 2, numWells - 1); -// if(verbose){System.out.println("Sequences removed");} + if(verbose){System.out.println("Removing singleton sequences and sequences present in all wells.");} + filterByOccupancyThresholds(allAlphas, 2, numWells - 1); + filterByOccupancyThresholds(allBetas, 2, numWells - 1); + if(verbose){System.out.println("Sequences removed");} int pairableAlphaCount = allAlphas.size(); if(verbose){System.out.println("Remaining alphas count: " + pairableAlphaCount);} int pairableBetaCount = allBetas.size();