efischer/BiGpairSEQ
1
0
Fork 0
Code Issues Pull Requests Projects Releases 13 Wiki Activity

Expanding simulation to attempt CDR3/CDR1 matching - doubles size of cells

Browse Source
This commit is contained in:
efischer
2021-11-15 15:00:11 -06:00
parent c425b05d8c
commit 13b58e3204
5 changed files with 445 additions and 56 deletions
Download Patch File Download Diff File Expand all files Collapse all files

219
src/main/java/Simulator.java
View File

@@ -71,7 +71,8 @@ public class Simulator {
return new CellSample(distinctCells);
}
public static void matchCDR3s(String filename, List<Integer[]> distinctCells, Plate samplePlate, Integer lowThreshold, Integer highThreshold){
public static MatchingResult matchCDR3s(List<Integer[]> distinctCells,
Plate samplePlate, Integer lowThreshold, Integer highThreshold){
System.out.println("Cells: " + distinctCells.size());
int numWells = samplePlate.getSize();
System.out.println("Making cell maps");
@@ -88,8 +89,8 @@ public class Simulator {
System.out.println("Cell maps made");
System.out.println("Making well maps");
Map<Integer, Integer> allAlphas = samplePlate.assayWellsAlpha();
Map<Integer, Integer> allBetas = samplePlate.assayWellsBeta();
Map<Integer, Integer> allAlphas = samplePlate.assayWellsCDR3Alpha();
Map<Integer, Integer> allBetas = samplePlate.assayWellsCDR3Beta();
int alphaCount = allAlphas.size();
System.out.println("all alphas count: " + alphaCount);
int betaCount = allBetas.size();
@@ -127,11 +128,11 @@ public class Simulator {
new SimpleWeightedGraph<>(DefaultWeightedEdge.class);
double[][] weights = new double[plateVtoAMap.size()][plateVtoBMap.size()];
for (int n = 0; n < numWells; n++) {
wellNAlphas = samplePlate.assayWellsAlpha(n);
wellNAlphas = samplePlate.assayWellsCDR3Alpha(n);
for (Integer a : wellNAlphas.keySet()) {
alphaWellCounts.merge(a, 1, (oldValue, newValue) -> oldValue + newValue);
}
wellNBetas = samplePlate.assayWellsBeta(n);
wellNBetas = samplePlate.assayWellsCDR3Beta(n);
for (Integer b : wellNBetas.keySet()) {
betaWellCounts.merge(b, 1, (oldValue, newValue) -> oldValue + newValue);
}
@@ -178,6 +179,7 @@ public class Simulator {
int trueCount = 0;
int falseCount = 0;
boolean check = false;
Map<Integer, Integer> matchMap = new HashMap<>();
while(weightIter.hasNext()) {
e = weightIter.next();
if(graph.getEdgeWeight(e) < lowThreshold || graph.getEdgeWeight(e) > highThreshold) {
@@ -185,6 +187,8 @@ public class Simulator {
}
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)));
if(check) {
trueCount++;
@@ -224,11 +228,210 @@ public class Simulator {
comments.add("Incorrect pairings: " + falseCount);
comments.add("Pairing error rate: " + pairingErrorRateTrunc);
//result writer
MatchingFileWriter writer = new MatchingFileWriter(filename, comments, header, allResults);
writer.writeResultsToFile();
return new MatchingResult(comments, header, allResults, matchMap);
}
public static MatchingResult matchCDR1s(List<Integer[]> distinctCells,
Plate samplePlate, Integer lowThreshold,
Integer highThreshold, Map<Integer, Integer> previousMatches){
int numWells = samplePlate.getSize();
System.out.println("Making previous match maps");
Map<Integer, Integer> CDR3AtoBMap = previousMatches;
Map<Integer, Integer> CDR3BtoAMap = invertVertexMap(CDR3AtoBMap);
System.out.println("Previous match maps made");
System.out.println("Making cell maps");
Map<Integer, Integer> alphaCDR3toCDR1Map = new HashMap<>();
for (Integer[] cell : distinctCells) {
alphaCDR3toCDR1Map.put(cell[0], cell[2]);
}
//HashMap keyed to Betas, values Alphas
Map<Integer, Integer> betaCDR3toCDR1Map = new HashMap<>();
for (Integer[] cell : distinctCells) {
betaCDR3toCDR1Map.put(cell[1], cell[3]);
}
System.out.println("Cell maps made");
System.out.println("Making well maps");
Map<Integer, Integer> allCDR3s = samplePlate.assayWellsCDR3();
Map<Integer, Integer> allCDR1s = samplePlate.assayWellsCDR1();
int CDR3Count = allCDR3s.size();
System.out.println("all alphas count: " + CDR3Count);
int CDR1Count = allCDR1s.size();
System.out.println("all betas count: " + CDR1Count);
System.out.println("Well maps made");
System.out.println("Making vertex maps");
//Using Integers instead of Strings to label vertices so I can do clever stuff with indices if I need to
// when I refactor to use a 2d array to make the graph
//For the autogenerator, all vertices must have distinct numbers associated with them
Integer vertexStartValue = 0;
//keys are sequential integer vertices, values are CDR3s
Map<Integer, Integer> plateVtoCDR3Map = getVertexToPeptideMap(allCDR3s, vertexStartValue);
//New start value for vertex to beta map should be one more than final vertex value in alpha map
vertexStartValue += plateVtoCDR3Map.size();
//keys are sequential integers vertices, values are CDR1s
Map<Integer, Integer> plateVtoCDR1Map = getVertexToPeptideMap(allCDR1s, vertexStartValue);
//keys are CDR3s, values are sequential integer vertices from previous map
Map<Integer, Integer> plateCDR3toVMap = invertVertexMap(plateVtoCDR3Map);
//keys are CDR1s, values are sequential integer vertices from previous map
Map<Integer, Integer> plateCDR1toVMap = invertVertexMap(plateVtoCDR1Map);
System.out.println("Vertex maps made");
System.out.println("Creating Graph");
//Count how many wells each alpha appears in
Map<Integer, Integer> CDR3WellCounts = new HashMap<>();
//count how many wells each beta appears in
Map<Integer, Integer> CDR1WellCounts = new HashMap<>();
//add edges, where weights are number of wells the peptides share in common.
//If this is too slow, can make a 2d array and use the SimpleWeightedGraphMatrixGenerator class
Map<Integer, Integer> wellNCDR3s = null;
Map<Integer, Integer> wellNCDR1s = null;
SimpleWeightedGraph<Integer, DefaultWeightedEdge> graph =
new SimpleWeightedGraph<>(DefaultWeightedEdge.class);
double[][] weights = new double[plateVtoCDR3Map.size()][plateVtoCDR1Map.size()];
for (int n = 0; n < numWells; n++) {
wellNCDR3s = samplePlate.assayWellsCDR3(n);
for (Integer a : wellNCDR3s.keySet()) {
CDR3WellCounts.merge(a, 1, (oldValue, newValue) -> oldValue + newValue);
}
wellNCDR1s = samplePlate.assayWellsCDR1(n);
for (Integer b : wellNCDR1s.keySet()) {
CDR1WellCounts.merge(b, 1, (oldValue, newValue) -> oldValue + newValue);
}
for (Integer i : wellNCDR3s.keySet()) {
for (Integer j : wellNCDR1s.keySet()) {
weights[plateCDR3toVMap.get(i)][plateCDR1toVMap.get(j) - vertexStartValue] += 1.0;
}
}
}
SimpleWeightedBipartiteGraphMatrixGenerator graphGenerator = new SimpleWeightedBipartiteGraphMatrixGenerator();
List<Integer> CDR3Vertices = new ArrayList<>();
CDR3Vertices.addAll(plateVtoCDR3Map.keySet()); //This will work because LinkedHashMap preserves order of entry
graphGenerator.first(CDR3Vertices);
List<Integer> CDR1Vertices = new ArrayList<>();
CDR1Vertices.addAll(plateVtoCDR1Map.keySet());
graphGenerator.second(CDR1Vertices); //This will work because LinkedHashMap preserves order of entry
graphGenerator.weights(weights);
graphGenerator.generateGraph(graph);
System.out.println("Graph created");
System.out.println("Finding first maximum weighted matching");
MaximumWeightBipartiteMatching firstMaxWeightMatching =
new MaximumWeightBipartiteMatching(graph, plateVtoCDR3Map.keySet(), plateVtoCDR1Map.keySet());
MatchingAlgorithm.Matching<String, DefaultWeightedEdge> graphMatching = firstMaxWeightMatching.getMatching();
System.out.println("First maximum weighted matching found");
//first processing run
Map<Integer, Integer> firstMatchCDR3toCDR1Map = new HashMap<>();
Iterator<DefaultWeightedEdge> weightIter = graphMatching.iterator();
DefaultWeightedEdge e = null;
while(weightIter.hasNext()){
e = weightIter.next();
if(graph.getEdgeWeight(e) < lowThreshold || graph.getEdgeWeight(e) > highThreshold) {
continue;
}
Integer source = graph.getEdgeSource(e);
if(!(CDR3AtoBMap.containsKey(source) || CDR3BtoAMap.containsKey(source))){
continue;
}
Integer target = graph.getEdgeTarget(e);
firstMatchCDR3toCDR1Map.put(plateVtoCDR3Map.get(source), plateVtoCDR1Map.get(target));
}
//zero out the edge weights in the matching
weightIter = graphMatching.iterator();
while(weightIter.hasNext()){
graph.removeEdge(weightIter.next());
}
//Generate a new matching
System.out.println("Finding second maximum weighted matching");
MaximumWeightBipartiteMatching secondMaxWeightMatching =
new MaximumWeightBipartiteMatching(graph, plateVtoCDR3Map.keySet(), plateVtoCDR1Map.keySet());
graphMatching = secondMaxWeightMatching.getMatching();
System.out.println("Second maximum weighted matching found");
//second processing run
Map<Integer, Integer> secondMatchCDR3toCDR1Map = new HashMap<>();
weightIter = graphMatching.iterator();
while(weightIter.hasNext()){
e = weightIter.next();
if(graph.getEdgeWeight(e) < lowThreshold || graph.getEdgeWeight(e) > highThreshold) {
continue;
}
Integer source = graph.getEdgeSource(e);
if(!(CDR3AtoBMap.containsKey(source) || CDR3BtoAMap.containsKey(source))){
continue;
}
Integer target = graph.getEdgeTarget(e);
secondMatchCDR3toCDR1Map.put(plateVtoCDR3Map.get(source), plateVtoCDR1Map.get(target));
}
//Look for matches that simply swapped already-matched alpha and beta CDR3s
Map<Integer, Integer> dualMatchesMap = new LinkedHashMap<>();
for(Integer alphaCDR3: CDR3AtoBMap.keySet()) {
if (!(firstMatchCDR3toCDR1Map.containsKey(alphaCDR3) && secondMatchCDR3toCDR1Map.containsKey(alphaCDR3))) {
continue;
}
Integer betaCDR3 = CDR3AtoBMap.get(alphaCDR3);
if (!(firstMatchCDR3toCDR1Map.containsKey(betaCDR3) && secondMatchCDR3toCDR1Map.containsKey(betaCDR3))) {
continue;
}
if(firstMatchCDR3toCDR1Map.get(alphaCDR3).equals(secondMatchCDR3toCDR1Map.get(betaCDR3))){
if(firstMatchCDR3toCDR1Map.get(betaCDR3).equals(secondMatchCDR3toCDR1Map.get(alphaCDR3))){
dualMatchesMap.put(alphaCDR3, firstMatchCDR3toCDR1Map.get(alphaCDR3));
dualMatchesMap.put(betaCDR3, firstMatchCDR3toCDR1Map.get(betaCDR3));
}
}
}
List<List<String>> allResults = new ArrayList<>();
Integer trueCount = 0;
Iterator iter = dualMatchesMap.keySet().iterator();
while(iter.hasNext()){
Boolean proven = false;
List<String> tmp = new ArrayList<>();
tmp.add(iter.next().toString());
tmp.add(iter.next().toString());
tmp.add(dualMatchesMap.get(Integer.valueOf(tmp.get(0))).toString());
tmp.add(dualMatchesMap.get(Integer.valueOf(tmp.get(1))).toString());
if(alphaCDR3toCDR1Map.get(Integer.valueOf(tmp.get(0))).equals(Integer.valueOf(tmp.get(2)))){
if(betaCDR3toCDR1Map.get(Integer.valueOf(tmp.get(1))).equals(Integer.valueOf(tmp.get(3)))){
proven = true;
}
}
else if(alphaCDR3toCDR1Map.get(Integer.valueOf(tmp.get(0))).equals(Integer.valueOf(tmp.get(3)))){
if(betaCDR3toCDR1Map.get(Integer.valueOf(tmp.get(1))).equals(Integer.valueOf(tmp.get(2)))){
proven = true;
}
}
tmp.add(proven.toString());
allResults.add(tmp);
if(proven){
trueCount++;
}
}
List<String> comments = new ArrayList<>();
comments.add("Previous pairs found: " + previousMatches.size());
comments.add("CDR1 matches attempted: " + allResults.size());
double attemptRate = (double) allResults.size() / previousMatches.size();
comments.add("Matching attempt rate: " + attemptRate);
comments.add("Number of correct matches: " + trueCount);
double correctRate = (double) trueCount / allResults.size();
comments.add("Correct matching rate: " + correctRate);
List<String> headers = new ArrayList<>();
headers.add("CDR3 alpha");
headers.add("CDR3 beta");
headers.add("first matched CDR1");
headers.add("second matched CDR1");
headers.add("Correct match?");
return new MatchingResult(comments, headers, allResults, dualMatchesMap);
}
/*