diff --git a/readme.md b/readme.md index fa8e894..b6f2251 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. @@ -43,13 +43,18 @@ 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` -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.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: ``` --------BiGPairSEQ SIMULATOR-------- @@ -66,6 +71,19 @@ 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 +93,25 @@ To run the simulation, the program reads and writes 4 kinds of files: * Matching Results files in CSV format These files are often generated in sequence. When entering filenames, it is not necessary to include the file extension -(.csv or .ser). When reading or writing files, the program will automatically add the correct extension to any filename without one. +(.csv or .ser). When reading or writing files, the program will automatically add the correct extension to any filename +without one. To save file I/O time, the most recent instance of each of these four -files either generated or read from disk can be cached in program memory. This is 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 @@ -107,7 +129,6 @@ Comments are preceded by `#` Structure: ---- # Sample contains 1 unique CDR1 for every 4 unique CDR3s. | Alpha CDR3 | Beta CDR3 | Alpha CDR1 | Beta CDR1 | |---|---|---|---| @@ -152,7 +173,6 @@ Dropout sequences are replaced with the value `-1`. Comments are preceded by `#` Structure: ---- ``` # Cell source file name: # Each row represents one well on the plate @@ -181,14 +201,24 @@ 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. + +*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. --- #### 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: @@ -203,7 +233,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 @@ -235,48 +264,105 @@ 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 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, +10 sample plate files were created. All of these sample plates had 96 wells, used an exponential distribution with a lambda of 0.6, and +had a sequence dropout rate of 10%. + +The well populations of the plates were: +* One sample plate with 1000 T cells/well +* One sample plate with 2000 T cells/well +* One sample plate with 3000 T cells/well +* One sample plate with 4000 T cells/well +* One sample plate with 5000 T cells/well +* Five sample plates with each individual well's population randomized, from 1000 to 5000 T cells. (Average population ~3000 T cells/well.) + +All BiGpairSEQ simulations were run with a low overlap threshold of 3 and a high overlap threshold of 94. +No optional filters were used, so pairing was attempted for all sequences with overlaps within the threshold values. + +Constant well population plate results: + +| |1000 Cell/Well Plate|2000 Cell/Well Plate|3000 Cell/Well Plate|4000 Cell/Well Plate|5000 Cell/Well Plate +|---|---|---|---|---|---| +|Total Alphas Found|6407|7330|7936|8278|8553| +|Total Betas Found|6405|7333|7968|8269|8582| +|Pairing Attempt Rate|0.661|0.653|0.600|0.579|0.559| +|Correct Pairing Count|4231|4749|4723|4761|4750| +|Incorrect Pairing Count|3|34|40|26|29| +|Pairing Error Rate|0.000709|0.00711|0.00840|0.00543|0.00607| +|Simulation Time (Seconds)|500|643|700|589|598| + +Randomized well population plate results: + +| |Random Plate 1 | Random Plate 2|Random Plate 3|Random Plate 4|Random Plate 5|Average| +|---|---|---|---|---|---|---| +Total Alphas Found|7853|7904|7964|7898|7917|7907| +Total Betas Found|7851|7891|7920|7910|7894|7893| +Pairing Attempt Rate|0.607|0.610|0.601|0.605|0.603|0.605| +Correct Pairing Count|4718|4782|4721|4755|4731|4741| +Incorrect Pairing Count|51|35|42|27|29|37| +Pairing Error Rate|0.0107|0.00727|0.00882|0.00565|0.00609|0.00771| +Simulation Time (Seconds)|590|677|730|618|615|646| + +The average results for the randomized plates are closest to the constant plate with 3000 T cells/well. +This and several other tests indicate that BiGpairSEQ treats a sample plate with a highly variable number of T cells/well +roughly as though it had a constant well population equal to the plate's average well population. + ## TODO * ~~Try invoking GC at end of workloads to reduce paging to disk~~ DONE * ~~Hold graph data in memory until another graph is read-in? ABANDONED UNABANDONED~~ DONE * ~~*No, this won't work, because BiGpairSEQ simulations alter the underlying graph based on filtering constraints. Changes would cascade with multiple experiments.*~~ * Might have figured out a way to do it, by taking edges out and then putting them back into the graph. This may actually be possible. - * It is possible, though the modifications to the graph incur their own performance penalties. Need testing to see which option is best. -* 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? + * It is possible, though the modifications to the graph incur their own performance penalties. Need testing to see which option is best. It may be computer-specific. * ~~Test whether pairing heap (currently used) or Fibonacci heap is more efficient for priority queue in current matching algorithm~~ DONE * ~~in theory Fibonacci heap should be more efficient, but complexity overhead may eliminate theoretical advantage~~ * ~~Add controllable heap-type parameter?~~ - * Parameter implemented. 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?~~ 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.~~ +* ~~Re-implement command line arguments, to enable scripting and statistical simulation studies~~ DONE +* Re-implement CDR1 matching method +* Implement Duan and Su's maximum weight matching algorithm + * Add controllable algorithm-type parameter? + * This would be fun and valuable, but probably take more time than I have for a hobby project. +* Implement an algorithm for approximating a maximum weight matching + * Some of these run in linear or near-linear time + * given that the underlying biological samples have many, many sources of error, this would probably be the most useful option in practice. It seems less mathematically elegant, though, and so less fun for me. +* Implement Vose's alias method for arbitrary statistical distributions of cells + * Should probably refactor to use apache commons rng for this +* Use commons JCS for caching +* Enable post-filtering instead of pre-filtering. Pre-filtering of things like singleton sequences or saturating-occupancy sequences reduces graph size, but could conceivably reduce pairing accuracy by throwing away data. While these sequences have very little signal, it would be interesting to compare unfiltered results to filtered results. This would require a much, much faster MWM algorithm, though, to handle the much larger graphs. Possible one of the linear-time approximation algorithms. + ## CITATIONS * 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) @@ -288,7 +374,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 diff --git a/src/main/java/BiGpairSEQ.java b/src/main/java/BiGpairSEQ.java index 23798b7..360e1f5 100644 --- a/src/main/java/BiGpairSEQ.java +++ b/src/main/java/BiGpairSEQ.java @@ -16,6 +16,7 @@ public class BiGpairSEQ { private static HeapType priorityQueueHeapType = HeapType.FIBONACCI; private static boolean outputBinary = true; private static boolean outputGraphML = false; + private static final String version = "version 3.0"; public static void main(String[] args) { if (args.length == 0) { @@ -23,8 +24,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."); } } @@ -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 0c527b1..1bf4dbc 100644 --- a/src/main/java/CommandLineInterface.java +++ b/src/main/java/CommandLineInterface.java @@ -1,5 +1,9 @@ import org.apache.commons.cli.*; +import java.io.IOException; +import java.util.Arrays; +import java.util.stream.Stream; + /* * Class for parsing options passed to program from command line * @@ -29,6 +33,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 @@ -43,242 +49,172 @@ import org.apache.commons.cli.*; 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 - Options mainOptions = new Options(); - Option makeCells = Option.builder("cells") - .longOpt("make-cells") - .desc("Makes a file of distinct cells") - .build(); - Option makePlate = Option.builder("plates") - .longOpt("make-plates") - .desc("Makes a sample plate file") - .build(); - Option makeGraph = Option.builder("graph") - .longOpt("make-graph") - .desc("Makes a graph and data file") - .build(); - Option matchCDR3 = Option.builder("match") - .longOpt("match-cdr3") - .desc("Match CDR3s. Requires a cell sample file and any number of plate files.") - .build(); - OptionGroup mainGroup = new OptionGroup(); - mainGroup.addOption(makeCells); - mainGroup.addOption(makePlate); - mainGroup.addOption(makeGraph); - mainGroup.addOption(matchCDR3); - mainGroup.setRequired(true); - mainOptions.addOptionGroup(mainGroup); - - //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") - .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") - .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); - - //Options plateOptions = new Options(); - Option inputCells = Option.builder("c") - .longOpt("cell-file") - .hasArg() - .argName("file") - .desc("The cell sample file used for filling wells") - .build(); - mainOptions.addOption(inputCells); - Option numWells = Option.builder("w") - .longOpt("num-wells") - .hasArg() - .argName("number") - .desc("The number of wells on each plate") - .build(); - mainOptions.addOption(numWells); - Option numPlates = Option.builder("np") - .longOpt("num-plates") - .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") - .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") - .hasArgs() - .argName("number 1, number 2, ...") - .desc("Number of T cells per well for each plate section") - .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") - .hasArgs() - .desc("Plate files to match") - .build(); - mainOptions.addOption(inputPlates); - - + //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, 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); - } + try{ + CommandLine line = parser.parse(mainOptions, Arrays.copyOfRange(args, 0, 1)); + + if (line.hasOption("help")) { + HelpFormatter formatter = new HelpFormatter(); + formatter.printHelp("BiGpairSEQ_Sim.jar", mainOptions); + System.out.println(); + formatter.printHelp("BiGpairSEQ_Sim.jar -cells", cellOptions); + System.out.println(); + formatter.printHelp("BiGpairSEQ_Sim.jar -plate", plateOptions); + System.out.println(); + formatter.printHelp("BiGpairSEQ_Sim.jar -graph", graphOptions); + System.out.println(); + formatter.printHelp("BiGpairSEQ_Sim.jar -match", matchOptions); } - else if(line.hasOption("cells")){ - //line = parser.parse(mainOptions, args); + 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")); + Integer diversity = Integer.valueOf(line.getOptionValue("d")); String filename = line.getOptionValue("o"); - Integer numDistCells = Integer.valueOf(line.getOptionValue("nc")); - Integer freq = Integer.valueOf(line.getOptionValue("d")); - makeCells(filename, numDistCells, freq); + makeCells(filename, number, diversity); } - 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; - Integer[] concentrationsToUse = new Integer[numSections]; - for(int i = 0; i graphAttributes = new HashMap<>(); - - 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.addGraphAttributeConsumer((str, attribute) -> { - graphAttributes.put(str, attribute.getValue()); - }); - importer.addVertexWithAttributesConsumer((vertex, attributes) -> { - vertex.setType(attributes.get("type").getValue()); - vertex.setSequence(attributes.get("sequence").getValue()); - vertex.setOccupancy((attributes.get("occupancy").getValue())); - }); - importer.importGraph(graph, reader); - } - catch (IOException ex) { - System.out.println("Graph file " + filename + " not found."); - System.err.println(ex); - } - } - - public SimpleWeightedGraph getGraph() { return graph; } - - public Map getGraphAttributes() { return graphAttributes; } - - public String getFilename() {return filename;} - -} 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; diff --git a/src/main/java/InteractiveInterface.java b/src/main/java/InteractiveInterface.java index 776f433..5055e00 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"); @@ -260,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) { @@ -292,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); } @@ -306,8 +304,7 @@ public class InteractiveInterface { System.out.println("Returning to main menu."); } else{ - List cells = cellSample.getCells(); - GraphWithMapData data = Simulator.makeCDR3Graph(cells, plate, true); + GraphWithMapData data = Simulator.makeGraph(cellSample, plate, true); assert filename != null; if(BiGpairSEQ.outputBinary()) { GraphDataObjectWriter dataWriter = new GraphDataObjectWriter(filename, data); @@ -378,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); @@ -507,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 { @@ -573,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(); 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)"); } } diff --git a/src/main/java/Plate.java b/src/main/java/Plate.java index 7f70c43..83d1c3b 100644 --- a/src/main/java/Plate.java +++ b/src/main/java/Plate.java @@ -8,7 +8,9 @@ TODO: Implement discrete frequency distributions using Vose's Alias Method import java.util.*; public class Plate { + private CellSample cells; private String sourceFile; + private String filename; private List> wells; private final Random rand = BiGpairSEQ.getRand(); private int size; @@ -18,6 +20,25 @@ public class Plate { private double lambda; boolean exponential = false; + public Plate(CellSample cells, String cellFilename, int numWells, Integer[] populations, + double dropoutRate, double stdDev_or_lambda, boolean exponential){ + this.cells = cells; + this.sourceFile = cellFilename; + this.size = numWells; + this.wells = new ArrayList<>(); + 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; @@ -26,8 +47,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(); @@ -43,10 +65,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 +95,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; @@ -160,4 +180,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 b673ba3..6ba2cbf 100644 --- a/src/main/java/Simulator.java +++ b/src/main/java/Simulator.java @@ -19,25 +19,19 @@ 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 { - //Replaced with SequenceType ordinals -// //These are the indices of the different sequences within a cell array -// 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 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 makeCDR3Graph(List distinctCells, Plate samplePlate, boolean verbose) { + public static GraphWithMapData makeGraph(CellSample cellSample, Plate samplePlate, boolean verbose) { Instant start = Instant.now(); - - int numWells = samplePlate.getSize(); - //The ordinal value of the sequence type enum is also that sequence's index in a cell array + List distinctCells = cellSample.getCells(); 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, 0, 1); @@ -53,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 sequences present in all wells.");} -// filterByOccupancyThresholds(allAlphas, 1, numWells - 1); -// filterByOccupancyThresholds(allBetas, 1, 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(); @@ -130,7 +125,7 @@ public class Simulator implements GraphModificationFunctions { //create GraphWithMapData object GraphWithMapData output = new GraphWithMapData(graph, numWells, samplePlate.getPopulations(), distCellsMapAlphaKey, 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; } @@ -266,17 +261,19 @@ 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); //rate of pairing error double pairingErrorRate = (double) falseCount / (trueCount + falseCount); BigDecimal pairingErrorRateTrunc; - if(pairingErrorRate == NaN || pairingErrorRate == POSITIVE_INFINITY || pairingErrorRate == NEGATIVE_INFINITY) { - pairingErrorRateTrunc = new BigDecimal(-1, mc); + if(Double.isFinite(pairingErrorRate)) { + pairingErrorRateTrunc = new BigDecimal(pairingErrorRate, mc); } else{ - pairingErrorRateTrunc = new BigDecimal(pairingErrorRate, mc); + pairingErrorRateTrunc = new BigDecimal(-1, mc); } //get list of well populations Integer[] wellPopulations = data.getWellPopulations(); @@ -296,6 +293,7 @@ public class Simulator implements GraphModificationFunctions { metadata.put("sample plate filename", data.getSourceFilename()); metadata.put("graph filename", dataFilename); metadata.put("algorithm type", algoType); + metadata.put("matching weight", totalMatchingWeight.toString()); metadata.put("well populations", wellPopulationsString); metadata.put("total alphas found", alphaCount.toString()); metadata.put("total betas found", betaCount.toString()); @@ -307,7 +305,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){ @@ -681,6 +679,7 @@ public class Simulator implements GraphModificationFunctions { } } } + } } @@ -695,7 +694,7 @@ public class Simulator implements GraphModificationFunctions { private static Map makeVertexToSequenceMap(Map sequences, Integer startValue) { Map map = new LinkedHashMap<>(); //LinkedHashMap to preserve order of entry - Integer index = startValue; + Integer index = startValue; //is this necessary? I don't think I use this. for (Integer k: sequences.keySet()) { map.put(index, k); index++;