add support for growth without subsampling (#60)

Project.toml

@@ -1,7 +1,7 @@
 name = "Crispulator"
 uuid = "39046a72-5a6f-4d18-ae88-b64c04dba158"
 authors = ["Tamas Nagy <[email protected]>"]
-version = "0.4.1"
+version = "0.5.0"
 
 [deps]
 ArgParse = "c7e460c6-2fb9-53a9-8c5b-16f535851c63"
@@ -18,6 +18,7 @@ HypothesisTests = "09f84164-cd44-5f33-b23f-e6b0d136a0d5"
 IterTools = "c8e1da08-722c-5040-9ed9-7db0dc04731e"
 Measures = "442fdcdd-2543-5da2-b0f3-8c86c306513e"
 ProgressMeter = "92933f4c-e287-5a05-a399-4b506db050ca"
+Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 StatsBase = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91"
 YAML = "ddb6d928-2868-570f-bddf-ab3f9cf99eb6"
 

src/Crispulator.jl

@@ -8,6 +8,7 @@ using IterTools
 using DocStringExtensions
 using Combinatorics
 using DataStructures
+using Random
 
 include("common.jl")
 include("utils.jl")

src/selection.jl

@@ -40,7 +40,7 @@ function grow!(cells::AbstractArray{Int},
                setup::GrowthScreen)
     num_inserted::Int = 0
     noise_dist = Normal(0, setup.noise)
-    @inbounds for i in 1:length(cells)
+    for i in 1:length(cells)
         ρ::Float64 = cell_phenotypes[i]
         ρ_noisy = ρ + rand(noise_dist)
         decision = abs(ρ_noisy) < rand() ? 2 : 2^trunc(Int, 1 + sign(ρ_noisy))
@@ -55,31 +55,80 @@ end
 """
 $(SIGNATURES)
 
-Growth Screen selection
+Given a population `initial_cells` with `initial_cell_phenotypes`, this function
+models the selection step of a [GrowthScreen](@ref). Cells are allowed to double
+according to their growth phenotype and then the population is subsampled to
+simulate passaging to maintain the `bottleneck_representation` as defined in
+`setup`. The number of doublings + subsamplings are controlled by the
+`num_bottlenecks` parameter in `setup`.
+
+So for a screen where there are three bottlenecks and the
+`length(initial_cells) == bottleneck_representation` then we would expect the
+following behavior:
+
+               Cell population
+  ┌────────────────────────────────────────┐
+2 │           .-           .r|          ..|│
+  │         r*`|         _-' |        .r/ |│
+  │       ."   |       .*    |      .r'   |│
+  │    .-"      |    ,"`     |    ,/'     |│
+  │  _-`        | .-/        |  .*`       |│
+1 │_*           Lr'          |r/          |│
+  └────────────────────────────────────────┘
+   0                                      3
+                 Bottleneck #
+
+If the coverage of the library by `initial_cells` is less than the
+`bottleneck_representation` than the population will be allowed to double until
+passaging is needed or the `num_bottlenecks` is hit, which ever comes first.
 """
 function select(setup::GrowthScreen,
                 initial_cells::AbstractArray{Int},
                 initial_cell_phenotypes::AbstractArray{Float64},
-                guides::Vector{Barcode};
-                debug=false)
+                guides::Vector{Barcode})
 
     bottleneck_representation = setup.bottleneck_representation
     num_bottlenecks = setup.num_bottlenecks
-    # all cells at all timepoints
+
+    # pre-allocated vectors
     cellmat = zeros(Int, length(guides)*bottleneck_representation)
     cpmat = zeros(Float64, size(cellmat))
     output_c = Array{Int}(undef, length(initial_cells)*4);
     output_p = Array{Float64}(undef, size(output_c));
-    cells = initial_cells # 1st timepoint slice
-    picked = Array{Int}(undef, size(cellmat, 1))
+
+    # 1st timepoint slice
+    cells = initial_cells
     cell_phenotypes = initial_cell_phenotypes
 
+    num_cells = 0
+    picked = Array{Int}(undef, size(cellmat, 1))
+
+    @debug "Initial cell population: $(length(initial_cells))"
+
     for k in 1:num_bottlenecks
-        num_inserted = grow!(cells, cell_phenotypes, output_c, output_p, setup)
-        cells, cell_phenotypes = view(cellmat, :), view(cpmat, :)
-        sample!(1:num_inserted, picked, replace=false)
-        copy!(view(cellmat, :), view(output_c, picked))
-        copy!(view(cpmat, :), view(output_p, picked))
+        num_cells = grow!(cells, cell_phenotypes, output_c, output_p, setup)
+        @debug "Doubling $k: cell population grew to $num_cells"
+
+        # if we have more cells than we are planning to keep, select a random
+        # subpopulation without replacement
+        if num_cells > length(picked)
+            sample!(1:num_cells, picked, replace=false)
+            @debug "Selecting $(length(picked))"
+        else
+            # if we have room to grow, don't subsample, just shuffle and grow
+            # the available space
+            picked .= 0
+            picked[1:num_cells] .= shuffle(1:num_cells)
+
+            resize!(output_c, num_cells * 4)
+            resize!(output_p, num_cells * 4)
+        end
+
+        ids = view(picked, picked .> 0)
+        num_cells = length(ids)
+        copy!(view(cellmat, 1:num_cells), view(output_c, ids))
+        copy!(view(cpmat, 1:num_cells), view(output_p, ids))
+        cells, cell_phenotypes = view(cellmat, 1:length(ids)), view(cpmat, 1:length(ids))
     end
-    return OrderedDict(:bin1 => initial_cells, :bin2 => cellmat)
+    return OrderedDict(:bin1 => initial_cells, :bin2 => cellmat[1:num_cells])
 end

test/selectionmethods.jl

@@ -63,6 +63,29 @@ function testgrowth(ideal, num)
     all(Bool[isapprox(x[1], x[2], atol=0.1) for x in tocompare])
 end
 
+function testloggrowth()
+    s = GrowthScreen()
+    s.num_genes = 3
+    s.coverage = 2
+    s.num_bottlenecks = 5
+    s.bottleneck_representation = 10
+
+    # true log growth where the cells double exactly once in one doubling period
+    max_phenotype_dists = Dict{Symbol, Tuple{Float64, Sampleable}}(
+           :inactive => (1, Delta(0.0)),
+       );
+    lib = Library(max_phenotype_dists, CRISPRi())
+    guides, guide_freqs_dist = construct_library(s, lib);
+    cells, cell_phenotypes = transfect(s, lib, guides, guide_freqs_dist);
+
+    s.bottleneck_representation = 1000 # no subsampling needed
+
+    bin_cells = select(s, cells, cell_phenotypes, guides)
+
+    num_doublings = log2(length(bin_cells[:bin2]) / length(bin_cells[:bin1]))
+    isapprox(num_doublings, 5, atol = 0.1)
+end
+
 begin
     Random.seed!(1)
     @test testfacs(0.5, -0.50138209, 1/3)
@@ -71,5 +94,8 @@ begin
 
     @test testgrowth([1/7, 2/7, 4/7]./(1/3), 1)
     @test testgrowth([1/21, 4/21,16/21]./(1/3), 2)
+
+    @test testloggrowth()
+
     Random.seed!()
 end