diff --git a/src/main/java/GraphModificationFunctions.java b/src/main/java/GraphModificationFunctions.java index 4b5591a..5a0f2cf 100644 --- a/src/main/java/GraphModificationFunctions.java +++ b/src/main/java/GraphModificationFunctions.java @@ -12,7 +12,6 @@ public interface GraphModificationFunctions { static Map filterByOverlapThresholds(SimpleWeightedGraph graph, int low, int high, boolean saveEdges) { Map removedEdges = new HashMap<>(); - //List removedEdges = new ArrayList<>(); for (DefaultWeightedEdge e : graph.edgeSet()) { if ((graph.getEdgeWeight(e) > high) || (graph.getEdgeWeight(e) < low)) { if(saveEdges) { @@ -94,6 +93,37 @@ public interface GraphModificationFunctions { return removedEdges; } + static Map filterByRelativeReadCount (SimpleWeightedGraph graph, Integer threshold, boolean saveEdges) { + Map removedEdges = new HashMap<>(); + Boolean passes; + for (DefaultWeightedEdge e : graph.edgeSet()) { + Integer alphaReadCount = graph.getEdgeSource(e).getReadCount(); + Integer betaReadCount = graph.getEdgeTarget(e).getReadCount(); + passes = RelativeReadCountFilterFunction(threshold, alphaReadCount, betaReadCount); + if (!passes) { + if (saveEdges) { + Vertex source = graph.getEdgeSource(e); + Vertex target = graph.getEdgeTarget(e); + Integer intWeight = (int) graph.getEdgeWeight(e); + Vertex[] edge = {source, target}; + removedEdges.put(edge, intWeight); + } + else { + graph.setEdgeWeight(e, 0.0); + } + } + } + if(saveEdges) { + for (Vertex[] edge : removedEdges.keySet()) { + graph.removeEdge(edge[0], edge[1]); + } + } + return removedEdges; + } + + static Boolean RelativeReadCountFilterFunction(Integer threshold, Integer alphaReadCount, Integer betaReadCount) { + return Math.abs(alphaReadCount - betaReadCount) < threshold; + } static void addRemovedEdges(SimpleWeightedGraph graph, Map removedEdges) { for (Vertex[] edge : removedEdges.keySet()) { @@ -102,4 +132,6 @@ public interface GraphModificationFunctions { } } + + } diff --git a/src/main/java/GraphWithMapData.java b/src/main/java/GraphWithMapData.java index 8588088..2d30519 100644 --- a/src/main/java/GraphWithMapData.java +++ b/src/main/java/GraphWithMapData.java @@ -15,6 +15,7 @@ public class GraphWithMapData implements java.io.Serializable { private Integer[] wellPopulations; private Integer alphaCount; private Integer betaCount; + private int readDepth; private final Map distCellsMapAlphaKey; // private final Map plateVtoAMap; // private final Map plateVtoBMap; @@ -25,7 +26,7 @@ public class GraphWithMapData implements java.io.Serializable { private final Duration time; public GraphWithMapData(SimpleWeightedGraph graph, Integer numWells, Integer[] wellConcentrations, - Map distCellsMapAlphaKey, Integer alphaCount, Integer betaCount, Duration time){ + Map distCellsMapAlphaKey, Integer alphaCount, Integer betaCount, int readDepth, Duration time){ // Map plateVtoAMap, // Map plateVtoBMap, Map plateAtoVMap, @@ -94,6 +95,8 @@ public class GraphWithMapData implements java.io.Serializable { // return betaWellCounts; // } + public int getReadDepth() { return readDepth; } + public Duration getTime() { return time; } diff --git a/src/main/java/InteractiveInterface.java b/src/main/java/InteractiveInterface.java index 8982d2d..25304c4 100644 --- a/src/main/java/InteractiveInterface.java +++ b/src/main/java/InteractiveInterface.java @@ -273,6 +273,9 @@ public class InteractiveInterface { if (simulateReadDepth) { System.out.print("\nPlease enter read depth (the integer number of reads per sequence): "); readDepth = sc.nextInt(); + if(readDepth < 1) { + throw new InputMismatchException("The read depth must be an integer >= 1"); + } System.out.print("\nPlease enter probability of a sequence read error (0.0 to 1.0): "); readErrorRate = sc.nextDouble(); if(readErrorRate < 0.0 || readErrorRate > 1.0) { @@ -332,7 +335,7 @@ public class InteractiveInterface { System.out.println("Returning to main menu."); } else{ - GraphWithMapData data = Simulator.makeGraph(cellSample, plate, true); + GraphWithMapData data = Simulator.makeGraph(cellSample, plate, readDepth, readErrorRate, errorCollisionRate, true); assert filename != null; if(BiGpairSEQ.outputBinary()) { GraphDataObjectWriter dataWriter = new GraphDataObjectWriter(filename, data); diff --git a/src/main/java/Plate.java b/src/main/java/Plate.java index 75ca2c8..de7d872 100644 --- a/src/main/java/Plate.java +++ b/src/main/java/Plate.java @@ -178,17 +178,30 @@ public class Plate { } } + //returns a map of the counts of the sequence at cell index sIndex, in a specific of wells + //Simulates read depth and read errors, counts the number of reads of a unique sequence into the given map. + public void assayWellsSequenceSWithReadDepth(Map occupancyMap, Map readCountMap, + int readDepth, double readErrorProb, double errorCollisionProb, int... sIndices) { + this.assayWellsSequenceSWithReadDepth(occupancyMap, readCountMap, readDepth, readErrorProb, errorCollisionProb, 0, size, sIndices); + } + //returns a map of the counts of the sequence at cell index sIndex, in a specific of wells + //Simulates read depth and read errors, counts the number of reads of a unique sequence into the given map. + public void assayWellsSequenceSWithReadDepth(Map occupancyMap, Map readCountMap, + int readDepth, double readErrorProb, double errorCollisionProb, + int n, int... sIndices) { + this.assayWellsSequenceSWithReadDepth(occupancyMap, readCountMap, readDepth, readErrorProb, errorCollisionProb, n, n+1, sIndices); + } //returns a map of the counts of the sequence at cell index sIndex, in a range of wells //Simulates read depth and read errors, counts the number of reads of a unique sequence into the given map. - public Map assayWellsSequenceSWithReadDepth(int start, int end, int... sIndices) { - Map assay = new HashMap<>(); + public void assayWellsSequenceSWithReadDepth(Map occupancyMap, Map readCountMap, + int readDepth, double readErrorProb, double errorCollisionProb, + int start, int end, int... sIndices) { for(int sIndex: sIndices){ for(int i = start; i < end; i++){ - countSequences(assay, wells.get(i), sIndex); + countSequencesWithReadDepth(occupancyMap, readCountMap, readDepth, readErrorProb, errorCollisionProb, wells.get(i), sIndex); } } - return assay; } //For the sequences at cell indices sIndices, counts number of unique sequences in the given well into the given map //Simulates read depth and read errors, counts the number of reads of a unique sequence into the given map. diff --git a/src/main/java/Simulator.java b/src/main/java/Simulator.java index 1db1853..9b30d7c 100644 --- a/src/main/java/Simulator.java +++ b/src/main/java/Simulator.java @@ -21,8 +21,12 @@ public class Simulator implements GraphModificationFunctions { //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) { + public static GraphWithMapData makeGraph(CellSample cellSample, Plate samplePlate, int readDepth, double readErrorRate, double errorCollisionRate, boolean verbose) { Instant start = Instant.now(); + Map allAlphas; + Map allBetas; + Map alphaReadCounts = null; + Map betaReadCounts = null; List distinctCells = cellSample.getCells(); int[] alphaIndices = {SequenceType.CDR3_ALPHA.ordinal()}; int[] betaIndices = {SequenceType.CDR3_BETA.ordinal()}; @@ -35,11 +39,21 @@ public class Simulator implements GraphModificationFunctions { if(verbose){System.out.println("Cell maps made");} if(verbose){System.out.println("Making well maps");} + if(readDepth == 1) { + allAlphas = new HashMap<>(); + samplePlate.assayWellsSequenceS(allAlphas, alphaIndices); + allBetas = new HashMap<>(); + samplePlate.assayWellsSequenceS(allBetas, betaIndices); + } + else { + allAlphas = new HashMap<>(); + alphaReadCounts = new HashMap<>(); + samplePlate.assayWellsSequenceSWithReadDepth(allAlphas, alphaReadCounts, readDepth, readErrorRate, errorCollisionRate, alphaIndices); + allBetas = new HashMap<>(); + betaReadCounts = new HashMap<>(); + samplePlate.assayWellsSequenceSWithReadDepth(allBetas, betaReadCounts, readDepth, readErrorRate, errorCollisionRate, betaIndices); + } - Map allAlphas = new HashMap<>(); - samplePlate.assayWellsSequenceS(allAlphas, alphaIndices); - Map allBetas = new HashMap<>(); - samplePlate.assayWellsSequenceS(allBetas, betaIndices); int alphaCount = allAlphas.size(); if(verbose){System.out.println("All alphas count: " + alphaCount);} int betaCount = allBetas.size(); @@ -101,7 +115,13 @@ public class Simulator implements GraphModificationFunctions { List alphaVertices = new ArrayList<>(); //start with map of all alphas mapped to vertex values, get occupancy from the alphaWellCounts map for (String seq : plateAtoVMap.keySet()) { - Vertex alphaVertex = new Vertex(SequenceType.CDR3_ALPHA, seq, alphaWellCounts.get(seq), plateAtoVMap.get(seq)); + Vertex alphaVertex; + if (readDepth == 1) { + alphaVertex = new Vertex(SequenceType.CDR3_ALPHA, seq, alphaWellCounts.get(seq), plateAtoVMap.get(seq)); + } + else { + alphaVertex = new Vertex(SequenceType.CDR3_ALPHA, seq, alphaWellCounts.get(seq), plateAtoVMap.get(seq), alphaReadCounts.get(seq)); + } alphaVertices.add(alphaVertex); } //Sort to make sure the order of vertices in list matches the order of the adjacency matrix @@ -112,7 +132,13 @@ public class Simulator implements GraphModificationFunctions { //List betaVertices = new ArrayList<>(plateVtoBMap.keySet());//This will work because LinkedHashMap preserves order of entry List betaVertices = new ArrayList<>(); for (String seq : plateBtoVMap.keySet()) { - Vertex betaVertex = new Vertex(SequenceType.CDR3_BETA, seq, betaWellCounts.get(seq), plateBtoVMap.get(seq)); + Vertex betaVertex; + if (readDepth == 1) { + betaVertex = new Vertex(SequenceType.CDR3_BETA, seq, betaWellCounts.get(seq), plateBtoVMap.get(seq)); + } + else { + betaVertex = new Vertex(SequenceType.CDR3_BETA, seq, betaWellCounts.get(seq), plateBtoVMap.get(seq), betaReadCounts.get(seq)); + } betaVertices.add(betaVertex); } //Sort to make sure the order of vertices in list matches the order of the adjacency matrix @@ -128,7 +154,7 @@ public class Simulator implements GraphModificationFunctions { Duration time = Duration.between(start, stop); //create GraphWithMapData object - GraphWithMapData output = new GraphWithMapData(graph, numWells, samplePlate.getPopulations(), distCellsMapAlphaKey, alphaCount, betaCount, time); + GraphWithMapData output = new GraphWithMapData(graph, numWells, samplePlate.getPopulations(), distCellsMapAlphaKey, alphaCount, betaCount, readDepth, time); //Set source file name in graph to name of sample plate output.setSourceFilename(samplePlate.getFilename()); //return GraphWithMapData object diff --git a/src/main/java/Vertex.java b/src/main/java/Vertex.java index 13a8dfc..90deacf 100644 --- a/src/main/java/Vertex.java +++ b/src/main/java/Vertex.java @@ -23,6 +23,7 @@ public class Vertex implements Serializable, Comparable { this.vertexLabel = vertexLabel; this.sequence = sequence; this.occupancy = occupancy; + this.readCount = 1; } public Vertex(SequenceType type, String sequence, Integer occupancy, Integer vertexLabel, Integer readCount) {