Compare commits
74 Commits
| Author | SHA1 | Date | |
|---|---|---|---|
|
|
e4e5a1f979 | ||
|
|
73c83bf35d | ||
|
|
06e72314b0 | ||
|
|
63317f2aa0 | ||
|
|
a054c0c20a | ||
|
|
29b844afd2 | ||
|
|
3ba305abdb | ||
|
|
3707923398 | ||
|
|
cf771ce574 | ||
| f980722b56 | |||
| 1df86f01df | |||
| 96ba57d653 | |||
| b602fb02f1 | |||
| 325e1ebe2b | |||
| df047267ee | |||
| 03e8d31210 | |||
| 582dc3ef40 | |||
| 4c872ed48e | |||
| 3fc39302c7 | |||
| 578bdc0fbf | |||
| 8275cf7740 | |||
| 64209691f0 | |||
| 1886800873 | |||
| bedf0894bc | |||
| 2ac3451842 | |||
| 67ec3f3764 | |||
| b5a8b7e2d5 | |||
| 9fb3095f0f | |||
| 25acf920c2 | |||
| f301327693 | |||
| e04d2d6777 | |||
| 3e41afaa64 | |||
| bc5d67680d | |||
| f2347e8fc2 | |||
| c8364d8a6e | |||
| 6f5afbc6ec | |||
| fb4d22e7a4 | |||
| e10350c214 | |||
| b1155f8100 | |||
| 12b003a69f | |||
| 32c5bcaaff | |||
| 2485ac4cf6 | |||
| 05556bce0c | |||
| a822f69ea4 | |||
| 3d1f8668ee | |||
| 40c743308b | |||
| 5246cc4a0c | |||
| a5f7c0641d | |||
| 8ebfc1469f | |||
| b53f5f1cc0 | |||
| 974d2d650c | |||
| 6b5837e6ce | |||
| b4cc240048 | |||
| ff72c9b359 | |||
| 88eb8aca50 | |||
| 98bf452891 | |||
| c2db4f87c1 | |||
| 8935407ade | |||
| 9fcc20343d | |||
| e4d094d796 | |||
| f385ebc31f | |||
| 8745550e11 | |||
| 41805135b3 | |||
| 373a5e02f9 | |||
| 7f18311054 | |||
| bcb816c3e6 | |||
| dad0fd35fd | |||
| 35d580cfcf | |||
| ab8d98ed81 | |||
| 3d9890e16a | |||
| dd64ac2731 | |||
| a5238624f1 | |||
| d8ba42b801 | |||
| 8edd89d784 |
188
readme.md
188
readme.md
@@ -12,7 +12,7 @@ Unlike pairSEQ, which calculates p-values for every TCR alpha/beta overlap and c
|
|||||||
against a null distribution, BiGpairSEQ does not do any statistical calculations
|
against a null distribution, BiGpairSEQ does not do any statistical calculations
|
||||||
directly.
|
directly.
|
||||||
|
|
||||||
BiGpairSEQ creates a [weightd bipartite graph](https://en.wikipedia.org/wiki/Bipartite_graph) representing the sample plate.
|
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
|
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.
|
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.)
|
(Sequences present in *all* wells are filtered out prior to creating the graph, as there is no signal in their occupancy pattern.)
|
||||||
@@ -20,8 +20,8 @@ The problem of pairing TCRA/TCRB sequences thus reduces to the "assignment probl
|
|||||||
matching on a bipartite graph--the subset of vertex-disjoint edges whose weights sum to the maximum possible value.
|
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.
|
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
|
The most efficient algorithm known to the author for maximum weight matching of a bipartite graph with strictly integral
|
||||||
is from Duan and Su (2012). For a graph with m edges, n vertices per side, and maximum integer edge weight N,
|
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
|
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.
|
bipartite with integer weights, this algorithm is ideal for BiGpairSEQ.
|
||||||
|
|
||||||
@@ -29,17 +29,13 @@ Unfortunately, it's a fairly new algorithm, and not yet implemented by the graph
|
|||||||
So this program instead uses the Fibonacci heap-based algorithm of Fredman and Tarjan (1987), which has a worst-case
|
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).
|
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
|
## USAGE
|
||||||
|
|
||||||
### RUNNING THE PROGRAM
|
### RUNNING THE PROGRAM
|
||||||
|
|
||||||
[Download the current version of BiGpairSEQ_Sim.](https://gitea.ejsf.synology.me/efischer/BiGpairSEQ/releases)
|
[Download the current version of BiGpairSEQ_Sim.](https://gitea.ejsf.synology.me/efischer/BiGpairSEQ/releases)
|
||||||
|
|
||||||
BiGpairSEQ_Sim is an executable .jar file. Requires Java 11 or higher. [OpenJDK 17](https://jdk.java.net/17/)
|
BiGpairSEQ_Sim is an executable .jar file. Requires Java 14 or higher. [OpenJDK 17](https://jdk.java.net/17/)
|
||||||
recommended.
|
recommended.
|
||||||
|
|
||||||
Run with the command:
|
Run with the command:
|
||||||
@@ -47,13 +43,18 @@ Run with the command:
|
|||||||
`java -jar BiGpairSEQ_Sim.jar`
|
`java -jar BiGpairSEQ_Sim.jar`
|
||||||
|
|
||||||
Processing sample plates with tens of thousands of sequences may require large amounts
|
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.
|
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:
|
For example, to run the program with 32 gigabytes of memory, use the command:
|
||||||
|
|
||||||
`java -Xmx32G -jar BiGpairSEQ_Sim.jar`
|
`java -Xmx32G -jar BiGpairSEQ_Sim.jar`
|
||||||
|
|
||||||
Once running, BiGpairSEQ_Sim has an interactive, menu-driven CLI for generating files and simulating TCR pairing. The
|
There are a number of command line options, to allow the program to be used in shell scripts. For a full list,
|
||||||
main menu looks like this:
|
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--------
|
--------BiGPairSEQ SIMULATOR--------
|
||||||
@@ -70,6 +71,19 @@ Please select an option:
|
|||||||
0) Exit
|
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) Maximum weight matching algorithm options
|
||||||
|
0) Return to main menu
|
||||||
|
```
|
||||||
|
|
||||||
|
|
||||||
### INPUT/OUTPUT
|
### INPUT/OUTPUT
|
||||||
|
|
||||||
To run the simulation, the program reads and writes 4 kinds of files:
|
To run the simulation, the program reads and writes 4 kinds of files:
|
||||||
@@ -79,20 +93,25 @@ To run the simulation, the program reads and writes 4 kinds of files:
|
|||||||
* Matching Results files in CSV 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
|
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.
|
(.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, the most recent instance of each of these four
|
To save file I/O time, the most recent instance of each of these four
|
||||||
files either generated or read from disk can be cached in program memory. This is especially important for Graph/Data files,
|
files either generated or read from disk can be cached in program memory. When caching is active, subsequent uses of the
|
||||||
which can be several gigabytes in size. Since some simulations may require running multiple,
|
same data file won't need to be read in again until another file of that type is used or generated,
|
||||||
differently-configured BiGpairSEQ matchings on the same graph, keeping the most recent graph cached can reduce execution time
|
|
||||||
|
|
||||||
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
|
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.
|
filenames as entered by the user. On encountering a new filename, the program flushes its cache and reads in the new file.
|
||||||
|
|
||||||
The program's caching behavior can be controlled in the Options menu. By default, caching for cell sample and
|
(Note that cached Graph/Data files must be transformed back into their original state after a matching experiment, which
|
||||||
sample plate files is OFF, and caching for graph/data files is ON.
|
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
|
||||||
Cell Sample files consist of any number of distinct "T cells." Every cell contains
|
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
|
four sequences: Alpha CDR3, Beta CDR3, Alpha CDR1, Beta CDR1. The sequences are represented by
|
||||||
@@ -110,7 +129,6 @@ Comments are preceded by `#`
|
|||||||
|
|
||||||
Structure:
|
Structure:
|
||||||
|
|
||||||
---
|
|
||||||
# Sample contains 1 unique CDR1 for every 4 unique CDR3s.
|
# Sample contains 1 unique CDR1 for every 4 unique CDR3s.
|
||||||
| Alpha CDR3 | Beta CDR3 | Alpha CDR1 | Beta CDR1 |
|
| Alpha CDR3 | Beta CDR3 | Alpha CDR1 | Beta CDR1 |
|
||||||
|---|---|---|---|
|
|---|---|---|---|
|
||||||
@@ -134,11 +152,14 @@ Options when making a Sample Plate file:
|
|||||||
* Standard deviation size
|
* Standard deviation size
|
||||||
* Exponential
|
* Exponential
|
||||||
* Lambda value
|
* 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))*
|
* *(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
|
* Total number of wells on the plate
|
||||||
* Number of sections on plate
|
* Well populations random or fixed
|
||||||
* Number of T cells per well
|
* If random, minimum and maximum population sizes
|
||||||
* per section, if more than one section
|
* If fixed
|
||||||
|
* Number of sections on plate
|
||||||
|
* Number of T cells per well
|
||||||
|
* per section, if more than one section
|
||||||
* Dropout rate
|
* Dropout rate
|
||||||
|
|
||||||
Files are in CSV format. There are no header labels. Every row represents a well.
|
Files are in CSV format. There are no header labels. Every row represents a well.
|
||||||
@@ -152,7 +173,6 @@ Dropout sequences are replaced with the value `-1`. Comments are preceded by `#`
|
|||||||
|
|
||||||
Structure:
|
Structure:
|
||||||
|
|
||||||
---
|
|
||||||
```
|
```
|
||||||
# Cell source file name:
|
# Cell source file name:
|
||||||
# Each row represents one well on the plate
|
# Each row represents one well on the plate
|
||||||
@@ -181,14 +201,24 @@ Options for creating a Graph/Data file:
|
|||||||
* The Cell Sample file to use
|
* The Cell Sample file to use
|
||||||
* The Sample Plate file to use. (This must have been generated from the selected Cell Sample file.)
|
* The Sample Plate file to use. (This must have been generated from the selected Cell Sample file.)
|
||||||
|
|
||||||
These files do not have a human-readable structure, and are not portable to other programs. (Export of graphs in a
|
These files do not have a human-readable structure, and are not portable to other programs.
|
||||||
portable data format may be implemented in the future. The tricky part is encoding the necessary metadata.)
|
|
||||||
|
*Optional GraphML output*
|
||||||
|
|
||||||
|
For portability of graph data to other software, turn on [GraphML](http://graphml.graphdrawing.org/index.html) output
|
||||||
|
in the Options menu in interactive mode, or use the `-graphml`command line argument. This will produce a .graphml file
|
||||||
|
for the weighted graph, with vertex attributes for sequence, type, and occupancy data. 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
|
||||||
Matching results files consist of the results of a BiGpairSEQ matching simulation. Making them requires a Graph and
|
Matching results files consist of the results of a BiGpairSEQ matching simulation. Making them requires a serialized
|
||||||
Data file. Matching results files are in CSV format. Rows are sequence pairings with extra relevant data. Columns are pairing-specific details.
|
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.)
|
||||||
|
|
||||||
|
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 `#`.
|
Metadata about the matching simulation is included as comments. Comments are preceded by `#`.
|
||||||
|
|
||||||
Options when running a BiGpairSEQ simulation of CDR3 alpha/beta matching:
|
Options when running a BiGpairSEQ simulation of CDR3 alpha/beta matching:
|
||||||
@@ -203,7 +233,6 @@ Options when running a BiGpairSEQ simulation of CDR3 alpha/beta matching:
|
|||||||
|
|
||||||
Example output:
|
Example output:
|
||||||
|
|
||||||
---
|
|
||||||
```
|
```
|
||||||
# Source Sample Plate file: 4MilCellsPlate.csv
|
# Source Sample Plate file: 4MilCellsPlate.csv
|
||||||
# Source Graph and Data file: 4MilCellsPlateGraph.ser
|
# Source Graph and Data file: 4MilCellsPlateGraph.ser
|
||||||
@@ -235,46 +264,105 @@ Example output:
|
|||||||
P-values are calculated *after* BiGpairSEQ matching is completed, for purposes of comparison only,
|
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)
|
using the (2021 corrected) formula from the original pairSEQ paper. (Howie, et al. 2015)
|
||||||
|
|
||||||
### PERFORMANCE
|
|
||||||
Performance details of the example excerpted above:
|
## 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),
|
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,
|
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.
|
taken from a sample of 4,000,000 distinct cells with an exponential frequency distribution (lambda 0.6).
|
||||||
|
|
||||||
With min/max occupancy threshold of 3 and 94 wells for matching, and no other pre-filtering, BiGpairSEQ identified 5,151
|
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%.
|
correct pairings and 18 incorrect pairings, for an accuracy of 99.652%.
|
||||||
|
|
||||||
The simulation time was 14'22". If intermediate results were held in memory, this would be equivalent to the total elapsed time.
|
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
|
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
|
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.
|
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 weighted matching.
|
||||||
|
|
||||||
|
## BEHAVIOR WITH RANDOMIZED WELL POPULATIONS
|
||||||
|
|
||||||
|
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
|
## TODO
|
||||||
|
|
||||||
* ~~Try invoking GC at end of workloads to reduce paging to disk~~ DONE
|
* ~~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
|
* ~~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.*~~
|
* ~~*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. If so, awesome.
|
* 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.
|
||||||
|
* 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.
|
* 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~~
|
* ~~Problem is variable number of cells in a well~~
|
||||||
* ~~Apache Commons CSV library writes entries a row at a time~~
|
* ~~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._
|
* _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.
|
||||||
* Re-implement command line arguments, to enable scripting and statistical simulation studies
|
* ~~Enable GraphML output in addition to serialized object binaries, for data portability~~ DONE
|
||||||
* Implement sample plates with random numbers of T cells per well.
|
* ~~Custom vertex type with attribute for sequence occupancy?~~ DONE
|
||||||
* Possible BiGpairSEQ advantage over pairSEQ: BiGpairSEQ is resilient to variations in well population sizes on a sample plate; pairSEQ is not.
|
* Advantage: would eliminate the need to use maps to associate vertices with sequences, which would make the code easier to understand.
|
||||||
* preliminary data suggests that BiGpairSEQ behaves roughly as though the whole plate had whatever the *average* well concentration is, but that's still speculative.
|
* ~~Have a branch where this is implemented, but there's a bug that broke matching. Don't currently have time to fix.~~
|
||||||
* Enable GraphML output in addition to serialized object binaries, for data portability
|
* ~~Re-implement command line arguments, to enable scripting and statistical simulation studies~~ DONE
|
||||||
* Custom vertex type with attribute for sequence occupancy?
|
|
||||||
* Re-implement CDR1 matching method
|
* Re-implement CDR1 matching method
|
||||||
* Implement Duan and Su's maximum weight matching algorithm
|
* Implement Duan and Su's maximum weight matching algorithm
|
||||||
* Add controllable algorithm-type parameter?
|
* Add controllable algorithm-type parameter?
|
||||||
* Test whether pairing heap (currently used) or Fibonacci heap is more efficient for priority queue in current matching algorithm
|
* This would be fun and valuable, but probably take more time than I have for a hobby project.
|
||||||
* in theory Fibonacci heap should be more efficient, but complexity overhead may eliminate theoretical advantage
|
* Implement an algorithm for approximating a maximum weight matching
|
||||||
* Add controllable heap-type parameter?
|
* 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
|
||||||
|
* Enable post-filtering instead of pre-filtering. Pre-filtering of things like singleton sequences or saturating-occupancy sequences reduces graph size, but could conceivably reduce pairing accuracy by throwing away data. While these sequences have very little signal, it would be interesting to compare unfiltered results to filtered results. This would require a much, much faster MWM algorithm, though, to handle the much larger graphs. Possible one of the linear-time approximation algorithms.
|
||||||
|
|
||||||
|
|
||||||
## CITATIONS
|
## 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)
|
* 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)
|
||||||
@@ -286,7 +374,7 @@ slightly less time than the simulation itself. Real elapsed time from start to f
|
|||||||
* [JGraphT](https://jgrapht.org) -- Graph theory data structures and algorithms
|
* [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
|
* [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 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**.)
|
* [Apache Commons CLI](https://commons.apache.org/proper/commons-cli/) -- To enable command line arguments for scripting.
|
||||||
|
|
||||||
## ACKNOWLEDGEMENTS
|
## ACKNOWLEDGEMENTS
|
||||||
BiGpairSEQ was conceived in collaboration with Dr. Alice MacQueen, who brought the original
|
BiGpairSEQ was conceived in collaboration with Dr. Alice MacQueen, who brought the original
|
||||||
|
|||||||
@@ -1,6 +1,6 @@
|
|||||||
import java.util.Random;
|
import java.util.Random;
|
||||||
|
|
||||||
//main class. For choosing interface type and caching file data
|
//main class. For choosing interface type and holding settings
|
||||||
public class BiGpairSEQ {
|
public class BiGpairSEQ {
|
||||||
|
|
||||||
private static final Random rand = new Random();
|
private static final Random rand = new Random();
|
||||||
@@ -12,7 +12,11 @@ public class BiGpairSEQ {
|
|||||||
private static String graphFilename = null;
|
private static String graphFilename = null;
|
||||||
private static boolean cacheCells = false;
|
private static boolean cacheCells = false;
|
||||||
private static boolean cachePlate = false;
|
private static boolean cachePlate = false;
|
||||||
private static boolean cacheGraph = true;
|
private static boolean cacheGraph = false;
|
||||||
|
private static String priorityQueueHeapType = "FIBONACCI";
|
||||||
|
private static boolean outputBinary = true;
|
||||||
|
private static boolean outputGraphML = false;
|
||||||
|
private static final String version = "version 3.0";
|
||||||
|
|
||||||
public static void main(String[] args) {
|
public static void main(String[] args) {
|
||||||
if (args.length == 0) {
|
if (args.length == 0) {
|
||||||
@@ -20,8 +24,8 @@ public class BiGpairSEQ {
|
|||||||
}
|
}
|
||||||
else {
|
else {
|
||||||
//This will be uncommented when command line arguments are re-implemented.
|
//This will be uncommented when command line arguments are re-implemented.
|
||||||
//CommandLineInterface.startCLI(args);
|
CommandLineInterface.startCLI(args);
|
||||||
System.out.println("Command line arguments are still being re-implemented.");
|
//System.out.println("Command line arguments are still being re-implemented.");
|
||||||
}
|
}
|
||||||
}
|
}
|
||||||
|
|
||||||
@@ -151,4 +155,23 @@ public class BiGpairSEQ {
|
|||||||
}
|
}
|
||||||
BiGpairSEQ.cacheGraph = cacheGraph;
|
BiGpairSEQ.cacheGraph = cacheGraph;
|
||||||
}
|
}
|
||||||
|
|
||||||
|
public static String getPriorityQueueHeapType() {
|
||||||
|
return priorityQueueHeapType;
|
||||||
|
}
|
||||||
|
|
||||||
|
public static void setPairingHeap() {
|
||||||
|
priorityQueueHeapType = "PAIRING";
|
||||||
|
}
|
||||||
|
|
||||||
|
public static void setFibonacciHeap() {
|
||||||
|
priorityQueueHeapType = "FIBONACCI";
|
||||||
|
}
|
||||||
|
|
||||||
|
public static boolean outputBinary() {return outputBinary;}
|
||||||
|
public static void setOutputBinary(boolean b) {outputBinary = b;}
|
||||||
|
|
||||||
|
public static boolean outputGraphML() {return outputGraphML;}
|
||||||
|
public static void setOutputGraphML(boolean b) {outputGraphML = b;}
|
||||||
|
public static String getVersion() { return version; }
|
||||||
}
|
}
|
||||||
|
|||||||
@@ -1,10 +1,37 @@
|
|||||||
|
import java.util.ArrayList;
|
||||||
|
import java.util.Collections;
|
||||||
import java.util.List;
|
import java.util.List;
|
||||||
|
import java.util.stream.IntStream;
|
||||||
|
|
||||||
public class CellSample {
|
public class CellSample {
|
||||||
|
|
||||||
private List<Integer[]> cells;
|
private List<Integer[]> cells;
|
||||||
private Integer cdr1Freq;
|
private Integer cdr1Freq;
|
||||||
|
|
||||||
|
public CellSample(Integer numDistinctCells, Integer cdr1Freq){
|
||||||
|
this.cdr1Freq = cdr1Freq;
|
||||||
|
List<Integer> numbersCDR3 = new ArrayList<>();
|
||||||
|
List<Integer> numbersCDR1 = new ArrayList<>();
|
||||||
|
Integer numDistCDR3s = 2 * numDistinctCells + 1;
|
||||||
|
IntStream.range(1, numDistCDR3s + 1).forEach(i -> numbersCDR3.add(i));
|
||||||
|
IntStream.range(numDistCDR3s + 1, numDistCDR3s + 1 + (numDistCDR3s / cdr1Freq) + 1).forEach(i -> numbersCDR1.add(i));
|
||||||
|
Collections.shuffle(numbersCDR3);
|
||||||
|
Collections.shuffle(numbersCDR1);
|
||||||
|
|
||||||
|
//Each cell represented by 4 values
|
||||||
|
//two CDR3s, and two CDR1s. First two values are CDR3s (alpha, beta), second two are CDR1s (alpha, beta)
|
||||||
|
List<Integer[]> distinctCells = new ArrayList<>();
|
||||||
|
for(int i = 0; i < numbersCDR3.size() - 1; i = i + 2){
|
||||||
|
Integer tmpCDR3a = numbersCDR3.get(i);
|
||||||
|
Integer tmpCDR3b = numbersCDR3.get(i+1);
|
||||||
|
Integer tmpCDR1a = numbersCDR1.get(i % numbersCDR1.size());
|
||||||
|
Integer tmpCDR1b = numbersCDR1.get((i+1) % numbersCDR1.size());
|
||||||
|
Integer[] tmp = {tmpCDR3a, tmpCDR3b, tmpCDR1a, tmpCDR1b};
|
||||||
|
distinctCells.add(tmp);
|
||||||
|
}
|
||||||
|
this.cells = distinctCells;
|
||||||
|
}
|
||||||
|
|
||||||
public CellSample(List<Integer[]> cells, Integer cdr1Freq){
|
public CellSample(List<Integer[]> cells, Integer cdr1Freq){
|
||||||
this.cells = cells;
|
this.cells = cells;
|
||||||
this.cdr1Freq = cdr1Freq;
|
this.cdr1Freq = cdr1Freq;
|
||||||
|
|||||||
@@ -1,5 +1,9 @@
|
|||||||
import org.apache.commons.cli.*;
|
import org.apache.commons.cli.*;
|
||||||
|
|
||||||
|
import java.io.IOException;
|
||||||
|
import java.util.Arrays;
|
||||||
|
import java.util.stream.Stream;
|
||||||
|
|
||||||
/*
|
/*
|
||||||
* Class for parsing options passed to program from command line
|
* Class for parsing options passed to program from command line
|
||||||
*
|
*
|
||||||
@@ -29,6 +33,8 @@ import org.apache.commons.cli.*;
|
|||||||
* cellfile : name of the cell sample file to use as input
|
* cellfile : name of the cell sample file to use as input
|
||||||
* platefile : name of the sample plate file to use as input
|
* platefile : name of the sample plate file to use as input
|
||||||
* output : name of the output file
|
* output : name of the output file
|
||||||
|
* graphml : output a graphml file
|
||||||
|
* binary : output a serialized binary object file
|
||||||
*
|
*
|
||||||
* Match flags:
|
* Match flags:
|
||||||
* graphFile : name of graph and data file to use as input
|
* graphFile : name of graph and data file to use as input
|
||||||
@@ -43,242 +49,172 @@ import org.apache.commons.cli.*;
|
|||||||
public class CommandLineInterface {
|
public class CommandLineInterface {
|
||||||
|
|
||||||
public static void startCLI(String[] args) {
|
public static void startCLI(String[] args) {
|
||||||
//These command line options are a big mess
|
//Options sets for the different modes
|
||||||
//Really, I don't think command line tools are expected to work in this many different modes
|
Options mainOptions = buildMainOptions();
|
||||||
//making cells, making plates, and matching are the sort of thing that UNIX philosophy would say
|
Options cellOptions = buildCellOptions();
|
||||||
//should be three separate programs.
|
Options plateOptions = buildPlateOptions();
|
||||||
//There might be a way to do it with option parameters?
|
Options graphOptions = buildGraphOptions();
|
||||||
|
Options matchOptions = buildMatchCDR3options();
|
||||||
//main options set
|
|
||||||
Options mainOptions = new Options();
|
|
||||||
Option makeCells = Option.builder("cells")
|
|
||||||
.longOpt("make-cells")
|
|
||||||
.desc("Makes a file of distinct cells")
|
|
||||||
.build();
|
|
||||||
Option makePlate = Option.builder("plates")
|
|
||||||
.longOpt("make-plates")
|
|
||||||
.desc("Makes a sample plate file")
|
|
||||||
.build();
|
|
||||||
Option makeGraph = Option.builder("graph")
|
|
||||||
.longOpt("make-graph")
|
|
||||||
.desc("Makes a graph and data file")
|
|
||||||
.build();
|
|
||||||
Option matchCDR3 = Option.builder("match")
|
|
||||||
.longOpt("match-cdr3")
|
|
||||||
.desc("Match CDR3s. Requires a cell sample file and any number of plate files.")
|
|
||||||
.build();
|
|
||||||
OptionGroup mainGroup = new OptionGroup();
|
|
||||||
mainGroup.addOption(makeCells);
|
|
||||||
mainGroup.addOption(makePlate);
|
|
||||||
mainGroup.addOption(makeGraph);
|
|
||||||
mainGroup.addOption(matchCDR3);
|
|
||||||
mainGroup.setRequired(true);
|
|
||||||
mainOptions.addOptionGroup(mainGroup);
|
|
||||||
|
|
||||||
//Reuse clones of this for other options groups, rather than making it lots of times
|
|
||||||
Option outputFile = Option.builder("o")
|
|
||||||
.longOpt("output-file")
|
|
||||||
.hasArg()
|
|
||||||
.argName("filename")
|
|
||||||
.desc("Name of output file")
|
|
||||||
.build();
|
|
||||||
mainOptions.addOption(outputFile);
|
|
||||||
|
|
||||||
//Options cellOptions = new Options();
|
|
||||||
Option numCells = Option.builder("nc")
|
|
||||||
.longOpt("num-cells")
|
|
||||||
.desc("The number of distinct cells to generate")
|
|
||||||
.hasArg()
|
|
||||||
.argName("number")
|
|
||||||
.build();
|
|
||||||
mainOptions.addOption(numCells);
|
|
||||||
Option cdr1Freq = Option.builder("d")
|
|
||||||
.longOpt("peptide-diversity-factor")
|
|
||||||
.hasArg()
|
|
||||||
.argName("number")
|
|
||||||
.desc("Number of distinct CDR3s for every CDR1")
|
|
||||||
.build();
|
|
||||||
mainOptions.addOption(cdr1Freq);
|
|
||||||
//Option cellOutput = (Option) outputFile.clone();
|
|
||||||
//cellOutput.setRequired(true);
|
|
||||||
//mainOptions.addOption(cellOutput);
|
|
||||||
|
|
||||||
//Options plateOptions = new Options();
|
|
||||||
Option inputCells = Option.builder("c")
|
|
||||||
.longOpt("cell-file")
|
|
||||||
.hasArg()
|
|
||||||
.argName("file")
|
|
||||||
.desc("The cell sample file used for filling wells")
|
|
||||||
.build();
|
|
||||||
mainOptions.addOption(inputCells);
|
|
||||||
Option numWells = Option.builder("w")
|
|
||||||
.longOpt("num-wells")
|
|
||||||
.hasArg()
|
|
||||||
.argName("number")
|
|
||||||
.desc("The number of wells on each plate")
|
|
||||||
.build();
|
|
||||||
mainOptions.addOption(numWells);
|
|
||||||
Option numPlates = Option.builder("np")
|
|
||||||
.longOpt("num-plates")
|
|
||||||
.hasArg()
|
|
||||||
.argName("number")
|
|
||||||
.desc("The number of plate files to output")
|
|
||||||
.build();
|
|
||||||
mainOptions.addOption(numPlates);
|
|
||||||
//Option plateOutput = (Option) outputFile.clone();
|
|
||||||
//plateOutput.setRequired(true);
|
|
||||||
//plateOutput.setDescription("Prefix for plate output filenames");
|
|
||||||
//mainOptions.addOption(plateOutput);
|
|
||||||
Option plateErr = Option.builder("err")
|
|
||||||
.longOpt("drop-out-rate")
|
|
||||||
.hasArg()
|
|
||||||
.argName("number")
|
|
||||||
.desc("Well drop-out rate. (Probability between 0 and 1)")
|
|
||||||
.build();
|
|
||||||
mainOptions.addOption(plateErr);
|
|
||||||
Option plateConcentrations = Option.builder("t")
|
|
||||||
.longOpt("t-cells-per-well")
|
|
||||||
.hasArgs()
|
|
||||||
.argName("number 1, number 2, ...")
|
|
||||||
.desc("Number of T cells per well for each plate section")
|
|
||||||
.build();
|
|
||||||
mainOptions.addOption(plateConcentrations);
|
|
||||||
|
|
||||||
//different distributions, mutually exclusive
|
|
||||||
OptionGroup plateDistributions = new OptionGroup();
|
|
||||||
Option plateExp = Option.builder("exponential")
|
|
||||||
.desc("Sample from distinct cells with exponential frequency distribution")
|
|
||||||
.build();
|
|
||||||
plateDistributions.addOption(plateExp);
|
|
||||||
Option plateGaussian = Option.builder("gaussian")
|
|
||||||
.desc("Sample from distinct cells with gaussain frequency distribution")
|
|
||||||
.build();
|
|
||||||
plateDistributions.addOption(plateGaussian);
|
|
||||||
Option platePoisson = Option.builder("poisson")
|
|
||||||
.desc("Sample from distinct cells with poisson frequency distribution")
|
|
||||||
.build();
|
|
||||||
plateDistributions.addOption(platePoisson);
|
|
||||||
mainOptions.addOptionGroup(plateDistributions);
|
|
||||||
|
|
||||||
Option plateStdDev = Option.builder("stddev")
|
|
||||||
.desc("Standard deviation for gaussian distribution")
|
|
||||||
.hasArg()
|
|
||||||
.argName("number")
|
|
||||||
.build();
|
|
||||||
mainOptions.addOption(plateStdDev);
|
|
||||||
|
|
||||||
Option plateLambda = Option.builder("lambda")
|
|
||||||
.desc("Lambda for exponential distribution")
|
|
||||||
.hasArg()
|
|
||||||
.argName("number")
|
|
||||||
.build();
|
|
||||||
mainOptions.addOption(plateLambda);
|
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
//
|
|
||||||
// String cellFile, String filename, Double stdDev,
|
|
||||||
// Integer numWells, Integer numSections,
|
|
||||||
// Integer[] concentrations, Double dropOutRate
|
|
||||||
//
|
|
||||||
|
|
||||||
//Options matchOptions = new Options();
|
|
||||||
inputCells.setDescription("The cell sample file to be used for matching.");
|
|
||||||
mainOptions.addOption(inputCells);
|
|
||||||
Option lowThresh = Option.builder("low")
|
|
||||||
.longOpt("low-threshold")
|
|
||||||
.hasArg()
|
|
||||||
.argName("number")
|
|
||||||
.desc("Sets the minimum occupancy overlap to attempt matching")
|
|
||||||
.build();
|
|
||||||
mainOptions.addOption(lowThresh);
|
|
||||||
Option highThresh = Option.builder("high")
|
|
||||||
.longOpt("high-threshold")
|
|
||||||
.hasArg()
|
|
||||||
.argName("number")
|
|
||||||
.desc("Sets the maximum occupancy overlap to attempt matching")
|
|
||||||
.build();
|
|
||||||
mainOptions.addOption(highThresh);
|
|
||||||
Option occDiff = Option.builder("occdiff")
|
|
||||||
.longOpt("occupancy-difference")
|
|
||||||
.hasArg()
|
|
||||||
.argName("Number")
|
|
||||||
.desc("Maximum difference in alpha/beta occupancy to attempt matching")
|
|
||||||
.build();
|
|
||||||
mainOptions.addOption(occDiff);
|
|
||||||
Option overlapPer = Option.builder("ovper")
|
|
||||||
.longOpt("overlap-percent")
|
|
||||||
.hasArg()
|
|
||||||
.argName("Percent")
|
|
||||||
.desc("Minimum overlap percent to attempt matching (0 -100)")
|
|
||||||
.build();
|
|
||||||
mainOptions.addOption(overlapPer);
|
|
||||||
Option inputPlates = Option.builder("p")
|
|
||||||
.longOpt("plate-files")
|
|
||||||
.hasArgs()
|
|
||||||
.desc("Plate files to match")
|
|
||||||
.build();
|
|
||||||
mainOptions.addOption(inputPlates);
|
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
CommandLineParser parser = new DefaultParser();
|
CommandLineParser parser = new DefaultParser();
|
||||||
try {
|
try{
|
||||||
CommandLine line = parser.parse(mainOptions, args);
|
CommandLine line = parser.parse(mainOptions, Arrays.copyOfRange(args, 0, 1));
|
||||||
if(line.hasOption("match")){
|
|
||||||
//line = parser.parse(mainOptions, args);
|
if (line.hasOption("help")) {
|
||||||
//String cellFile = line.getOptionValue("c");
|
HelpFormatter formatter = new HelpFormatter();
|
||||||
String graphFile = line.getOptionValue("g");
|
formatter.printHelp("BiGpairSEQ_Sim.jar", mainOptions);
|
||||||
Integer lowThreshold = Integer.valueOf(line.getOptionValue(lowThresh));
|
System.out.println();
|
||||||
Integer highThreshold = Integer.valueOf(line.getOptionValue(highThresh));
|
formatter.printHelp("BiGpairSEQ_Sim.jar -cells", cellOptions);
|
||||||
Integer occupancyDifference = Integer.valueOf(line.getOptionValue(occDiff));
|
System.out.println();
|
||||||
Integer overlapPercent = Integer.valueOf(line.getOptionValue(overlapPer));
|
formatter.printHelp("BiGpairSEQ_Sim.jar -plate", plateOptions);
|
||||||
for(String plate: line.getOptionValues("p")) {
|
System.out.println();
|
||||||
matchCDR3s(graphFile, lowThreshold, highThreshold, occupancyDifference, overlapPercent);
|
formatter.printHelp("BiGpairSEQ_Sim.jar -graph", graphOptions);
|
||||||
}
|
System.out.println();
|
||||||
|
formatter.printHelp("BiGpairSEQ_Sim.jar -match", matchOptions);
|
||||||
}
|
}
|
||||||
else if(line.hasOption("cells")){
|
else if (line.hasOption("version")) {
|
||||||
//line = parser.parse(mainOptions, args);
|
System.out.println("BiGpairSEQ_Sim " + BiGpairSEQ.getVersion());
|
||||||
|
}
|
||||||
|
else if (line.hasOption("cells")) {
|
||||||
|
line = parser.parse(cellOptions, Arrays.copyOfRange(args, 1, args.length));
|
||||||
|
Integer number = Integer.valueOf(line.getOptionValue("n"));
|
||||||
|
Integer diversity = Integer.valueOf(line.getOptionValue("d"));
|
||||||
String filename = line.getOptionValue("o");
|
String filename = line.getOptionValue("o");
|
||||||
Integer numDistCells = Integer.valueOf(line.getOptionValue("nc"));
|
makeCells(filename, number, diversity);
|
||||||
Integer freq = Integer.valueOf(line.getOptionValue("d"));
|
|
||||||
makeCells(filename, numDistCells, freq);
|
|
||||||
}
|
}
|
||||||
else if(line.hasOption("plates")){
|
|
||||||
//line = parser.parse(mainOptions, args);
|
|
||||||
String cellFile = line.getOptionValue("c");
|
|
||||||
String filenamePrefix = line.getOptionValue("o");
|
|
||||||
Integer numWellsOnPlate = Integer.valueOf(line.getOptionValue("w"));
|
|
||||||
Integer numPlatesToMake = Integer.valueOf(line.getOptionValue("np"));
|
|
||||||
String[] concentrationsToUseString = line.getOptionValues("t");
|
|
||||||
Integer numSections = concentrationsToUseString.length;
|
|
||||||
|
|
||||||
Integer[] concentrationsToUse = new Integer[numSections];
|
else if (line.hasOption("plate")) {
|
||||||
for(int i = 0; i <numSections; i++){
|
line = parser.parse(plateOptions, Arrays.copyOfRange(args, 1, args.length));
|
||||||
concentrationsToUse[i] = Integer.valueOf(concentrationsToUseString[i]);
|
//get the cells
|
||||||
|
String cellFilename = line.getOptionValue("c");
|
||||||
|
CellSample cells = getCells(cellFilename);
|
||||||
|
//get the rest of the parameters
|
||||||
|
Integer[] populations;
|
||||||
|
String outputFilename = line.getOptionValue("o");
|
||||||
|
Integer numWells = Integer.parseInt(line.getOptionValue("w"));
|
||||||
|
Double dropoutRate = Double.parseDouble(line.getOptionValue("err"));
|
||||||
|
if (line.hasOption("random")) {
|
||||||
|
//Array holding values of minimum and maximum populations
|
||||||
|
Integer[] min_max = Stream.of(line.getOptionValues("random"))
|
||||||
|
.mapToInt(Integer::parseInt)
|
||||||
|
.boxed()
|
||||||
|
.toArray(Integer[]::new);
|
||||||
|
populations = BiGpairSEQ.getRand().ints(min_max[0], min_max[1] + 1)
|
||||||
|
.limit(numWells)
|
||||||
|
.boxed()
|
||||||
|
.toArray(Integer[]::new);
|
||||||
}
|
}
|
||||||
Double dropOutRate = Double.valueOf(line.getOptionValue("err"));
|
else if (line.hasOption("pop")) {
|
||||||
if(line.hasOption("exponential")){
|
populations = Stream.of(line.getOptionValues("pop"))
|
||||||
Double lambda = Double.valueOf(line.getOptionValue("lambda"));
|
.mapToInt(Integer::parseInt)
|
||||||
for(int i = 1; i <= numPlatesToMake; i++){
|
.boxed()
|
||||||
makePlateExp(cellFile, filenamePrefix + i, lambda, numWellsOnPlate,
|
.toArray(Integer[]::new);
|
||||||
concentrationsToUse,dropOutRate);
|
|
||||||
}
|
|
||||||
}
|
}
|
||||||
else if(line.hasOption("gaussian")){
|
else{
|
||||||
Double stdDev = Double.valueOf(line.getOptionValue("std-dev"));
|
populations = new Integer[1];
|
||||||
for(int i = 1; i <= numPlatesToMake; i++){
|
populations[0] = 1;
|
||||||
makePlate(cellFile, filenamePrefix + i, stdDev, numWellsOnPlate,
|
}
|
||||||
concentrationsToUse,dropOutRate);
|
//make the plate
|
||||||
}
|
Plate plate;
|
||||||
|
if (line.hasOption("poisson")) {
|
||||||
|
Double stdDev = Math.sqrt(numWells);
|
||||||
|
plate = new Plate(cells, cellFilename, numWells, populations, dropoutRate, stdDev, false);
|
||||||
|
}
|
||||||
|
else if (line.hasOption("gaussian")) {
|
||||||
|
Double stdDev = Double.parseDouble(line.getOptionValue("stddev"));
|
||||||
|
plate = new Plate(cells, cellFilename, numWells, populations, dropoutRate, stdDev, false);
|
||||||
|
}
|
||||||
|
else {
|
||||||
|
assert line.hasOption("exponential");
|
||||||
|
Double lambda = Double.parseDouble(line.getOptionValue("lambda"));
|
||||||
|
plate = new Plate(cells, cellFilename, numWells, populations, dropoutRate, lambda, true);
|
||||||
|
}
|
||||||
|
PlateFileWriter writer = new PlateFileWriter(outputFilename, plate);
|
||||||
|
writer.writePlateFile();
|
||||||
|
}
|
||||||
|
|
||||||
|
else if (line.hasOption("graph")) { //Making a graph
|
||||||
|
line = parser.parse(graphOptions, Arrays.copyOfRange(args, 1, args.length));
|
||||||
|
String cellFilename = line.getOptionValue("c");
|
||||||
|
String plateFilename = line.getOptionValue("p");
|
||||||
|
String outputFilename = line.getOptionValue("o");
|
||||||
|
//get cells
|
||||||
|
CellSample cells = getCells(cellFilename);
|
||||||
|
//get plate
|
||||||
|
Plate plate = getPlate(plateFilename);
|
||||||
|
GraphWithMapData graph = Simulator.makeGraph(cells, plate, false);
|
||||||
|
if (!line.hasOption("no-binary")) { //output binary file unless told not to
|
||||||
|
GraphDataObjectWriter writer = new GraphDataObjectWriter(outputFilename, graph, false);
|
||||||
|
writer.writeDataToFile();
|
||||||
}
|
}
|
||||||
else if(line.hasOption("poisson")){
|
if (line.hasOption("graphml")) { //if told to, output graphml file
|
||||||
for(int i = 1; i <= numPlatesToMake; i++){
|
GraphMLFileWriter gmlwriter = new GraphMLFileWriter(outputFilename, graph);
|
||||||
makePlatePoisson(cellFile, filenamePrefix + i, numWellsOnPlate,
|
gmlwriter.writeGraphToFile();
|
||||||
concentrationsToUse,dropOutRate);
|
}
|
||||||
|
}
|
||||||
|
|
||||||
|
else if (line.hasOption("match")) { //can add a flag for which match type in future, spit this in two
|
||||||
|
line = parser.parse(matchOptions, Arrays.copyOfRange(args, 1, args.length));
|
||||||
|
String graphFilename = line.getOptionValue("g");
|
||||||
|
|
||||||
|
String outputFilename;
|
||||||
|
if(line.hasOption("o")) {
|
||||||
|
outputFilename = line.getOptionValue("o");
|
||||||
|
}
|
||||||
|
else {
|
||||||
|
outputFilename = null;
|
||||||
|
}
|
||||||
|
Integer minThreshold = Integer.parseInt(line.getOptionValue("min"));
|
||||||
|
Integer maxThreshold = Integer.parseInt(line.getOptionValue("max"));
|
||||||
|
int minOverlapPct;
|
||||||
|
if (line.hasOption("minpct")) { //see if this filter is being used
|
||||||
|
minOverlapPct = Integer.parseInt(line.getOptionValue("minpct"));
|
||||||
|
}
|
||||||
|
else {
|
||||||
|
minOverlapPct = 0;
|
||||||
|
}
|
||||||
|
int maxOccupancyDiff;
|
||||||
|
if (line.hasOption("maxdiff")) { //see if this filter is being used
|
||||||
|
maxOccupancyDiff = Integer.parseInt(line.getOptionValue("maxdiff"));
|
||||||
|
}
|
||||||
|
else {
|
||||||
|
maxOccupancyDiff = Integer.MAX_VALUE;
|
||||||
|
}
|
||||||
|
GraphWithMapData graph = getGraph(graphFilename);
|
||||||
|
MatchingResult result = Simulator.matchCDR3s(graph, graphFilename, minThreshold, maxThreshold,
|
||||||
|
maxOccupancyDiff, minOverlapPct, false);
|
||||||
|
if(outputFilename != null){
|
||||||
|
MatchingFileWriter writer = new MatchingFileWriter(outputFilename, result);
|
||||||
|
writer.writeResultsToFile();
|
||||||
|
}
|
||||||
|
//can put a bunch of ifs for outputting various things from the MatchingResult to System.out here
|
||||||
|
//after I put those flags in the matchOptions
|
||||||
|
if(line.hasOption("print-metadata")) {
|
||||||
|
for (String k : result.getMetadata().keySet()) {
|
||||||
|
System.out.println(k + ": " + result.getMetadata().get(k));
|
||||||
}
|
}
|
||||||
}
|
}
|
||||||
|
if(line.hasOption("print-error")) {
|
||||||
|
System.out.println("pairing error rate: " + result.getPairingErrorRate());
|
||||||
|
}
|
||||||
|
if(line.hasOption("print-attempt")) {
|
||||||
|
System.out.println("pairing attempt rate: " +result.getPairingAttemptRate());
|
||||||
|
}
|
||||||
|
if(line.hasOption("print-correct")) {
|
||||||
|
System.out.println("correct pairings: " + result.getCorrectPairingCount());
|
||||||
|
}
|
||||||
|
if(line.hasOption("print-incorrect")) {
|
||||||
|
System.out.println("incorrect pairings: " + result.getIncorrectPairingCount());
|
||||||
|
}
|
||||||
|
if(line.hasOption("print-alphas")) {
|
||||||
|
System.out.println("total alphas found: " + result.getAlphaCount());
|
||||||
|
}
|
||||||
|
if(line.hasOption("print-betas")) {
|
||||||
|
System.out.println("total betas found: " + result.getBetaCount());
|
||||||
|
}
|
||||||
|
if(line.hasOption("print-time")) {
|
||||||
|
System.out.println("simulation time (seconds): " + result.getSimulationTime());
|
||||||
|
}
|
||||||
}
|
}
|
||||||
}
|
}
|
||||||
catch (ParseException exp) {
|
catch (ParseException exp) {
|
||||||
@@ -286,43 +222,278 @@ public class CommandLineInterface {
|
|||||||
}
|
}
|
||||||
}
|
}
|
||||||
|
|
||||||
|
private static Option outputFileOption() {
|
||||||
|
Option outputFile = Option.builder("o")
|
||||||
|
.longOpt("output-file")
|
||||||
|
.hasArg()
|
||||||
|
.argName("filename")
|
||||||
|
.desc("Name of output file")
|
||||||
|
.required()
|
||||||
|
.build();
|
||||||
|
return outputFile;
|
||||||
|
}
|
||||||
|
|
||||||
|
private static Options buildMainOptions() {
|
||||||
|
Options mainOptions = new Options();
|
||||||
|
Option help = Option.builder("help")
|
||||||
|
.desc("Displays this help menu")
|
||||||
|
.build();
|
||||||
|
Option makeCells = Option.builder("cells")
|
||||||
|
.longOpt("make-cells")
|
||||||
|
.desc("Makes a cell sample file of distinct T cells")
|
||||||
|
.build();
|
||||||
|
Option makePlate = Option.builder("plate")
|
||||||
|
.longOpt("make-plate")
|
||||||
|
.desc("Makes a sample plate file. Requires a cell sample file.")
|
||||||
|
.build();
|
||||||
|
Option makeGraph = Option.builder("graph")
|
||||||
|
.longOpt("make-graph")
|
||||||
|
.desc("Makes a graph/data file. Requires a cell sample file and a sample plate file")
|
||||||
|
.build();
|
||||||
|
Option matchCDR3 = Option.builder("match")
|
||||||
|
.longOpt("match-cdr3")
|
||||||
|
.desc("Matches CDR3s. Requires a graph/data file.")
|
||||||
|
.build();
|
||||||
|
Option printVersion = Option.builder("version")
|
||||||
|
.desc("Prints the program version number to stdout").build();
|
||||||
|
OptionGroup mainGroup = new OptionGroup();
|
||||||
|
mainGroup.addOption(help);
|
||||||
|
mainGroup.addOption(printVersion);
|
||||||
|
mainGroup.addOption(makeCells);
|
||||||
|
mainGroup.addOption(makePlate);
|
||||||
|
mainGroup.addOption(makeGraph);
|
||||||
|
mainGroup.addOption(matchCDR3);
|
||||||
|
mainGroup.setRequired(true);
|
||||||
|
mainOptions.addOptionGroup(mainGroup);
|
||||||
|
return mainOptions;
|
||||||
|
}
|
||||||
|
|
||||||
|
private static Options buildCellOptions() {
|
||||||
|
Options cellOptions = new Options();
|
||||||
|
Option numCells = Option.builder("n")
|
||||||
|
.longOpt("num-cells")
|
||||||
|
.desc("The number of distinct cells to generate")
|
||||||
|
.hasArg()
|
||||||
|
.argName("number")
|
||||||
|
.required().build();
|
||||||
|
Option cdr3Diversity = Option.builder("d")
|
||||||
|
.longOpt("diversity-factor")
|
||||||
|
.desc("The factor by which unique CDR3s outnumber unique CDR1s")
|
||||||
|
.hasArg()
|
||||||
|
.argName("factor")
|
||||||
|
.required().build();
|
||||||
|
cellOptions.addOption(numCells);
|
||||||
|
cellOptions.addOption(cdr3Diversity);
|
||||||
|
cellOptions.addOption(outputFileOption());
|
||||||
|
return cellOptions;
|
||||||
|
}
|
||||||
|
|
||||||
|
private static Options buildPlateOptions() {
|
||||||
|
Options plateOptions = new Options();
|
||||||
|
Option cellFile = Option.builder("c") // add this to plate options
|
||||||
|
.longOpt("cell-file")
|
||||||
|
.desc("The cell sample file to use")
|
||||||
|
.hasArg()
|
||||||
|
.argName("filename")
|
||||||
|
.required().build();
|
||||||
|
Option numWells = Option.builder("w")// add this to plate options
|
||||||
|
.longOpt("wells")
|
||||||
|
.desc("The number of wells on the sample plate")
|
||||||
|
.hasArg()
|
||||||
|
.argName("number")
|
||||||
|
.required().build();
|
||||||
|
//options group for choosing with distribution to use
|
||||||
|
OptionGroup distributions = new OptionGroup();// add this to plate options
|
||||||
|
distributions.setRequired(true);
|
||||||
|
Option poisson = Option.builder("poisson")
|
||||||
|
.desc("Use a Poisson distribution for cell sample")
|
||||||
|
.build();
|
||||||
|
Option gaussian = Option.builder("gaussian")
|
||||||
|
.desc("Use a Gaussian distribution for cell sample")
|
||||||
|
.build();
|
||||||
|
Option exponential = Option.builder("exponential")
|
||||||
|
.desc("Use an exponential distribution for cell sample")
|
||||||
|
.build();
|
||||||
|
distributions.addOption(poisson);
|
||||||
|
distributions.addOption(gaussian);
|
||||||
|
distributions.addOption(exponential);
|
||||||
|
//options group for statistical distribution parameters
|
||||||
|
OptionGroup statParams = new OptionGroup();// add this to plate options
|
||||||
|
Option stdDev = Option.builder("stddev")
|
||||||
|
.desc("If using -gaussian flag, standard deviation for distrbution")
|
||||||
|
.hasArg()
|
||||||
|
.argName("value")
|
||||||
|
.build();
|
||||||
|
Option lambda = Option.builder("lambda")
|
||||||
|
.desc("If using -exponential flag, lambda value for distribution")
|
||||||
|
.hasArg()
|
||||||
|
.argName("value")
|
||||||
|
.build();
|
||||||
|
statParams.addOption(stdDev);
|
||||||
|
statParams.addOption(lambda);
|
||||||
|
//Option group for random plate or set populations
|
||||||
|
OptionGroup wellPopOptions = new OptionGroup(); // add this to plate options
|
||||||
|
wellPopOptions.setRequired(true);
|
||||||
|
Option randomWellPopulations = Option.builder("random")
|
||||||
|
.desc("Randomize well populations on sample plate. Takes two arguments: the minimum possible population and the maximum possible population.")
|
||||||
|
.hasArgs()
|
||||||
|
.numberOfArgs(2)
|
||||||
|
.argName("min> <max")
|
||||||
|
.build();
|
||||||
|
Option specificWellPopulations = Option.builder("pop")
|
||||||
|
.desc("The well populations for each section of the sample plate. There will be as many sections as there are populations given.")
|
||||||
|
.hasArgs()
|
||||||
|
.argName("number [number]...")
|
||||||
|
.build();
|
||||||
|
Option dropoutRate = Option.builder("err") //add this to plate options
|
||||||
|
.hasArg()
|
||||||
|
.desc("The sequence dropout rate due to amplification error. (0.0 - 1.0)")
|
||||||
|
.argName("rate")
|
||||||
|
.required()
|
||||||
|
.build();
|
||||||
|
wellPopOptions.addOption(randomWellPopulations);
|
||||||
|
wellPopOptions.addOption(specificWellPopulations);
|
||||||
|
plateOptions.addOption(cellFile);
|
||||||
|
plateOptions.addOption(numWells);
|
||||||
|
plateOptions.addOptionGroup(distributions);
|
||||||
|
plateOptions.addOptionGroup(statParams);
|
||||||
|
plateOptions.addOptionGroup(wellPopOptions);
|
||||||
|
plateOptions.addOption(dropoutRate);
|
||||||
|
plateOptions.addOption(outputFileOption());
|
||||||
|
return plateOptions;
|
||||||
|
}
|
||||||
|
|
||||||
|
private static Options buildGraphOptions() {
|
||||||
|
Options graphOptions = new Options();
|
||||||
|
Option cellFilename = Option.builder("c")
|
||||||
|
.longOpt("cell-file")
|
||||||
|
.desc("Cell sample file to use for checking pairing accuracy")
|
||||||
|
.hasArg()
|
||||||
|
.argName("filename")
|
||||||
|
.required().build();
|
||||||
|
Option plateFilename = Option.builder("p")
|
||||||
|
.longOpt("plate-filename")
|
||||||
|
.desc("Sample plate file from which to construct graph")
|
||||||
|
.hasArg()
|
||||||
|
.argName("filename")
|
||||||
|
.required().build();
|
||||||
|
Option outputGraphML = Option.builder("graphml")
|
||||||
|
.desc("(Optional) Output GraphML file")
|
||||||
|
.build();
|
||||||
|
Option outputSerializedBinary = Option.builder("nb")
|
||||||
|
.longOpt("no-binary")
|
||||||
|
.desc("(Optional) Don't output serialized binary file")
|
||||||
|
.build();
|
||||||
|
graphOptions.addOption(cellFilename);
|
||||||
|
graphOptions.addOption(plateFilename);
|
||||||
|
graphOptions.addOption(outputFileOption());
|
||||||
|
graphOptions.addOption(outputGraphML);
|
||||||
|
graphOptions.addOption(outputSerializedBinary);
|
||||||
|
return graphOptions;
|
||||||
|
}
|
||||||
|
|
||||||
|
private static Options buildMatchCDR3options() {
|
||||||
|
Options matchCDR3options = new Options();
|
||||||
|
Option graphFilename = Option.builder("g")
|
||||||
|
.longOpt("graph-file")
|
||||||
|
.desc("The graph/data file to use")
|
||||||
|
.hasArg()
|
||||||
|
.argName("filename")
|
||||||
|
.required().build();
|
||||||
|
Option minOccupancyOverlap = Option.builder("min")
|
||||||
|
.desc("The minimum number of shared wells to attempt to match a sequence pair")
|
||||||
|
.hasArg()
|
||||||
|
.argName("number")
|
||||||
|
.required().build();
|
||||||
|
Option maxOccupancyOverlap = Option.builder("max")
|
||||||
|
.desc("The maximum number of shared wells to attempt to match a sequence pair")
|
||||||
|
.hasArg()
|
||||||
|
.argName("number")
|
||||||
|
.required().build();
|
||||||
|
Option minOverlapPercent = Option.builder("minpct")
|
||||||
|
.desc("(Optional) The minimum percentage of a sequence's total occupancy shared by another sequence to attempt matching. (0 - 100) ")
|
||||||
|
.hasArg()
|
||||||
|
.argName("percent")
|
||||||
|
.build();
|
||||||
|
Option maxOccupancyDifference = Option.builder("maxdiff")
|
||||||
|
.desc("(Optional) The maximum difference in total occupancy between two sequences to attempt matching.")
|
||||||
|
.hasArg()
|
||||||
|
.argName("number")
|
||||||
|
.build();
|
||||||
|
Option outputFile = Option.builder("o") //can't call the method this time, because this one's optional
|
||||||
|
.longOpt("output-file")
|
||||||
|
.hasArg()
|
||||||
|
.argName("filename")
|
||||||
|
.desc("(Optional) Name of output the output file. If not present, no file will be written.")
|
||||||
|
.build();
|
||||||
|
matchCDR3options.addOption(graphFilename)
|
||||||
|
.addOption(minOccupancyOverlap)
|
||||||
|
.addOption(maxOccupancyOverlap)
|
||||||
|
.addOption(minOverlapPercent)
|
||||||
|
.addOption(maxOccupancyDifference)
|
||||||
|
.addOption(outputFile);
|
||||||
|
|
||||||
|
//options for output to System.out
|
||||||
|
Option printAlphaCount = Option.builder().longOpt("print-alphas")
|
||||||
|
.desc("(Optional) Print the number of distinct alpha sequences to stdout.").build();
|
||||||
|
Option printBetaCount = Option.builder().longOpt("print-betas")
|
||||||
|
.desc("(Optional) Print the number of distinct beta sequences to stdout.").build();
|
||||||
|
Option printTime = Option.builder().longOpt("print-time")
|
||||||
|
.desc("(Optional) Print the total simulation time to stdout.").build();
|
||||||
|
Option printErrorRate = Option.builder().longOpt("print-error")
|
||||||
|
.desc("(Optional) Print the pairing error rate to stdout").build();
|
||||||
|
Option printAttempt = Option.builder().longOpt("print-attempt")
|
||||||
|
.desc("(Optional) Print the pairing attempt rate to stdout").build();
|
||||||
|
Option printCorrect = Option.builder().longOpt("print-correct")
|
||||||
|
.desc("(Optional) Print the number of correct pairs to stdout").build();
|
||||||
|
Option printIncorrect = Option.builder().longOpt("print-incorrect")
|
||||||
|
.desc("(Optional) Print the number of incorrect pairs to stdout").build();
|
||||||
|
Option printMetadata = Option.builder().longOpt("print-metadata")
|
||||||
|
.desc("(Optional) Print a full summary of the matching results to stdout.").build();
|
||||||
|
|
||||||
|
matchCDR3options
|
||||||
|
.addOption(printErrorRate)
|
||||||
|
.addOption(printAttempt)
|
||||||
|
.addOption(printCorrect)
|
||||||
|
.addOption(printIncorrect)
|
||||||
|
.addOption(printMetadata)
|
||||||
|
.addOption(printAlphaCount)
|
||||||
|
.addOption(printBetaCount)
|
||||||
|
.addOption(printTime);
|
||||||
|
return matchCDR3options;
|
||||||
|
}
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
private static CellSample getCells(String cellFilename) {
|
||||||
|
assert cellFilename != null;
|
||||||
|
CellFileReader reader = new CellFileReader(cellFilename);
|
||||||
|
return reader.getCellSample();
|
||||||
|
}
|
||||||
|
|
||||||
|
private static Plate getPlate(String plateFilename) {
|
||||||
|
assert plateFilename != null;
|
||||||
|
PlateFileReader reader = new PlateFileReader(plateFilename);
|
||||||
|
return reader.getSamplePlate();
|
||||||
|
}
|
||||||
|
|
||||||
|
private static GraphWithMapData getGraph(String graphFilename) {
|
||||||
|
assert graphFilename != null;
|
||||||
|
try{
|
||||||
|
GraphDataObjectReader reader = new GraphDataObjectReader(graphFilename, false);
|
||||||
|
return reader.getData();
|
||||||
|
|
||||||
|
}
|
||||||
|
catch (IOException ex) {
|
||||||
|
ex.printStackTrace();
|
||||||
|
return null;
|
||||||
|
}
|
||||||
|
}
|
||||||
|
|
||||||
//for calling from command line
|
//for calling from command line
|
||||||
public static void makeCells(String filename, Integer numCells, Integer cdr1Freq){
|
public static void makeCells(String filename, Integer numCells, Integer cdr1Freq) {
|
||||||
CellSample sample = Simulator.generateCellSample(numCells, cdr1Freq);
|
CellSample sample = new CellSample(numCells, cdr1Freq);
|
||||||
CellFileWriter writer = new CellFileWriter(filename, sample);
|
CellFileWriter writer = new CellFileWriter(filename, sample);
|
||||||
writer.writeCellsToFile();
|
writer.writeCellsToFile();
|
||||||
}
|
}
|
||||||
|
|
||||||
public static void makePlateExp(String cellFile, String filename, Double lambda,
|
|
||||||
Integer numWells, Integer[] concentrations, Double dropOutRate){
|
|
||||||
CellFileReader cellReader = new CellFileReader(cellFile);
|
|
||||||
Plate samplePlate = new Plate(numWells, dropOutRate, concentrations);
|
|
||||||
samplePlate.fillWellsExponential(cellReader.getFilename(), cellReader.getListOfDistinctCellsDEPRECATED(), lambda);
|
|
||||||
PlateFileWriter writer = new PlateFileWriter(filename, samplePlate);
|
|
||||||
writer.writePlateFile();
|
|
||||||
}
|
|
||||||
|
|
||||||
private static void makePlatePoisson(String cellFile, String filename, Integer numWells,
|
|
||||||
Integer[] concentrations, Double dropOutRate){
|
|
||||||
CellFileReader cellReader = new CellFileReader(cellFile);
|
|
||||||
Double stdDev = Math.sqrt(cellReader.getCellCountDEPRECATED());
|
|
||||||
Plate samplePlate = new Plate(numWells, dropOutRate, concentrations);
|
|
||||||
samplePlate.fillWells(cellReader.getFilename(), cellReader.getListOfDistinctCellsDEPRECATED(), stdDev);
|
|
||||||
PlateFileWriter writer = new PlateFileWriter(filename, samplePlate);
|
|
||||||
writer.writePlateFile();
|
|
||||||
}
|
|
||||||
|
|
||||||
private static void makePlate(String cellFile, String filename, Double stdDev,
|
|
||||||
Integer numWells, Integer[] concentrations, Double dropOutRate){
|
|
||||||
CellFileReader cellReader = new CellFileReader(cellFile);
|
|
||||||
Plate samplePlate = new Plate(numWells, dropOutRate, concentrations);
|
|
||||||
samplePlate.fillWells(cellReader.getFilename(), cellReader.getListOfDistinctCellsDEPRECATED(), stdDev);
|
|
||||||
PlateFileWriter writer = new PlateFileWriter(filename, samplePlate);
|
|
||||||
writer.writePlateFile();
|
|
||||||
}
|
|
||||||
|
|
||||||
private static void matchCDR3s(String graphFile, Integer lowThreshold, Integer highThreshold,
|
|
||||||
Integer occupancyDifference, Integer overlapPercent) {
|
|
||||||
|
|
||||||
}
|
|
||||||
}
|
}
|
||||||
|
|||||||
@@ -1,10 +1,12 @@
|
|||||||
import java.io.*;
|
import java.io.*;
|
||||||
|
|
||||||
public class GraphDataObjectReader {
|
public class GraphDataObjectReader {
|
||||||
|
|
||||||
private GraphWithMapData data;
|
private GraphWithMapData data;
|
||||||
private String filename;
|
private String filename;
|
||||||
|
|
||||||
public GraphDataObjectReader(String filename) throws IOException {
|
|
||||||
|
public GraphDataObjectReader(String filename, boolean verbose) throws IOException {
|
||||||
if(!filename.matches(".*\\.ser")){
|
if(!filename.matches(".*\\.ser")){
|
||||||
filename = filename + ".ser";
|
filename = filename + ".ser";
|
||||||
}
|
}
|
||||||
@@ -13,10 +15,13 @@ public class GraphDataObjectReader {
|
|||||||
BufferedInputStream fileIn = new BufferedInputStream(new FileInputStream(filename));
|
BufferedInputStream fileIn = new BufferedInputStream(new FileInputStream(filename));
|
||||||
ObjectInputStream in = new ObjectInputStream(fileIn))
|
ObjectInputStream in = new ObjectInputStream(fileIn))
|
||||||
{
|
{
|
||||||
System.out.println("Reading graph data from file. This may take some time");
|
if (verbose) {
|
||||||
System.out.println("File I/O time is not included in results");
|
System.out.println("Reading graph data from file. This may take some time");
|
||||||
|
System.out.println("File I/O time is not included in results");
|
||||||
|
}
|
||||||
data = (GraphWithMapData) in.readObject();
|
data = (GraphWithMapData) in.readObject();
|
||||||
} catch (FileNotFoundException | ClassNotFoundException ex) {
|
} catch (FileNotFoundException | ClassNotFoundException ex) {
|
||||||
|
System.out.println("Graph/data file " + filename + " not found.");
|
||||||
ex.printStackTrace();
|
ex.printStackTrace();
|
||||||
}
|
}
|
||||||
}
|
}
|
||||||
|
|||||||
@@ -1,3 +1,5 @@
|
|||||||
|
import org.jgrapht.Graph;
|
||||||
|
|
||||||
import java.io.BufferedOutputStream;
|
import java.io.BufferedOutputStream;
|
||||||
import java.io.FileOutputStream;
|
import java.io.FileOutputStream;
|
||||||
import java.io.IOException;
|
import java.io.IOException;
|
||||||
@@ -7,6 +9,7 @@ public class GraphDataObjectWriter {
|
|||||||
|
|
||||||
private GraphWithMapData data;
|
private GraphWithMapData data;
|
||||||
private String filename;
|
private String filename;
|
||||||
|
private boolean verbose = true;
|
||||||
|
|
||||||
public GraphDataObjectWriter(String filename, GraphWithMapData data) {
|
public GraphDataObjectWriter(String filename, GraphWithMapData data) {
|
||||||
if(!filename.matches(".*\\.ser")){
|
if(!filename.matches(".*\\.ser")){
|
||||||
@@ -16,13 +19,24 @@ public class GraphDataObjectWriter {
|
|||||||
this.data = data;
|
this.data = data;
|
||||||
}
|
}
|
||||||
|
|
||||||
|
public GraphDataObjectWriter(String filename, GraphWithMapData data, boolean verbose) {
|
||||||
|
this.verbose = verbose;
|
||||||
|
if(!filename.matches(".*\\.ser")){
|
||||||
|
filename = filename + ".ser";
|
||||||
|
}
|
||||||
|
this.filename = filename;
|
||||||
|
this.data = data;
|
||||||
|
}
|
||||||
|
|
||||||
public void writeDataToFile() {
|
public void writeDataToFile() {
|
||||||
try (BufferedOutputStream bufferedOut = new BufferedOutputStream(new FileOutputStream(filename));
|
try (BufferedOutputStream bufferedOut = new BufferedOutputStream(new FileOutputStream(filename));
|
||||||
|
|
||||||
ObjectOutputStream out = new ObjectOutputStream(bufferedOut);
|
ObjectOutputStream out = new ObjectOutputStream(bufferedOut);
|
||||||
){
|
){
|
||||||
System.out.println("Writing graph and occupancy data to file. This may take some time.");
|
if(verbose) {
|
||||||
System.out.println("File I/O time is not included in results.");
|
System.out.println("Writing graph and occupancy data to file. This may take some time.");
|
||||||
|
System.out.println("File I/O time is not included in results.");
|
||||||
|
}
|
||||||
out.writeObject(data);
|
out.writeObject(data);
|
||||||
} catch (IOException ex) {
|
} catch (IOException ex) {
|
||||||
ex.printStackTrace();
|
ex.printStackTrace();
|
||||||
|
|||||||
@@ -1,35 +0,0 @@
|
|||||||
import org.jgrapht.graph.SimpleWeightedGraph;
|
|
||||||
import org.jgrapht.nio.graphml.GraphMLImporter;
|
|
||||||
|
|
||||||
import java.io.BufferedReader;
|
|
||||||
import java.io.IOException;
|
|
||||||
import java.nio.file.Files;
|
|
||||||
import java.nio.file.Path;
|
|
||||||
|
|
||||||
public class GraphMLFileReader {
|
|
||||||
|
|
||||||
private String filename;
|
|
||||||
private SimpleWeightedGraph graph;
|
|
||||||
|
|
||||||
public GraphMLFileReader(String filename, SimpleWeightedGraph graph) {
|
|
||||||
if(!filename.matches(".*\\.graphml")){
|
|
||||||
filename = filename + ".graphml";
|
|
||||||
}
|
|
||||||
this.filename = filename;
|
|
||||||
this.graph = graph;
|
|
||||||
|
|
||||||
try(//don't need to close reader bc of try-with-resources auto-closing
|
|
||||||
BufferedReader reader = Files.newBufferedReader(Path.of(filename));
|
|
||||||
){
|
|
||||||
GraphMLImporter<SimpleWeightedGraph, BufferedReader> importer = new GraphMLImporter<>();
|
|
||||||
importer.importGraph(graph, reader);
|
|
||||||
}
|
|
||||||
catch (IOException ex) {
|
|
||||||
System.out.println("Graph file " + filename + " not found.");
|
|
||||||
System.err.println(ex);
|
|
||||||
}
|
|
||||||
}
|
|
||||||
|
|
||||||
public SimpleWeightedGraph getGraph() { return graph; }
|
|
||||||
|
|
||||||
}
|
|
||||||
@@ -1,20 +1,38 @@
|
|||||||
|
import org.jgrapht.graph.DefaultWeightedEdge;
|
||||||
import org.jgrapht.graph.SimpleWeightedGraph;
|
import org.jgrapht.graph.SimpleWeightedGraph;
|
||||||
import org.jgrapht.nio.dot.DOTExporter;
|
import org.jgrapht.nio.Attribute;
|
||||||
|
import org.jgrapht.nio.AttributeType;
|
||||||
|
import org.jgrapht.nio.DefaultAttribute;
|
||||||
import org.jgrapht.nio.graphml.GraphMLExporter;
|
import org.jgrapht.nio.graphml.GraphMLExporter;
|
||||||
|
import org.jgrapht.nio.graphml.GraphMLExporter.AttributeCategory;
|
||||||
|
import org.w3c.dom.Attr;
|
||||||
|
|
||||||
import java.io.BufferedWriter;
|
import java.io.BufferedWriter;
|
||||||
import java.io.IOException;
|
import java.io.IOException;
|
||||||
import java.nio.file.Files;
|
import java.nio.file.Files;
|
||||||
import java.nio.file.Path;
|
import java.nio.file.Path;
|
||||||
import java.nio.file.StandardOpenOption;
|
import java.nio.file.StandardOpenOption;
|
||||||
|
import java.util.HashMap;
|
||||||
|
import java.util.Map;
|
||||||
|
|
||||||
public class GraphMLFileWriter {
|
public class GraphMLFileWriter {
|
||||||
|
|
||||||
String filename;
|
String filename;
|
||||||
SimpleWeightedGraph graph;
|
SimpleWeightedGraph graph;
|
||||||
|
GraphWithMapData data;
|
||||||
|
Map<String, Attribute> graphAttributes;
|
||||||
|
|
||||||
|
public GraphMLFileWriter(String filename, GraphWithMapData data) {
|
||||||
|
if(!filename.matches(".*\\.graphml")){
|
||||||
|
filename = filename + ".graphml";
|
||||||
|
}
|
||||||
|
this.filename = filename;
|
||||||
|
this.data = data;
|
||||||
|
this.graph = data.getGraph();
|
||||||
|
graphAttributes = createGraphAttributes();
|
||||||
|
}
|
||||||
|
|
||||||
public GraphMLFileWriter(String filename, SimpleWeightedGraph graph) {
|
public GraphMLFileWriter(String filename, SimpleWeightedGraph<Vertex, DefaultWeightedEdge> graph) {
|
||||||
if(!filename.matches(".*\\.graphml")){
|
if(!filename.matches(".*\\.graphml")){
|
||||||
filename = filename + ".graphml";
|
filename = filename + ".graphml";
|
||||||
}
|
}
|
||||||
@@ -22,10 +40,50 @@ public class GraphMLFileWriter {
|
|||||||
this.graph = graph;
|
this.graph = graph;
|
||||||
}
|
}
|
||||||
|
|
||||||
|
private Map<String, Attribute> createGraphAttributes(){
|
||||||
|
Map<String, Attribute> ga = new HashMap<>();
|
||||||
|
//Sample plate filename
|
||||||
|
ga.put("sample plate filename", DefaultAttribute.createAttribute(data.getSourceFilename()));
|
||||||
|
// Number of wells
|
||||||
|
ga.put("well count", DefaultAttribute.createAttribute(data.getNumWells().toString()));
|
||||||
|
//Well populations
|
||||||
|
Integer[] wellPopulations = data.getWellPopulations();
|
||||||
|
StringBuilder populationsStringBuilder = new StringBuilder();
|
||||||
|
populationsStringBuilder.append(wellPopulations[0].toString());
|
||||||
|
for(int i = 1; i < wellPopulations.length; i++){
|
||||||
|
populationsStringBuilder.append(", ");
|
||||||
|
populationsStringBuilder.append(wellPopulations[i].toString());
|
||||||
|
}
|
||||||
|
String wellPopulationsString = populationsStringBuilder.toString();
|
||||||
|
ga.put("well populations", DefaultAttribute.createAttribute(wellPopulationsString));
|
||||||
|
return ga;
|
||||||
|
}
|
||||||
|
|
||||||
public void writeGraphToFile() {
|
public void writeGraphToFile() {
|
||||||
try(BufferedWriter writer = Files.newBufferedWriter(Path.of(filename), StandardOpenOption.CREATE_NEW);
|
try(BufferedWriter writer = Files.newBufferedWriter(Path.of(filename), StandardOpenOption.CREATE_NEW);
|
||||||
){
|
){
|
||||||
GraphMLExporter<SimpleWeightedGraph, BufferedWriter> exporter = new GraphMLExporter<>();
|
//create exporter. Let the vertex labels be the unique ids for the vertices
|
||||||
|
GraphMLExporter<Vertex, SimpleWeightedGraph<Vertex, DefaultWeightedEdge>> exporter = new GraphMLExporter<>(v -> v.getVertexLabel().toString());
|
||||||
|
//set to export weights
|
||||||
|
exporter.setExportEdgeWeights(true);
|
||||||
|
//Set graph attributes
|
||||||
|
exporter.setGraphAttributeProvider( () -> graphAttributes);
|
||||||
|
//set type, sequence, and occupancy attributes for each vertex
|
||||||
|
exporter.setVertexAttributeProvider( v -> {
|
||||||
|
Map<String, Attribute> attributes = new HashMap<>();
|
||||||
|
attributes.put("type", DefaultAttribute.createAttribute(v.getType().name()));
|
||||||
|
attributes.put("sequence", DefaultAttribute.createAttribute(v.getSequence()));
|
||||||
|
attributes.put("occupancy", DefaultAttribute.createAttribute(v.getOccupancy()));
|
||||||
|
return attributes;
|
||||||
|
});
|
||||||
|
//register the attributes
|
||||||
|
for(String s : graphAttributes.keySet()) {
|
||||||
|
exporter.registerAttribute(s, AttributeCategory.GRAPH, AttributeType.STRING);
|
||||||
|
}
|
||||||
|
exporter.registerAttribute("type", AttributeCategory.NODE, AttributeType.STRING);
|
||||||
|
exporter.registerAttribute("sequence", AttributeCategory.NODE, AttributeType.STRING);
|
||||||
|
exporter.registerAttribute("occupancy", AttributeCategory.NODE, AttributeType.STRING);
|
||||||
|
//export the graph
|
||||||
exporter.exportGraph(graph, writer);
|
exporter.exportGraph(graph, writer);
|
||||||
} catch(IOException ex){
|
} catch(IOException ex){
|
||||||
System.out.println("Could not make new file named "+filename);
|
System.out.println("Could not make new file named "+filename);
|
||||||
|
|||||||
@@ -2,88 +2,103 @@ import org.jgrapht.graph.DefaultWeightedEdge;
|
|||||||
import org.jgrapht.graph.SimpleWeightedGraph;
|
import org.jgrapht.graph.SimpleWeightedGraph;
|
||||||
|
|
||||||
import java.util.ArrayList;
|
import java.util.ArrayList;
|
||||||
|
import java.util.HashMap;
|
||||||
import java.util.List;
|
import java.util.List;
|
||||||
import java.util.Map;
|
import java.util.Map;
|
||||||
import java.util.Set;
|
|
||||||
|
|
||||||
public abstract class GraphModificationFunctions {
|
public interface GraphModificationFunctions {
|
||||||
|
|
||||||
//remove over- and under-weight edges
|
//remove over- and under-weight edges
|
||||||
public static List<Integer[]> filterByOverlapThresholds(SimpleWeightedGraph<Integer, DefaultWeightedEdge> graph,
|
static Map<Vertex[], Integer> filterByOverlapThresholds(SimpleWeightedGraph<Vertex, DefaultWeightedEdge> graph,
|
||||||
int low, int high) {
|
int low, int high, boolean saveEdges) {
|
||||||
List<Integer[]> removedEdges = new ArrayList<>();
|
Map<Vertex[], Integer> removedEdges = new HashMap<>();
|
||||||
for(DefaultWeightedEdge e: graph.edgeSet()){
|
//List<Integer[]> removedEdges = new ArrayList<>();
|
||||||
if ((graph.getEdgeWeight(e) > high) || (graph.getEdgeWeight(e) < low)){
|
for (DefaultWeightedEdge e : graph.edgeSet()) {
|
||||||
Integer source = graph.getEdgeSource(e);
|
if ((graph.getEdgeWeight(e) > high) || (graph.getEdgeWeight(e) < low)) {
|
||||||
Integer target = graph.getEdgeTarget(e);
|
if(saveEdges) {
|
||||||
Integer weight = (int) graph.getEdgeWeight(e);
|
Vertex source = graph.getEdgeSource(e);
|
||||||
Integer[] edge = {source, target, weight};
|
Vertex target = graph.getEdgeTarget(e);
|
||||||
removedEdges.add(edge);
|
Integer weight = (int) graph.getEdgeWeight(e);
|
||||||
|
Vertex[] edge = {source, target};
|
||||||
|
removedEdges.put(edge, weight);
|
||||||
|
}
|
||||||
|
else {
|
||||||
|
graph.setEdgeWeight(e, 0.0);
|
||||||
|
}
|
||||||
}
|
}
|
||||||
}
|
}
|
||||||
for (Integer[] edge : removedEdges) {
|
if(saveEdges) {
|
||||||
graph.removeEdge(edge[0], edge[1]);
|
for (Vertex[] edge : removedEdges.keySet()) {
|
||||||
|
graph.removeEdge(edge[0], edge[1]);
|
||||||
|
}
|
||||||
}
|
}
|
||||||
return removedEdges;
|
return removedEdges;
|
||||||
}
|
}
|
||||||
|
|
||||||
//Remove edges for pairs with large occupancy discrepancy
|
//Remove edges for pairs with large occupancy discrepancy
|
||||||
public static List<Integer[]> filterByRelativeOccupancy(SimpleWeightedGraph<Integer, DefaultWeightedEdge> graph,
|
static Map<Vertex[], Integer> filterByRelativeOccupancy(SimpleWeightedGraph<Vertex, DefaultWeightedEdge> graph,
|
||||||
Map<Integer, Integer> alphaWellCounts,
|
Integer maxOccupancyDifference, boolean saveEdges) {
|
||||||
Map<Integer, Integer> betaWellCounts,
|
Map<Vertex[], Integer> removedEdges = new HashMap<>();
|
||||||
Map<Integer, Integer> plateVtoAMap,
|
|
||||||
Map<Integer, Integer> plateVtoBMap,
|
|
||||||
Integer maxOccupancyDifference) {
|
|
||||||
List<Integer[]> removedEdges = new ArrayList<>();
|
|
||||||
for (DefaultWeightedEdge e : graph.edgeSet()) {
|
for (DefaultWeightedEdge e : graph.edgeSet()) {
|
||||||
Integer alphaOcc = alphaWellCounts.get(plateVtoAMap.get(graph.getEdgeSource(e)));
|
Integer alphaOcc = graph.getEdgeSource(e).getOccupancy();
|
||||||
Integer betaOcc = betaWellCounts.get(plateVtoBMap.get(graph.getEdgeTarget(e)));
|
Integer betaOcc = graph.getEdgeTarget(e).getOccupancy();
|
||||||
if (Math.abs(alphaOcc - betaOcc) >= maxOccupancyDifference) {
|
if (Math.abs(alphaOcc - betaOcc) >= maxOccupancyDifference) {
|
||||||
Integer source = graph.getEdgeSource(e);
|
if (saveEdges) {
|
||||||
Integer target = graph.getEdgeTarget(e);
|
Vertex source = graph.getEdgeSource(e);
|
||||||
Integer weight = (int) graph.getEdgeWeight(e);
|
Vertex target = graph.getEdgeTarget(e);
|
||||||
Integer[] edge = {source, target, weight};
|
Integer weight = (int) graph.getEdgeWeight(e);
|
||||||
removedEdges.add(edge);
|
Vertex[] edge = {source, target};
|
||||||
|
removedEdges.put(edge, weight);
|
||||||
|
}
|
||||||
|
else {
|
||||||
|
graph.setEdgeWeight(e, 0.0);
|
||||||
|
}
|
||||||
}
|
}
|
||||||
}
|
}
|
||||||
for (Integer[] edge : removedEdges) {
|
if(saveEdges) {
|
||||||
graph.removeEdge(edge[0], edge[1]);
|
for (Vertex[] edge : removedEdges.keySet()) {
|
||||||
|
graph.removeEdge(edge[0], edge[1]);
|
||||||
|
}
|
||||||
}
|
}
|
||||||
return removedEdges;
|
return removedEdges;
|
||||||
}
|
}
|
||||||
|
|
||||||
//Remove edges for pairs where overlap size is significantly lower than the well occupancy
|
//Remove edges for pairs where overlap size is significantly lower than the well occupancy
|
||||||
public static List<Integer[]> filterByOverlapPercent(SimpleWeightedGraph<Integer, DefaultWeightedEdge> graph,
|
static Map<Vertex[], Integer> filterByOverlapPercent(SimpleWeightedGraph<Vertex, DefaultWeightedEdge> graph,
|
||||||
Map<Integer, Integer> alphaWellCounts,
|
Integer minOverlapPercent,
|
||||||
Map<Integer, Integer> betaWellCounts,
|
boolean saveEdges) {
|
||||||
Map<Integer, Integer> plateVtoAMap,
|
Map<Vertex[], Integer> removedEdges = new HashMap<>();
|
||||||
Map<Integer, Integer> plateVtoBMap,
|
|
||||||
Integer minOverlapPercent) {
|
|
||||||
List<Integer[]> removedEdges = new ArrayList<>();
|
|
||||||
for (DefaultWeightedEdge e : graph.edgeSet()) {
|
for (DefaultWeightedEdge e : graph.edgeSet()) {
|
||||||
Integer alphaOcc = alphaWellCounts.get(plateVtoAMap.get(graph.getEdgeSource(e)));
|
Integer alphaOcc = graph.getEdgeSource(e).getOccupancy();
|
||||||
Integer betaOcc = betaWellCounts.get(plateVtoBMap.get(graph.getEdgeTarget(e)));
|
Integer betaOcc = graph.getEdgeTarget(e).getOccupancy();
|
||||||
double weight = graph.getEdgeWeight(e);
|
double weight = graph.getEdgeWeight(e);
|
||||||
double min = minOverlapPercent / 100.0;
|
double min = minOverlapPercent / 100.0;
|
||||||
if ((weight / alphaOcc < min) || (weight / betaOcc < min)) {
|
if ((weight / alphaOcc < min) || (weight / betaOcc < min)) {
|
||||||
Integer source = graph.getEdgeSource(e);
|
if (saveEdges) {
|
||||||
Integer target = graph.getEdgeTarget(e);
|
Vertex source = graph.getEdgeSource(e);
|
||||||
Integer intWeight = (int) graph.getEdgeWeight(e);
|
Vertex target = graph.getEdgeTarget(e);
|
||||||
Integer[] edge = {source, target, intWeight};
|
Integer intWeight = (int) graph.getEdgeWeight(e);
|
||||||
removedEdges.add(edge);
|
Vertex[] edge = {source, target};
|
||||||
|
removedEdges.put(edge, intWeight);
|
||||||
|
}
|
||||||
|
else {
|
||||||
|
graph.setEdgeWeight(e, 0.0);
|
||||||
|
}
|
||||||
}
|
}
|
||||||
}
|
}
|
||||||
for (Integer[] edge : removedEdges) {
|
if(saveEdges) {
|
||||||
graph.removeEdge(edge[0], edge[1]);
|
for (Vertex[] edge : removedEdges.keySet()) {
|
||||||
|
graph.removeEdge(edge[0], edge[1]);
|
||||||
|
}
|
||||||
}
|
}
|
||||||
return removedEdges;
|
return removedEdges;
|
||||||
}
|
}
|
||||||
|
|
||||||
public static void addRemovedEdges(SimpleWeightedGraph<Integer, DefaultWeightedEdge> graph,
|
static void addRemovedEdges(SimpleWeightedGraph<Vertex, DefaultWeightedEdge> graph,
|
||||||
List<Integer[]> removedEdges) {
|
Map<Vertex[], Integer> removedEdges) {
|
||||||
for (Integer[] edge : removedEdges) {
|
for (Vertex[] edge : removedEdges.keySet()) {
|
||||||
DefaultWeightedEdge e = graph.addEdge(edge[0], edge[1]);
|
DefaultWeightedEdge e = graph.addEdge(edge[0], edge[1]);
|
||||||
graph.setEdgeWeight(e, (double) edge[2]);
|
graph.setEdgeWeight(e, removedEdges.get(edge));
|
||||||
}
|
}
|
||||||
}
|
}
|
||||||
|
|
||||||
|
|||||||
@@ -6,6 +6,7 @@ import java.util.Map;
|
|||||||
//Can't just write the graph, because I need the occupancy data too.
|
//Can't just write the graph, because I need the occupancy data too.
|
||||||
//Makes most sense to serialize object and write that to a file.
|
//Makes most sense to serialize object and write that to a file.
|
||||||
//Which means there's no reason to split map data and graph data up.
|
//Which means there's no reason to split map data and graph data up.
|
||||||
|
//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 {
|
public class GraphWithMapData implements java.io.Serializable {
|
||||||
|
|
||||||
private String sourceFilename;
|
private String sourceFilename;
|
||||||
@@ -15,32 +16,33 @@ public class GraphWithMapData implements java.io.Serializable {
|
|||||||
private Integer alphaCount;
|
private Integer alphaCount;
|
||||||
private Integer betaCount;
|
private Integer betaCount;
|
||||||
private final Map<Integer, Integer> distCellsMapAlphaKey;
|
private final Map<Integer, Integer> distCellsMapAlphaKey;
|
||||||
private final Map<Integer, Integer> plateVtoAMap;
|
// private final Map<Integer, Integer> plateVtoAMap;
|
||||||
private final Map<Integer, Integer> plateVtoBMap;
|
// private final Map<Integer, Integer> plateVtoBMap;
|
||||||
private final Map<Integer, Integer> plateAtoVMap;
|
// private final Map<Integer, Integer> plateAtoVMap;
|
||||||
private final Map<Integer, Integer> plateBtoVMap;
|
// private final Map<Integer, Integer> plateBtoVMap;
|
||||||
private final Map<Integer, Integer> alphaWellCounts;
|
// private final Map<Integer, Integer> alphaWellCounts;
|
||||||
private final Map<Integer, Integer> betaWellCounts;
|
// private final Map<Integer, Integer> betaWellCounts;
|
||||||
private final Duration time;
|
private final Duration time;
|
||||||
|
|
||||||
public GraphWithMapData(SimpleWeightedGraph graph, Integer numWells, Integer[] wellConcentrations,
|
public GraphWithMapData(SimpleWeightedGraph graph, Integer numWells, Integer[] wellConcentrations,
|
||||||
Integer alphaCount, Integer betaCount,
|
Map<Integer, Integer> distCellsMapAlphaKey, Duration time){
|
||||||
Map<Integer, Integer> distCellsMapAlphaKey, Map<Integer, Integer> plateVtoAMap,
|
|
||||||
Map<Integer,Integer> plateVtoBMap, Map<Integer, Integer> plateAtoVMap,
|
// Map<Integer, Integer> plateVtoAMap, Integer alphaCount, Integer betaCount,
|
||||||
Map<Integer, Integer> plateBtoVMap, Map<Integer, Integer> alphaWellCounts,
|
// Map<Integer,Integer> plateVtoBMap, Map<Integer, Integer> plateAtoVMap,
|
||||||
Map<Integer, Integer> betaWellCounts, Duration time) {
|
// Map<Integer, Integer> plateBtoVMap, Map<Integer, Integer> alphaWellCounts,
|
||||||
|
// Map<Integer, Integer> betaWellCounts,) {
|
||||||
this.graph = graph;
|
this.graph = graph;
|
||||||
this.numWells = numWells;
|
this.numWells = numWells;
|
||||||
this.wellPopulations = wellConcentrations;
|
this.wellPopulations = wellConcentrations;
|
||||||
this.alphaCount = alphaCount;
|
this.alphaCount = alphaCount;
|
||||||
this.betaCount = betaCount;
|
this.betaCount = betaCount;
|
||||||
this.distCellsMapAlphaKey = distCellsMapAlphaKey;
|
this.distCellsMapAlphaKey = distCellsMapAlphaKey;
|
||||||
this.plateVtoAMap = plateVtoAMap;
|
// this.plateVtoAMap = plateVtoAMap;
|
||||||
this.plateVtoBMap = plateVtoBMap;
|
// this.plateVtoBMap = plateVtoBMap;
|
||||||
this.plateAtoVMap = plateAtoVMap;
|
// this.plateAtoVMap = plateAtoVMap;
|
||||||
this.plateBtoVMap = plateBtoVMap;
|
// this.plateBtoVMap = plateBtoVMap;
|
||||||
this.alphaWellCounts = alphaWellCounts;
|
// this.alphaWellCounts = alphaWellCounts;
|
||||||
this.betaWellCounts = betaWellCounts;
|
// this.betaWellCounts = betaWellCounts;
|
||||||
this.time = time;
|
this.time = time;
|
||||||
}
|
}
|
||||||
|
|
||||||
@@ -56,41 +58,41 @@ public class GraphWithMapData implements java.io.Serializable {
|
|||||||
return wellPopulations;
|
return wellPopulations;
|
||||||
}
|
}
|
||||||
|
|
||||||
public Integer getAlphaCount() {
|
// public Integer getAlphaCount() {
|
||||||
return alphaCount;
|
// return alphaCount;
|
||||||
}
|
// }
|
||||||
|
//
|
||||||
public Integer getBetaCount() {
|
// public Integer getBetaCount() {
|
||||||
return betaCount;
|
// return betaCount;
|
||||||
}
|
// }
|
||||||
|
|
||||||
public Map<Integer, Integer> getDistCellsMapAlphaKey() {
|
public Map<Integer, Integer> getDistCellsMapAlphaKey() {
|
||||||
return distCellsMapAlphaKey;
|
return distCellsMapAlphaKey;
|
||||||
}
|
}
|
||||||
|
|
||||||
public Map<Integer, Integer> getPlateVtoAMap() {
|
// public Map<Integer, Integer> getPlateVtoAMap() {
|
||||||
return plateVtoAMap;
|
// return plateVtoAMap;
|
||||||
}
|
// }
|
||||||
|
//
|
||||||
public Map<Integer, Integer> getPlateVtoBMap() {
|
// public Map<Integer, Integer> getPlateVtoBMap() {
|
||||||
return plateVtoBMap;
|
// return plateVtoBMap;
|
||||||
}
|
// }
|
||||||
|
//
|
||||||
public Map<Integer, Integer> getPlateAtoVMap() {
|
// public Map<Integer, Integer> getPlateAtoVMap() {
|
||||||
return plateAtoVMap;
|
// return plateAtoVMap;
|
||||||
}
|
// }
|
||||||
|
//
|
||||||
public Map<Integer, Integer> getPlateBtoVMap() {
|
// public Map<Integer, Integer> getPlateBtoVMap() {
|
||||||
return plateBtoVMap;
|
// return plateBtoVMap;
|
||||||
}
|
// }
|
||||||
|
//
|
||||||
public Map<Integer, Integer> getAlphaWellCounts() {
|
// public Map<Integer, Integer> getAlphaWellCounts() {
|
||||||
return alphaWellCounts;
|
// return alphaWellCounts;
|
||||||
}
|
// }
|
||||||
|
//
|
||||||
public Map<Integer, Integer> getBetaWellCounts() {
|
// public Map<Integer, Integer> getBetaWellCounts() {
|
||||||
return betaWellCounts;
|
// return betaWellCounts;
|
||||||
}
|
// }
|
||||||
|
|
||||||
public Duration getTime() {
|
public Duration getTime() {
|
||||||
return time;
|
return time;
|
||||||
|
|||||||
@@ -38,10 +38,10 @@ public class InteractiveInterface {
|
|||||||
case 3 -> makeCDR3Graph();
|
case 3 -> makeCDR3Graph();
|
||||||
case 4 -> matchCDR3s();
|
case 4 -> matchCDR3s();
|
||||||
//case 6 -> matchCellsCDR1();
|
//case 6 -> matchCellsCDR1();
|
||||||
case 8 -> options();
|
case 8 -> mainOptions();
|
||||||
case 9 -> acknowledge();
|
case 9 -> acknowledge();
|
||||||
case 0 -> quit = true;
|
case 0 -> quit = true;
|
||||||
default -> throw new InputMismatchException("Invalid input.");
|
default -> System.out.println("Invalid input.");
|
||||||
}
|
}
|
||||||
} catch (InputMismatchException | IOException ex) {
|
} catch (InputMismatchException | IOException ex) {
|
||||||
System.out.println(ex);
|
System.out.println(ex);
|
||||||
@@ -74,7 +74,7 @@ public class InteractiveInterface {
|
|||||||
System.out.println(ex);
|
System.out.println(ex);
|
||||||
sc.next();
|
sc.next();
|
||||||
}
|
}
|
||||||
CellSample sample = Simulator.generateCellSample(numCells, cdr1Freq);
|
CellSample sample = new CellSample(numCells, cdr1Freq);
|
||||||
assert filename != null;
|
assert filename != null;
|
||||||
System.out.println("Writing cells to file");
|
System.out.println("Writing cells to file");
|
||||||
CellFileWriter writer = new CellFileWriter(filename, sample);
|
CellFileWriter writer = new CellFileWriter(filename, sample);
|
||||||
@@ -227,16 +227,14 @@ public class InteractiveInterface {
|
|||||||
Plate samplePlate;
|
Plate samplePlate;
|
||||||
PlateFileWriter writer;
|
PlateFileWriter writer;
|
||||||
if(exponential){
|
if(exponential){
|
||||||
samplePlate = new Plate(numWells, dropOutRate, populations);
|
samplePlate = new Plate(cells, cellFile, numWells, populations, dropOutRate, lambda, true);
|
||||||
samplePlate.fillWellsExponential(cellFile, cells.getCells(), lambda);
|
|
||||||
writer = new PlateFileWriter(filename, samplePlate);
|
writer = new PlateFileWriter(filename, samplePlate);
|
||||||
}
|
}
|
||||||
else {
|
else {
|
||||||
if (poisson) {
|
if (poisson) {
|
||||||
stdDev = Math.sqrt(cells.getCellCount()); //gaussian with square root of elements approximates poisson
|
stdDev = Math.sqrt(cells.getCellCount()); //gaussian with square root of elements approximates poisson
|
||||||
}
|
}
|
||||||
samplePlate = new Plate(numWells, dropOutRate, populations);
|
samplePlate = new Plate(cells, cellFile, numWells, populations, dropOutRate, stdDev, false);
|
||||||
samplePlate.fillWells(cellFile, cells.getCells(), stdDev);
|
|
||||||
writer = new PlateFileWriter(filename, samplePlate);
|
writer = new PlateFileWriter(filename, samplePlate);
|
||||||
}
|
}
|
||||||
System.out.println("Writing Sample Plate to file");
|
System.out.println("Writing Sample Plate to file");
|
||||||
@@ -252,7 +250,6 @@ public class InteractiveInterface {
|
|||||||
String filename = null;
|
String filename = null;
|
||||||
String cellFile = null;
|
String cellFile = null;
|
||||||
String plateFile = null;
|
String plateFile = null;
|
||||||
|
|
||||||
try {
|
try {
|
||||||
String str = "\nGenerating bipartite weighted graph encoding occupancy overlap data ";
|
String str = "\nGenerating bipartite weighted graph encoding occupancy overlap data ";
|
||||||
str = str.concat("\nrequires a cell sample file and a sample plate file.");
|
str = str.concat("\nrequires a cell sample file and a sample plate file.");
|
||||||
@@ -261,7 +258,7 @@ public class InteractiveInterface {
|
|||||||
cellFile = sc.next();
|
cellFile = sc.next();
|
||||||
System.out.print("\nPlease enter name of an existing sample plate file: ");
|
System.out.print("\nPlease enter name of an existing sample plate file: ");
|
||||||
plateFile = sc.next();
|
plateFile = sc.next();
|
||||||
System.out.println("\nThe graph and occupancy data will be written to a serialized binary file.");
|
System.out.println("\nThe graph and occupancy data will be written to a file.");
|
||||||
System.out.print("Please enter a name for the output file: ");
|
System.out.print("Please enter a name for the output file: ");
|
||||||
filename = sc.next();
|
filename = sc.next();
|
||||||
} catch (InputMismatchException ex) {
|
} catch (InputMismatchException ex) {
|
||||||
@@ -293,7 +290,7 @@ public class InteractiveInterface {
|
|||||||
else {
|
else {
|
||||||
System.out.println("Reading Sample Plate file: " + plateFile);
|
System.out.println("Reading Sample Plate file: " + plateFile);
|
||||||
PlateFileReader plateReader = new PlateFileReader(plateFile);
|
PlateFileReader plateReader = new PlateFileReader(plateFile);
|
||||||
plate = new Plate(plateReader.getFilename(), plateReader.getWells());
|
plate = plateReader.getSamplePlate();
|
||||||
if(BiGpairSEQ.cachePlate()) {
|
if(BiGpairSEQ.cachePlate()) {
|
||||||
BiGpairSEQ.setPlateInMemory(plate, plateFile);
|
BiGpairSEQ.setPlateInMemory(plate, plateFile);
|
||||||
}
|
}
|
||||||
@@ -307,12 +304,18 @@ public class InteractiveInterface {
|
|||||||
System.out.println("Returning to main menu.");
|
System.out.println("Returning to main menu.");
|
||||||
}
|
}
|
||||||
else{
|
else{
|
||||||
List<Integer[]> cells = cellSample.getCells();
|
GraphWithMapData data = Simulator.makeGraph(cellSample, plate, true);
|
||||||
GraphWithMapData data = Simulator.makeGraph(cells, plate, true);
|
|
||||||
assert filename != null;
|
assert filename != null;
|
||||||
GraphDataObjectWriter dataWriter = new GraphDataObjectWriter(filename, data);
|
if(BiGpairSEQ.outputBinary()) {
|
||||||
dataWriter.writeDataToFile();
|
GraphDataObjectWriter dataWriter = new GraphDataObjectWriter(filename, data);
|
||||||
System.out.println("Graph and Data file written to: " + filename);
|
dataWriter.writeDataToFile();
|
||||||
|
System.out.println("Serialized binary graph/data file written to: " + filename);
|
||||||
|
}
|
||||||
|
if(BiGpairSEQ.outputGraphML()) {
|
||||||
|
GraphMLFileWriter graphMLWriter = new GraphMLFileWriter(filename, data);
|
||||||
|
graphMLWriter.writeGraphToFile();
|
||||||
|
System.out.println("GraphML file written to: " + filename);
|
||||||
|
}
|
||||||
if(BiGpairSEQ.cacheGraph()) {
|
if(BiGpairSEQ.cacheGraph()) {
|
||||||
BiGpairSEQ.setGraphInMemory(data, filename);
|
BiGpairSEQ.setGraphInMemory(data, filename);
|
||||||
|
|
||||||
@@ -372,7 +375,7 @@ public class InteractiveInterface {
|
|||||||
data = BiGpairSEQ.getGraphInMemory();
|
data = BiGpairSEQ.getGraphInMemory();
|
||||||
}
|
}
|
||||||
else {
|
else {
|
||||||
GraphDataObjectReader dataReader = new GraphDataObjectReader(graphFilename);
|
GraphDataObjectReader dataReader = new GraphDataObjectReader(graphFilename, true);
|
||||||
data = dataReader.getData();
|
data = dataReader.getData();
|
||||||
if(BiGpairSEQ.cacheGraph()) {
|
if(BiGpairSEQ.cacheGraph()) {
|
||||||
BiGpairSEQ.setGraphInMemory(data, graphFilename);
|
BiGpairSEQ.setGraphInMemory(data, graphFilename);
|
||||||
@@ -493,13 +496,16 @@ public class InteractiveInterface {
|
|||||||
// }
|
// }
|
||||||
// }
|
// }
|
||||||
|
|
||||||
private static void options(){
|
private static void mainOptions(){
|
||||||
boolean backToMain = false;
|
boolean backToMain = false;
|
||||||
while(!backToMain) {
|
while(!backToMain) {
|
||||||
System.out.println("\n--------------OPTIONS---------------");
|
System.out.println("\n--------------OPTIONS---------------");
|
||||||
System.out.println("1) Turn " + getOnOff(!BiGpairSEQ.cacheCells()) + " cell sample file caching");
|
System.out.println("1) Turn " + getOnOff(!BiGpairSEQ.cacheCells()) + " cell sample file caching");
|
||||||
System.out.println("2) Turn " + getOnOff(!BiGpairSEQ.cachePlate()) + " plate file caching");
|
System.out.println("2) Turn " + getOnOff(!BiGpairSEQ.cachePlate()) + " plate file caching");
|
||||||
System.out.println("3) Turn " + getOnOff(!BiGpairSEQ.cacheGraph()) + " graph/data file caching");
|
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("0) Return to main menu");
|
System.out.println("0) Return to main menu");
|
||||||
try {
|
try {
|
||||||
input = sc.nextInt();
|
input = sc.nextInt();
|
||||||
@@ -507,8 +513,11 @@ public class InteractiveInterface {
|
|||||||
case 1 -> BiGpairSEQ.setCacheCells(!BiGpairSEQ.cacheCells());
|
case 1 -> BiGpairSEQ.setCacheCells(!BiGpairSEQ.cacheCells());
|
||||||
case 2 -> BiGpairSEQ.setCachePlate(!BiGpairSEQ.cachePlate());
|
case 2 -> BiGpairSEQ.setCachePlate(!BiGpairSEQ.cachePlate());
|
||||||
case 3 -> BiGpairSEQ.setCacheGraph(!BiGpairSEQ.cacheGraph());
|
case 3 -> BiGpairSEQ.setCacheGraph(!BiGpairSEQ.cacheGraph());
|
||||||
|
case 4 -> BiGpairSEQ.setOutputBinary(!BiGpairSEQ.outputBinary());
|
||||||
|
case 5 -> BiGpairSEQ.setOutputGraphML(!BiGpairSEQ.outputGraphML());
|
||||||
|
case 6 -> algorithmOptions();
|
||||||
case 0 -> backToMain = true;
|
case 0 -> backToMain = true;
|
||||||
default -> throw new InputMismatchException("Invalid input.");
|
default -> System.out.println("Invalid input");
|
||||||
}
|
}
|
||||||
} catch (InputMismatchException ex) {
|
} catch (InputMismatchException ex) {
|
||||||
System.out.println(ex);
|
System.out.println(ex);
|
||||||
@@ -517,12 +526,52 @@ public class InteractiveInterface {
|
|||||||
}
|
}
|
||||||
}
|
}
|
||||||
|
|
||||||
|
/**
|
||||||
|
* Helper function for printing menu items in mainOptions(). Returns a string based on the value of parameter.
|
||||||
|
*
|
||||||
|
* @param b - a boolean value
|
||||||
|
* @return String "on" if b is true, "off" if b is false
|
||||||
|
*/
|
||||||
private static String getOnOff(boolean b) {
|
private static String getOnOff(boolean b) {
|
||||||
if (b) { return "on";}
|
if (b) { return "on";}
|
||||||
else { return "off"; }
|
else { return "off"; }
|
||||||
}
|
}
|
||||||
|
|
||||||
|
private static void algorithmOptions(){
|
||||||
|
boolean backToOptions = false;
|
||||||
|
while(!backToOptions) {
|
||||||
|
System.out.println("\n---------ALGORITHM OPTIONS----------");
|
||||||
|
System.out.println("1) Use scaling algorithm by Duan and Su.");
|
||||||
|
System.out.println("2) Use LEDA book algorithm with Fibonacci heap priority queue");
|
||||||
|
System.out.println("3) Use LEDA book algorithm with pairing heap priority queue");
|
||||||
|
System.out.println("0) Return to Options menu");
|
||||||
|
try {
|
||||||
|
input = sc.nextInt();
|
||||||
|
switch (input) {
|
||||||
|
case 1 -> System.out.println("This option is not yet implemented. Choose another.");
|
||||||
|
case 2 -> {
|
||||||
|
BiGpairSEQ.setFibonacciHeap();
|
||||||
|
System.out.println("MWM algorithm set to LEDA with Fibonacci heap");
|
||||||
|
backToOptions = true;
|
||||||
|
}
|
||||||
|
case 3 -> {
|
||||||
|
BiGpairSEQ.setPairingHeap();
|
||||||
|
System.out.println("MWM algorithm set to LEDA with pairing heap");
|
||||||
|
backToOptions = true;
|
||||||
|
}
|
||||||
|
case 0 -> backToOptions = true;
|
||||||
|
default -> System.out.println("Invalid input");
|
||||||
|
}
|
||||||
|
} catch (InputMismatchException ex) {
|
||||||
|
System.out.println(ex);
|
||||||
|
sc.next();
|
||||||
|
}
|
||||||
|
}
|
||||||
|
}
|
||||||
|
|
||||||
private static void acknowledge(){
|
private static void acknowledge(){
|
||||||
|
System.out.println("BiGpairSEQ_Sim " + BiGpairSEQ.getVersion());
|
||||||
|
System.out.println();
|
||||||
System.out.println("This program simulates BiGpairSEQ, a graph theory based adaptation");
|
System.out.println("This program simulates BiGpairSEQ, a graph theory based adaptation");
|
||||||
System.out.println("of the pairSEQ algorithm for pairing T cell receptor sequences.");
|
System.out.println("of the pairSEQ algorithm for pairing T cell receptor sequences.");
|
||||||
System.out.println();
|
System.out.println();
|
||||||
|
|||||||
@@ -21,15 +21,15 @@ public class MatchingResult {
|
|||||||
* well populations *
|
* well populations *
|
||||||
* total alphas found *
|
* total alphas found *
|
||||||
* total betas found *
|
* total betas found *
|
||||||
* high overlap threshold
|
* high overlap threshold *
|
||||||
* low overlap threshold
|
* low overlap threshold *
|
||||||
* maximum occupancy difference
|
* maximum occupancy difference *
|
||||||
* minimum overlap percent
|
* minimum overlap percent *
|
||||||
* pairing attempt rate
|
* pairing attempt rate *
|
||||||
* correct pairing count
|
* correct pairing count *
|
||||||
* incorrect pairing count
|
* incorrect pairing count *
|
||||||
* pairing error rate
|
* pairing error rate *
|
||||||
* simulation time
|
* simulation time (seconds)
|
||||||
*/
|
*/
|
||||||
this.metadata = metadata;
|
this.metadata = metadata;
|
||||||
this.comments = new ArrayList<>();
|
this.comments = new ArrayList<>();
|
||||||
@@ -91,6 +91,22 @@ public class MatchingResult {
|
|||||||
return Integer.parseInt(metadata.get("total beta count"));
|
return Integer.parseInt(metadata.get("total beta count"));
|
||||||
}
|
}
|
||||||
|
|
||||||
//put in the rest of these methods following the same pattern
|
public Integer getHighOverlapThreshold() { return Integer.parseInt(metadata.get("high overlap threshold"));}
|
||||||
|
|
||||||
|
public Integer getLowOverlapThreshold() { return Integer.parseInt(metadata.get("low overlap threshold"));}
|
||||||
|
|
||||||
|
public Integer getMaxOccupancyDifference() { return Integer.parseInt(metadata.get("maximum occupancy difference"));}
|
||||||
|
|
||||||
|
public Integer getMinOverlapPercent() { return Integer.parseInt(metadata.get("minimum overlap percent"));}
|
||||||
|
|
||||||
|
public Double getPairingAttemptRate() { return Double.parseDouble(metadata.get("pairing attempt rate"));}
|
||||||
|
|
||||||
|
public Integer getCorrectPairingCount() { return Integer.parseInt(metadata.get("correct pairing count"));}
|
||||||
|
|
||||||
|
public Integer getIncorrectPairingCount() { return Integer.parseInt(metadata.get("incorrect pairing count"));}
|
||||||
|
|
||||||
|
public Double getPairingErrorRate() { return Double.parseDouble(metadata.get("pairing error rate"));}
|
||||||
|
|
||||||
|
public String getSimulationTime() { return metadata.get("simulation time (seconds)"); }
|
||||||
|
|
||||||
}
|
}
|
||||||
|
|||||||
@@ -8,7 +8,9 @@ TODO: Implement discrete frequency distributions using Vose's Alias Method
|
|||||||
import java.util.*;
|
import java.util.*;
|
||||||
|
|
||||||
public class Plate {
|
public class Plate {
|
||||||
|
private CellSample cells;
|
||||||
private String sourceFile;
|
private String sourceFile;
|
||||||
|
private String filename;
|
||||||
private List<List<Integer[]>> wells;
|
private List<List<Integer[]>> wells;
|
||||||
private final Random rand = BiGpairSEQ.getRand();
|
private final Random rand = BiGpairSEQ.getRand();
|
||||||
private int size;
|
private int size;
|
||||||
@@ -18,6 +20,25 @@ public class Plate {
|
|||||||
private double lambda;
|
private double lambda;
|
||||||
boolean exponential = false;
|
boolean exponential = false;
|
||||||
|
|
||||||
|
public Plate(CellSample cells, String cellFilename, int numWells, Integer[] populations,
|
||||||
|
double dropoutRate, double stdDev_or_lambda, boolean exponential){
|
||||||
|
this.cells = cells;
|
||||||
|
this.sourceFile = cellFilename;
|
||||||
|
this.size = numWells;
|
||||||
|
this.wells = new ArrayList<>();
|
||||||
|
this.error = dropoutRate;
|
||||||
|
this.populations = populations;
|
||||||
|
this.exponential = exponential;
|
||||||
|
if (this.exponential) {
|
||||||
|
this.lambda = stdDev_or_lambda;
|
||||||
|
fillWellsExponential(cells.getCells(), this.lambda);
|
||||||
|
}
|
||||||
|
else {
|
||||||
|
this.stdDev = stdDev_or_lambda;
|
||||||
|
fillWells(cells.getCells(), this.stdDev);
|
||||||
|
}
|
||||||
|
}
|
||||||
|
|
||||||
|
|
||||||
public Plate(int size, double error, Integer[] populations) {
|
public Plate(int size, double error, Integer[] populations) {
|
||||||
this.size = size;
|
this.size = size;
|
||||||
@@ -26,8 +47,9 @@ public class Plate {
|
|||||||
wells = new ArrayList<>();
|
wells = new ArrayList<>();
|
||||||
}
|
}
|
||||||
|
|
||||||
public Plate(String sourceFileName, List<List<Integer[]>> wells) {
|
//constructor for returning a Plate from a PlateFileReader
|
||||||
this.sourceFile = sourceFileName;
|
public Plate(String filename, List<List<Integer[]>> wells) {
|
||||||
|
this.filename = filename;
|
||||||
this.wells = wells;
|
this.wells = wells;
|
||||||
this.size = wells.size();
|
this.size = wells.size();
|
||||||
|
|
||||||
@@ -43,10 +65,9 @@ public class Plate {
|
|||||||
}
|
}
|
||||||
}
|
}
|
||||||
|
|
||||||
public void fillWellsExponential(String sourceFileName, List<Integer[]> cells, double lambda){
|
private void fillWellsExponential(List<Integer[]> cells, double lambda){
|
||||||
this.lambda = lambda;
|
this.lambda = lambda;
|
||||||
exponential = true;
|
exponential = true;
|
||||||
sourceFile = sourceFileName;
|
|
||||||
int numSections = populations.length;
|
int numSections = populations.length;
|
||||||
int section = 0;
|
int section = 0;
|
||||||
double m;
|
double m;
|
||||||
@@ -74,9 +95,8 @@ public class Plate {
|
|||||||
}
|
}
|
||||||
}
|
}
|
||||||
|
|
||||||
public void fillWells(String sourceFileName, List<Integer[]> cells, double stdDev) {
|
private void fillWells( List<Integer[]> cells, double stdDev) {
|
||||||
this.stdDev = stdDev;
|
this.stdDev = stdDev;
|
||||||
sourceFile = sourceFileName;
|
|
||||||
int numSections = populations.length;
|
int numSections = populations.length;
|
||||||
int section = 0;
|
int section = 0;
|
||||||
double m;
|
double m;
|
||||||
@@ -159,4 +179,6 @@ public class Plate {
|
|||||||
public String getSourceFileName() {
|
public String getSourceFileName() {
|
||||||
return sourceFile;
|
return sourceFile;
|
||||||
}
|
}
|
||||||
|
|
||||||
|
public String getFilename() { return filename; }
|
||||||
}
|
}
|
||||||
|
|||||||
@@ -56,11 +56,8 @@ public class PlateFileReader {
|
|||||||
|
|
||||||
}
|
}
|
||||||
|
|
||||||
public List<List<Integer[]>> getWells() {
|
public Plate getSamplePlate() {
|
||||||
return wells;
|
return new Plate(filename, wells);
|
||||||
}
|
}
|
||||||
|
|
||||||
public String getFilename() {
|
|
||||||
return filename;
|
|
||||||
}
|
|
||||||
}
|
}
|
||||||
8
src/main/java/SequenceType.java
Normal file
8
src/main/java/SequenceType.java
Normal file
@@ -0,0 +1,8 @@
|
|||||||
|
//enum for tagging types of sequences
|
||||||
|
//Listed in order that they appear in a cell array, so ordinal() method will return correct index
|
||||||
|
public enum SequenceType {
|
||||||
|
CDR3_ALPHA,
|
||||||
|
CDR3_BETA,
|
||||||
|
CDR1_ALPHA,
|
||||||
|
CDR1_BETA
|
||||||
|
}
|
||||||
@@ -3,6 +3,7 @@ import org.jgrapht.alg.matching.MaximumWeightBipartiteMatching;
|
|||||||
import org.jgrapht.generate.SimpleWeightedBipartiteGraphMatrixGenerator;
|
import org.jgrapht.generate.SimpleWeightedBipartiteGraphMatrixGenerator;
|
||||||
import org.jgrapht.graph.DefaultWeightedEdge;
|
import org.jgrapht.graph.DefaultWeightedEdge;
|
||||||
import org.jgrapht.graph.SimpleWeightedGraph;
|
import org.jgrapht.graph.SimpleWeightedGraph;
|
||||||
|
import org.jheaps.tree.FibonacciHeap;
|
||||||
import org.jheaps.tree.PairingHeap;
|
import org.jheaps.tree.PairingHeap;
|
||||||
|
|
||||||
import java.math.BigDecimal;
|
import java.math.BigDecimal;
|
||||||
@@ -16,41 +17,18 @@ import java.util.stream.IntStream;
|
|||||||
import static java.lang.Float.*;
|
import static java.lang.Float.*;
|
||||||
|
|
||||||
//NOTE: "sequence" in method and variable names refers to a peptide sequence from a simulated T cell
|
//NOTE: "sequence" in method and variable names refers to a peptide sequence from a simulated T cell
|
||||||
public class Simulator {
|
public class Simulator implements GraphModificationFunctions {
|
||||||
private static final int cdr3AlphaIndex = 0;
|
|
||||||
private static final int cdr3BetaIndex = 1;
|
|
||||||
private static final int cdr1AlphaIndex = 2;
|
|
||||||
private static final int cdr1BetaIndex = 3;
|
|
||||||
|
|
||||||
public static CellSample generateCellSample(Integer numDistinctCells, Integer cdr1Freq) {
|
|
||||||
//In real T cells, CDR1s have about one third the diversity of CDR3s
|
|
||||||
List<Integer> numbersCDR3 = new ArrayList<>();
|
|
||||||
List<Integer> numbersCDR1 = new ArrayList<>();
|
|
||||||
Integer numDistCDR3s = 2 * numDistinctCells + 1;
|
|
||||||
IntStream.range(1, numDistCDR3s + 1).forEach(i -> numbersCDR3.add(i));
|
|
||||||
IntStream.range(numDistCDR3s + 1, numDistCDR3s + 1 + (numDistCDR3s / cdr1Freq) + 1).forEach(i -> numbersCDR1.add(i));
|
|
||||||
Collections.shuffle(numbersCDR3);
|
|
||||||
Collections.shuffle(numbersCDR1);
|
|
||||||
|
|
||||||
//Each cell represented by 4 values
|
//Make the graph needed for matching sequences.
|
||||||
//two CDR3s, and two CDR1s. First two values are CDR3s (alpha, beta), second two are CDR1s (alpha, beta)
|
//sourceVertexIndices and targetVertexIndices are indices within the cell to use as for the two sets of vertices
|
||||||
List<Integer[]> distinctCells = new ArrayList<>();
|
//in the bipartite graph. "Source" and "target" are JGraphT terms for the two vertices an edge touches,
|
||||||
for(int i = 0; i < numbersCDR3.size() - 1; i = i + 2){
|
//even if not directed.
|
||||||
Integer tmpCDR3a = numbersCDR3.get(i);
|
public static GraphWithMapData makeGraph(CellSample cellSample, Plate samplePlate, boolean verbose) {
|
||||||
Integer tmpCDR3b = numbersCDR3.get(i+1);
|
|
||||||
Integer tmpCDR1a = numbersCDR1.get(i % numbersCDR1.size());
|
|
||||||
Integer tmpCDR1b = numbersCDR1.get((i+1) % numbersCDR1.size());
|
|
||||||
Integer[] tmp = {tmpCDR3a, tmpCDR3b, tmpCDR1a, tmpCDR1b};
|
|
||||||
distinctCells.add(tmp);
|
|
||||||
}
|
|
||||||
return new CellSample(distinctCells, cdr1Freq);
|
|
||||||
}
|
|
||||||
|
|
||||||
//Make the graph needed for matching CDR3s
|
|
||||||
public static GraphWithMapData makeGraph(List<Integer[]> distinctCells, Plate samplePlate, boolean verbose) {
|
|
||||||
Instant start = Instant.now();
|
Instant start = Instant.now();
|
||||||
int[] alphaIndex = {cdr3AlphaIndex};
|
List<Integer[]> distinctCells = cellSample.getCells();
|
||||||
int[] betaIndex = {cdr3BetaIndex};
|
int[] alphaIndices = {SequenceType.CDR3_ALPHA.ordinal()};
|
||||||
|
int[] betaIndices = {SequenceType.CDR3_BETA.ordinal()};
|
||||||
|
|
||||||
int numWells = samplePlate.getSize();
|
int numWells = samplePlate.getSize();
|
||||||
|
|
||||||
@@ -60,18 +38,20 @@ public class Simulator {
|
|||||||
if(verbose){System.out.println("Cell maps made");}
|
if(verbose){System.out.println("Cell maps made");}
|
||||||
|
|
||||||
if(verbose){System.out.println("Making well maps");}
|
if(verbose){System.out.println("Making well maps");}
|
||||||
Map<Integer, Integer> allAlphas = samplePlate.assayWellsSequenceS(alphaIndex);
|
|
||||||
Map<Integer, Integer> allBetas = samplePlate.assayWellsSequenceS(betaIndex);
|
Map<Integer, Integer> allAlphas = samplePlate.assayWellsSequenceS(alphaIndices);
|
||||||
|
Map<Integer, Integer> allBetas = samplePlate.assayWellsSequenceS(betaIndices);
|
||||||
int alphaCount = allAlphas.size();
|
int alphaCount = allAlphas.size();
|
||||||
if(verbose){System.out.println("All alphas count: " + alphaCount);}
|
if(verbose){System.out.println("All alphas count: " + alphaCount);}
|
||||||
int betaCount = allBetas.size();
|
int betaCount = allBetas.size();
|
||||||
if(verbose){System.out.println("All betas count: " + betaCount);}
|
if(verbose){System.out.println("All betas count: " + betaCount);}
|
||||||
if(verbose){System.out.println("Well maps made");}
|
if(verbose){System.out.println("Well maps made");}
|
||||||
|
|
||||||
if(verbose){System.out.println("Removing sequences present in all wells.");}
|
|
||||||
filterByOccupancyThresholds(allAlphas, 1, numWells - 1);
|
// if(verbose){System.out.println("Removing singleton sequences and sequences present in all wells.");}
|
||||||
filterByOccupancyThresholds(allBetas, 1, numWells - 1);
|
// filterByOccupancyThresholds(allAlphas, 2, numWells - 1);
|
||||||
if(verbose){System.out.println("Sequences removed");}
|
// filterByOccupancyThresholds(allBetas, 2, numWells - 1);
|
||||||
|
// if(verbose){System.out.println("Sequences removed");}
|
||||||
int pairableAlphaCount = allAlphas.size();
|
int pairableAlphaCount = allAlphas.size();
|
||||||
if(verbose){System.out.println("Remaining alphas count: " + pairableAlphaCount);}
|
if(verbose){System.out.println("Remaining alphas count: " + pairableAlphaCount);}
|
||||||
int pairableBetaCount = allBetas.size();
|
int pairableBetaCount = allBetas.size();
|
||||||
@@ -100,29 +80,40 @@ public class Simulator {
|
|||||||
//(technically this is only 1/4 of an adjacency matrix, but that's all you need
|
//(technically this is only 1/4 of an adjacency matrix, but that's all you need
|
||||||
//for a bipartite graph, and all the SimpleWeightedBipartiteGraphMatrixGenerator class expects.)
|
//for a bipartite graph, and all the SimpleWeightedBipartiteGraphMatrixGenerator class expects.)
|
||||||
if(verbose){System.out.println("Creating adjacency matrix");}
|
if(verbose){System.out.println("Creating adjacency matrix");}
|
||||||
//Count how many wells each alpha appears in
|
//Count how many wells each alpha sequence appears in
|
||||||
Map<Integer, Integer> alphaWellCounts = new HashMap<>();
|
Map<Integer, Integer> alphaWellCounts = new HashMap<>();
|
||||||
//count how many wells each beta appears in
|
//count how many wells each beta sequence appears in
|
||||||
Map<Integer, Integer> betaWellCounts = new HashMap<>();
|
Map<Integer, Integer> betaWellCounts = new HashMap<>();
|
||||||
//the adjacency matrix to be used by the graph generator
|
//the adjacency matrix to be used by the graph generator
|
||||||
double[][] weights = new double[plateVtoAMap.size()][plateVtoBMap.size()];
|
double[][] weights = new double[plateVtoAMap.size()][plateVtoBMap.size()];
|
||||||
countSequencesAndFillMatrix(samplePlate, allAlphas, allBetas, plateAtoVMap,
|
countSequencesAndFillMatrix(samplePlate, allAlphas, allBetas, plateAtoVMap,
|
||||||
plateBtoVMap, alphaIndex, betaIndex, alphaWellCounts, betaWellCounts, weights);
|
plateBtoVMap, alphaIndices, betaIndices, alphaWellCounts, betaWellCounts, weights);
|
||||||
if(verbose){System.out.println("Matrix created");}
|
if(verbose){System.out.println("Matrix created");}
|
||||||
|
|
||||||
//create bipartite graph
|
//create bipartite graph
|
||||||
if(verbose){System.out.println("Creating graph");}
|
if(verbose){System.out.println("Creating graph");}
|
||||||
//the graph object
|
//the graph object
|
||||||
SimpleWeightedGraph<Integer, DefaultWeightedEdge> graph =
|
SimpleWeightedGraph<Vertex, DefaultWeightedEdge> graph =
|
||||||
new SimpleWeightedGraph<>(DefaultWeightedEdge.class);
|
new SimpleWeightedGraph<>(DefaultWeightedEdge.class);
|
||||||
//the graph generator
|
//the graph generator
|
||||||
SimpleWeightedBipartiteGraphMatrixGenerator graphGenerator = new SimpleWeightedBipartiteGraphMatrixGenerator();
|
SimpleWeightedBipartiteGraphMatrixGenerator graphGenerator = new SimpleWeightedBipartiteGraphMatrixGenerator();
|
||||||
//the list of alpha vertices
|
//the list of alpha vertices
|
||||||
List<Integer> alphaVertices = new ArrayList<>(plateVtoAMap.keySet()); //This will work because LinkedHashMap preserves order of entry
|
//List<Integer> alphaVertices = new ArrayList<>(plateVtoAMap.keySet()); //This will work because LinkedHashMap preserves order of entry
|
||||||
|
List<Vertex> alphaVertices = new ArrayList<>();
|
||||||
|
//start with map of all alphas mapped to vertex values, get occupancy from the alphaWellCounts map
|
||||||
|
for (Integer seq : plateAtoVMap.keySet()) {
|
||||||
|
Vertex alphaVertex = new Vertex(SequenceType.CDR3_ALPHA, seq, alphaWellCounts.get(seq), plateAtoVMap.get(seq));
|
||||||
|
alphaVertices.add(alphaVertex);
|
||||||
|
}
|
||||||
graphGenerator.first(alphaVertices);
|
graphGenerator.first(alphaVertices);
|
||||||
//the list of beta vertices
|
//the list of beta vertices
|
||||||
List<Integer> betaVertices = new ArrayList<>(plateVtoBMap.keySet());
|
//List<Integer> betaVertices = new ArrayList<>(plateVtoBMap.keySet());//This will work because LinkedHashMap preserves order of entry
|
||||||
graphGenerator.second(betaVertices); //This will work because LinkedHashMap preserves order of entry
|
List<Vertex> betaVertices = new ArrayList<>();
|
||||||
|
for (Integer seq : plateBtoVMap.keySet()) {
|
||||||
|
Vertex betaVertex = new Vertex(SequenceType.CDR3_BETA, seq, betaWellCounts.get(seq), plateBtoVMap.get(seq));
|
||||||
|
betaVertices.add(betaVertex);
|
||||||
|
}
|
||||||
|
graphGenerator.second(betaVertices);
|
||||||
//use adjacency matrix of weight created previously
|
//use adjacency matrix of weight created previously
|
||||||
graphGenerator.weights(weights);
|
graphGenerator.weights(weights);
|
||||||
graphGenerator.generateGraph(graph);
|
graphGenerator.generateGraph(graph);
|
||||||
@@ -132,11 +123,9 @@ public class Simulator {
|
|||||||
Duration time = Duration.between(start, stop);
|
Duration time = Duration.between(start, stop);
|
||||||
|
|
||||||
//create GraphWithMapData object
|
//create GraphWithMapData object
|
||||||
GraphWithMapData output = new GraphWithMapData(graph, numWells, samplePlate.getPopulations(), alphaCount, betaCount,
|
GraphWithMapData output = new GraphWithMapData(graph, numWells, samplePlate.getPopulations(), distCellsMapAlphaKey, time);
|
||||||
distCellsMapAlphaKey, plateVtoAMap, plateVtoBMap, plateAtoVMap,
|
|
||||||
plateBtoVMap, alphaWellCounts, betaWellCounts, time);
|
|
||||||
//Set source file name in graph to name of sample plate
|
//Set source file name in graph to name of sample plate
|
||||||
output.setSourceFilename(samplePlate.getSourceFileName());
|
output.setSourceFilename(samplePlate.getFilename());
|
||||||
//return GraphWithMapData object
|
//return GraphWithMapData object
|
||||||
return output;
|
return output;
|
||||||
}
|
}
|
||||||
@@ -146,47 +135,70 @@ public class Simulator {
|
|||||||
Integer highThreshold, Integer maxOccupancyDifference,
|
Integer highThreshold, Integer maxOccupancyDifference,
|
||||||
Integer minOverlapPercent, boolean verbose) {
|
Integer minOverlapPercent, boolean verbose) {
|
||||||
Instant start = Instant.now();
|
Instant start = Instant.now();
|
||||||
//Integer arrays will contain TO VERTEX, FROM VERTEX, and WEIGHT (which I'll need to cast to double)
|
SimpleWeightedGraph<Vertex, DefaultWeightedEdge> graph = data.getGraph();
|
||||||
List<Integer[]> removedEdges = new ArrayList<>();
|
Map<Vertex[], Integer> removedEdges = new HashMap<>();
|
||||||
|
boolean saveEdges = BiGpairSEQ.cacheGraph();
|
||||||
int numWells = data.getNumWells();
|
int numWells = data.getNumWells();
|
||||||
Integer alphaCount = data.getAlphaCount();
|
//Integer alphaCount = data.getAlphaCount();
|
||||||
Integer betaCount = data.getBetaCount();
|
//Integer betaCount = data.getBetaCount();
|
||||||
Map<Integer, Integer> distCellsMapAlphaKey = data.getDistCellsMapAlphaKey();
|
Map<Integer, Integer> distCellsMapAlphaKey = data.getDistCellsMapAlphaKey();
|
||||||
Map<Integer, Integer> plateVtoAMap = data.getPlateVtoAMap();
|
Set<Vertex> alphas = new HashSet<>();
|
||||||
Map<Integer, Integer> plateVtoBMap = data.getPlateVtoBMap();
|
Set<Vertex> betas = new HashSet<>();
|
||||||
Map<Integer, Integer> alphaWellCounts = data.getAlphaWellCounts();
|
for(Vertex v: graph.vertexSet()) {
|
||||||
Map<Integer, Integer> betaWellCounts = data.getBetaWellCounts();
|
if (SequenceType.CDR3_ALPHA.equals(v.getType())){
|
||||||
SimpleWeightedGraph<Integer, DefaultWeightedEdge> graph = data.getGraph();
|
alphas.add(v);
|
||||||
|
}
|
||||||
|
else {
|
||||||
|
betas.add(v);
|
||||||
|
}
|
||||||
|
}
|
||||||
|
Integer alphaCount = alphas.size();
|
||||||
|
Integer betaCount = betas.size();
|
||||||
|
|
||||||
//remove edges with weights outside given overlap thresholds, add those to removed edge list
|
//remove edges with weights outside given overlap thresholds, add those to removed edge list
|
||||||
if(verbose){System.out.println("Eliminating edges with weights outside overlap threshold values");}
|
if(verbose){System.out.println("Eliminating edges with weights outside overlap threshold values");}
|
||||||
removedEdges.addAll(GraphModificationFunctions.filterByOverlapThresholds(graph, lowThreshold, highThreshold));
|
removedEdges.putAll(GraphModificationFunctions.filterByOverlapThresholds(graph, lowThreshold, highThreshold, saveEdges));
|
||||||
if(verbose){System.out.println("Over- and under-weight edges removed");}
|
if(verbose){System.out.println("Over- and under-weight edges removed");}
|
||||||
|
|
||||||
//remove edges between vertices with too small an overlap size, add those to removed edge list
|
//remove edges between vertices with too small an overlap size, add those to removed edge list
|
||||||
if(verbose){System.out.println("Eliminating edges with weights less than " + minOverlapPercent.toString() +
|
if(verbose){System.out.println("Eliminating edges with weights less than " + minOverlapPercent.toString() +
|
||||||
" percent of vertex occupancy value.");}
|
" percent of vertex occupancy value.");}
|
||||||
removedEdges.addAll(GraphModificationFunctions.filterByOverlapPercent(graph, alphaWellCounts, betaWellCounts,
|
removedEdges.putAll(GraphModificationFunctions.filterByOverlapPercent(graph, minOverlapPercent, saveEdges));
|
||||||
plateVtoAMap, plateVtoBMap, minOverlapPercent));
|
|
||||||
if(verbose){System.out.println("Edges with weights too far below a vertex occupancy value removed");}
|
if(verbose){System.out.println("Edges with weights too far below a vertex occupancy value removed");}
|
||||||
|
|
||||||
//Filter by relative occupancy
|
//Filter by relative occupancy
|
||||||
if(verbose){System.out.println("Eliminating edges between vertices with occupancy difference > "
|
if(verbose){System.out.println("Eliminating edges between vertices with occupancy difference > "
|
||||||
+ maxOccupancyDifference);}
|
+ maxOccupancyDifference);}
|
||||||
removedEdges.addAll(GraphModificationFunctions.filterByRelativeOccupancy(graph, alphaWellCounts, betaWellCounts,
|
removedEdges.putAll(GraphModificationFunctions.filterByRelativeOccupancy(graph, maxOccupancyDifference, saveEdges));
|
||||||
plateVtoAMap, plateVtoBMap, maxOccupancyDifference));
|
|
||||||
if(verbose){System.out.println("Edges between vertices of with excessively different occupancy values " +
|
if(verbose){System.out.println("Edges between vertices of with excessively different occupancy values " +
|
||||||
"removed");}
|
"removed");}
|
||||||
|
|
||||||
//Find Maximum Weighted Matching
|
//Find Maximum Weighted Matching
|
||||||
//using jheaps library class PairingHeap for improved efficiency
|
//using jheaps library class PairingHeap for improved efficiency
|
||||||
if(verbose){System.out.println("Finding maximum weighted matching");}
|
if(verbose){System.out.println("Finding maximum weighted matching");}
|
||||||
//Attempting to use addressable heap to improve performance
|
MaximumWeightBipartiteMatching maxWeightMatching;
|
||||||
MaximumWeightBipartiteMatching maxWeightMatching =
|
//Use correct heap type for priority queue
|
||||||
new MaximumWeightBipartiteMatching(graph,
|
String heapType = BiGpairSEQ.getPriorityQueueHeapType();
|
||||||
plateVtoAMap.keySet(),
|
switch (heapType) {
|
||||||
plateVtoBMap.keySet(),
|
case "PAIRING" -> {
|
||||||
|
maxWeightMatching = new MaximumWeightBipartiteMatching(graph,
|
||||||
|
alphas,
|
||||||
|
betas,
|
||||||
i -> new PairingHeap(Comparator.naturalOrder()));
|
i -> new PairingHeap(Comparator.naturalOrder()));
|
||||||
|
}
|
||||||
|
case "FIBONACCI" -> {
|
||||||
|
maxWeightMatching = new MaximumWeightBipartiteMatching(graph,
|
||||||
|
alphas,
|
||||||
|
betas,
|
||||||
|
i -> new FibonacciHeap(Comparator.naturalOrder()));
|
||||||
|
}
|
||||||
|
default -> {
|
||||||
|
maxWeightMatching = new MaximumWeightBipartiteMatching(graph,
|
||||||
|
alphas,
|
||||||
|
betas);
|
||||||
|
}
|
||||||
|
}
|
||||||
|
//get the matching
|
||||||
MatchingAlgorithm.Matching<String, DefaultWeightedEdge> graphMatching = maxWeightMatching.getMatching();
|
MatchingAlgorithm.Matching<String, DefaultWeightedEdge> graphMatching = maxWeightMatching.getMatching();
|
||||||
if(verbose){System.out.println("Matching completed");}
|
if(verbose){System.out.println("Matching completed");}
|
||||||
Instant stop = Instant.now();
|
Instant stop = Instant.now();
|
||||||
@@ -213,11 +225,14 @@ public class Simulator {
|
|||||||
Map<Integer, Integer> matchMap = new HashMap<>();
|
Map<Integer, Integer> matchMap = new HashMap<>();
|
||||||
while(weightIter.hasNext()) {
|
while(weightIter.hasNext()) {
|
||||||
e = weightIter.next();
|
e = weightIter.next();
|
||||||
Integer source = graph.getEdgeSource(e);
|
Vertex source = graph.getEdgeSource(e);
|
||||||
Integer target = graph.getEdgeTarget(e);
|
Vertex target = graph.getEdgeTarget(e);
|
||||||
|
//Integer source = graph.getEdgeSource(e);
|
||||||
|
//Integer target = graph.getEdgeTarget(e);
|
||||||
//The match map is all matches found, not just true matches!
|
//The match map is all matches found, not just true matches!
|
||||||
matchMap.put(plateVtoAMap.get(source), plateVtoBMap.get(target));
|
matchMap.put(source.getSequence(), target.getSequence());
|
||||||
check = plateVtoBMap.get(target).equals(distCellsMapAlphaKey.get(plateVtoAMap.get(source)));
|
check = target.getSequence().equals(distCellsMapAlphaKey.get(source.getSequence()));
|
||||||
|
//check = plateVtoBMap.get(target).equals(distCellsMapAlphaKey.get(plateVtoAMap.get(source)));
|
||||||
if(check) {
|
if(check) {
|
||||||
trueCount++;
|
trueCount++;
|
||||||
}
|
}
|
||||||
@@ -225,35 +240,40 @@ public class Simulator {
|
|||||||
falseCount++;
|
falseCount++;
|
||||||
}
|
}
|
||||||
List<String> result = new ArrayList<>();
|
List<String> result = new ArrayList<>();
|
||||||
result.add(plateVtoAMap.get(source).toString());
|
//alpha sequence
|
||||||
|
result.add(source.getSequence().toString());
|
||||||
//alpha well count
|
//alpha well count
|
||||||
result.add(alphaWellCounts.get(plateVtoAMap.get(source)).toString());
|
result.add(source.getOccupancy().toString());
|
||||||
result.add(plateVtoBMap.get(target).toString());
|
//beta sequence
|
||||||
|
result.add(target.getSequence().toString());
|
||||||
//beta well count
|
//beta well count
|
||||||
result.add(betaWellCounts.get(plateVtoBMap.get(target)).toString());
|
result.add(target.getOccupancy().toString());
|
||||||
//overlap count
|
//overlap count
|
||||||
result.add(Double.toString(graph.getEdgeWeight(e)));
|
result.add(Double.toString(graph.getEdgeWeight(e)));
|
||||||
result.add(Boolean.toString(check));
|
result.add(Boolean.toString(check));
|
||||||
double pValue = Equations.pValue(numWells, alphaWellCounts.get(plateVtoAMap.get(source)),
|
double pValue = Equations.pValue(numWells, source.getOccupancy(),
|
||||||
betaWellCounts.get(plateVtoBMap.get(target)), graph.getEdgeWeight(e));
|
target.getOccupancy(), graph.getEdgeWeight(e));
|
||||||
BigDecimal pValueTrunc = new BigDecimal(pValue, mc);
|
BigDecimal pValueTrunc = new BigDecimal(pValue, mc);
|
||||||
result.add(pValueTrunc.toString());
|
result.add(pValueTrunc.toString());
|
||||||
allResults.add(result);
|
allResults.add(result);
|
||||||
}
|
}
|
||||||
|
|
||||||
//Metadata comments for CSV file
|
//Metadata comments for CSV file
|
||||||
|
String algoType = "LEDA book with heap: " + heapType;
|
||||||
int min = Math.min(alphaCount, betaCount);
|
int min = Math.min(alphaCount, betaCount);
|
||||||
|
//matching weight
|
||||||
|
BigDecimal totalMatchingWeight = maxWeightMatching.getMatchingWeight();
|
||||||
//rate of attempted matching
|
//rate of attempted matching
|
||||||
double attemptRate = (double) (trueCount + falseCount) / min;
|
double attemptRate = (double) (trueCount + falseCount) / min;
|
||||||
BigDecimal attemptRateTrunc = new BigDecimal(attemptRate, mc);
|
BigDecimal attemptRateTrunc = new BigDecimal(attemptRate, mc);
|
||||||
//rate of pairing error
|
//rate of pairing error
|
||||||
double pairingErrorRate = (double) falseCount / (trueCount + falseCount);
|
double pairingErrorRate = (double) falseCount / (trueCount + falseCount);
|
||||||
BigDecimal pairingErrorRateTrunc;
|
BigDecimal pairingErrorRateTrunc;
|
||||||
if(pairingErrorRate == NaN || pairingErrorRate == POSITIVE_INFINITY || pairingErrorRate == NEGATIVE_INFINITY) {
|
if(Double.isFinite(pairingErrorRate)) {
|
||||||
pairingErrorRateTrunc = new BigDecimal(-1, mc);
|
pairingErrorRateTrunc = new BigDecimal(pairingErrorRate, mc);
|
||||||
}
|
}
|
||||||
else{
|
else{
|
||||||
pairingErrorRateTrunc = new BigDecimal(pairingErrorRate, mc);
|
pairingErrorRateTrunc = new BigDecimal(-1, mc);
|
||||||
}
|
}
|
||||||
//get list of well populations
|
//get list of well populations
|
||||||
Integer[] wellPopulations = data.getWellPopulations();
|
Integer[] wellPopulations = data.getWellPopulations();
|
||||||
@@ -272,6 +292,8 @@ public class Simulator {
|
|||||||
Map<String, String> metadata = new LinkedHashMap<>();
|
Map<String, String> metadata = new LinkedHashMap<>();
|
||||||
metadata.put("sample plate filename", data.getSourceFilename());
|
metadata.put("sample plate filename", data.getSourceFilename());
|
||||||
metadata.put("graph filename", dataFilename);
|
metadata.put("graph filename", dataFilename);
|
||||||
|
metadata.put("algorithm type", algoType);
|
||||||
|
metadata.put("matching weight", totalMatchingWeight.toString());
|
||||||
metadata.put("well populations", wellPopulationsString);
|
metadata.put("well populations", wellPopulationsString);
|
||||||
metadata.put("total alphas found", alphaCount.toString());
|
metadata.put("total alphas found", alphaCount.toString());
|
||||||
metadata.put("total betas found", betaCount.toString());
|
metadata.put("total betas found", betaCount.toString());
|
||||||
@@ -283,7 +305,7 @@ public class Simulator {
|
|||||||
metadata.put("correct pairing count", Integer.toString(trueCount));
|
metadata.put("correct pairing count", Integer.toString(trueCount));
|
||||||
metadata.put("incorrect pairing count", Integer.toString(falseCount));
|
metadata.put("incorrect pairing count", Integer.toString(falseCount));
|
||||||
metadata.put("pairing error rate", pairingErrorRateTrunc.toString());
|
metadata.put("pairing error rate", pairingErrorRateTrunc.toString());
|
||||||
metadata.put("simulation time", nf.format(time.toSeconds()));
|
metadata.put("simulation time (seconds)", nf.format(time.toSeconds()));
|
||||||
//create MatchingResult object
|
//create MatchingResult object
|
||||||
MatchingResult output = new MatchingResult(metadata, header, allResults, matchMap, time);
|
MatchingResult output = new MatchingResult(metadata, header, allResults, matchMap, time);
|
||||||
if(verbose){
|
if(verbose){
|
||||||
@@ -292,10 +314,11 @@ public class Simulator {
|
|||||||
}
|
}
|
||||||
}
|
}
|
||||||
|
|
||||||
//put the removed edges back on the graph
|
if(saveEdges) {
|
||||||
System.out.println("Restoring removed edges to graph.");
|
//put the removed edges back on the graph
|
||||||
GraphModificationFunctions.addRemovedEdges(graph, removedEdges);
|
System.out.println("Restoring removed edges to graph.");
|
||||||
|
GraphModificationFunctions.addRemovedEdges(graph, removedEdges);
|
||||||
|
}
|
||||||
//return MatchingResult object
|
//return MatchingResult object
|
||||||
return output;
|
return output;
|
||||||
}
|
}
|
||||||
@@ -671,7 +694,7 @@ public class Simulator {
|
|||||||
|
|
||||||
private static Map<Integer, Integer> makeVertexToSequenceMap(Map<Integer, Integer> sequences, Integer startValue) {
|
private static Map<Integer, Integer> makeVertexToSequenceMap(Map<Integer, Integer> sequences, Integer startValue) {
|
||||||
Map<Integer, Integer> map = new LinkedHashMap<>(); //LinkedHashMap to preserve order of entry
|
Map<Integer, Integer> map = new LinkedHashMap<>(); //LinkedHashMap to preserve order of entry
|
||||||
Integer index = startValue;
|
Integer index = startValue; //is this necessary? I don't think I use this.
|
||||||
for (Integer k: sequences.keySet()) {
|
for (Integer k: sequences.keySet()) {
|
||||||
map.put(index, k);
|
map.put(index, k);
|
||||||
index++;
|
index++;
|
||||||
|
|||||||
@@ -1,17 +1,92 @@
|
|||||||
public class Vertex {
|
import java.io.Serializable;
|
||||||
private final Integer peptide;
|
|
||||||
private final Integer occupancy;
|
|
||||||
|
|
||||||
public Vertex(Integer peptide, Integer occupancy) {
|
public class Vertex implements Serializable {
|
||||||
this.peptide = peptide;
|
private SequenceType type;
|
||||||
|
private Integer vertexLabel;
|
||||||
|
private Integer sequence;
|
||||||
|
private Integer occupancy;
|
||||||
|
|
||||||
|
public Vertex(Integer vertexLabel) {
|
||||||
|
this.vertexLabel = vertexLabel;
|
||||||
|
}
|
||||||
|
public Vertex(String vertexLabel) {
|
||||||
|
this.vertexLabel = Integer.parseInt((vertexLabel));
|
||||||
|
}
|
||||||
|
|
||||||
|
public Vertex(SequenceType type, Integer sequence, Integer occupancy, Integer vertexLabel) {
|
||||||
|
this.type = type;
|
||||||
|
this.vertexLabel = vertexLabel;
|
||||||
|
this.sequence = sequence;
|
||||||
this.occupancy = occupancy;
|
this.occupancy = occupancy;
|
||||||
}
|
}
|
||||||
|
|
||||||
public Integer getPeptide() {
|
|
||||||
return peptide;
|
public SequenceType getType() {
|
||||||
|
return type;
|
||||||
|
}
|
||||||
|
|
||||||
|
public void setType(String type) {
|
||||||
|
this.type = SequenceType.valueOf(type);
|
||||||
|
}
|
||||||
|
|
||||||
|
public Integer getVertexLabel() {
|
||||||
|
return vertexLabel;
|
||||||
|
}
|
||||||
|
|
||||||
|
public void setVertexLabel(String label) {
|
||||||
|
this.vertexLabel = Integer.parseInt(label);
|
||||||
|
}
|
||||||
|
|
||||||
|
public Integer getSequence() {
|
||||||
|
|
||||||
|
return sequence;
|
||||||
|
}
|
||||||
|
|
||||||
|
public void setSequence(String sequence) {
|
||||||
|
this.sequence = Integer.parseInt(sequence);
|
||||||
}
|
}
|
||||||
|
|
||||||
public Integer getOccupancy() {
|
public Integer getOccupancy() {
|
||||||
return occupancy;
|
return occupancy;
|
||||||
}
|
}
|
||||||
|
|
||||||
|
public void setOccupancy(String occupancy) {
|
||||||
|
this.occupancy = Integer.parseInt(occupancy);
|
||||||
|
}
|
||||||
|
|
||||||
|
@Override //adapted from JGraphT example code
|
||||||
|
public int hashCode()
|
||||||
|
{
|
||||||
|
return (sequence == null) ? 0 : sequence.hashCode();
|
||||||
|
}
|
||||||
|
|
||||||
|
@Override //adapted from JGraphT example code
|
||||||
|
public boolean equals(Object obj)
|
||||||
|
{
|
||||||
|
if (this == obj)
|
||||||
|
return true;
|
||||||
|
if (obj == null)
|
||||||
|
return false;
|
||||||
|
if (getClass() != obj.getClass())
|
||||||
|
return false;
|
||||||
|
Vertex other = (Vertex) obj;
|
||||||
|
if (sequence == null) {
|
||||||
|
return other.sequence == null;
|
||||||
|
} else {
|
||||||
|
return sequence.equals(other.sequence);
|
||||||
|
}
|
||||||
|
}
|
||||||
|
|
||||||
|
|
||||||
|
@Override //adapted from JGraphT example code
|
||||||
|
public String toString()
|
||||||
|
{
|
||||||
|
StringBuilder sb = new StringBuilder();
|
||||||
|
sb.append("(").append(vertexLabel)
|
||||||
|
.append(", Type: ").append(type.name())
|
||||||
|
.append(", Sequence: ").append(sequence)
|
||||||
|
.append(", Occupancy: ").append(occupancy).append(")");
|
||||||
|
return sb.toString();
|
||||||
|
}
|
||||||
|
|
||||||
}
|
}
|
||||||
|
|||||||
Reference in New Issue
Block a user