commit 82d5500ae849869be79815f82bcb4a8c4adf5be7 Author: efischer Date: Tue Nov 9 20:24:22 2021 -0600 Initial Commit - works well, but too memory intensive diff --git a/src/main/java/Equations.java b/src/main/java/Equations.java new file mode 100644 index 0000000..908133c --- /dev/null +++ b/src/main/java/Equations.java @@ -0,0 +1,6 @@ +public abstract class Equations { + + public static double pValue(Integer w, Integer w_a, Integer w_b, Integer w_ab) { + return 1.0; + } +} diff --git a/src/main/java/Plate.java b/src/main/java/Plate.java new file mode 100644 index 0000000..5b5d81a --- /dev/null +++ b/src/main/java/Plate.java @@ -0,0 +1,102 @@ +import java.sql.Array; +import java.util.*; + +public class Plate { + private List> wells; + private Random rand = new Random(); + private int size; + private double error; + + public Plate (int size, double error) { + this.size = size; + this.error = error; + wells = new ArrayList<>(); + } + + public void fillWells(List cells, int[] concentrations, double stdDev) { + int numSections = concentrations.length; + int section = 0; + double m; + int n; + boolean drop; + while (section < numSections){ + for (int i = 0; i < (size / numSections); i++) { + List well = new ArrayList<>(); + for (int j = 0; j < concentrations[section]; j++) { + do { + m = Math.abs(rand.nextGaussian()) * stdDev; + } while (m >= cells.size()); + n = (int) Math.floor(m); + Integer[] cellToAdd = cells.get(n).clone(); + drop = Math.abs(rand.nextDouble()) < error; + if (drop) { + if (rand.nextBoolean()) { + cellToAdd[0] = -1; + } else { + cellToAdd[1] = -1; + } + } + well.add(cellToAdd); + } + wells.add(well); + } + section++; + } + } + + public Map assayWellsAlpha() { + return this.assayWellsAlpha(0, size); + } + + public Map assayWellsAlpha(int n) { + return this.assayWellsAlpha(n, n+1); + } + + public Map assayWellsAlpha(int start, int end) { + Map assay = new HashMap<>(); + for(int i = start; i < end; i++){ + countAlphas(assay, wells.get(i)); + } + return assay; + } + + //given a map, counts distinct alphas in a well + private void countAlphas(Map wellMap, List well){ + for(Integer[] cell : well) { + if(cell[0] != -1){ + //keys are alphas, value is how many of them have been assayed + wellMap.merge(cell[0], 1, (oldValue, newValue) -> oldValue + newValue); + } + } + } + + //assays all wells + public Map assayWellsBeta() { + return this.assayWellsBeta(0, size); + } + + //assays a specific well + public Map assayWellsBeta(int n) { + return this.assayWellsBeta(n, n+1); + } + + //assays a range of wells + public Map assayWellsBeta(int start, int end) { + Map assay = new HashMap<>(); + for(int i = start; i < end; i++){ + countBetas(assay, wells.get(i)); + } + return assay; + } + + //given a map, counts distinct betas in a well + private void countBetas(Map wellMap, List well){ + for(Integer[] cell : well) { + if(cell[1] != -1){ + wellMap.merge(cell[1], 1, (oldValue, newValue) -> oldValue + newValue); + } + } + } + + +} diff --git a/src/main/java/Simulation.java b/src/main/java/Simulation.java new file mode 100644 index 0000000..0b9bf4e --- /dev/null +++ b/src/main/java/Simulation.java @@ -0,0 +1,231 @@ +import org.jgrapht.Graph; +import org.jgrapht.alg.interfaces.MatchingAlgorithm; +import org.jgrapht.alg.matching.MaximumWeightBipartiteMatching; +import org.jgrapht.graph.DefaultWeightedEdge; +import org.jgrapht.graph.SimpleWeightedGraph; + +import java.math.BigDecimal; +import java.math.MathContext; +import java.text.NumberFormat; +import java.time.Duration; +import java.time.Instant; +import java.util.*; +import java.util.stream.IntStream; + +public class Simulation { + private static Integer numDistinctCells = 4_000_000; + private static double stdDeviation; + private static int numWells = 96; + private static int numConcentrations = 1; + private static double errorRate = 0.1; + private static int[] concentrations = {50, 500, 5000}; + private static int lowThreshold = 2; //min number of shared wells to attempt pairing + private static int highThreshold = 96; //max number of shared wells to attempt pairing + + public static void main(String[] args) { + Instant start = Instant.now(); + + stdDeviation = Math.sqrt(numDistinctCells.doubleValue()); + //concentrations = new int[numConcentrations]; + + //OUTPUT + System.out.println("Generating cells"); + + List numbers = new ArrayList<>(); + IntStream.range(1, (2 * numDistinctCells) + 1).forEach(i -> numbers.add(i)); + Collections.shuffle(numbers); + + //Each cell represented by two numbers from the random permutation + //These represent unique alpha and beta peptides + List distinctCells = new ArrayList<>(); + for(int i = 0; i < numbers.size() - 1; i = i + 2) { + Integer tmp1 = numbers.get(i); + Integer tmp2 = numbers.get(i+1); + Integer[] tmp = {tmp1, tmp2}; + distinctCells.add(tmp); + } + + //OUTPUT + System.out.println("Cells generated"); + System.out.println("Filling wells"); + + //HashMap keyed to Alphas, values Betas + Map distCellsMapAlphaKey = new HashMap<>(); + for (Integer[] cell : distinctCells) { + distCellsMapAlphaKey.put(cell[0], cell[1]); + } + + //HashMap keyed to Betas, values Alphas + Map distCellsMapBetaKey = new HashMap<>(); + for (Integer[] cell : distinctCells) { + distCellsMapBetaKey.put(cell[1], cell[0]); + } + + Plate samplePlate = new Plate(numWells, errorRate); + samplePlate.fillWells(distinctCells, concentrations, stdDeviation); + + //OUTPUT + System.out.println("Wells filled"); + System.out.println("Making maps"); + + //All alphas on plate, with + Map allAlphas = samplePlate.assayWellsAlpha(); + Map allBetas = samplePlate.assayWellsBeta(); + int alphaCount = allAlphas.size(); + int betaCount = allBetas.size(); + + //keys are sequential integer vertices, values are alphas + Map plateVtoAMap = getVertexToPeptideMap(allAlphas, "A"); + //keys are sequential integers vertices, values are betas + Map plateVtoBMap = getVertexToPeptideMap(allBetas, "B"); + //keys are alphas, values are sequential integer vertices from previous map + Map plateAtoVMap = invertVertexMap(plateVtoAMap); + //keys are betas, values are sequential integer vertices from previous map + Map plateBtoVMap = invertVertexMap(plateVtoBMap); + + //OUTPUT + System.out.println("Maps made"); + System.out.println("Creating Graph"); + + //construct graph + SimpleWeightedGraph graph = new SimpleWeightedGraph(DefaultWeightedEdge.class); + //add vertices + for(String v: plateVtoAMap.keySet()) { + graph.addVertex(v); + } + for(String v: plateVtoBMap.keySet()) { + graph.addVertex(v); + } + //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 wellNAlphas = null; + Map wellNBetas = null; + //Count how many wells each alpha appears in + Map alphaWellCounts = new HashMap<>(); + //count how many wells each beta appears in + Map betaWellCounts = new HashMap<>(); + for (int n = 0; n < numWells; n++) { + wellNAlphas = samplePlate.assayWellsAlpha(n); + for(Integer a: wellNAlphas.keySet()){ + alphaWellCounts.merge(a, 1, (oldValue, newValue) -> oldValue + newValue); + } + wellNBetas = samplePlate.assayWellsBeta(n); + for(Integer b: wellNBetas.keySet()){ + betaWellCounts.merge(b, 1, (oldValue, newValue) -> oldValue + newValue); + } + for(Integer i: wellNAlphas.keySet()) { + for(Integer j: wellNBetas.keySet()){ + addEdgeOrWeight(graph, plateAtoVMap.get(i), plateBtoVMap.get(j)); + } + } + } + + //OUTPUT + System.out.println("Graph created"); + System.out.println("Finding maximum weighted matching"); + + MaximumWeightBipartiteMatching maxWeightMatching = new MaximumWeightBipartiteMatching(graph, plateVtoAMap.keySet(), plateVtoBMap.keySet()); + MatchingAlgorithm.Matching graphMatching = maxWeightMatching.getMatching(); + + Instant stop = Instant.now(); + + //OUTPUT + System.out.println("Matching completed"); + + + Iterator weightIter = graphMatching.iterator(); + DefaultWeightedEdge e = null; + int trueCount = 0; + int falseCount = 0; + boolean check = false; + while(weightIter.hasNext()) { + e = weightIter.next(); + if(graph.getEdgeWeight(e) <= 2.0) { + continue; + } + String source = graph.getEdgeSource(e); + String target = graph.getEdgeTarget(e); + //Need map of each peptide and number of wells it appears in. + check =printData(graph.getEdgeWeight(e), source, alphaWellCounts.get(plateVtoAMap.get(source)), target, betaWellCounts.get(plateVtoBMap.get(target)), distCellsMapAlphaKey, plateVtoAMap, plateVtoBMap); + if(check) { + trueCount++; + } + else { + falseCount++; + } + } + + //RESULTS SUMMARY + NumberFormat nf = NumberFormat.getInstance(Locale.US); + System.out.println(""); + System.out.println("Distinct cell population: " + nf.format(numDistinctCells)); + System.out.print("Sample plate: " + numWells +" wells, "); + System.out.print("divided into " + concentrations.length + " sections of "); + System.out.println(Arrays.toString(concentrations) + " cells per well."); + System.out.print("Attempted to match alpha/beta overlaps of at least " + lowThreshold + " wells "); + System.out.println("and no more than " + highThreshold + " wells."); + System.out.println(nf.format(alphaCount) + " alpha peptides found"); + System.out.println(nf.format(betaCount) + " beta peptides found"); + MathContext mc = new MathContext(3); + int min = alphaCount > betaCount ? betaCount : alphaCount; + double attemptRate = (double) (trueCount + falseCount) / min; + BigDecimal attR = new BigDecimal(attemptRate, mc); + System.out.println("Pairing attempt rate: " + attR); + System.out.println("Correct pairings: " + nf.format(trueCount)); + System.out.println("Incorrect pairings: " + nf.format(falseCount)); + double experimentalError = (double) falseCount / (trueCount + falseCount); + BigDecimal expErr = new BigDecimal(experimentalError, mc); + System.out.println("Pairing error rate: " + expErr); + Duration time = Duration.between(start, stop); + System.out.println("Simulation time: " + nf.format(time.toSeconds()) + " seconds"); + System.out.println(""); + + + } + + private static Map getVertexToPeptideMap(Map peptides, String prefix) { + Map map = new LinkedHashMap<>(); + Integer index = 0; + for (Integer k: peptides.keySet()) { + map.put(prefix + index, k); + index++; + } + return map; + } + + private static Map invertVertexMap(Map map) { + Map inverse = new HashMap<>(); + for(String k: map.keySet()) { + inverse.put(map.get(k), k); + } + return inverse; + } + + private static void addEdgeOrWeight(Graph g, String v1, String v2) { + if (g.containsEdge(v1, v2)) { + g.setEdgeWeight(v1, v2, g.getEdgeWeight(g.getEdge(v1, v2)) + 1.0); + } + else { + g.addEdge(v1, v2); + } + //System.out.println("added!"); + } + + private static boolean printData (Double edgeWeight, String source, Integer sourceCount, String target, Integer targetCount, Map AtoBMap, Map VtoAMap, Map VtoBMap) { + //placeholder for real calculations + int trueCount = 0; + int falseCount = 0; + System.out.println("Alpha peptide " + source + " present in " + sourceCount + " wells"); + System.out.println("Beta peptide " + target + " present in " + targetCount + " wells"); + System.out.println("Both alpha and beta present in " + edgeWeight + " wells"); + if(VtoBMap.get(target).equals(AtoBMap.get(VtoAMap.get(source)))){ + System.out.println("True"); + return true; + } + else { + System.out.println("False"); + System.out.println("The real pair is " + source + " and " + AtoBMap.get(VtoAMap.get(source))); + return false; + } + } +}