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 cell (of which there may be many clones).

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 is 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 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 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 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. 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

Sequence a sample plate of T cells as in pairSEQ.
Pre-filter the sequence data to minimize the size of the necessary graph.
- Saturating sequence filter: remove any sequences present in all of the wells on the sample plate, as there is no signal in the occupancy data of saturating sequences.
- 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. If a sequence is read correctly at least half the time, then the total read count of a sequence (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 C < (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 causes a sequence to seem to be in a well where it is not, in fact, present. This can be detected by looking at 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 3 or more wells--that is, edges with weight >= 3.0.)
The resultant matching represents the likeliest TCRA/TCRB sequence pairs based on the occupancy pattern of the sample plate.

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

There are a number of command line options, to allow the program to be used in shell scripts. For a full list, use the -help flag:

java -jar BiGpairSEQ_Sim.jar -help

If no command line arguments are given, BiGpairSEQ_Sim will launch with an interactive, menu-driven CLI for generating files and simulating TCR pairing. The main menu looks like this:

--------BiGPairSEQ SIMULATOR--------
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
    - (Based on the slope of the graph in Figure 4C of the pairSEQ paper, the distribution of the original experiment was very roughly exponential with a lambda ~0.6. (Howie, et al. 2015) The actual distribution was certainly quite different.)
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:

# sample plate filename: 8MilCells_Plate.csv
# 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
# pre-filter sequences present in all wells: true
# pre-filter sequences based on occupancy/read count discrepancy: true
# 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 abpve 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 in 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

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)

Simulation 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 comprising ~11,900 TCRAs and TCRBs, taken from a sample of 8,400,000 distinct cells with an exponential frequency distribution (lambda 0.6). The sequence dropout rate was 10%, as the analysis from the original paper concluded that most TCR sequences "have less than a 10% chance of going unobserved." (Howie, et al. 2015)

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 any sequencing errors, to reduce the processing time.

With min/max occupancy thresholds of 3 and 95 wells respectively for matching, BiGpairSEQ identified:

8,495 correct pairings
5 incorrect pairings

for an overall pairing accuracy of 99.9992%.

The total simulation time (excluding file I/O) was 28m52. The total elapsed time with file I/O was 41m23s. Calculation of p-values was enabled for this simulation, increasing the overall processing time.

BEHAVIOR WITH RANDOMIZED WELL POPULATIONS (old results, need updating for new version of the simulator (though resilience to varying well populations is unchanged))

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 may be computer-specific.
~~Test whether pairing heap (currently used) or Fibonacci heap is more efficient for priority queue in current matching algorithm~~ DONE
- ~~in theory Fibonacci heap should be more efficient, but complexity overhead may eliminate theoretical advantage~~
- ~~Add controllable heap-type parameter?~~
  - Parameter implemented. 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.
~~See if there's a reasonable way to reformat Sample Plate files so that wells are columns instead of rows.~~
- ~~Problem is variable number of cells in a well~~
- ~~Apache Commons CSV library writes entries a row at a time~~
  - Got this working, but at the cost of a profoundly strange bug in graph occupancy filtering. Have reverted the repo until I can figure out what caused that. Given how easily Thingiverse transposes CSV matrices in R, might not even be worth fixing.
~~Enable GraphML output in addition to serialized object binaries, for data portability~~ DONE
- ~~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
~~Implement custom Vertex class to simplify code and make it easier to implement different MWM algorithms~~ DONE
- 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 CLI~~ DONE
~~Update graphml output to reflect current Vertex class attributes~~ DONE
- 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
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 plate
- Currently 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 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. 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.

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 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.

readme.md Unescape Escape

BiGpairSEQ SIMULATOR

CONTENTS

ABOUT