BiGpairSEQ SIMULATOR
CONTENTS
- ABOUT
- THEORY
- THE BiGpairSEQ ALGORITHM
- USAGE
- RUNNING THE PROGRAM
- COMMAND LINE OPTIONS
- INTERACTIVE INTERFACE
- INPUT/OUTPUT
- Cell Sample Files
- Sample Plate Files
- Graph/Data Files
- Matching Results Files
- RESULTS
- SAMPLE PLATES WITH VARYING NUMBERS OF CELLS PER WELL
- SIMULATING EXPERIMENTS FROM pairSEQ PAPER
- TODO
- CITATIONS
- ACKNOWLEDGEMENTS
- AUTHOR
- DISCLOSURE
ABOUT
This program simulates BiGpairSEQ (Bipartite Graph pairSEQ), a graph theory-based adaptation of the pairSEQ algorithm (Howie, et al. 2015) for pairing T cell receptor sequences.
THEORY
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 representing the sample plate. 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.
The problem of pairing TCRA/TCRB sequences thus reduces to the 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.
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.
So this program instead uses the Fibonacci heap based algorithm of Fredman and Tarjan (1987) (essentially the 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 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.
There have been some studies which show that auction algorithms 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. The author has implemented a basic forward auction algorithm, which produces optimal assignment for unbalanced bipartite graphs with integer weights. To allow for unbalanced assignment, this algorithim eschews epsilon-scaling, and as a result is prone to "bidding-wars" which increase run time, making it less efficient than the implementation of the Fredman-Tarjan algorithm in JGraphT. A forward/reverse auction algorithm as developed by Bertsekas and Castañon should be able to handle unbalanced (or, as they call it, asymmetric) assignment much more efficiently, but has yet to be implemented.
The relative time/space efficiencies of BiGpairSEQ when backed by different MWM algorithms remains an open problem.
THE BiGpairSEQ ALGORITHM
- Sequence a sample plate of T cells as in pairSEQ.
- Pre-filter the sequence data to reduce error and minimize the size of the necessary graph.
- 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).
- 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.
- 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).
- 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.
- Find a maximum weight matching of the bipartite graph, using any MWM algorithm 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.
- 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
RUNNING THE PROGRAM
Download the current version of BiGpairSEQ_Sim.
BiGpairSEQ_Sim is an executable .jar file. Requires Java 14 or higher. OpenJDK 17 recommended.
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.
For example, to run the program with 32 gigabytes of memory, use the command:
java -Xmx32G -jar BiGpairSEQ_Sim.jar
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:
--------BiGPairSEQ SIMULATOR--------
ALPHA/BETA T CELL RECEPTOR MATCHING
USING WEIGHTED BIPARTITE GRAPHS
------------------------------------
Please select an option:
1) Generate a population of distinct cells
2) Generate a sample plate of T cells
3) Generate CDR3 alpha/beta occupancy data and overlap graph
4) Simulate bipartite graph CDR3 alpha/beta matching (BiGpairSEQ)
8) Options
9) About/Acknowledgments
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) Turn on calculation of p-values
7) 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:
- Cell Sample files in CSV format
- Sample Plate files in CSV format
- Graph/Data files in binary object serialization format
- 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.
To save file I/O time when using the interactive interface, the most recent instance of each of these four 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 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 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.)
Options when making a Cell Sample file:
- Number of T cells to generate
- Factor by which CDR3s are more diverse than CDR1s
Files are in CSV format. Rows are distinct T cells, columns are sequences within the cells.
Comments are preceded by #
Structure:
# Sample contains 1 unique CDR1 for every 4 unique CDR3s.
| Alpha CDR3 | Beta CDR3 | Alpha CDR1 | Beta CDR1 |
|---|---|---|---|
| unique number | unique number | number | number |
| ... | ... | ... | ... |
Sample Plate Files
Sample Plate files consist of any number of "wells" containing any number of T cells (as described above). The wells are filled randomly from a Cell Sample file, according to a selected 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. 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
- Statistical distribution to apply to Cell Sample file
- Poisson
- Gaussian
- Standard deviation size
- Exponential
- Lambda value
- Total number of wells on the plate
- Well populations random or fixed
- If random, minimum and maximum population sizes
- If fixed
- Number of sections on plate
- Number of T cells per well
- per section, if more than one section
- Sequence dropout rate
Files are in CSV format. There are no header labels. Every row represents a well.
Every value represents an individual cell, containing four sequences, depicted as an array string:
[CDR3A, CDR3B, CDR1A, CDR1B]. So a representative cell might look like this:
[525902, 791533, -1, 866282]
Notice that the CDR1 Alpha is missing in the cell above--sequence dropout from simulated amplification error.
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
# Plate size:
# Concentrations:
# Lambda -or- StdDev:
| Well 1, cell 1 | Well 1, cell 2 | Well 1, cell 3 | ... |
|---|---|---|---|
| Well 2, cell 1 | Well 2, cell 2 | Well 2, cell 3 | ... |
| Well 3, cell 1 | Well 3, cell 2 | Well 3, cell 3 | ... |
| ... | ... | ... | ... |
Graph/Data Files
Graph/Data files are serialized binaries of a Java object containing the weigthed bipartite graph representation of a Sample Plate, along with the necessary metadata for matching and results output. Making them requires a Cell Sample file (to construct a list of correct sequence pairs for checking the accuracy of BiGpairSEQ simulations) and a Sample Plate file (to construct the associated occupancy graph).
These files can be several gigabytes in size. Writing them to a file lets us generate a graph and its metadata once, 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 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.
Optional GraphML output
For portability of graph data to other software, turn on GraphML output
in the Options menu in interactive mode, or use the -graphmlcommand line argument. This will produce a .graphml file
for the weighted graph, with vertex attributes for sequence, type, total occupancy, total read count, and the read count for every individual occupied well.
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 serialized binary Graph/Data file (.ser). (Because .graphML files are larger than .ser files, BiGpairSEQ_Sim supports .graphML 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 #.
Options when running a BiGpairSEQ simulation of CDR3 alpha/beta matching:
- The minimum number of alpha/beta overlap wells to attempt to match
- (must be >= 1)
- The maximum number of alpha/beta overlap wells to attempt to match
- (must be <= the number of wells on the plate - 1)
- The maximum difference in alpha/beta occupancy to attempt to match
- (Optional. To skip using this filter, enter a value >= the number of wells on the plate)
- The minimum overlap percentage--the percentage of a sequence's occupied wells shared by another sequence--to attempt to match. Given as value in range 0 - 100.
- (Optional. To skip using this filter, enter 0)
Example output:
# 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 |
|---|---|---|---|---|---|---|
| 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 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.
RESULTS
Several BiGpairSEQ simulations were performed on a home computer with the following specs:
- Ryzen 5600X CPU
- 128GB of 3200MHz DDR4 RAM
- 2TB PCIe 3.0 SSD
- Linux Mint 21 (5.15 kernel)
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 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.
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 diskDONEHold graph data in memory until another graph is read-in? ABANDONED UNABANDONEDDONENo, 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 may be computer-specific.
Test whether pairing heap (currently used) or Fibonacci heap is more efficient for priority queue in current matching algorithmDONEin theory Fibonacci heap should be more efficient, but complexity overhead may eliminate theoretical advantageAdd controllable heap-type parameter?- Parameter implemented. Fibonacci heap the current default.
Implement sample plates with random numbers of T cells per well.DONE- Possible BiGpairSEQ advantage over pairSEQ: BiGpairSEQ is resilient to variations in well population sizes on a sample plate; pairSEQ is not due to nature of probability calculations.
- preliminary data suggests that BiGpairSEQ behaves roughly as though the whole plate had whatever the average well concentration is, but that's still speculative.
- Possible BiGpairSEQ advantage over pairSEQ: BiGpairSEQ is resilient to variations in well population sizes on a sample plate; pairSEQ is not due to nature of probability calculations.
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 wellApache 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 portabilityDONEHave 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 studiesDONEImplement custom Vertex class to simplify code and make it easier to implement different MWM algorithmsDONE- Advantage: would eliminate the need to use maps to associate vertices with sequences, which would make the code easier to understand.
- This also seems to be faster when using the same algorithm than the version with lots of maps, which is a nice bonus!
Implement simulation of read depth, and of read errors. Pre-filter graph for difference in read count to eliminate spurious sequences.DONE- Pre-filtering based on comparing (read depth) * (occupancy) to (read count) for each sequence works extremely well
Add read depth simulation options to CLIDONEUpdate graphml output to reflect current Vertex class attributesDONE- 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 sequenceDONE- Implement redistributive heap for LEDA matching algorithm to achieve theoretical worst case of O(n(m + n log C)) where C is highest edge weight.
- 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
Refactor simulator code to collect all needed data in a single scan of the plateDONECurrently it scans once for the vertices and then again for the edge weights. This made simulating read depth awkward, and incompatible with caching of plate files.This would be a fairly major rewrite of the simulator code, but could make things faster, and would definitely make them cleaner.
- Implement Duan and Su's maximum weight matching algorithm
- Add controllable algorithm-type parameter?
- This would be fun and valuable, but probably take more time than I have for a hobby project.
Implement an auction algorithm for maximum weight matchingDONE- Implement a forward/reverse auction algorithm for maximum weight matching
- Implement an algorithm for approximating a maximum weight matching
- Some of these run in linear or near-linear time
- given that the underlying biological samples have many, many sources of error, this would probably be the most useful option in practice. It seems less mathematically elegant, though, and so less fun for me.
- Implement Vose's alias method for arbitrary statistical distributions of cells
- Should probably refactor to use apache commons rng for this
- Use commons JCS for caching
- Parameterize pre-filtering options
CITATIONS
- Howie, B., Sherwood, A. M., et al. "High-throughput pairing of T cell receptor alpha and beta sequences." Sci. Transl. Med. 7, 301ra131 (2015)
- Duan, R., Su H. "A Scaling Algorithm for Maximum Weight Matching in Bipartite Graphs." Proceedings of the Twenty-Third Annual ACM-SIAM Symposium on Discrete Algorithms, p. 1413-1424. (2012)
- Melhorn, K., Näher, St. The LEDA Platform of Combinatorial and Geometric Computing. Cambridge University Press. Chapter 7, Graph Algorithms; p. 132-162 (1999)
- Fredman, M., Tarjan, R. "Fibonacci heaps and their uses in improved network optimization algorithms." J. ACM, 34(3):596–615 (1987))
EXTERNAL LIBRARIES USED
- JGraphT -- Graph theory data structures and algorithms
- JHeaps -- For pairing heap priority queue used in maximum weight matching algorithm
- Apache Commons CSV -- For CSV file output
- Apache Commons CLI -- To enable command line arguments for scripting.
ACKNOWLEDGEMENTS
BiGpairSEQ was conceived in collaboration with the author's spouse, Dr. Alice MacQueen, who brought the original 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. Improvements and documentation, 2022–2023.
DISCLOSURE
The earliest versions of the BiGpairSEQ simulator were written in 2021 to let Dr. MacQueen test hypothetical extensions of the published pairSEQ protocol while she was interviewing for a position at Adaptive Biotechnologies. She has been employed at Adaptive Biotechnologies since 2022.
The author has worked on this BiGpairSEQ simulator since 2021 without Dr. MacQueen's involvement, since she has had access to related, proprietary technologies. The author has had no such access, relying exclusively on the 2015 pairSEQ paper and other academic publications. He continues to work on BiGpairSEQ recreationally, as it involves some very beautiful math.