Remove singleton sequences from wells in initial filtering

readme.md

@@ -20,8 +20,8 @@ The problem of pairing TCRA/TCRB sequences thus reduces to the "assignment probl
 matching on a bipartite graph--the subset of vertex-disjoint edges whose weights sum to the maximum possible value.
 
 This is a well-studied combinatorial optimization problem, with many known solutions.
-The most efficient algorithm known to the author for maximum weight matching of a bipartite graph with strictly integral weights
-is from Duan and Su (2012). For a graph with m edges, n vertices per side, and maximum integer edge weight N, 
+The most efficient algorithm known to the author for maximum weight matching of a bipartite graph with strictly integral
+weights is from Duan and Su (2012). For a graph with m edges, n vertices per side, and maximum integer edge weight N, 
 their algorithm runs in **O(m sqrt(n) log(N))** time. As the graph representation of a pairSEQ experiment is 
 bipartite with integer weights, this algorithm is ideal for BiGpairSEQ.
 
@@ -43,13 +43,13 @@ Run with the command:
 `java -jar BiGpairSEQ_Sim.jar`
 
 Processing sample plates with tens of thousands of sequences may require large amounts 
-of RAM. It is often desirable to increase the JVM maximum heap allocation with the -Xmx flag.
+of RAM. It is often desirable to increase the JVM maximum heap allocation with the `-Xmx` flag.
 For example, to run the program with 32 gigabytes of memory, use the command:
 
 `java -Xmx32G -jar BiGpairSEQ_Sim.jar`
 
 There are a number of command line options, to allow the program to be used in shell scripts. For a full list,
-use the -help flag:
+use the `-help` flag:
 
 `java -jar BiGpairSEQ_Sim.jar -help`
 
@@ -108,7 +108,7 @@ device-specific.)
 
 The program's caching behavior can be controlled in the Options menu. By default, all caching is OFF.
 
-The program can optionally output Graph/Data files in .GraphML format (.graphml) for data portability. This can be 
+The program can optionally output Graph/Data files in GraphML format (.graphml) for data portability. This can be 
 turned on in the Options menu. By default, GraphML output is OFF.
 
 ---
