82 Commits

Author SHA1 Message Date
eugenefischer
d1a56c3578 Hand-merge of some things from Dev_Vertex branch that didn't make it in for some reason 2022-09-25 19:07:25 -05:00
eugenefischer
16daf02dd6 Merge branch 'Dev_Vertex'
# Conflicts:
#	src/main/java/GraphModificationFunctions.java
#	src/main/java/GraphWithMapData.java
#	src/main/java/Simulator.java
#	src/main/java/Vertex.java
2022-09-25 18:33:26 -05:00
eugenefischer
04a077da2e update Readme 2022-09-25 18:24:12 -05:00
eugenefischer
740835f814 fix typo 2022-09-25 17:47:07 -05:00
eugenefischer
8a77d53f1f Output sequence counts before and after pre-filtering (currently pre-filtering only sequences present in all wells) 2022-09-25 17:20:50 -05:00
eugenefischer
58fa140ee5 add comments 2022-09-25 16:10:17 -05:00
eugenefischer
475bbf3107 Sort vertex lists by vertex label before making adjacency matrix 2022-09-25 15:54:28 -05:00
eugenefischer
4f2fa4cbbe Pre-filter saturating sequences only. Retaining singletons seems to improve matching accuracy in high sample rate test (well populations 10% of total cell sample size) 2022-09-25 15:19:56 -05:00
eugenefischer
58d418e44b Pre-filter saturating sequences only. Retaining singletons seems to improve matching accuracy in high sample rate test (well populations 10% of total cell sample size) 2022-09-25 15:06:46 -05:00
eugenefischer
1971a96467 Remove pre-filtering of singleton and saturating sequences 2022-09-25 14:55:43 -05:00
eugenefischer
e699795521 Revert "by-hand merge of needed code from custom vertex branch"
This reverts commit 29b844afd2.
2022-09-25 14:34:31 -05:00
eugenefischer
bd6d010b0b Revert "update TODO"
This reverts commit a054c0c20a.
2022-09-25 14:34:31 -05:00
eugenefischer
61d1eb3eb1 Revert "Reword output message"
This reverts commit 63317f2aa0.
2022-09-25 14:34:31 -05:00
eugenefischer
cb41b45204 Revert "Reword option menu item"
This reverts commit 06e72314b0.
2022-09-25 14:34:31 -05:00
eugenefischer
a84d2e1bfe Revert "Add comment on map data encodng"
This reverts commit 73c83bf35d.
2022-09-25 14:34:31 -05:00
eugenefischer
7b61d2c0d7 Revert "update version number"
This reverts commit e4e5a1f979.
2022-09-25 14:34:31 -05:00
eugenefischer
56454417c0 Revert "Restore pre-filtering of singleton and saturating sequences"
This reverts commit 5c03909a11.
2022-09-25 14:34:31 -05:00
eugenefischer
8ee1c5903e Merge branch 'master' into Dev_Vertex
# Conflicts:
#	src/main/java/GraphMLFileReader.java
#	src/main/java/InteractiveInterface.java
#	src/main/java/Simulator.java
2022-09-25 14:18:56 -05:00
eugenefischer
5c03909a11 Restore pre-filtering of singleton and saturating sequences 2022-09-22 01:39:13 -05:00
eugenefischer
e4e5a1f979 update version number 2022-09-22 00:00:02 -05:00
eugenefischer
73c83bf35d Add comment on map data encodng 2022-09-21 21:46:00 -05:00
eugenefischer
06e72314b0 Reword option menu item 2022-09-21 21:43:47 -05:00
eugenefischer
63317f2aa0 Reword output message 2022-09-21 18:08:52 -05:00
eugenefischer
a054c0c20a update TODO 2022-09-21 16:50:00 -05:00
eugenefischer
29b844afd2 by-hand merge of needed code from custom vertex branch 2022-09-21 16:48:26 -05:00
eugenefischer
dea4972927 remove prefiltering of singletons and saturating sequences 2022-09-21 16:09:08 -05:00
eugenefischer
9ae38bf247 Fix bug in correct match counter 2022-09-21 15:59:23 -05:00
eugenefischer
3ba305abdb Update ToDo 2022-09-21 13:30:30 -05:00
eugenefischer
3707923398 Merge remote-tracking branch 'origin/master' 2022-09-21 13:16:52 -05:00
eugenefischer
cf771ce574 parameterized sequence indices 2022-09-21 13:15:49 -05:00
f980722b56 update TODO 2022-09-21 18:09:37 +00:00
1df86f01df parameterized sequence indices 2022-03-05 12:03:31 -06:00
96ba57d653 Remove singleton sequences from wells in initial filtering 2022-03-04 16:14:17 -06:00
b602fb02f1 Remove obsolete comments 2022-03-02 23:35:24 -06:00
325e1ebe2b Add data on randomized well population behavior 2022-03-02 23:21:56 -06:00
df047267ee Add data on randomized well population behavior 2022-03-02 22:54:17 -06:00
03e8d31210 Add data on randomized well population behavior 2022-03-02 18:55:19 -06:00
582dc3ef40 Update readme 2022-03-02 12:39:40 -06:00
4c872ed48e Add optional stdout print flags 2022-03-01 15:27:04 -06:00
3fc39302c7 Add detail to error message 2022-03-01 15:24:14 -06:00
578bdc0fbf clarify help menu text 2022-03-01 15:08:43 -06:00
8275cf7740 Check for finite pairing error rate 2022-03-01 09:01:53 -06:00
64209691f0 Check for finite pairing error rate 2022-03-01 09:00:58 -06:00
1886800873 update readme 2022-03-01 08:54:32 -06:00
bedf0894bc update readme 2022-03-01 08:45:40 -06:00
2ac3451842 update readme 2022-03-01 08:43:48 -06:00
67ec3f3764 update readme 2022-03-01 08:43:18 -06:00
b5a8b7e2d5 update readme 2022-03-01 08:41:57 -06:00
9fb3095f0f Clarify help text 2022-03-01 08:40:34 -06:00
25acf920c2 Add version information 2022-03-01 08:34:35 -06:00
f301327693 Update readme with -graphml flag 2022-03-01 08:24:43 -06:00
e04d2d6777 Fix typos in help menu 2022-03-01 08:16:06 -06:00
3e41afaa64 bugfix 2022-02-27 19:08:29 -06:00
bc5d67680d Add flag to print metadata to stdout 2022-02-27 17:36:23 -06:00
f2347e8fc2 check verbose flag 2022-02-27 17:35:50 -06:00
c8364d8a6e check verbose flag 2022-02-27 17:34:20 -06:00
6f5afbc6ec Update readme with CLI arguments 2022-02-27 17:01:12 -06:00
fb4d22e7a4 Update readme with CLI arguments 2022-02-27 17:00:54 -06:00
e10350c214 Update readme with CLI arguments 2022-02-27 16:56:58 -06:00
b1155f8100 Format -help CLI option 2022-02-27 16:53:46 -06:00
12b003a69f Add -help CLI option 2022-02-27 16:45:30 -06:00
32c5bcaaff Deactivate file I/O announcement for CLI 2022-02-27 16:16:24 -06:00
2485ac4cf6 Add getters to MatchingResult 2022-02-27 16:15:26 -06:00
05556bce0c Add units to metadata 2022-02-27 16:08:59 -06:00
a822f69ea4 Control verbose output 2022-02-27 16:07:17 -06:00
3d1f8668ee Control verbose output 2022-02-27 16:03:57 -06:00
40c743308b Initialize wells 2022-02-27 15:54:47 -06:00
5246cc4a0c Re-implement command line options 2022-02-27 15:35:07 -06:00
a5f7c0641d Refactor for better encapsulation with CellSamples 2022-02-27 14:51:53 -06:00
8ebfc1469f Refactor plate to fill its own wells in its constructor 2022-02-27 14:25:53 -06:00
b53f5f1cc0 Refactor plate to fill its own wells in its constructor 2022-02-27 14:17:16 -06:00
974d2d650c Refactor plate to fill its own wells in its constructor 2022-02-27 14:17:11 -06:00
6b5837e6ce Add Vose's alias method to to-dos 2022-02-27 11:46:11 -06:00
817fe51708 Code cleanup 2022-02-26 09:56:46 -06:00
1ea68045ce Refactor cdr3 matching to use new Vertex class 2022-02-26 09:49:16 -06:00
75b2aa9553 testing graph attributes 2022-02-26 08:58:52 -06:00
b3dc10f287 add graph attributes to graphml writer 2022-02-26 08:15:48 -06:00
fb8d8d8785 make heap type an enum 2022-02-26 08:15:31 -06:00
ab437512e9 make Vertex serializable 2022-02-26 07:45:36 -06:00
7b03a3cce8 bugfix 2022-02-26 07:35:34 -06:00
f032d3e852 rewrite GraphML importer/exporter 2022-02-26 07:34:07 -06:00
b604b1d3cd Changing graph to use Vertex class 2022-02-26 06:19:08 -06:00
17 changed files with 929 additions and 501 deletions

