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

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.

BiGpairSEQ creates a simple bipartite weighted 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 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.

This is a well-studied combinatorial optimization problem, with many known solutions. The best currently-known algorithm for bipartite graphs with integer weights--which is what BiGpairSEQ uses-- 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. This is the best known efficiency for finding a maximum weight matching on a bipartite graph, and the integer edge weight requirement makes it ideal for BiGpairSEQ.

Unfortunately, it's a fairly new algorithm, and the integer edge weight requirement makes it less generically useful. It is not 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).

The current version of the program uses a pairing heap instead of a Fibonacci heap for its priority queue, which has lower theoretical efficiency but also lower complexity overhead, and is often equivalently performant in practice.

USAGE

RUNNING THE PROGRAM

BiGpairSEQ_Sim is an executable .jar file. Requires Java 11 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

Once running, BiGpairSEQ_Sim has 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)
9) About/Acknowledgments
0) Exit

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 and Data files in binary object serialization format
Matching Results files in CSV format

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.

Cell Sample Files

Cell Sample files consist of any number of distinct "T cells." Every cell contains four sequences: Alpha CDR3, Beta CDR, Alpha CDR1, Beta CDR1. The sequences are represented by random integers. CDR3 Alpha and Beta sequences are all unique. CDR1 Alpha and Beta sequences are not necessarily unique; the relative diversity can be set when making a Cell Sample 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 example:

# 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 (approximately) evenly-sized sections, each of which can have a different number of T cells per well.

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 exponential with a lambda of approximately 0.6. (Howie et al. 2015)
Total number of wells on the plate
Number of sections on plate
Number of T cells per well
- per section, if more than one section
Dropout rate

Files are in CSV format. There are no header labels. Every row represents a well. Every column represents an individual cell, containing four sequences, represented by an array string: [CDR3A, CDR3B, CDR1A, CDR1B]. So a representative cell might look like this:

[525902, 791533, -1, 866282]

Notice that the Alpha CDR1 is missing in the cell above, due to sequence dropout. Dropouts are represented by replacing sequences with the value -1. Comments are preceded by #

Structure Example:

# Cell source file name: 4MilCells.csv
# Plate size: 96
# Error rate: 0.1
# Concentrations: 10000 5000 500
# Lambda: 0.6

well 1	well 2	well 3	...
[105383, 786528, 959247, 925928]	[525902, 791533, -1, 866282]	[409236, 132303, 804465, 942261]	...
[249930, 301502, 970003, 881099]	[523787, 552952, 997194, 970507]	[425363, 417411, 845399, -1]	...
...	...	...	...

Graph and Data Files

Graph and Data files are serialized binaries of a Java object containing the graph representation of a Sample Plate and 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 and Data file:

The Cell Sample file to use
The Sample Plate file (generated from the given Cell Sample file) to use.

These files do not have a human-readable structure, and are not portable to other programs. (Export of graphs in a portable data format may be implemented in the future. The tricky part is encoding the necessary metadata.)

Matching Results Files

Matching results files consist of the results of a BiGpairSEQ matching simulation. 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
- (To skip using this filter, enter a value >= the number of wells on the plate)
The minimum percentage of a sequence's occupied wells shared by another sequence to attempt to match
- given value from 0 to 100
- (To skip using this filter, enter 0)

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

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

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, using the (2021 corrected) formula from the original pairSEQ paper. (Howie, et al. 2015)

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?~~
- No, this won't work, because BiGpairSEQ simulations alter the underlying graph based on filtering constraints. Changes would cascade with multiple experiments.
~~See if there's a reasonable way to reformat Sample Plate files so that wells are columns instead of rows~~ DONE
Enable GraphML output in addition to serialized object binaries, for data portability
- Custom vertex type with attribute for sequence occupancy?
Re-implement CDR1 matching method
Re-implement command line arguments, to enable scripting and statistical simulation studies
Implement Duan and Su's maximum weight matching algorithms
- Add controllable algorithm-type parameter?
Test whether pairing heap (currently used) or Fibonacci heap is more efficient for current matching algorithm
- in theory Fibonacci heap should be more efficient, but complexity overhead may eliminate theoretical advantage
- Add controllable heap-type parameter?
Implement sample plates with random numbers of T cells per well
- Possible BiGpairSEQ advantage over pairSEQ: BiGpairSEQ is resilient to variations in well populations; pairSEQ is not.
  - preliminary data suggests that BiGpairSEQ behaves roughly as though the whole plate had whatever the average well concentration is, but that's still speculative.

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)
K. Melhorn, St. Näher. The LEDA Platform of Combinatorial and Geometric Computing. Cambridge University Press. Chapter 7, Graph Algorithms; p. 132-162 (1999)
M. Fredman, R. Tarjan. "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. (Awaiting re-implementation.)

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

Eugene Fischer, 2021. UI improvements and documentation, 2022.

readme.md Unescape Escape

BiGpairSEQ SIMULATOR

ABOUT