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

Initial Commit - works well, but too memory intensive

Browse Source
This commit is contained in:
efischer
2021-11-09 20:24:22 -06:00
commit 82d5500ae8
3 changed files with 339 additions and 0 deletions
Download Patch File Download Diff File Expand all files Collapse all files

6
src/main/java/Equations.java Normal file
View File

@@ -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;
}
}

102
src/main/java/Plate.java Normal file
View File

@@ -0,0 +1,102 @@
import java.sql.Array;
import java.util.*;
public class Plate {
private List<List<Integer[]>> 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<Integer[]> 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<Integer[]> 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<Integer, Integer> assayWellsAlpha() {
return this.assayWellsAlpha(0, size);
}
public Map<Integer, Integer> assayWellsAlpha(int n) {
return this.assayWellsAlpha(n, n+1);
}
public Map<Integer, Integer> assayWellsAlpha(int start, int end) {
Map<Integer, Integer> 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<Integer, Integer> wellMap, List<Integer[]> 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<Integer, Integer> assayWellsBeta() {
return this.assayWellsBeta(0, size);
}
//assays a specific well
public Map<Integer, Integer> assayWellsBeta(int n) {
return this.assayWellsBeta(n, n+1);
}
//assays a range of wells
public Map<Integer, Integer> assayWellsBeta(int start, int end) {
Map<Integer, Integer> 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<Integer, Integer> wellMap, List<Integer[]> well){
for(Integer[] cell : well) {
if(cell[1] != -1){
wellMap.merge(cell[1], 1, (oldValue, newValue) -> oldValue + newValue);
}
}
}
}

231
src/main/java/Simulation.java Normal file
View File

@@ -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<Integer> 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<Integer[]> 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<Integer, Integer> distCellsMapAlphaKey = new HashMap<>();
for (Integer[] cell : distinctCells) {
distCellsMapAlphaKey.put(cell[0], cell[1]);
}
//HashMap keyed to Betas, values Alphas
Map<Integer, Integer> 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<Integer, Integer> allAlphas = samplePlate.assayWellsAlpha();
Map<Integer, Integer> allBetas = samplePlate.assayWellsBeta();
int alphaCount = allAlphas.size();
int betaCount = allBetas.size();
//keys are sequential integer vertices, values are alphas
Map<String, Integer> plateVtoAMap = getVertexToPeptideMap(allAlphas, "A");
//keys are sequential integers vertices, values are betas
Map<String, Integer> plateVtoBMap = getVertexToPeptideMap(allBetas, "B");
//keys are alphas, values are sequential integer vertices from previous map
Map<Integer, String> plateAtoVMap = invertVertexMap(plateVtoAMap);
//keys are betas, values are sequential integer vertices from previous map
Map<Integer, String> plateBtoVMap = invertVertexMap(plateVtoBMap);
//OUTPUT
System.out.println("Maps made");
System.out.println("Creating Graph");
//construct graph
SimpleWeightedGraph<String, DefaultWeightedEdge> graph = new SimpleWeightedGraph<String, DefaultWeightedEdge>(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<Integer, Integer> wellNAlphas = null;
Map<Integer, Integer> wellNBetas = null;
//Count how many wells each alpha appears in
Map<Integer, Integer> alphaWellCounts = new HashMap<>();
//count how many wells each beta appears in
Map<Integer, Integer> 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<String, DefaultWeightedEdge> graphMatching = maxWeightMatching.getMatching();
Instant stop = Instant.now();
//OUTPUT
System.out.println("Matching completed");
Iterator<DefaultWeightedEdge> 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<String, Integer> getVertexToPeptideMap(Map<Integer, Integer> peptides, String prefix) {
Map<String, Integer> map = new LinkedHashMap<>();
Integer index = 0;
for (Integer k: peptides.keySet()) {
map.put(prefix + index, k);
index++;
}
return map;
}
private static Map<Integer, String> invertVertexMap(Map<String, Integer> map) {
Map<Integer, String> 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<Integer,Integer> AtoBMap, Map<String,Integer> VtoAMap, Map<String, Integer> 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;
}
}
}