103
readme.md
View File

@@ -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.
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
is from Duan and Su (2012). For a graph with m edges, n vertices per side, and maximum integer edge weight N,
The most efficient algorithm known to the author for maximum weight matching of a bipartite graph with strictly integral
weights is from Duan and Su (2012). For a graph with m edges, n vertices per side, and maximum integer edge weight N,
their algorithm runs in **O(m sqrt(n) log(N))** time. As the graph representation of a pairSEQ experiment is
bipartite with integer weights, this algorithm is ideal for BiGpairSEQ.
@@ -43,13 +43,18 @@ Run with the command:
`java -jar BiGpairSEQ_Sim.jar`
Processing sample plates with tens of thousands of sequences may require large amounts
of RAM. It is often desirable to increase the JVM maximum heap allocation with the -Xmx flag.
of RAM. It is often desirable to increase the JVM maximum heap allocation with the `-Xmx` flag.
For example, to run the program with 32 gigabytes of memory, use the command:
`java -Xmx32G -jar BiGpairSEQ_Sim.jar`
Once running, BiGpairSEQ_Sim has an interactive, menu-driven CLI for generating files and simulating TCR pairing. The
main menu looks like this:
There are a number of command line options, to allow the program to be used in shell scripts. For a full list,
use the `-help` flag:
`java -jar BiGpairSEQ_Sim.jar -help`
If no command line arguments are given, BiGpairSEQ_Sim will launch with an interactive, menu-driven CLI for
generating files and simulating TCR pairing. The main menu looks like this:
```
--------BiGPairSEQ SIMULATOR--------
@@ -78,6 +83,7 @@ By default, the Options menu looks like this:
0) Return to main menu
```
### INPUT/OUTPUT
To run the simulation, the program reads and writes 4 kinds of files:
@@ -102,7 +108,7 @@ 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
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.
---
@@ -197,8 +203,13 @@ Options for creating a Graph/Data file:
These files do not have a human-readable structure, and are not portable to other programs.
(For portability to other software, turn on GraphML output in the Options menu. This will produce a .graphml file
for the weighted graph, with vertex attributes sequence, type, and occupancy data.)
*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.
---
@@ -253,29 +264,77 @@ Example output:
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)
### 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),
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
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
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.
As mentioned in the theory section, performance could be improved by implementing a more efficient algorithm for finding
the maximum weight 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
* ~~Try invoking GC at end of workloads to reduce paging to disk~~ DONE
* ~~Hold graph data in memory until another graph is read-in? ABANDONED UNABANDONED~~ DONE
* ~~*No, this won't work, because BiGpairSEQ simulations alter the underlying graph based on filtering constraints. Changes would cascade with multiple experiments.*~~
* Might have figured out a way to do it, by taking edges out and then putting them back into the graph. This may actually be possible.
* It is possible, though the modifications to the graph incur their own performance penalties. Need testing to see which option is best.
* It 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?~~
@@ -288,13 +347,23 @@ slightly less time than the simulation itself. Real elapsed time from start to f
* ~~Apache Commons CSV library writes entries a row at a time~~
* _Got this working, but at the cost of a profoundly strange bug in graph occupancy filtering. Have reverted the repo until I can figure out what caused that. Given how easily Thingiverse transposes CSV matrices in R, might not even be worth fixing.
* ~~Enable GraphML output in addition to serialized object binaries, for data portability~~ DONE
* ~~Custom vertex type with attribute for sequence occupancy?~~ ABANDONED
* Have a branch where this is implemented, but there's a bug that broke matching. Don't currently have time to fix.
* Re-implement command line arguments, to enable scripting and statistical simulation studies
* ~~Custom vertex type with attribute for sequence occupancy?~~ DONE
* Advantage: would eliminate the need to use maps to associate vertices with sequences, which would make the code easier to understand.
* ~~Have a branch where this is implemented, but there's a bug that broke matching. Don't currently have time to fix.~~
* ~~Re-implement command line arguments, to enable scripting and statistical simulation studies~~ DONE
* ~~Implement custom Vertex class to simplify code and make it easier to implement different MWM algorithms~~ DONE
* This also seems to be faster when using the same algorithm than the version with lots of maps, which is a nice bonus!
* Re-implement CDR1 matching method
* Implement Duan and Su's maximum weight matching algorithm
* Add controllable algorithm-type parameter?
* This would be fun and valuable, but probably take more time than I have for a hobby project.
* Implement an algorithm for approximating a maximum weight matching
* Some of these run in linear or near-linear time
* given that the underlying biological samples have many, many sources of error, this would probably be the most useful option in practice. It seems less mathematically elegant, though, and so less fun for me.
* Implement Vose's alias method for arbitrary statistical distributions of cells
* Should probably refactor to use apache commons rng for this
* Use commons JCS for caching
* Parameterize pre-filtering. Currently, sequences present in all wells are filtered out before constructing the graph, which massively reduces graph size. But, ideally, no pre-filtering would be necessary.
## CITATIONS
@@ -307,7 +376,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
* [JHeaps](https://www.jheaps.org) -- For pairing heap priority queue used in maximum weight matching algorithm
* [Apache Commons CSV](https://commons.apache.org/proper/commons-csv/) -- For CSV file output
* [Apache Commons CLI](https://commons.apache.org/proper/commons-cli/) -- To enable command line arguments for scripting. (**Awaiting re-implementation**.)
* [Apache Commons CLI](https://commons.apache.org/proper/commons-cli/) -- To enable command line arguments for scripting.
## ACKNOWLEDGEMENTS
BiGpairSEQ was conceived in collaboration with Dr. Alice MacQueen, who brought the original

View File

@@ -13,9 +13,10 @@ public class BiGpairSEQ {
private static boolean cacheCells = false;
private static boolean cachePlate = false;
private static boolean cacheGraph = false;
private static String priorityQueueHeapType = "FIBONACCI";
private static HeapType priorityQueueHeapType = HeapType.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) {
if (args.length == 0) {
@@ -23,8 +24,8 @@ public class BiGpairSEQ {
}
else {
//This will be uncommented when command line arguments are re-implemented.
//CommandLineInterface.startCLI(args);
System.out.println("Command line arguments are still being re-implemented.");
CommandLineInterface.startCLI(args);
//System.out.println("Command line arguments are still being re-implemented.");
}
}
@@ -156,15 +157,15 @@ public class BiGpairSEQ {
}
public static String getPriorityQueueHeapType() {
return priorityQueueHeapType;
return priorityQueueHeapType.name();
}
public static void setPairingHeap() {
priorityQueueHeapType = "PAIRING";
priorityQueueHeapType = HeapType.PAIRING;
}
public static void setFibonacciHeap() {
priorityQueueHeapType = "FIBONACCI";
priorityQueueHeapType = HeapType.FIBONACCI;
}
public static boolean outputBinary() {return outputBinary;}
@@ -172,5 +173,5 @@ public class BiGpairSEQ {
public static boolean outputGraphML() {return outputGraphML;}
public static void setOutputGraphML(boolean b) {outputGraphML = b;}
public static String getVersion() { return version; }
}

View File

@@ -13,8 +13,12 @@ public class CellSample {
List<Integer> numbersCDR3 = new ArrayList<>();
List<Integer> numbersCDR1 = new ArrayList<>();
Integer numDistCDR3s = 2 * numDistinctCells + 1;
//Assign consecutive integers for each CDR3. This ensures they are all unique.
IntStream.range(1, numDistCDR3s + 1).forEach(i -> numbersCDR3.add(i));
//After all CDR3s are assigned, start assigning consecutive integers to CDR1s
//There will usually be fewer integers in the CDR1 list, which will allow repeats below
IntStream.range(numDistCDR3s + 1, numDistCDR3s + 1 + (numDistCDR3s / cdr1Freq) + 1).forEach(i -> numbersCDR1.add(i));
//randomize the order of the numbers in the lists
Collections.shuffle(numbersCDR3);
Collections.shuffle(numbersCDR1);
@@ -22,11 +26,15 @@ public class CellSample {
//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){
//Go through entire CDR3 list once, make pairs of alphas and betas
Integer tmpCDR3a = numbersCDR3.get(i);
Integer tmpCDR3b = numbersCDR3.get(i+1);
//Go through (likely shorter) CDR1 list as many times as necessary, make pairs of alphas and betas
Integer tmpCDR1a = numbersCDR1.get(i % numbersCDR1.size());
Integer tmpCDR1b = numbersCDR1.get((i+1) % numbersCDR1.size());
//Make the array representing the cell
Integer[] tmp = {tmpCDR3a, tmpCDR3b, tmpCDR1a, tmpCDR1b};
//Add the cell to the list of distinct cells
distinctCells.add(tmp);
}
this.cells = distinctCells;

View File

@@ -1,5 +1,9 @@
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
*
@@ -29,6 +33,8 @@ import org.apache.commons.cli.*;
* cellfile : name of the cell sample file to use as input
* platefile : name of the sample plate file to use as input
* output : name of the output file
* graphml : output a graphml file
* binary : output a serialized binary object file
*
* Match flags:
* graphFile : name of graph and data file to use as input
@@ -43,246 +49,444 @@ import org.apache.commons.cli.*;
public class CommandLineInterface {
public static void startCLI(String[] args) {
//These command line options are a big mess
//Really, I don't think command line tools are expected to work in this many different modes
//making cells, making plates, and matching are the sort of thing that UNIX philosophy would say
//should be three separate programs.
//There might be a way to do it with option parameters?
//Options sets for the different modes
Options mainOptions = buildMainOptions();
Options cellOptions = buildCellOptions();
Options plateOptions = buildPlateOptions();
Options graphOptions = buildGraphOptions();
Options matchOptions = buildMatchCDR3options();
//main options set
CommandLineParser parser = new DefaultParser();
try{
CommandLine line = parser.parse(mainOptions, Arrays.copyOfRange(args, 0, 1));
if (line.hasOption("help")) {
HelpFormatter formatter = new HelpFormatter();
formatter.printHelp("BiGpairSEQ_Sim.jar", mainOptions);
System.out.println();
formatter.printHelp("BiGpairSEQ_Sim.jar -cells", cellOptions);
System.out.println();
formatter.printHelp("BiGpairSEQ_Sim.jar -plate", plateOptions);
System.out.println();
formatter.printHelp("BiGpairSEQ_Sim.jar -graph", graphOptions);
System.out.println();
formatter.printHelp("BiGpairSEQ_Sim.jar -match", matchOptions);
}
else if (line.hasOption("version")) {
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");
makeCells(filename, number, diversity);
}
else if (line.hasOption("plate")) {
line = parser.parse(plateOptions, Arrays.copyOfRange(args, 1, args.length));
//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);
}
else if (line.hasOption("pop")) {
populations = Stream.of(line.getOptionValues("pop"))
.mapToInt(Integer::parseInt)
.boxed()
.toArray(Integer[]::new);
}
else{
populations = new Integer[1];
populations[0] = 1;
}
//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();
}
if (line.hasOption("graphml")) { //if told to, output graphml file
GraphMLFileWriter gmlwriter = new GraphMLFileWriter(outputFilename, graph);
gmlwriter.writeGraphToFile();
}
}
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) {
System.err.println("Parsing failed. Reason: " + exp.getMessage());
}
}
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 file of distinct cells")
.desc("Makes a cell sample file of distinct T cells")
.build();
Option makePlate = Option.builder("plates")
.longOpt("make-plates")
.desc("Makes a sample plate file")
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 and data file")
.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("Match CDR3s. Requires a cell sample file and any number of plate files.")
.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;
}
//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")
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")
.build();
mainOptions.addOption(numCells);
Option cdr1Freq = Option.builder("d")
.longOpt("peptide-diversity-factor")
.required().build();
Option cdr3Diversity = Option.builder("d")
.longOpt("diversity-factor")
.desc("The factor by which unique CDR3s outnumber unique CDR1s")
.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);
.argName("factor")
.required().build();
cellOptions.addOption(numCells);
cellOptions.addOption(cdr3Diversity);
cellOptions.addOption(outputFileOption());
return cellOptions;
}
//Options plateOptions = new Options();
Option inputCells = Option.builder("c")
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("file")
.desc("The cell sample file used for filling wells")
.build();
mainOptions.addOption(inputCells);
Option numWells = Option.builder("w")
.longOpt("num-wells")
.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")
.desc("The number of wells on each plate")
.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();
mainOptions.addOption(numWells);
Option numPlates = Option.builder("np")
.longOpt("num-plates")
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("number")
.desc("The number of plate files to output")
.argName("value")
.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")
Option lambda = Option.builder("lambda")
.desc("If using -exponential flag, lambda value for distribution")
.hasArg()
.argName("number")
.desc("Well drop-out rate. (Probability between 0 and 1)")
.argName("value")
.build();
mainOptions.addOption(plateErr);
Option plateConcentrations = Option.builder("t")
.longOpt("t-cells-per-well")
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()
.argName("number 1, number 2, ...")
.desc("Number of T cells per well for each plate section")
.numberOfArgs(2)
.argName("min> <max")
.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")
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()
.desc("Plate files to match")
.argName("number [number]...")
.build();
mainOptions.addOption(inputPlates);
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;
}
CommandLineParser parser = new DefaultParser();
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{
CommandLine line = parser.parse(mainOptions, args);
if(line.hasOption("match")){
//line = parser.parse(mainOptions, args);
//String cellFile = line.getOptionValue("c");
String graphFile = line.getOptionValue("g");
Integer lowThreshold = Integer.valueOf(line.getOptionValue(lowThresh));
Integer highThreshold = Integer.valueOf(line.getOptionValue(highThresh));
Integer occupancyDifference = Integer.valueOf(line.getOptionValue(occDiff));
Integer overlapPercent = Integer.valueOf(line.getOptionValue(overlapPer));
for(String plate: line.getOptionValues("p")) {
matchCDR3s(graphFile, lowThreshold, highThreshold, occupancyDifference, overlapPercent);
}
}
else if(line.hasOption("cells")){
//line = parser.parse(mainOptions, args);
String filename = line.getOptionValue("o");
Integer numDistCells = Integer.valueOf(line.getOptionValue("nc"));
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];
for(int i = 0; i <numSections; i++){
concentrationsToUse[i] = Integer.valueOf(concentrationsToUseString[i]);
}
Double dropOutRate = Double.valueOf(line.getOptionValue("err"));
if(line.hasOption("exponential")){
Double lambda = Double.valueOf(line.getOptionValue("lambda"));
for(int i = 1; i <= numPlatesToMake; i++){
makePlateExp(cellFile, filenamePrefix + i, lambda, numWellsOnPlate,
concentrationsToUse,dropOutRate);
}
}
else if(line.hasOption("gaussian")){
Double stdDev = Double.valueOf(line.getOptionValue("std-dev"));
for(int i = 1; i <= numPlatesToMake; i++){
makePlate(cellFile, filenamePrefix + i, stdDev, numWellsOnPlate,
concentrationsToUse,dropOutRate);
}
GraphDataObjectReader reader = new GraphDataObjectReader(graphFilename, false);
return reader.getData();
}
else if(line.hasOption("poisson")){
for(int i = 1; i <= numPlatesToMake; i++){
makePlatePoisson(cellFile, filenamePrefix + i, numWellsOnPlate,
concentrationsToUse,dropOutRate);
}
}
}
}
catch (ParseException exp) {
System.err.println("Parsing failed. Reason: " + exp.getMessage());
catch (IOException ex) {
ex.printStackTrace();
return null;
}
}
@@ -292,37 +496,4 @@ public class CommandLineInterface {
CellFileWriter writer = new CellFileWriter(filename, sample);
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) {
}
}

