262 lines
13 KiB
Markdown
262 lines
13 KiB
Markdown
# BiGpairSEQ SIMULATOR
|
||
|
||
|
||
## 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
|
||
|
||
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](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 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](https://jdk.java.net/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:
|
||
|
||
---
|
||
# 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:
|
||
|
||
---
|
||
```
|
||
# Cell source file name:
|
||
# Each row represents one well on the plate
|
||
# Plate size:
|
||
# Concentrations:
|
||
# Lambda:
|
||
```
|
||
| 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 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
|
||
* ~~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.
|
||
* 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.
|
||
* 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
|
||
* Can possibly sort the wells by length first, then construct entries
|
||
|
||
## CITATIONS
|
||
* Howie, B., Sherwood, A. M., et al. ["High-throughput pairing of T cell receptor alpha and beta sequences."](https://pubmed.ncbi.nlm.nih.gov/26290413/) Sci. Transl. Med. 7, 301ra131 (2015)
|
||
* Duan, R., Su H. ["A Scaling Algorithm for Maximum Weight Matching in Bipartite Graphs."](https://web.eecs.umich.edu/~pettie/matching/Duan-Su-scaling-bipartite-matching.pdf) 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.](https://people.mpi-inf.mpg.de/~mehlhorn/LEDAbook.html) 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."](https://www.cl.cam.ac.uk/teaching/1011/AlgorithII/1987-FredmanTar-fibonacci.pdf) J. ACM, 34(3):596–615 (1987))
|
||
|
||
## EXTERNAL LIBRARIES USED
|
||
* [JGraphT](https://jgrapht.org) -- Graph theory data structures and algorithms
|
||
* [JHeaps](https://www.jheaps.org) -- For pairing heap priority queue used in maximum weight matching algorithm
|
||
* [Apache Commons CSV](https://commons.apache.org/proper/commons-csv/) -- For CSV file output
|
||
* [Apache Commons CLI](https://commons.apache.org/proper/commons-cli/) -- To enable command line arguments for scripting. (**Awaiting re-implementation**.)
|
||
|
||
## 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. |