StemCellTheoryPaperCode/plot_two_compartment_clone_sizes.py at master · jvzonlab/StemCellTheoryPaperCode · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
from typing import List, Tuple

import numpy as np
import matplotlib.pyplot as plt
import pickle as pickle
import matplotlib as mpl

from stem_cell_model import sweeper, two_compartment_model_space, clone_size_distributions
from stem_cell_model.clone_size_distributions import CloneSizeDistribution
from stem_cell_model.parameters import SimulationParameters, SimulationConfig
from stem_cell_model.results import MultiRunStats

mpl.rcParams['pdf.fonttype'] = 42

def load_data() -> List[Tuple[SimulationParameters, MultiRunStats]]:
    return list(sweeper.load_sweep_results("two_comp_sweep_data_fixed_D_aT1"))

sim_data = load_data()
Np=40
start = .025
end = 1

alpha_n_range = np.linspace(start,end,Np)
alpha_m_range = np.linspace(-end,-start,Np)
phi_range = np.linspace(start,end,Np)

# parameters (alpha_n, alpha_m and phi) for which trajectories will be plotted
plot_param = [ [0.2,-0.9,1.0], [0.2,-0.2,1.0], [0.9,-0.9,1.0], [0.9, -0.2, 1.0],
                [0.5,-0.5,0.6], [0.2,-0.2,0.3] ]

plot_run_ind_list=[]
# for each parameter set in plot_param
for s in range(0,len(plot_param)):
    # determine index for run with simulation parameters closest to plot_param[s]
    d_min=6e66
    ind_min=-1
    for i in range(0,len(sim_data)):
        # get parameters for this run
        sweep_param=sim_data[i][0]
        # look at difference in alpha_n,m and phi
        da0 = sweep_param.alpha[0]-plot_param[s][0]
        da1 = sweep_param.alpha[1]-plot_param[s][1]
        dph = sweep_param.phi[0]-plot_param[s][2]
        # if smaller than current minimum distance
        if (da0**2+da1**2+dph**2)<d_min:
            # set new minimum to this distance
            d_min=da0**2+da1**2+dph**2
            # and save index of run
            ind_min=i
    # save index of overal minimum
    plot_run_ind_list.append( ind_min )

fig = plt.figure(figsize=(4, 8))

cmap = plt.cm.jet
cmaplist = [cmap(i) for i in range(cmap.N)]
cmaplist[0] = (.7, .7, .7, 1.0) # force the first color entry to be grey
cmap = mpl.colors.LinearSegmentedColormap.from_list('Custom cmap', cmaplist, cmap.N)


# plot trajectories

random = np.random.Generator(np.random.MT19937(seed=15))
t_lineage_range = (0, 160)
crypts = 50

for n in range(0,len(plot_run_ind_list)):
    i=plot_run_ind_list[n]

    # set model parameters
    params = sim_data[i][0]
    S = int(np.round(params.S))
    alpha = params.alpha
    a = params.a
    phi = params.phi[0]
    N_0 = params.n0[0]
    M_0 = params.n0[1]

    print("%d/%d, a_n:%1.1f, a_m:%1.1f, phi:%1.1f, S:%1.1f, N:%1.1f" % (i,len(sim_data),alpha[0],alpha[1],phi,S,N_0) )

    # run simulations
    t_sim=t_lineage_range[1] + 1
    clone_sizes = CloneSizeDistribution()
    clone_size_duration = 40
    min_clone_count_time, max_clone_count_time = t_lineage_range
    config = SimulationConfig(t_sim=t_sim, track_n_vs_t=True, track_lineage_time_interval=t_lineage_range, random=random)

    for i in range(crypts):  # simulate a crypt 50 times, so that the uncertainty in the clone size is small
        res = two_compartment_model_space.run_simulation_niche(config, params)
        lineages = res.lineages
        clone_sizes.merge(clone_size_distributions.get_clone_size_distributions_with_duration(lineages, min_clone_count_time, max_clone_count_time, clone_size_duration))

    # plot clone sizes
    max_clone_size = max(12, clone_sizes.max())  # Show at least 12 bins
    plt.subplot2grid((6,1),(n,0), rowspan=1)
    plt.text(1,1,"#" + str(n + 1), transform=plt.gca().transAxes, verticalalignment="top")
    plt.hist(clone_sizes.to_flat_array(), bins=np.arange(2, max_clone_size + 1) - 0.5, color="black")
    if n == 5:
        plt.xlabel(f'clone size for {clone_size_duration}h')
        plt.xticks(range(2, max_clone_size + 2, 2))
    else:
        plt.xticks([])
    if n == 3:
        plt.ylabel("count")

    # print sister and cousin statistics
    stats = lineages.count_divisions()
    print(f"STATS FOR PARAMETER SET {n}: {stats}")

plt.show()