View File

@@ -1,10 +1,12 @@
import java.io.*;
public class GraphDataObjectReader {
private GraphWithMapData data;
private String filename;
public GraphDataObjectReader(String filename) throws IOException {
public GraphDataObjectReader(String filename, boolean verbose) throws IOException {
if(!filename.matches(".*\\.ser")){
filename = filename + ".ser";
}
@@ -13,10 +15,13 @@ public class GraphDataObjectReader {
BufferedInputStream fileIn = new BufferedInputStream(new FileInputStream(filename));
ObjectInputStream in = new ObjectInputStream(fileIn))
{
if (verbose) {
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();
} catch (FileNotFoundException | ClassNotFoundException ex) {
System.out.println("Graph/data file " + filename + " not found.");
ex.printStackTrace();
}
}

View File

@@ -1,3 +1,5 @@
import org.jgrapht.Graph;
import java.io.BufferedOutputStream;
import java.io.FileOutputStream;
import java.io.IOException;
@@ -7,6 +9,7 @@ public class GraphDataObjectWriter {
private GraphWithMapData data;
private String filename;
private boolean verbose = true;
public GraphDataObjectWriter(String filename, GraphWithMapData data) {
if(!filename.matches(".*\\.ser")){
@@ -16,13 +19,24 @@ public class GraphDataObjectWriter {
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() {
try (BufferedOutputStream bufferedOut = new BufferedOutputStream(new FileOutputStream(filename));
ObjectOutputStream out = new ObjectOutputStream(bufferedOut);
){
if(verbose) {
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);
} catch (IOException ex) {
ex.printStackTrace();

View File

@@ -3,8 +3,9 @@ import org.jgrapht.graph.SimpleWeightedGraph;
import org.jgrapht.nio.Attribute;
import org.jgrapht.nio.AttributeType;
import org.jgrapht.nio.DefaultAttribute;
import org.jgrapht.nio.dot.DOTExporter;
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.IOException;
@@ -12,14 +13,14 @@ import java.nio.file.Files;
import java.nio.file.Path;
import java.nio.file.StandardOpenOption;
import java.util.HashMap;
import java.util.LinkedHashMap;
import java.util.Map;
public class GraphMLFileWriter {
String filename;
SimpleWeightedGraph graph;
GraphWithMapData data;
Map<String, Attribute> graphAttributes;
public GraphMLFileWriter(String filename, GraphWithMapData data) {
if(!filename.matches(".*\\.graphml")){
@@ -27,52 +28,61 @@ public class GraphMLFileWriter {
}
this.filename = filename;
this.data = data;
this.graph = data.getGraph();
graphAttributes = createGraphAttributes();
}
// public void writeGraphToFile() {
// try(BufferedWriter writer = Files.newBufferedWriter(Path.of(filename), StandardOpenOption.CREATE_NEW);
// ){
// GraphMLExporter<SimpleWeightedGraph, BufferedWriter> exporter = new GraphMLExporter<>();
// exporter.exportGraph(graph, writer);
// } catch(IOException ex){
// System.out.println("Could not make new file named "+filename);
// System.err.println(ex);
// }
// }
public GraphMLFileWriter(String filename, SimpleWeightedGraph<Vertex, DefaultWeightedEdge> graph) {
if(!filename.matches(".*\\.graphml")){
filename = filename + ".graphml";
}
this.filename = filename;
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() {
SimpleWeightedGraph graph = data.getGraph();
Map<Integer, Integer> vertexToAlphaMap = data.getPlateVtoAMap();
Map<Integer, Integer> vertexToBetaMap = data.getPlateVtoBMap();
Map<Integer, Integer> alphaOccs = data.getAlphaWellCounts();
Map<Integer, Integer> betaOccs = data.getBetaWellCounts();
try(BufferedWriter writer = Files.newBufferedWriter(Path.of(filename), StandardOpenOption.CREATE_NEW);
){
//create exporter. Let the vertex labels be the unique ids for the vertices
GraphMLExporter<Integer, SimpleWeightedGraph<Vertex, DefaultWeightedEdge>> exporter = new GraphMLExporter<>(v -> v.toString());
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<>();
if(vertexToAlphaMap.containsKey(v)) {
attributes.put("type", DefaultAttribute.createAttribute("CDR3 Alpha"));
attributes.put("sequence", DefaultAttribute.createAttribute(vertexToAlphaMap.get(v)));
attributes.put("occupancy", DefaultAttribute.createAttribute(
alphaOccs.get(vertexToAlphaMap.get(v))));
}
else if(vertexToBetaMap.containsKey(v)) {
attributes.put("type", DefaultAttribute.createAttribute("CDR3 Beta"));
attributes.put("sequence", DefaultAttribute.createAttribute(vertexToBetaMap.get(v)));
attributes.put("occupancy", DefaultAttribute.createAttribute(
betaOccs.get(vertexToBetaMap.get(v))));
}
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
exporter.registerAttribute("type", GraphMLExporter.AttributeCategory.NODE, AttributeType.STRING);
exporter.registerAttribute("sequence", GraphMLExporter.AttributeCategory.NODE, AttributeType.STRING);
exporter.registerAttribute("occupancy", GraphMLExporter.AttributeCategory.NODE, AttributeType.STRING);
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);
} catch(IOException ex){
@@ -81,4 +91,3 @@ public class GraphMLFileWriter {
}
}
}

View File

@@ -2,23 +2,25 @@ import org.jgrapht.graph.DefaultWeightedEdge;
import org.jgrapht.graph.SimpleWeightedGraph;
import java.util.ArrayList;
import java.util.HashMap;
import java.util.List;
import java.util.Map;
public interface GraphModificationFunctions {
//remove over- and under-weight edges
static List<Integer[]> filterByOverlapThresholds(SimpleWeightedGraph<Integer, DefaultWeightedEdge> graph,
//remove over- and under-weight edges, return removed edges
static Map<Vertex[], Integer> filterByOverlapThresholds(SimpleWeightedGraph<Vertex, DefaultWeightedEdge> graph,
int low, int high, boolean saveEdges) {
List<Integer[]> removedEdges = new ArrayList<>();
Map<Vertex[], Integer> removedEdges = new HashMap<>();
//List<Integer[]> removedEdges = new ArrayList<>();
for (DefaultWeightedEdge e : graph.edgeSet()) {
if ((graph.getEdgeWeight(e) > high) || (graph.getEdgeWeight(e) < low)) {
if(saveEdges) {
Integer source = graph.getEdgeSource(e);
Integer target = graph.getEdgeTarget(e);
Vertex source = graph.getEdgeSource(e);
Vertex target = graph.getEdgeTarget(e);
Integer weight = (int) graph.getEdgeWeight(e);
Integer[] edge = {source, target, weight};
removedEdges.add(edge);
Vertex[] edge = {source, target};
removedEdges.put(edge, weight);
}
else {
graph.setEdgeWeight(e, 0.0);
@@ -26,31 +28,27 @@ public interface GraphModificationFunctions {
}
}
if(saveEdges) {
for (Integer[] edge : removedEdges) {
for (Vertex[] edge : removedEdges.keySet()) {
graph.removeEdge(edge[0], edge[1]);
}
}
return removedEdges;
}
//Remove edges for pairs with large occupancy discrepancy
static List<Integer[]> filterByRelativeOccupancy(SimpleWeightedGraph<Integer, DefaultWeightedEdge> graph,
Map<Integer, Integer> alphaWellCounts,
Map<Integer, Integer> betaWellCounts,
Map<Integer, Integer> plateVtoAMap,
Map<Integer, Integer> plateVtoBMap,
//Remove edges for pairs with large occupancy discrepancy, return removed edges
static Map<Vertex[], Integer> filterByRelativeOccupancy(SimpleWeightedGraph<Vertex, DefaultWeightedEdge> graph,
Integer maxOccupancyDifference, boolean saveEdges) {
List<Integer[]> removedEdges = new ArrayList<>();
Map<Vertex[], Integer> removedEdges = new HashMap<>();
for (DefaultWeightedEdge e : graph.edgeSet()) {
Integer alphaOcc = alphaWellCounts.get(plateVtoAMap.get(graph.getEdgeSource(e)));
Integer betaOcc = betaWellCounts.get(plateVtoBMap.get(graph.getEdgeTarget(e)));
Integer alphaOcc = graph.getEdgeSource(e).getOccupancy();
Integer betaOcc = graph.getEdgeTarget(e).getOccupancy();
if (Math.abs(alphaOcc - betaOcc) >= maxOccupancyDifference) {
if (saveEdges) {
Integer source = graph.getEdgeSource(e);
Integer target = graph.getEdgeTarget(e);
Vertex source = graph.getEdgeSource(e);
Vertex target = graph.getEdgeTarget(e);
Integer weight = (int) graph.getEdgeWeight(e);
Integer[] edge = {source, target, weight};
removedEdges.add(edge);
Vertex[] edge = {source, target};
removedEdges.put(edge, weight);
}
else {
graph.setEdgeWeight(e, 0.0);
@@ -58,34 +56,30 @@ public interface GraphModificationFunctions {
}
}
if(saveEdges) {
for (Integer[] edge : removedEdges) {
for (Vertex[] edge : removedEdges.keySet()) {
graph.removeEdge(edge[0], edge[1]);
}
}
return removedEdges;
}
//Remove edges for pairs where overlap size is significantly lower than the well occupancy
static List<Integer[]> filterByOverlapPercent(SimpleWeightedGraph<Integer, DefaultWeightedEdge> graph,
Map<Integer, Integer> alphaWellCounts,
Map<Integer, Integer> betaWellCounts,
Map<Integer, Integer> plateVtoAMap,
Map<Integer, Integer> plateVtoBMap,
//Remove edges for pairs where overlap size is significantly lower than the well occupancy, return removed edges
static Map<Vertex[], Integer> filterByOverlapPercent(SimpleWeightedGraph<Vertex, DefaultWeightedEdge> graph,
Integer minOverlapPercent,
boolean saveEdges) {
List<Integer[]> removedEdges = new ArrayList<>();
Map<Vertex[], Integer> removedEdges = new HashMap<>();
for (DefaultWeightedEdge e : graph.edgeSet()) {
Integer alphaOcc = alphaWellCounts.get(plateVtoAMap.get(graph.getEdgeSource(e)));
Integer betaOcc = betaWellCounts.get(plateVtoBMap.get(graph.getEdgeTarget(e)));
Integer alphaOcc = graph.getEdgeSource(e).getOccupancy();
Integer betaOcc = graph.getEdgeTarget(e).getOccupancy();
double weight = graph.getEdgeWeight(e);
double min = minOverlapPercent / 100.0;
if ((weight / alphaOcc < min) || (weight / betaOcc < min)) {
if (saveEdges) {
Integer source = graph.getEdgeSource(e);
Integer target = graph.getEdgeTarget(e);
Vertex source = graph.getEdgeSource(e);
Vertex target = graph.getEdgeTarget(e);
Integer intWeight = (int) graph.getEdgeWeight(e);
Integer[] edge = {source, target, intWeight};
removedEdges.add(edge);
Vertex[] edge = {source, target};
removedEdges.put(edge, intWeight);
}
else {
graph.setEdgeWeight(e, 0.0);
@@ -93,18 +87,18 @@ public interface GraphModificationFunctions {
}
}
if(saveEdges) {
for (Integer[] edge : removedEdges) {
for (Vertex[] edge : removedEdges.keySet()) {
graph.removeEdge(edge[0], edge[1]);
}
}
return removedEdges;
}
static void addRemovedEdges(SimpleWeightedGraph<Integer, DefaultWeightedEdge> graph,
List<Integer[]> removedEdges) {
for (Integer[] edge : removedEdges) {
static void addRemovedEdges(SimpleWeightedGraph<Vertex, DefaultWeightedEdge> graph,
Map<Vertex[], Integer> removedEdges) {
for (Vertex[] edge : removedEdges.keySet()) {
DefaultWeightedEdge e = graph.addEdge(edge[0], edge[1]);
graph.setEdgeWeight(e, (double) edge[2]);
graph.setEdgeWeight(e, removedEdges.get(edge));
}
}

View File

@@ -6,6 +6,7 @@ import java.util.Map;
//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.
//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 {
private String sourceFilename;
@@ -15,32 +16,33 @@ public class GraphWithMapData implements java.io.Serializable {
private Integer alphaCount;
private Integer betaCount;
private final Map<Integer, Integer> distCellsMapAlphaKey;
private final Map<Integer, Integer> plateVtoAMap;
private final Map<Integer, Integer> plateVtoBMap;
private final Map<Integer, Integer> plateAtoVMap;
private final Map<Integer, Integer> plateBtoVMap;
private final Map<Integer, Integer> alphaWellCounts;
private final Map<Integer, Integer> betaWellCounts;
// private final Map<Integer, Integer> plateVtoAMap;
// private final Map<Integer, Integer> plateVtoBMap;
// private final Map<Integer, Integer> plateAtoVMap;
// private final Map<Integer, Integer> plateBtoVMap;
// private final Map<Integer, Integer> alphaWellCounts;
// private final Map<Integer, Integer> betaWellCounts;
private final Duration time;
public GraphWithMapData(SimpleWeightedGraph graph, Integer numWells, Integer[] wellConcentrations,
Integer alphaCount, Integer betaCount,
Map<Integer, Integer> distCellsMapAlphaKey, Map<Integer, Integer> plateVtoAMap,
Map<Integer,Integer> plateVtoBMap, Map<Integer, Integer> plateAtoVMap,
Map<Integer, Integer> plateBtoVMap, Map<Integer, Integer> alphaWellCounts,
Map<Integer, Integer> betaWellCounts, Duration time) {
Map<Integer, Integer> distCellsMapAlphaKey, Integer alphaCount, Integer betaCount, Duration time){
// Map<Integer, Integer> plateVtoAMap,
// Map<Integer,Integer> plateVtoBMap, Map<Integer, Integer> plateAtoVMap,
// Map<Integer, Integer> plateBtoVMap, Map<Integer, Integer> alphaWellCounts,
// Map<Integer, Integer> betaWellCounts,) {
this.graph = graph;
this.numWells = numWells;
this.wellPopulations = wellConcentrations;
this.alphaCount = alphaCount;
this.betaCount = betaCount;
this.distCellsMapAlphaKey = distCellsMapAlphaKey;
this.plateVtoAMap = plateVtoAMap;
this.plateVtoBMap = plateVtoBMap;
this.plateAtoVMap = plateAtoVMap;
this.plateBtoVMap = plateBtoVMap;
this.alphaWellCounts = alphaWellCounts;
this.betaWellCounts = betaWellCounts;
// this.plateVtoAMap = plateVtoAMap;
// this.plateVtoBMap = plateVtoBMap;
// this.plateAtoVMap = plateAtoVMap;
// this.plateBtoVMap = plateBtoVMap;
// this.alphaWellCounts = alphaWellCounts;
// this.betaWellCounts = betaWellCounts;
this.time = time;
}
@@ -68,29 +70,29 @@ public class GraphWithMapData implements java.io.Serializable {
return distCellsMapAlphaKey;
}
public Map<Integer, Integer> getPlateVtoAMap() {
return plateVtoAMap;
}
public Map<Integer, Integer> getPlateVtoBMap() {
return plateVtoBMap;
}
public Map<Integer, Integer> getPlateAtoVMap() {
return plateAtoVMap;
}
public Map<Integer, Integer> getPlateBtoVMap() {
return plateBtoVMap;
}
public Map<Integer, Integer> getAlphaWellCounts() {
return alphaWellCounts;
}
public Map<Integer, Integer> getBetaWellCounts() {
return betaWellCounts;
}
// public Map<Integer, Integer> getPlateVtoAMap() {
// return plateVtoAMap;
// }
//
// public Map<Integer, Integer> getPlateVtoBMap() {
// return plateVtoBMap;
// }
//
// public Map<Integer, Integer> getPlateAtoVMap() {
// return plateAtoVMap;
// }
//
// public Map<Integer, Integer> getPlateBtoVMap() {
// return plateBtoVMap;
// }
//
// public Map<Integer, Integer> getAlphaWellCounts() {
// return alphaWellCounts;
// }
//
// public Map<Integer, Integer> getBetaWellCounts() {
// return betaWellCounts;
// }
public Duration getTime() {
return time;

View File

@@ -0,0 +1,4 @@
public enum HeapType {
FIBONACCI,
PAIRING
}

View File

@@ -227,16 +227,14 @@ public class InteractiveInterface {
Plate samplePlate;
PlateFileWriter writer;
if(exponential){
samplePlate = new Plate(numWells, dropOutRate, populations);
samplePlate.fillWellsExponential(cellFile, cells.getCells(), lambda);
samplePlate = new Plate(cells, cellFile, numWells, populations, dropOutRate, lambda, true);
writer = new PlateFileWriter(filename, samplePlate);
}
else {
if (poisson) {
stdDev = Math.sqrt(cells.getCellCount()); //gaussian with square root of elements approximates poisson
}
samplePlate = new Plate(numWells, dropOutRate, populations);
samplePlate.fillWells(cellFile, cells.getCells(), stdDev);
samplePlate = new Plate(cells, cellFile, numWells, populations, dropOutRate, stdDev, false);
writer = new PlateFileWriter(filename, samplePlate);
}
System.out.println("Writing Sample Plate to file");
@@ -260,7 +258,7 @@ public class InteractiveInterface {
cellFile = sc.next();
System.out.print("\nPlease enter name of an existing sample plate file: ");
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: ");
filename = sc.next();
} catch (InputMismatchException ex) {
@@ -292,7 +290,7 @@ public class InteractiveInterface {
else {
System.out.println("Reading Sample Plate file: " + plateFile);
PlateFileReader plateReader = new PlateFileReader(plateFile);
plate = new Plate(plateReader.getFilename(), plateReader.getWells());
plate = plateReader.getSamplePlate();
if(BiGpairSEQ.cachePlate()) {
BiGpairSEQ.setPlateInMemory(plate, plateFile);
}
@@ -306,8 +304,7 @@ public class InteractiveInterface {
System.out.println("Returning to main menu.");
}
else{
List<Integer[]> cells = cellSample.getCells();
GraphWithMapData data = Simulator.makeGraph(cells, plate, true);
GraphWithMapData data = Simulator.makeGraph(cellSample, plate, true);
assert filename != null;
if(BiGpairSEQ.outputBinary()) {
GraphDataObjectWriter dataWriter = new GraphDataObjectWriter(filename, data);
@@ -378,7 +375,7 @@ public class InteractiveInterface {
data = BiGpairSEQ.getGraphInMemory();
}
else {
GraphDataObjectReader dataReader = new GraphDataObjectReader(graphFilename);
GraphDataObjectReader dataReader = new GraphDataObjectReader(graphFilename, true);
data = dataReader.getData();
if(BiGpairSEQ.cacheGraph()) {
BiGpairSEQ.setGraphInMemory(data, graphFilename);
@@ -507,7 +504,7 @@ public class InteractiveInterface {
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("4) Turn " + getOnOff(!BiGpairSEQ.outputBinary()) + " serialized binary graph output");
System.out.println("5) Turn " + getOnOff(!BiGpairSEQ.outputGraphML()) + " GraphML 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");
try {
@@ -573,6 +570,8 @@ public class InteractiveInterface {
}
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("of the pairSEQ algorithm for pairing T cell receptor sequences.");
System.out.println();

View File

@@ -21,15 +21,15 @@ public class MatchingResult {
* well populations *
* total alphas found *
* total betas found *
* high overlap threshold
* low overlap threshold
* maximum occupancy difference
* minimum overlap percent
* pairing attempt rate
* correct pairing count
* incorrect pairing count
* pairing error rate
* simulation time
* high overlap threshold *
* low overlap threshold *
* maximum occupancy difference *
* minimum overlap percent *
* pairing attempt rate *
* correct pairing count *
* incorrect pairing count *
* pairing error rate *
* simulation time (seconds)
*/
this.metadata = metadata;
this.comments = new ArrayList<>();
@@ -91,6 +91,22 @@ public class MatchingResult {
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)"); }
}

View File

@@ -8,7 +8,9 @@ TODO: Implement discrete frequency distributions using Vose's Alias Method
import java.util.*;
public class Plate {
private CellSample cells;
private String sourceFile;
private String filename;
private List<List<Integer[]>> wells;
private final Random rand = BiGpairSEQ.getRand();
private int size;
@@ -18,6 +20,25 @@ public class Plate {
private double lambda;
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) {
this.size = size;
@@ -26,8 +47,9 @@ public class Plate {
wells = new ArrayList<>();
}
public Plate(String sourceFileName, List<List<Integer[]>> wells) {
this.sourceFile = sourceFileName;
//constructor for returning a Plate from a PlateFileReader
public Plate(String filename, List<List<Integer[]>> wells) {
this.filename = filename;
this.wells = wells;
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;
exponential = true;
sourceFile = sourceFileName;
int numSections = populations.length;
int section = 0;
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;
sourceFile = sourceFileName;
int numSections = populations.length;
int section = 0;
double m;
@@ -138,9 +158,9 @@ public class Plate {
//returns a map of the counts of the sequence at cell index sIndex, in a range of wells
public Map<Integer, Integer> assayWellsSequenceS(int start, int end, int... sIndices) {
Map<Integer,Integer> assay = new HashMap<>();
for(int pIndex: sIndices){
for(int sIndex: sIndices){
for(int i = start; i < end; i++){
countSequences(assay, wells.get(i), pIndex);
countSequences(assay, wells.get(i), sIndex);
}
}
return assay;
@@ -149,6 +169,7 @@ public class Plate {
private void countSequences(Map<Integer, Integer> wellMap, List<Integer[]> well, int... sIndices) {
for(Integer[] cell : well) {
for(int sIndex: sIndices){
//skip dropout sequences, which have value -1
if(cell[sIndex] != -1){
wellMap.merge(cell[sIndex], 1, (oldValue, newValue) -> oldValue + newValue);
}
@@ -159,4 +180,6 @@ public class Plate {
public String getSourceFileName() {
return sourceFile;
}
public String getFilename() { return filename; }
}

View File

@@ -56,11 +56,8 @@ public class PlateFileReader {
}
public List<List<Integer[]>> getWells() {
return wells;
public Plate getSamplePlate() {
return new Plate(filename, wells);
}
public String getFilename() {
return filename;
}
}

View 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
}

View File

@@ -18,16 +18,17 @@ import static java.lang.Float.*;
//NOTE: "sequence" in method and variable names refers to a peptide sequence from a simulated T cell
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;
//Make the graph needed for matching CDR3s
public static GraphWithMapData makeGraph(List<Integer[]> distinctCells, Plate samplePlate, boolean verbose) {
//Make the graph needed for matching sequences.
//sourceVertexIndices and targetVertexIndices are indices within the cell to use as for the two sets of vertices
//in the bipartite graph. "Source" and "target" are JGraphT terms for the two vertices an edge touches,
//even if not directed.
public static GraphWithMapData makeGraph(CellSample cellSample, Plate samplePlate, boolean verbose) {
Instant start = Instant.now();
int[] alphaIndex = {cdr3AlphaIndex};
int[] betaIndex = {cdr3BetaIndex};
List<Integer[]> distinctCells = cellSample.getCells();
int[] alphaIndices = {SequenceType.CDR3_ALPHA.ordinal()};
int[] betaIndices = {SequenceType.CDR3_BETA.ordinal()};
int numWells = samplePlate.getSize();
@@ -37,14 +38,16 @@ public class Simulator implements GraphModificationFunctions {
if(verbose){System.out.println("Cell maps made");}
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();
if(verbose){System.out.println("All alphas count: " + alphaCount);}
int betaCount = allBetas.size();
if(verbose){System.out.println("All betas count: " + betaCount);}
if(verbose){System.out.println("Well maps made");}
//ideally we wouldn't do any graph pre-filtering. But sequences present in all wells add a huge number of edges to the graph and don't carry any signal value
if(verbose){System.out.println("Removing sequences present in all wells.");}
filterByOccupancyThresholds(allAlphas, 1, numWells - 1);
filterByOccupancyThresholds(allBetas, 1, numWells - 1);
@@ -77,29 +80,46 @@ public class Simulator implements GraphModificationFunctions {
//(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.)
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<>();
//count how many wells each beta appears in
//count how many wells each beta sequence appears in
Map<Integer, Integer> betaWellCounts = new HashMap<>();
//the adjacency matrix to be used by the graph generator
double[][] weights = new double[plateVtoAMap.size()][plateVtoBMap.size()];
countSequencesAndFillMatrix(samplePlate, allAlphas, allBetas, plateAtoVMap,
plateBtoVMap, alphaIndex, betaIndex, alphaWellCounts, betaWellCounts, weights);
plateBtoVMap, alphaIndices, betaIndices, alphaWellCounts, betaWellCounts, weights);
if(verbose){System.out.println("Matrix created");}
//create bipartite graph
if(verbose){System.out.println("Creating graph");}
//the graph object
SimpleWeightedGraph<Integer, DefaultWeightedEdge> graph =
SimpleWeightedGraph<Vertex, DefaultWeightedEdge> graph =
new SimpleWeightedGraph<>(DefaultWeightedEdge.class);
//the graph generator
SimpleWeightedBipartiteGraphMatrixGenerator graphGenerator = new SimpleWeightedBipartiteGraphMatrixGenerator();
//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);
}
//Sort to make sure the order of vertices in list matches the order of the adjacency matrix
Collections.sort(alphaVertices);
//Add ordered list of vertices to the graph
graphGenerator.first(alphaVertices);
//the list of beta vertices
List<Integer> betaVertices = new ArrayList<>(plateVtoBMap.keySet());
graphGenerator.second(betaVertices); //This will work because LinkedHashMap preserves order of entry
//List<Integer> betaVertices = new ArrayList<>(plateVtoBMap.keySet());//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);
}
//Sort to make sure the order of vertices in list matches the order of the adjacency matrix
Collections.sort(betaVertices);
//Add ordered list of vertices to the graph
graphGenerator.second(betaVertices);
//use adjacency matrix of weight created previously
graphGenerator.weights(weights);
graphGenerator.generateGraph(graph);
@@ -109,11 +129,9 @@ public class Simulator implements GraphModificationFunctions {
Duration time = Duration.between(start, stop);
//create GraphWithMapData object
GraphWithMapData output = new GraphWithMapData(graph, numWells, samplePlate.getPopulations(), alphaCount, betaCount,
distCellsMapAlphaKey, plateVtoAMap, plateVtoBMap, plateAtoVMap,
plateBtoVMap, alphaWellCounts, betaWellCounts, time);
GraphWithMapData output = new GraphWithMapData(graph, numWells, samplePlate.getPopulations(), distCellsMapAlphaKey, alphaCount, betaCount, time);
//Set source file name in graph to name of sample plate
output.setSourceFilename(samplePlate.getSourceFileName());
output.setSourceFilename(samplePlate.getFilename());
//return GraphWithMapData object
return output;
}
@@ -123,61 +141,67 @@ public class Simulator implements GraphModificationFunctions {
Integer highThreshold, Integer maxOccupancyDifference,
Integer minOverlapPercent, boolean verbose) {
Instant start = Instant.now();
List<Integer[]> removedEdges = new ArrayList<>();
SimpleWeightedGraph<Vertex, DefaultWeightedEdge> graph = data.getGraph();
Map<Vertex[], Integer> removedEdges = new HashMap<>();
boolean saveEdges = BiGpairSEQ.cacheGraph();
int numWells = data.getNumWells();
Integer alphaCount = data.getAlphaCount();
Integer betaCount = data.getBetaCount();
//Integer alphaCount = data.getAlphaCount();
//Integer betaCount = data.getBetaCount();
Map<Integer, Integer> distCellsMapAlphaKey = data.getDistCellsMapAlphaKey();
Map<Integer, Integer> plateVtoAMap = data.getPlateVtoAMap();
Map<Integer, Integer> plateVtoBMap = data.getPlateVtoBMap();
Map<Integer, Integer> alphaWellCounts = data.getAlphaWellCounts();
Map<Integer, Integer> betaWellCounts = data.getBetaWellCounts();
SimpleWeightedGraph<Integer, DefaultWeightedEdge> graph = data.getGraph();
Set<Vertex> alphas = new HashSet<>();
Set<Vertex> betas = new HashSet<>();
for(Vertex v: graph.vertexSet()) {
if (SequenceType.CDR3_ALPHA.equals(v.getType())){
alphas.add(v);
}
else {
betas.add(v);
}
}
Integer graphAlphaCount = alphas.size();
Integer graphBetaCount = betas.size();
//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");}
removedEdges.addAll(GraphModificationFunctions.filterByOverlapThresholds(graph, lowThreshold, highThreshold, saveEdges));
removedEdges.putAll(GraphModificationFunctions.filterByOverlapThresholds(graph, lowThreshold, highThreshold, saveEdges));
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
if(verbose){System.out.println("Eliminating edges with weights less than " + minOverlapPercent.toString() +
" percent of vertex occupancy value.");}
removedEdges.addAll(GraphModificationFunctions.filterByOverlapPercent(graph, alphaWellCounts, betaWellCounts,
plateVtoAMap, plateVtoBMap, minOverlapPercent, saveEdges));
removedEdges.putAll(GraphModificationFunctions.filterByOverlapPercent(graph, minOverlapPercent, saveEdges));
if(verbose){System.out.println("Edges with weights too far below a vertex occupancy value removed");}
//Filter by relative occupancy
if(verbose){System.out.println("Eliminating edges between vertices with occupancy difference > "
+ maxOccupancyDifference);}
removedEdges.addAll(GraphModificationFunctions.filterByRelativeOccupancy(graph, alphaWellCounts, betaWellCounts,
plateVtoAMap, plateVtoBMap, maxOccupancyDifference, saveEdges));
removedEdges.putAll(GraphModificationFunctions.filterByRelativeOccupancy(graph, maxOccupancyDifference, saveEdges));
if(verbose){System.out.println("Edges between vertices of with excessively different occupancy values " +
"removed");}
//Find Maximum Weighted Matching
//Find Maximum Weight Matching
//using jheaps library class PairingHeap for improved efficiency
if(verbose){System.out.println("Finding maximum weighted matching");}
if(verbose){System.out.println("Finding maximum weight matching");}
MaximumWeightBipartiteMatching maxWeightMatching;
//Use correct heap type for priority queue
String heapType = BiGpairSEQ.getPriorityQueueHeapType();
switch (heapType) {
case "PAIRING" -> {
maxWeightMatching = new MaximumWeightBipartiteMatching(graph,
plateVtoAMap.keySet(),
plateVtoBMap.keySet(),
alphas,
betas,
i -> new PairingHeap(Comparator.naturalOrder()));
}
case "FIBONACCI" -> {
maxWeightMatching = new MaximumWeightBipartiteMatching(graph,
plateVtoAMap.keySet(),
plateVtoBMap.keySet(),
alphas,
betas,
i -> new FibonacciHeap(Comparator.naturalOrder()));
}
default -> {
maxWeightMatching = new MaximumWeightBipartiteMatching(graph,
plateVtoAMap.keySet(),
plateVtoBMap.keySet());
alphas,
betas);
}
}
//get the matching
@@ -207,11 +231,14 @@ public class Simulator implements GraphModificationFunctions {
Map<Integer, Integer> matchMap = new HashMap<>();
while(weightIter.hasNext()) {
e = weightIter.next();
Integer source = graph.getEdgeSource(e);
Integer target = graph.getEdgeTarget(e);
Vertex source = graph.getEdgeSource(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!
matchMap.put(plateVtoAMap.get(source), plateVtoBMap.get(target));
check = plateVtoBMap.get(target).equals(distCellsMapAlphaKey.get(plateVtoAMap.get(source)));
matchMap.put(source.getSequence(), target.getSequence());
check = target.getSequence().equals(distCellsMapAlphaKey.get(source.getSequence()));
//check = plateVtoBMap.get(target).equals(distCellsMapAlphaKey.get(plateVtoAMap.get(source)));
if(check) {
trueCount++;
}
@@ -219,17 +246,19 @@ public class Simulator implements GraphModificationFunctions {
falseCount++;
}
List<String> result = new ArrayList<>();
result.add(plateVtoAMap.get(source).toString());
//alpha sequence
result.add(source.getSequence().toString());
//alpha well count
result.add(alphaWellCounts.get(plateVtoAMap.get(source)).toString());
result.add(plateVtoBMap.get(target).toString());
result.add(source.getOccupancy().toString());
//beta sequence
result.add(target.getSequence().toString());
//beta well count
result.add(betaWellCounts.get(plateVtoBMap.get(target)).toString());
result.add(target.getOccupancy().toString());
//overlap count
result.add(Double.toString(graph.getEdgeWeight(e)));
result.add(Boolean.toString(check));
double pValue = Equations.pValue(numWells, alphaWellCounts.get(plateVtoAMap.get(source)),
betaWellCounts.get(plateVtoBMap.get(target)), graph.getEdgeWeight(e));
double pValue = Equations.pValue(numWells, source.getOccupancy(),
target.getOccupancy(), graph.getEdgeWeight(e));
BigDecimal pValueTrunc = new BigDecimal(pValue, mc);
result.add(pValueTrunc.toString());
allResults.add(result);
@@ -237,18 +266,20 @@ public class Simulator implements GraphModificationFunctions {
//Metadata comments for CSV file
String algoType = "LEDA book with heap: " + heapType;
int min = Math.min(alphaCount, betaCount);
int min = Math.min(graphAlphaCount, graphBetaCount);
//matching weight
BigDecimal totalMatchingWeight = maxWeightMatching.getMatchingWeight();
//rate of attempted matching
double attemptRate = (double) (trueCount + falseCount) / min;
BigDecimal attemptRateTrunc = new BigDecimal(attemptRate, mc);
//rate of pairing error
double pairingErrorRate = (double) falseCount / (trueCount + falseCount);
BigDecimal pairingErrorRateTrunc;
if(pairingErrorRate == NaN || pairingErrorRate == POSITIVE_INFINITY || pairingErrorRate == NEGATIVE_INFINITY) {
pairingErrorRateTrunc = new BigDecimal(-1, mc);
if(Double.isFinite(pairingErrorRate)) {
pairingErrorRateTrunc = new BigDecimal(pairingErrorRate, mc);
}
else{
pairingErrorRateTrunc = new BigDecimal(pairingErrorRate, mc);
pairingErrorRateTrunc = new BigDecimal(-1, mc);
}
//get list of well populations
Integer[] wellPopulations = data.getWellPopulations();
@@ -268,9 +299,12 @@ public class Simulator implements GraphModificationFunctions {
metadata.put("sample plate filename", data.getSourceFilename());
metadata.put("graph filename", dataFilename);
metadata.put("algorithm type", algoType);
metadata.put("matching weight", totalMatchingWeight.toString());
metadata.put("well populations", wellPopulationsString);
metadata.put("total alphas found", alphaCount.toString());
metadata.put("total betas found", betaCount.toString());
metadata.put("total alphas on plate", data.getAlphaCount().toString());
metadata.put("total betas on plate", data.getBetaCount().toString());
metadata.put("alphas in graph (after pre-filtering)", graphAlphaCount.toString());
metadata.put("betas in graph (after pre-filtering)", graphBetaCount.toString());
metadata.put("high overlap threshold", highThreshold.toString());
metadata.put("low overlap threshold", lowThreshold.toString());
metadata.put("minimum overlap percent", minOverlapPercent.toString());
@@ -279,7 +313,7 @@ public class Simulator implements GraphModificationFunctions {
metadata.put("correct pairing count", Integer.toString(trueCount));
metadata.put("incorrect pairing count", Integer.toString(falseCount));
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
MatchingResult output = new MatchingResult(metadata, header, allResults, matchMap, time);
if(verbose){

View File

@@ -1,23 +1,97 @@
import java.io.Serializable;
public class Vertex implements Serializable, Comparable<Vertex> {
private SequenceType type;
private Integer vertexLabel;
private Integer sequence;
private Integer occupancy;
public class Vertex {
private final Integer vertexLabel;
private final Integer sequence;
private final Integer occupancy;
public Vertex(Integer vertexLabel) {
this.vertexLabel = vertexLabel;
}
public Vertex(String vertexLabel) {
this.vertexLabel = Integer.parseInt((vertexLabel));
}
public Vertex(Integer vertexLabel, Integer sequence, Integer occupancy) {
public Vertex(SequenceType type, Integer sequence, Integer occupancy, Integer vertexLabel) {
this.type = type;
this.vertexLabel = vertexLabel;
this.sequence = sequence;
this.occupancy = occupancy;
}
public Integer getVertexLabel() { return vertexLabel; }
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() {
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();
}
@Override
public int compareTo(Vertex other) {
return this.vertexLabel - other.getVertexLabel();
}
}