@@ -203,8 +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. This will produce a .graphml file
-for the weighted graph, with vertex attributes sequence, type, and occupancy data.)
+*Optional GraphML output*
+
+For portability of graph data to other software, turn on [GraphML](http://graphml.graphdrawing.org/index.html) output 
+in the Options menu in interactive mode, or use the `-graphml`command line argument. This will produce a .graphml file 
+for the weighted graph, with vertex attributes for sequence, type, and occupancy data. This graph contains all the data 
+necessary for the BiGpairSEQ matching algorithm. It does not include the data to measure pairing accuracy; for that, 
+compare the matching results to the original Cell Sample .csv file.
 
 ---
 
@@ -259,29 +264,77 @@ 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.
+  * 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?~~
@@ -295,6 +348,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
@@ -314,7 +368,7 @@ slightly less time than the simulation itself. Real elapsed time from start to f
 * [JGraphT](https://jgrapht.org) -- Graph theory data structures and algorithms
 * [JHeaps](https://www.jheaps.org) -- For pairing heap priority queue used in maximum weight matching algorithm
 * [Apache Commons CSV](https://commons.apache.org/proper/commons-csv/) -- For CSV file output
-* [Apache Commons CLI](https://commons.apache.org/proper/commons-cli/) -- To enable command line arguments for scripting. (**Awaiting re-implementation**.)
+* [Apache Commons CLI](https://commons.apache.org/proper/commons-cli/) -- To enable command line arguments for scripting.
 
 ## ACKNOWLEDGEMENTS
 BiGpairSEQ was conceived in collaboration with Dr. Alice MacQueen, who brought the original 

src/main/java/BiGpairSEQ.java

@@ -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; }
 }

src/main/java/CommandLineInterface.java

@@ -62,15 +62,18 @@ public class CommandLineInterface {
 
             if (line.hasOption("help")) {
                 HelpFormatter formatter = new HelpFormatter();
-                formatter.printHelp("BiGpairSEQ_Sim", mainOptions);
+                formatter.printHelp("BiGpairSEQ_Sim.jar", mainOptions);
                 System.out.println();
-                formatter.printHelp("BiGpairSEQ_SIM -cells", cellOptions);
+                formatter.printHelp("BiGpairSEQ_Sim.jar -cells", cellOptions);
                 System.out.println();
-                formatter.printHelp("BiGpairSEQ_Sim -plate", plateOptions);
+                formatter.printHelp("BiGpairSEQ_Sim.jar -plate", plateOptions);
                 System.out.println();
-                formatter.printHelp("BiGpairSEQ_Sim -graph", graphOptions);
+                formatter.printHelp("BiGpairSEQ_Sim.jar -graph", graphOptions);
                 System.out.println();
-                formatter.printHelp("BiGpairSEQ_Sim -match", matchOptions);
+                formatter.printHelp("BiGpairSEQ_Sim.jar -match", matchOptions);
+            }
+            else if (line.hasOption("version")) {
+                System.out.println("BiGpairSEQ_Sim " + BiGpairSEQ.getVersion());
             }
             else if (line.hasOption("cells")) {
                 line = parser.parse(cellOptions, Arrays.copyOfRange(args, 1, args.length));
@@ -153,17 +156,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 +183,38 @@ public class CommandLineInterface {
                 GraphWithMapData graph = getGraph(graphFilename);
                 MatchingResult result = Simulator.matchCDR3s(graph, graphFilename, minThreshold, maxThreshold,
                         maxOccupancyDiff, minOverlapPct, false);
-                MatchingFileWriter writer = new MatchingFileWriter(outputFilename, result);
-                writer.writeResultsToFile();
+                if(outputFilename != null){
+                    MatchingFileWriter writer = new MatchingFileWriter(outputFilename, result);
+                    writer.writeResultsToFile();
+                }
                 //can put a bunch of ifs for outputting various things from the MatchingResult to System.out here
                 //after I put those flags in the matchOptions
+                if(line.hasOption("print-metadata")) {
+                    for (String k : result.getMetadata().keySet()) {
+                        System.out.println(k + ": " + result.getMetadata().get(k));
+                    }
+                }
+                if(line.hasOption("print-error")) {
+                    System.out.println("pairing error rate: " + result.getPairingErrorRate());
+                }
+                if(line.hasOption("print-attempt")) {
+                    System.out.println("pairing attempt rate: " +result.getPairingAttemptRate());
+                }
+                if(line.hasOption("print-correct")) {
+                    System.out.println("correct pairings: " + result.getCorrectPairingCount());
+                }
+                if(line.hasOption("print-incorrect")) {
+                    System.out.println("incorrect pairings: " + result.getIncorrectPairingCount());
+                }
+                if(line.hasOption("print-alphas")) {
+                    System.out.println("total alphas found: " + result.getAlphaCount());
+                }
+                if(line.hasOption("print-betas")) {
+                    System.out.println("total betas found: " + result.getBetaCount());
+                }
+                if(line.hasOption("print-time")) {
+                    System.out.println("simulation time (seconds): " + result.getSimulationTime());
+                }
             }
         }
         catch (ParseException exp) {
@@ -216,8 +254,11 @@ public class CommandLineInterface {
                 .longOpt("match-cdr3")
                 .desc("Matches CDR3s. Requires a graph/data file.")
                 .build();
+        Option printVersion = Option.builder("version")
+                .desc("Prints the program version number to stdout").build();
         OptionGroup mainGroup = new OptionGroup();
         mainGroup.addOption(help);
+        mainGroup.addOption(printVersion);
         mainGroup.addOption(makeCells);
         mainGroup.addOption(makePlate);
         mainGroup.addOption(makeGraph);
@@ -297,7 +338,7 @@ public class CommandLineInterface {
                 .desc("Randomize well populations on sample plate. Takes two arguments: the minimum possible population and the maximum possible population.")
                 .hasArgs()
                 .numberOfArgs(2)
-                .argName("minimum maximum")
+                .argName("min> <max")
                 .build();
         Option specificWellPopulations = Option.builder("pop")
                 .desc("The well populations for each section of the sample plate. There will be as many sections as there are populations given.")
@@ -326,22 +367,22 @@ public class CommandLineInterface {
         Options graphOptions = new Options();
         Option cellFilename = Option.builder("c")
                 .longOpt("cell-file")
-                .desc("Cell sample file to use for checking accuracy")
+                .desc("Cell sample file to use for checking pairing accuracy")
                 .hasArg()
                 .argName("filename")
                 .required().build();
         Option plateFilename = Option.builder("p")
                 .longOpt("plate-filename")
-                .desc("Sample plate file (made from given cell sample file) to construct graph from")
+                .desc("Sample plate file from which to construct graph")
                 .hasArg()
                 .argName("filename")
                 .required().build();
         Option outputGraphML = Option.builder("graphml")
-                .desc("Output GraphML file")
+                .desc("(Optional) Output GraphML file")
                 .build();
         Option outputSerializedBinary = Option.builder("nb")
                 .longOpt("no-binary")
-                .desc("Don't output serialized binary file")
+                .desc("(Optional) Don't output serialized binary file")
                 .build();
         graphOptions.addOption(cellFilename);
         graphOptions.addOption(plateFilename);
@@ -379,15 +420,46 @@ public class CommandLineInterface {
                 .hasArg()
                 .argName("number")
                 .build();
-        matchCDR3options.addOption(graphFilename);
-        matchCDR3options.addOption(minOccupancyOverlap);
-        matchCDR3options.addOption(maxOccupancyOverlap);
-        matchCDR3options.addOption(minOverlapPercent);
-        matchCDR3options.addOption(maxOccupancyDifference);
-        matchCDR3options.addOption(outputFileOption());
-        //options for output to System.out
-       //Option printPairingErrorRate = Option.builder()
+        Option outputFile = Option.builder("o") //can't call the method this time, because this one's optional
+                .longOpt("output-file")
+                .hasArg()
+                .argName("filename")
+                .desc("(Optional) Name of output the output file. If not present, no file will be written.")
+                .build();
+        matchCDR3options.addOption(graphFilename)
+                .addOption(minOccupancyOverlap)
+                .addOption(maxOccupancyOverlap)
+                .addOption(minOverlapPercent)
+                .addOption(maxOccupancyDifference)
+                .addOption(outputFile);
 
+        //options for output to System.out
+        Option printAlphaCount = Option.builder().longOpt("print-alphas")
+                .desc("(Optional) Print the number of distinct alpha sequences to stdout.").build();
+        Option printBetaCount = Option.builder().longOpt("print-betas")
+                .desc("(Optional) Print the number of distinct beta sequences to stdout.").build();
+        Option printTime = Option.builder().longOpt("print-time")
+                .desc("(Optional) Print the total simulation time to stdout.").build();
+        Option printErrorRate = Option.builder().longOpt("print-error")
+               .desc("(Optional) Print the pairing error rate to stdout").build();
+        Option printAttempt = Option.builder().longOpt("print-attempt")
+               .desc("(Optional) Print the pairing attempt rate to stdout").build();
+        Option printCorrect = Option.builder().longOpt("print-correct")
+               .desc("(Optional) Print the number of correct pairs to stdout").build();
+        Option printIncorrect = Option.builder().longOpt("print-incorrect")
+               .desc("(Optional) Print the number of incorrect pairs to stdout").build();
+        Option printMetadata = Option.builder().longOpt("print-metadata")
+                       .desc("(Optional) Print a full summary of the matching results to stdout.").build();
+
+       matchCDR3options
+               .addOption(printErrorRate)
+               .addOption(printAttempt)
+               .addOption(printCorrect)
+               .addOption(printIncorrect)
+               .addOption(printMetadata)
+               .addOption(printAlphaCount)
+               .addOption(printBetaCount)
+               .addOption(printTime);
         return matchCDR3options;
     }
 

src/main/java/GraphDataObjectReader.java

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

src/main/java/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();

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);}
@@ -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
@@ -245,11 +244,11 @@ 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) {
-            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();