Update performance section and TODO

readme.md

@@ -1,5 +1,25 @@
 # BiGpairSEQ SIMULATOR
 
+## CONTENTS
+1. ABOUT
+2. THEORY
+3. THE BiGpairSEQ ALGORITHM
+4. USAGE
+   1. RUNNING THE PROGRAM
+   2. COMMAND LINE OPTIONS
+   3. INTERACTIVE INTERFACE
+   4. INPUT/OUTPUT
+      1. Cell Sample Files
+      2. Sample Plate Files
+      3. Graph/Data Files
+      4. Matching Results Files
+5. PERFORMANCE (needs revision!)
+   1. SIMULATING EXPERIMENTS FROM pairSEQ PAPER
+   2. BEHAVIOR WITH RANDOMIZED WELL POPULATIONS
+6. TODO
+7. CITATIONS
+8. ACKNOWLEDGEMENTS
+9. AUTHOR
 
 ## ABOUT
 
@@ -8,26 +28,77 @@ of the pairSEQ algorithm (Howie, et al. 2015) for pairing T cell receptor sequen
 
 ## THEORY
 
-Unlike pairSEQ, which calculates p-values for every TCR alpha/beta overlap and compares
-against a null distribution, BiGpairSEQ does not do any statistical calculations
-directly.
+T cell receptors (TCRs) are encoded by pairs of sequences, alpha sequences (TCRAs) and beta sequences (TCRBs). These sequences
+are extremely diverse; to the first approximation, this pair of sequences uniquely identifies a line of T cells. 
+
+As described in the original 2015 paper, pairSEQ pairs TCRAs and TCRBs by distributing a 
+sample of T cells across a 96-well sample plate, then sequencing the contents of each well. It then calculates p-values for 
+every TCRA/TCRB sequence overlap and compares that against a null distribution, to find the most statistically probable pairings.
+
+BiGpairSEQ uses the same fundamental idea of using occupancy overlap to pair TCR sequences, but unlike pairSEQ it
+does not require performing any statistical calculations at all. Instead, BiGpairSEQ uses graph theory methods which
+produce provably optimal solutions.
 
 BiGpairSEQ creates a [weighted bipartite graph](https://en.wikipedia.org/wiki/Bipartite_graph) representing the sample plate.
-The distinct TCRA and TCRB sequences form the two sets of vertices. Every TCRA/TCRB pair that share a well
-are connected by an edge, with the edge weight set to the number of wells in which both sequences appear.
-(Sequences present in *all* wells are filtered out prior to creating the graph, as there is no signal in their occupancy pattern.)
-The problem of pairing TCRA/TCRB sequences thus reduces to the "assignment problem" of finding a maximum weight
-matching on a bipartite graph--the subset of vertex-disjoint edges whose weights sum to the maximum possible value.
+The distinct TCRA and TCRB sequences form the two sets of vertices. Every TCRA/TCRB pair that share a well on the sample plate
+are connected by an edge in the graph, with the edge weight set to the number of wells in which both sequences appear. The vertices
+themselves are labeled with the occupancy data for the individual sequences they represent, which is useful for pre-filtering
+before finding a maximum weight matching. Such a graph fully encodes the distribution data from the sample plate.
 
-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, 
-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.
+The problem of pairing TCRA/TCRB sequences thus reduces to the [assignment problem](https://en.wikipedia.org/wiki/Assignment_problem) of finding a maximum weight
+matching (MWM) on a bipartite graph--the subset of vertex-disjoint edges whose weights sum to the maximum possible value.
 
-Unfortunately, it's a fairly new algorithm, and not yet implemented by the graph theory library used in this simulator.
-So this program instead uses the Fibonacci heap-based algorithm of Fredman and Tarjan (1987), which has a worst-case
-runtime of **O(n (n log(n) + m))**. The algorithm is implemented as described in Melhorn and Näher (1999).
+This is a well-studied combinatorial optimization problem, with many known algorithms that produce 
+provably-optimal solutions. The most theoretically 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 seems ideal for BiGpairSEQ. Unfortunately, it's a 
+fairly new algorithm, and not yet implemented by the graph theory library used in this simulator (JGraphT), nor has the author had 
+time to implement it himself.
+
+There have been some studies which show that [auction algorithms](https://en.wikipedia.org/wiki/Auction_algorithm) for the assignment problem can have superior performance in 
+real-world implementations, due to their simplicity, than more complex algorithms with better theoretical asymptotic 
+performance. But, again, there is no such algorithms implemented by JGraphT, nor has the author yet had time to implement one.
+
+So this program instead uses the [Fibonacci heap](https://en.wikipedia.org/wiki/Fibonacci_heap) based algorithm of Fredman and Tarjan (1987) (essentially
+[the Hungarian algorithm](https://en.wikipedia.org/wiki/Hungarian_algorithm) augmented with a more efficeint priority queue) which has a worst-case
+runtime of **O(n (n log(n) + m))**. The algorithm is implemented as described in Melhorn and Näher (1999). (The simulator
+allows the substitution of a [pairing heap](https://en.wikipedia.org/wiki/Pairing_heap) for a Fibonacci heap, though the relative performance difference of the two
+has not yet been thoroughly tested.)
+
+One possible advantage of this less efficient algorithm is that the Hungarian algorithm and its variations work with both the balanced and the unbalanced assignment problem
+(that is, cases where both sides of the bipartite graph have the same number of vertices and those in which they don't.)
+Many other MWM algorithms only work for the balanced assignment problem. While pairSEQ-style experiments should theoretically
+be balanced assignment problems, in practice sequence dropout can cause them to be unbalanced. The unbalanced case 
+*can* be reduced to the balanced case, but doing so involves doubling the graph size. Since the current implementation uses only
+the Hungarian algorithm, graph doubling--which could be challenging with the computational resources available to the 
+author--has not yet been necessary.
+
+The relative time/space efficiencies of BiGpairSEQ when backed by different MWM algorithms remains an open problem.
+
+## THE BiGpairSEQ ALGORITHM
+
+1. Sequence a sample plate of T cells as in pairSEQ.
+2. Pre-filter the sequence data to reduce error and minimize the size of the necessary graph.
+   1. *Saturating sequence filter*: remove any sequences present in all wells on the sample plate, as there is no signal in the occupancy data of saturating sequences (and each saturating sequence will have an edge to every vertex on the opposite side of the graph, vastly increasing the total graph size).
+   2. *Non-existent sequence filter*: sequencing misreads can pollute the data from the sample plate with non-existent sequences. These can be identified by the discrepancy between their occupancy and their total read count. Assuming sequences are read correctly at least half the time, then a sequence's  total read count (R) should be at least half the well occupancy of that sequence (O) times the read depth of the sequencing run (D). Remove any sequences for which R < (O * D) / 2.
+   3. *Misidentified sequence filter*: sequencing misreads can cause one real sequence to be misidentified as a different real sequence. This should be fairly infrequent, but is a problem if it skews a sequence's overall occupancy pattern by causing the sequence to seem to be in a well where it's not. This can be detected by looking for discrepancies in a sequence's per-well read count. On average, the read count for a sequence in an individual well (r) should be equal to its total read count (R) divided by its total well occupancy (O). Remove from the list of wells occupied by a sequence any wells for which r < R / (2 * O).
+3. Encode the occupancy data from the sample plate as a weighted bipartite graph, where one set of vertices represent the distinct TCRAs and the other set represents distinct TCRBs. Between any TCRA and TCRB that share a well, draw an edge. Assign that edge a weight equal to the total number of wells shared by both sequences.
+4. Find a maximum weight matching of the bipartite graph, using any [MWM algorithm](https://en.wikipedia.org/wiki/Assignment_problem#Algorithms) that produces a provably optimal result.
+    * If desired, restrict the matching to a subset of the graph. (Example: restricting matching attempts to cases where the occupancy overlap is 4 or more wells--that is, edges with weight >= 4.0.) See below for discussion of why this might be desirable.
+5. The resultant matching represents the likeliest TCRA/TCRB sequence pairs based on the occupancy pattern of the sample plate.
+
+It is important to note that a maximum weight matching is not necessarily unique. If two different sets of vertex-disjoint edges
+sum to the same maximal weight, then a MWM algorithms might find either one of them.
+
+For example, consider a well that contains four rare sequences found only in that well, two TCRAs and two TCRBs. 
+In the graph, both of those TCRAs would have edges to both TCRBs (and to others of course, but since those edges will have a weight of 1.0, 
+they are unlikely be paired in a MWM to sequences with total occupancy of more than one well). If these four sequences 
+represent two unique T cells, then only one of the two possible pairings between these sequences is correct. But both 
+the correct and incorrect pairing will add 2.0 to the total graph weight, so either one could be part of a maximum weight matching.
+
+It is to minimize the number of possible equivalent-weight matchings that one might restrict the algorithm to examining
+only a subset of the graph, as described in step 4 above.
 
 ## USAGE
 
@@ -48,11 +119,117 @@ 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:
+### COMMAND LINE OPTIONS
+
+There are a number of command line options, to allow the program to be used in shell scripts. These can be viewed with
+ the `-help` flag:
 
 `java -jar BiGpairSEQ_Sim.jar -help`
 
+```
+usage: BiGpairSEQ_Sim.jar
+ -cells,--make-cells   Makes a cell sample file of distinct T cells
+ -graph,--make-graph   Makes a graph/data file. Requires a cell sample
+                       file and a sample plate file
+ -help                 Displays this help menu
+ -match,--match-cdr3   Matches CDR3s. Requires a graph/data file.
+ -plate,--make-plate   Makes a sample plate file. Requires a cell sample
+                       file.
+ -version              Prints the program version number to stdout
+
+usage: BiGpairSEQ_Sim.jar -cells
+ -d,--diversity-factor <factor>   The factor by which unique CDR3s
+                                  outnumber unique CDR1s
+ -n,--num-cells <number>          The number of distinct cells to generate
+ -o,--output-file <filename>      Name of output file
+
+usage: BiGpairSEQ_Sim.jar -plate
+ -c,--cell-file <filename>     The cell sample file to use
+ -d,--dropout-rate <rate>      The sequence dropout rate due to
+                               amplification error. (0.0 - 1.0)
+ -exponential                  Use an exponential distribution for cell
+                               sample
+ -gaussian                     Use a Gaussian distribution for cell sample
+ -lambda <value>               If using -exponential flag, lambda value
+                               for distribution
+ -o,--output-file <filename>   Name of output file
+ -poisson                      Use a Poisson distribution for cell sample
+ -pop <number [number]...>     The well populations for each section of
+                               the sample plate. There will be as many
+                               sections as there are populations given.
+ -random <min> <max>           Randomize well populations on sample plate.
+                               Takes two arguments: the minimum possible
+                               population and the maximum possible
+                               population.
+ -stddev <value>               If using -gaussian flag, standard deviation
+                               for distrbution
+ -w,--wells <number>           The number of wells on the sample plate
+
+usage: BiGpairSEQ_Sim.jar -graph
+ -c,--cell-file <filename>                Cell sample file to use for
+                                          checking pairing accuracy
+ -err,--read-error-prob <prob>            (Optional) The probability that
+                                          a sequence will be misread. (0.0
+                                          - 1.0)
+ -errcoll,--error-collision-prob <prob>   (Optional) The probability that
+                                          two misreads will produce the
+                                          same spurious sequence. (0.0 -
+                                          1.0)
+ -graphml                                 (Optional) Output GraphML file
+ -nb,--no-binary                          (Optional) Don't output
+                                          serialized binary file
+ -o,--output-file <filename>              Name of output file
+ -p,--plate-filename <filename>           Sample plate file from which to
+                                          construct graph
+ -rd,--read-depth <depth>                 (Optional) The number of times
+                                          to read each sequence.
+ -realcoll,--real-collision-prob <prob>   (Optional) The probability that
+                                          a sequence will be misread as
+                                          another real sequence. (Only
+                                          applies to unique misreads;
+                                          after this has happened once,
+                                          future error collisions could
+                                          produce the real sequence again)
+                                          (0.0 - 1.0)
+
+usage: BiGpairSEQ_Sim.jar -match
+ -g,--graph-file <filename>    The graph/data file to use
+ -max <number>                 The maximum number of shared wells to
+                               attempt to match a sequence pair
+ -maxdiff <number>             (Optional) The maximum difference in total
+                               occupancy between two sequences to attempt
+                               matching.
+ -min <number>                 The minimum number of shared wells to
+                               attempt to match a sequence pair
+ -minpct <percent>             (Optional) The minimum percentage of a
+                               sequence's total occupancy shared by
+                               another sequence to attempt matching. (0 -
+                               100)
+ -o,--output-file <filename>   (Optional) Name of output the output file.
+                               If not present, no file will be written.
+    --print-alphas             (Optional) Print the number of distinct
+                               alpha sequences to stdout.
+    --print-attempt            (Optional) Print the pairing attempt rate
+                               to stdout
+    --print-betas              (Optional) Print the number of distinct
+                               beta sequences to stdout.
+    --print-correct            (Optional) Print the number of correct
+                               pairs to stdout
+    --print-error              (Optional) Print the pairing error rate to
+                               stdout
+    --print-incorrect          (Optional) Print the number of incorrect
+                               pairs to stdout
+    --print-metadata           (Optional) Print a full summary of the
+                               matching results to stdout.
+    --print-time               (Optional) Print the total simulation time
+                               to stdout.
+ -pv,--p-value                 (Optional) Calculate p-values for sequence
+                               pairs.
+
+```
+
+### INTERACTIVE INTERFACE
+
 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:
 
@@ -79,7 +256,8 @@ By default, the Options menu looks like this:
 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
+6) Turn on calculation of p-values
+7) Maximum weight matching algorithm options
 0) Return to main menu
 ```
 
@@ -116,7 +294,7 @@ turned on in the Options menu. By default, GraphML output is OFF.
 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
 random integers. CDR3 Alpha and Beta sequences are all unique within a given Cell Sample file. CDR1 Alpha and Beta sequences
-are not necessarily unique; the relative diversity can be set when making the file.
+are not necessarily unique; the relative diversity of CRD1s with respect to CDR3s can be set when making the file.
 
 (Note: though cells still have CDR1 sequences, matching of CDR1s is currently awaiting re-implementation.)
 
@@ -133,7 +311,7 @@ Structure:
 | Alpha CDR3 | Beta CDR3 | Alpha CDR1 | Beta CDR1 |
 |---|---|---|---|
 |unique number|unique number|number|number|
-
+| ... | ...  |... | ... |
 ---
 
 #### Sample Plate Files
@@ -142,7 +320,8 @@ described above). The wells are filled randomly from a Cell Sample file, accordi
 frequency distribution. Additionally, every individual sequence within each cell may, with some
 given dropout probability, be omitted from the file; this simulates the effect of amplification errors
 prior to sequencing. Plates can also be partitioned into any number of sections, each of which can have a 
-different concentration of T cells per well.
+different concentration of T cells per well. Alternatively, the number of T cells in each well can be randomized between 
+given minimum and maximum population values.
 
 Options when making a Sample Plate file:
 * Cell Sample file to use
@@ -152,7 +331,6 @@ Options when making a Sample Plate file:
     * Standard deviation size 
   * Exponential
     * Lambda value
-      * *(Based on the slope of the graph in Figure 4C of the pairSEQ paper, the distribution of the original experiment was approximately exponential with a lambda ~0.6. (Howie, et al. 2015))*
 * Total number of wells on the plate
 * Well populations random or fixed
   * If random, minimum and maximum population sizes
@@ -199,12 +377,12 @@ then use it for multiple different BiGpairSEQ simulations.
 
 Options for creating a Graph/Data file:
 * The Cell Sample file to use
-* The Sample Plate file to use. (This must have been generated from the selected Cell Sample file.)
-* Whether to simulate sequence read depth. If simulated:
-  * The read depth (number of times each sequence is read)
-  * The read error rate (probability a sequence is misread)
-  * The error collision rate (probability two misreads produce the same spurious sequence)
-  * The real sequence collision rate (probability that a misread will produce a different, real sequence from the sample plate. Only applies to new misreads; once an error of this type has occurred, it's likelihood of ocurring again is dominated by the error collision probability.)
+* The Sample Plate file to use (This must have been generated from the selected Cell Sample file.)
+* Whether to simulate sequencing read depth. If simulated:
+  * The read depth (The number of times each sequence is read)
+  * The read error rate (The probability a sequence is misread)
+  * The error collision rate (The probability two misreads produce the same spurious sequence)
+  * The real sequence collision rate (The probability that a misread will produce a different, real sequence from the sample plate. Only applies to new misreads; once an error of this type has occurred, it's likelihood of occurring again is dominated by the error collision probability.)
 
 These files do not have a human-readable structure, and are not portable to other programs.
 
@@ -221,7 +399,7 @@ 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 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.)
+output only. Graph 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 `#`.
@@ -239,56 +417,65 @@ 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
-# T cell counts in sample plate wells: 30000
-# Total alphas found: 11813
-# Total betas found: 11808
-# High overlap threshold: 94
-# Low overlap threshold: 3
-# Minimum overlap percent: 0
-# Maximum occupancy difference: 96
-# Pairing attempt rate: 0.438
-# Correct pairings: 5151
-# Incorrect pairings: 18
-# Pairing error rate: 0.00348
-# Simulation time: 862 seconds
+# cell sample filename: 8MilCells.csv
+# cell sample size: 8000000
+# sample plate filename: 8MilCells_Plate.csv
+# sample plate well count: 96
+# sequence dropout rate: 0.1
+# graph filename: 8MilGraph_rd2
+# MWM algorithm type: LEDA book with heap: FIBONACCI
+# matching weight: 218017.0
+# well populations: 4000
+# sequence read depth: 100
+# sequence read error rate: 0.01
+# read error collision rate: 0.1
+# real sequence collision rate: 0.05
+# total alphas read from plate: 323711
+# total betas read from plate: 323981
+# alphas in graph (after pre-filtering): 11707
+# betas in graph (after pre-filtering): 11705
+# high overlap threshold for pairing: 95
+# low overlap threshold for pairing: 3
+# minimum overlap percent for pairing: 0
+# maximum occupancy difference for pairing: 2147483647
+# pairing attempt rate: 0.716
+# correct pairing count: 8373
+# incorrect pairing count: 7
+# pairing error rate: 0.000835
+# time to generate graph (seconds): 293
+# time to pair sequences (seconds): 1,416
+# total simulation time (seconds): 1,709
 ```
 
 | Alpha | Alpha well count | Beta | Beta well count | Overlap count | Matched Correctly? | P-value |
 |---|---|---|---|---|---|---|
-|5242972|17|1571520|18|17|true|1.41E-18|
-|5161027|18|2072219|18|18|true|7.31E-20|
-|4145198|33|1064455|30|29|true|2.65E-21|
-|7700582|18|112748|18|18|true|7.31E-20|
+|10258642|73|1172093|72|70.0|true|4.19E-21|
+|6186865|34|4290363|37|34.0|true|4.56E-26|
+|10222686|70|11044018|72|68.0|true|9.55E-25|
+|5338100|75|2422988|76|74.0|true|4.57E-22|
+|12363907|33|6569852|35|33.0|true|5.70E-26|
 |...|...|...|...|...|...|...|
 
 ---
 
-**NOTE: The p-values in the output are not used for matching**—they aren't part of the BiGpairSEQ algorithm at all. 
-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)
+**NOTE: The p-values in the sample output above are not used for matching**—they aren't part of the BiGpairSEQ algorithm at all. 
+P-values (if enabled in the interactive menu options or by using the -pv flag on the command line) are calculated *after* 
+BiGpairSEQ matching is completed, for purposes of comparison with pairSEQ only, using the (corrected) formula from the 
+original pairSEQ paper. (Howie, et al. 2015) Calculation of p-values is off by default to reduce processing time.
 
 
-## PERFORMANCE (old results; need updating to reflect current, improved simulator performance)
+## 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 (lambda 0.6).
+Several BiGpairSEQ simulations were performed on a home computer with the following specs:
 
-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%.
+* Ryzen 5600X CPU
+* 128GB of 3200MHz DDR4 RAM
+* 2TB PCIe 3.0 SSD
+* Linux Mint 21 (5.15 kernel)
 
-The total simulation time was 14'22". If intermediate results were held in memory, this would be equivalent to the total elapsed time.
-
-Since this implementation of BiGpairSEQ writes intermediate results to disk (to improve the efficiency of *repeated* simulations
-with different filtering options), the actual elapsed time was greater. File I/O time was not measured, but took 
-slightly less time than the simulation itself. Real elapsed time from start to finish was under 30 minutes.
-
-As mentioned in the theory section, performance could be improved by implementing a more efficient algorithm for finding
-the maximum weight matching.
-
-## BEHAVIOR WITH RANDOMIZED WELL POPULATIONS
+### SAMPLE PLATES WITH VARYING NUMBERS OF CELLS PER WELL
+NOTE: these results were obtained with an earlier version of BiGpairSEQ_Sim, and should be re-run with the current version.
+The observed behavior is not believed to be likely to change, however.
 
 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
@@ -333,6 +520,70 @@ The average results for the randomized plates are closest to the constant plate
 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.
 
+### SIMULATING EXPERIMENTS FROM THE 2015 pairSEQ PAPER
+#### Experiment 1
+This simulation was an attempt to replicate the conditions of experiment 1 from the 2015 pairSEQ paper: a matching was found for a 
+96-well sample plate with 4,000 T cells/well, taken from a sample of 8,400,000 
+distinct cells sampled with an exponential frequency distribution. Examination of Figure 4C from the paper seems to show the points
+(-5, 4) and (-4.5, 3.3) on the line at the boundary of the shaded region, so a lambda value of 1.4 was used for the 
+exponential distribution.
+
+The sequence dropout rate was 10%, as the analysis in the paper concluded that most TCR 
+sequences "have less than a 10% chance of going unobserved." (Howie, et al. 2015) Given this choice of 10%, the simulated
+sample plate is likely to have more sequence dropout, and thus greater error, than the real experiment.
+
+The original paper does not contain (or the author of this document failed to identify) information on sequencing depth, 
+read error probability, or the probabilities of different kinds of read error collisions. As the pre-filtering of BiGpairSEQ
+has successfully filtered out all such errors for any reasonable error rates the author has yet tested, this simulation was
+done without simulating any sequencing errors, to reduce the processing time.
+
+This simulation was performed 5 times with min/max occupancy thresholds of 3 and 95 wells respectively for matching.
+
+| |Run 1|Run 2|Run 3|Run 4|Run 5| Average|
+|---|---|---|---|---|---|---|
+|Total pairs|4398|4420|4404|4409|4414|4409|
+|Correct pairs|4322|4313|4337|4336|4339|4329.4|
+|Incorrect pairs|76|107|67|73|75|79.6|
+|Error rate|0.0173|0.0242|0.0152|0.0166|0.0170|0.018|
+|Simulation time (seconds)|697|484|466|473|463|516.6|
+
+The experiment in the original paper called 4143 pairs with a false discovery rate of 0.01.
+
+Given the roughness of the estimation for the cell frequency distribution of the original experiment and the likely 
+higher rate of sequence dropout in the simulation, these simulated results match the real experiment fairly well.
+
+#### Experiment 3
+To simulate experiment 3 from the original paper, a matching was made for a 96-well sample plate with 160,000 T cells/well,
+taken from a sample of 4.5 million distinct T cells sampled with an exponential frequency distribution (lambda 1.4). The 
+sequence dropout rate was again 10%, and no sequencing errors were simulated. Once again, deviation from the original 
+experiment is expected due to the roughness of the estimated frequency distribution, and due to the high sequence dropout 
+rate.
+
+Results metadata:
+```
+# total alphas read from plate: 6929
+# total betas read from plate: 6939
+# alphas in graph (after pre-filtering): 4452
+# betas in graph (after pre-filtering): 4461
+# high overlap threshold for pairing: 95
+# low overlap threshold for pairing: 3
+# minimum overlap percent for pairing: 0
+# maximum occupancy difference for pairing: 100
+# pairing attempt rate: 0.767
+# correct pairing count: 3233
+# incorrect pairing count: 182
+# pairing error rate: 0.0533
+# time to generate graph (seconds): 40
+# time to pair sequences (seconds): 230
+# total simulation time (seconds): 270
+```
+
+The simulation ony found 6929 distinct TCRAs and 6939 TCRBs on the sample plate, orders of magnitude fewer than the number of
+pairs called in the pairSEQ experiment. These results show that at very high sampling depths, the differences in the 
+underlying frequency distribution drastically affect the results. The real distribution clearly has a much longer "tail" 
+than the simulated exponential distribution. Implementing a way to exert finer control over the sampling distribution from 
+the file of distinct cells may enable better simulated replication of this experiment.
+
 ## TODO
 
 * ~~Try invoking GC at end of workloads to reduce paging to disk~~ DONE
@@ -364,6 +615,8 @@ roughly as though it had a constant well population equal to the plate's average
   * Individual well data from the SequenceRecords could be included, if there's ever a reason for it
 * ~~Implement simulation of sequences being misread as other real sequence~~ DONE
 * Update matching metadata output options in CLI
+* Add frequency distribution details to metadata output
+  * need to make an enum for the different distribution types and refactor the Plate class and user interfaces, also add the necessary fields to GraphWithMapData and then call if from Simulator
 * Update performance data in this readme
 * Add section to ReadMe describing data filtering methods.
 * Re-implement CDR1 matching method
@@ -380,7 +633,7 @@ roughly as though it had a constant well population equal to the plate's average
 * Implement Vose's alias method for arbitrary statistical distributions of cells
   * Should probably refactor to use apache commons rng for this
 * Use commons JCS for caching
-* Parameterize pre-filtering. Currently, sequences present in all wells are filtered out before constructing the graph, which massively reduces graph size. But, ideally, no pre-filtering would be necessary.
+* Parameterize pre-filtering options
   
   
 ## CITATIONS
@@ -400,4 +653,4 @@ BiGpairSEQ was conceived in collaboration with Dr. Alice MacQueen, who brought t
 pairSEQ paper to the author's attention and explained all the biology terms he didn't know.
 
 ## AUTHOR
-BiGpairSEQ algorithm and simulation by Eugene Fischer, 2021. UI improvements and documentation, 2022.
+BiGpairSEQ algorithm and simulation by Eugene Fischer, 2021. Improvements and documentation, 2022.

src/main/java/BiGpairSEQ.java

@@ -16,7 +16,8 @@ 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 4.1";
+    private static boolean calculatePValue = false;
+    private static final String version = "version 4.2";
 
     public static void main(String[] args) {
         if (args.length == 0) {
@@ -173,5 +174,9 @@ public class BiGpairSEQ {
 
     public static boolean outputGraphML() {return outputGraphML;}
     public static void setOutputGraphML(boolean b) {outputGraphML = b;}
+
+    public static boolean calculatePValue() {return calculatePValue; }
+    public static void setCalculatePValue(boolean b) {calculatePValue = b; }
+
     public static String getVersion() { return version; }
 }

src/main/java/CellFileReader.java

@@ -58,7 +58,9 @@ public class CellFileReader {
     }
 
     public CellSample getCellSample() {
-        return new CellSample(distinctCells, cdr1Freq);
+        CellSample sample = new CellSample(distinctCells, cdr1Freq);
+        sample.setFilename(filename);
+        return sample;
     }
 
     public String getFilename() { return filename;}

src/main/java/CellSample.java

@@ -7,6 +7,7 @@ public class CellSample {
 
     private List<String[]> cells;
     private Integer cdr1Freq;
+    private String filename;
 
     public CellSample(Integer numDistinctCells, Integer cdr1Freq){
         this.cdr1Freq = cdr1Freq;
@@ -38,6 +39,7 @@ public class CellSample {
             distinctCells.add(tmp);
         }
         this.cells = distinctCells;
+        this.filename = filename;
     }
 
     public CellSample(List<String[]> cells, Integer cdr1Freq){
@@ -57,4 +59,8 @@ public class CellSample {
         return cells.size();
     }
 
+    public String getFilename() { return filename; }
+
+    public void setFilename(String filename) { this.filename = filename; }
+
 }

src/main/java/CommandLineInterface.java

@@ -48,6 +48,7 @@ import java.util.stream.Stream;
  * minpercent : (optional) the minimum percent overlap to attempt a matching.
  * writefile : (optional) the filename to write results to
  * output : the values to print to System.out for piping
+ * pv : (optional) calculate p-values
  *
  */
 public class CommandLineInterface {
@@ -96,7 +97,7 @@ public class CommandLineInterface {
                 Integer[] populations;
                 String outputFilename = line.getOptionValue("o");
                 Integer numWells = Integer.parseInt(line.getOptionValue("w"));
-                Double dropoutRate = Double.parseDouble(line.getOptionValue("err"));
+                Double dropoutRate = Double.parseDouble(line.getOptionValue("d"));
                 if (line.hasOption("random")) {
                     //Array holding values of minimum and maximum populations
                    Integer[] min_max = Stream.of(line.getOptionValues("random"))
@@ -202,9 +203,12 @@ public class CommandLineInterface {
                 else {
                     maxOccupancyDiff = Integer.MAX_VALUE;
                 }
+                if (line.hasOption("pv")) {
+                    BiGpairSEQ.setCalculatePValue(true);
+                }
                 GraphWithMapData graph = getGraph(graphFilename);
                 MatchingResult result = Simulator.matchCDR3s(graph, graphFilename, minThreshold, maxThreshold,
-                        maxOccupancyDiff, minOverlapPct, false);
+                        maxOccupancyDiff, minOverlapPct, false, BiGpairSEQ.calculatePValue());
                 if(outputFilename != null){
                     MatchingFileWriter writer = new MatchingFileWriter(outputFilename, result);
                     writer.writeResultsToFile();
@@ -367,7 +371,8 @@ public class CommandLineInterface {
                 .hasArgs()
                 .argName("number [number]...")
                 .build();
-        Option dropoutRate = Option.builder("err") //add this to plate options
+        Option dropoutRate = Option.builder("d") //add this to plate options
+                .longOpt("dropout-rate")
                 .hasArg()
                 .desc("The sequence dropout rate due to amplification error. (0.0 - 1.0)")
                 .argName("rate")
@@ -478,12 +483,17 @@ public class CommandLineInterface {
                 .argName("filename")
                 .desc("(Optional) Name of output the output file. If not present, no file will be written.")
                 .build();
+        Option pValue = Option.builder("pv") //can't call the method this time, because this one's optional
+                .longOpt("p-value")
+                .desc("(Optional) Calculate p-values for sequence pairs.")
+                .build();
         matchCDR3options.addOption(graphFilename)
                 .addOption(minOccupancyOverlap)
                 .addOption(maxOccupancyOverlap)
                 .addOption(minOverlapPercent)
                 .addOption(maxOccupancyDifference)
-                .addOption(outputFile);
+                .addOption(outputFile)
+                .addOption(pValue);
 
         //options for output to System.out
         Option printAlphaCount = Option.builder().longOpt("print-alphas")

src/main/java/GraphMLFileWriter.java

@@ -43,7 +43,7 @@ public class GraphMLFileWriter {
     private Map<String, Attribute> createGraphAttributes(){
         Map<String, Attribute> attributes = new HashMap<>();
         //Sample plate filename
-        attributes.put("sample plate filename", DefaultAttribute.createAttribute(data.getSourceFilename()));
+        attributes.put("sample plate filename", DefaultAttribute.createAttribute(data.getPlateFilename()));
         // Number of wells
         attributes.put("well count", DefaultAttribute.createAttribute(data.getNumWells().toString()));
         //Well populations

src/main/java/GraphWithMapData.java

@@ -9,12 +9,15 @@ import java.util.Map;
 //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;
+    private String cellFilename;
+    private int cellSampleSize;
+    private String plateFilename;
     private final SimpleWeightedGraph graph;
     private final int numWells;
     private final Integer[] wellPopulations;
     private final int alphaCount;
     private final int betaCount;
+    private final double dropoutRate;
     private final int readDepth;
     private final double readErrorRate;
     private final double errorCollisionRate;
@@ -30,7 +33,7 @@ public class GraphWithMapData implements java.io.Serializable {
 
     public GraphWithMapData(SimpleWeightedGraph graph, Integer numWells, Integer[] wellConcentrations,
                             Map<String, String> distCellsMapAlphaKey, Integer alphaCount, Integer betaCount,
-                            Integer readDepth, Double readErrorRate, Double errorCollisionRate,
+                            Double dropoutRate, Integer readDepth, Double readErrorRate, Double errorCollisionRate,
                             Double realSequenceCollisionRate, Duration time){
 
 //                            Map<Integer, Integer> plateVtoAMap,
@@ -49,6 +52,7 @@ public class GraphWithMapData implements java.io.Serializable {
 //        this.plateBtoVMap = plateBtoVMap;
 //        this.alphaWellCounts = alphaWellCounts;
 //        this.betaWellCounts = betaWellCounts;
+        this.dropoutRate = dropoutRate;
         this.readDepth = readDepth;
         this.readErrorRate = readErrorRate;
         this.errorCollisionRate = errorCollisionRate;
@@ -110,12 +114,20 @@ public class GraphWithMapData implements java.io.Serializable {
         return time;
     }
 
-    public void setSourceFilename(String filename) {
-        this.sourceFilename = filename;
+    public void setCellFilename(String filename) { this.cellFilename = filename; }
+
+    public String getCellFilename() { return this.cellFilename; }
+
+    public Integer getCellSampleSize() { return this.cellSampleSize; }
+
+    public void setCellSampleSize(int size) { this.cellSampleSize = size;}
+
+    public void setPlateFilename(String filename) {
+        this.plateFilename = filename;
     }
 
-    public String getSourceFilename() {
-        return sourceFilename;
+    public String getPlateFilename() {
+        return plateFilename;
     }
 
     public Double getReadErrorRate() {
@@ -127,4 +139,6 @@ public class GraphWithMapData implements java.io.Serializable {
     }
 
     public Double getRealSequenceCollisionRate() { return realSequenceCollisionRate; }
+
+    public Double getDropoutRate() { return dropoutRate; }
 }

src/main/java/InteractiveInterface.java

@@ -114,8 +114,8 @@ public class InteractiveInterface {
             System.out.println("1) Poisson");
             System.out.println("2) Gaussian");
             System.out.println("3) Exponential");
-            System.out.println("(Note: approximate distribution in original paper is exponential, lambda = 0.6)");
-            System.out.println("(lambda value approximated from slope of log-log graph in figure 4c)");
+//            System.out.println("(Note: approximate distribution in original paper is exponential, lambda = 0.6)");
+//            System.out.println("(lambda value approximated from slope of log-log graph in figure 4c)");
             System.out.println("(Note: wider distributions are more memory intensive to match)");
             System.out.print("Enter selection value: ");
             input = sc.nextInt();
@@ -422,7 +422,7 @@ public class InteractiveInterface {
         }
         //simulate matching
         MatchingResult results = Simulator.matchCDR3s(data, graphFilename, lowThreshold, highThreshold, maxOccupancyDiff,
-                minOverlapPercent, true);
+                minOverlapPercent, true, BiGpairSEQ.calculatePValue());
         //write results to file
         assert filename != null;
         MatchingFileWriter writer = new MatchingFileWriter(filename, results);
@@ -544,7 +544,8 @@ public class InteractiveInterface {
             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 (for data portability to other programs)");
-            System.out.println("6) Maximum weight matching algorithm options");
+            System.out.println("6) Turn " + getOnOff(!BiGpairSEQ.calculatePValue()) + " calculation of p-values");
+            System.out.println("7) Maximum weight matching algorithm options");
             System.out.println("0) Return to main menu");
             try {
                 input = sc.nextInt();
@@ -554,7 +555,8 @@ public class InteractiveInterface {
                     case 3 -> BiGpairSEQ.setCacheGraph(!BiGpairSEQ.cacheGraph());
                     case 4 -> BiGpairSEQ.setOutputBinary(!BiGpairSEQ.outputBinary());
                     case 5 -> BiGpairSEQ.setOutputGraphML(!BiGpairSEQ.outputGraphML());
-                    case 6 -> algorithmOptions();
+                    case 6 -> BiGpairSEQ.setCalculatePValue(!BiGpairSEQ.calculatePValue());
+                    case 7 -> algorithmOptions();
                     case 0 -> backToMain = true;
                     default -> System.out.println("Invalid input");
                 }

src/main/java/SequenceRecord.java

@@ -39,6 +39,11 @@ public class SequenceRecord implements Serializable {
         wells.put(wellNumber, readCount);
     }
 
+    //Method to remove a well from the occupancy map.
+    //Useful for cases where one sequence is misread as another sequence that isn't actually present in the well
+    //This can reveal itself as an anomalously low read count in that well.
+    public void deleteWell(Integer wellNumber) { wells.remove(wellNumber); }
+
     public Set<Integer> getWells() {
         return wells.keySet();
     }

src/main/java/Simulator.java

@@ -12,14 +12,6 @@ import java.text.NumberFormat;
 import java.time.Instant;
 import java.time.Duration;
 import java.util.*;
-/*
- Refactor notes
- What would be necessary to do everything with only one scan through the sample plate?
- I would need to keep a list of sequences (real and spurious), and metadata about each sequence.
- I would need the data:
- * # of each well the sequence appears in
- * Read count in that well
- */
 
 
 //NOTE: "sequence" in method and variable names refers to a peptide sequence from a simulated T cell
@@ -69,6 +61,14 @@ public class Simulator implements GraphModificationFunctions {
             if(verbose){System.out.println("Remaining alpha sequence count: " + alphaSequences.size());}
             if(verbose){System.out.println("Remaining beta sequence count: " + betaSequences.size());}
         }
+        if (realSequenceCollisionRate > 0.0) {
+            if(verbose){System.out.println("Removing wells with anomalous read counts from sequence records");}
+            int alphaWellsRemoved = filterWellsByReadCount(alphaSequences);
+            int betaWellsRemoved = filterWellsByReadCount(betaSequences);
+            if(verbose){System.out.println("Wells with anomalous read counts removed from sequence records");}
+            if(verbose){System.out.println("Total alpha sequence wells removed: " + alphaWellsRemoved);}
+            if(verbose){System.out.println("Total beta sequence wells removed: " + betaWellsRemoved);}
+        }
 
         //construct the graph. For simplicity, going to make
         if(verbose){System.out.println("Making vertex maps");}
@@ -128,9 +128,13 @@ public class Simulator implements GraphModificationFunctions {
         Duration time = Duration.between(start, stop);
         //create GraphWithMapData object
         GraphWithMapData output = new GraphWithMapData(graph, numWells, samplePlate.getPopulations(), distCellsMapAlphaKey,
-                alphaCount, betaCount, readDepth, readErrorRate, errorCollisionRate, realSequenceCollisionRate, time);
-        //Set source file name in graph to name of sample plate
-        output.setSourceFilename(samplePlate.getFilename());
+                alphaCount, betaCount, samplePlate.getError(), readDepth, readErrorRate, errorCollisionRate, realSequenceCollisionRate, time);
+        //Set cell sample file name in graph to name of cell sample
+        output.setCellFilename(cellSample.getFilename());
+        //Set cell sample size in graph
+        output.setCellSampleSize(cellSample.getCellCount());
+        //Set sample plate file name in graph to name of sample plate
+        output.setPlateFilename(samplePlate.getFilename());
         //return GraphWithMapData object
         return output;
     }
@@ -138,7 +142,7 @@ public class Simulator implements GraphModificationFunctions {
     //match CDR3s.
     public static MatchingResult matchCDR3s(GraphWithMapData data, String dataFilename, Integer lowThreshold,
                                             Integer highThreshold, Integer maxOccupancyDifference,
-                                            Integer minOverlapPercent, boolean verbose) {
+                                            Integer minOverlapPercent, boolean verbose, boolean calculatePValue) {
         Instant start = Instant.now();
         SimpleWeightedGraph<Vertex, DefaultWeightedEdge> graph = data.getGraph();
         Map<Vertex[], Integer> removedEdges = new HashMap<>();
@@ -216,7 +220,7 @@ public class Simulator implements GraphModificationFunctions {
         header.add("Beta well count");
         header.add("Overlap well count");
         header.add("Matched correctly?");
-        header.add("P-value");
+        if(calculatePValue) { header.add("P-value"); }
 
         //Results for csv file
         List<List<String>> allResults = new ArrayList<>();
@@ -253,10 +257,12 @@ public class Simulator implements GraphModificationFunctions {
             //overlap count
             result.add(Double.toString(graph.getEdgeWeight(e)));
             result.add(Boolean.toString(check));
-            double pValue = Equations.pValue(numWells, source.getOccupancy(),
+            if (calculatePValue) {
+                double pValue = Equations.pValue(numWells, source.getOccupancy(),
                     target.getOccupancy(), graph.getEdgeWeight(e));
-            BigDecimal pValueTrunc = new BigDecimal(pValue, mc);
-            result.add(pValueTrunc.toString());
+                BigDecimal pValueTrunc = new BigDecimal(pValue, mc);
+                result.add(pValueTrunc.toString());
+            }
             allResults.add(result);
         }
 
@@ -296,7 +302,11 @@ public class Simulator implements GraphModificationFunctions {
 
 
         Map<String, String> metadata = new LinkedHashMap<>();
-        metadata.put("sample plate filename", data.getSourceFilename());
+        metadata.put("cell sample filename", data.getCellFilename());
+        metadata.put("cell sample size", data.getCellSampleSize().toString());
+        metadata.put("sample plate filename", data.getPlateFilename());
+        metadata.put("sample plate well count", data.getNumWells().toString());
+        metadata.put("sequence dropout rate", data.getDropoutRate().toString());
         metadata.put("graph filename", dataFilename);
         metadata.put("MWM algorithm type", algoType);
         metadata.put("matching weight", totalMatchingWeight.toString());
@@ -307,10 +317,6 @@ public class Simulator implements GraphModificationFunctions {
         metadata.put("real sequence collision rate", data.getRealSequenceCollisionRate().toString());
         metadata.put("total alphas read from plate", data.getAlphaCount().toString());
         metadata.put("total betas read from plate", data.getBetaCount().toString());
-        //HARD CODED, PARAMETERIZE LATER
-        metadata.put("pre-filter sequences present in all wells", "true");
-        //HARD CODED, PARAMETERIZE LATER
-        metadata.put("pre-filter sequences based on occupancy/read count discrepancy", "true");
         metadata.put("alphas in graph (after pre-filtering)", graphAlphaCount.toString());
         metadata.put("betas in graph (after pre-filtering)", graphBetaCount.toString());
         metadata.put("high overlap threshold for pairing", highThreshold.toString());
@@ -647,7 +653,7 @@ public class Simulator implements GraphModificationFunctions {
 //    }
 
     //Remove sequences based on occupancy
-    public static void filterByOccupancyThresholds(Map<String, SequenceRecord> wellMap, int low, int high){
+    private static void filterByOccupancyThresholds(Map<String, SequenceRecord> wellMap, int low, int high){
         List<String> noise = new ArrayList<>();
         for(String k: wellMap.keySet()){
             if((wellMap.get(k).getOccupancy() > high) || (wellMap.get(k).getOccupancy() < low)){
@@ -659,10 +665,10 @@ public class Simulator implements GraphModificationFunctions {
         }
     }
 
-    public static void filterByOccupancyAndReadCount(Map<String, SequenceRecord> sequences, int readDepth) {
+    private static void filterByOccupancyAndReadCount(Map<String, SequenceRecord> sequences, int readDepth) {
         List<String> noise = new ArrayList<>();
         for(String k : sequences.keySet()){
-            //occupancy times read depth should be more than half the sequence read count if the read error rate is low
+            //the sequence read count should be more than half the occupancy times read depth if the read error rate is low
             Integer threshold = (sequences.get(k).getOccupancy() * readDepth) / 2;
             if(sequences.get(k).getReadCount() < threshold) {
                 noise.add(k);
@@ -673,6 +679,26 @@ public class Simulator implements GraphModificationFunctions {
         }
     }
 
+    private static int filterWellsByReadCount(Map<String, SequenceRecord> sequences) {
+        int count = 0;
+        for (String k: sequences.keySet()) {
+            //If a sequence has read count R and appears in W wells, then on average its read count in each
+            //well should be R/W. Delete any wells where the read count is less than R/2W.
+            Integer threshold = sequences.get(k).getReadCount() / (2 * sequences.get(k).getOccupancy());
+            List<Integer> noise = new ArrayList<>();
+            for (Integer well: sequences.get(k).getWells()) {
+                if (sequences.get(k).getReadCount(well) < threshold) {
+                    noise.add(well);
+                    count++;
+                }
+            }
+            for (Integer well: noise) {
+                sequences.get(k).deleteWell(well);
+            }
+        }
+        return count;
+    }
+
     private static Map<String, String> makeSequenceToSequenceMap(List<String[]> cells, int keySequenceIndex,
                                                                  int valueSequenceIndex){
         Map<String, String> keySequenceToValueSequenceMap = new HashMap<>();