Simulate single cell lineage tracing data using scPhyloX

1. Simulation of tissue overshoot model

Import necessary packages

import scPhyloX as spx
from scipy.integrate import solve_ivp
import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import poisson
from Bio import Phylo
from io import StringIO
plt.rcParams['font.size'] = 12
plt.rcParams['pdf.fonttype'] = 42

WARNING (pytensor.tensor.blas): Using NumPy C-API based implementation for BLAS functions.

Specifies the phylodynamics parameters, derive analytical solution of developmental dynamics.

d, p, r, a, b, k, t0 = (0.01, 0.6, 0.4, 0.9, 0.1, 0.3, 12)
T = 35
t = range(T)
x0 = [100, 0]

sol = solve_ivp(spx.sim_tissue.cellnumber, t_span=(0, T), y0=x0, t_eval=range(T), method='RK45', args=(a, b, k, t0, p, r, d))
fig, ax = plt.subplots(1, 2, figsize=(6,3.2))
ax[0].plot(sol.t, sol.y[0], label='Stem cells', c='#9098d9', lw=4)
ax[0].plot(sol.t, sol.y[1], label='Non-stem cells', c='#ed9e44', lw=4)
ax[0].legend(loc=2)
ax[1].plot(spx.sim_tissue.bt(np.array(t), a, b, k, t0), lw=3, c='tab:gray')
ax[0].set_xlabel('Time')
ax[1].set_xlabel('Time')
ax[0].set_ylabel('Cell number')
ax[1].set_ylabel(r'$\beta(t)$')
ax[1].set_ylim([0,1.05])
# ax[0].ticklabel_format(style='sci', scilimits= (-1,2), axis='y', useMathText=True)
ax[1].spines['right'].set_visible(False)
ax[1].spines['top'].set_visible(False)
ax[0].spines['right'].set_visible(False)
ax[0].spines['top'].set_visible(False)
ax[0].set_title('Analytical\ndevelopmental dynamics')
ax[1].set_title('Stem cell\nself-renewal prob.')
plt.tight_layout()

Specify DNA mutation rates and perform simulation

mu = 2
system = spx.sim_tissue.simulation(x0, T, mu, a, b, p, r, k, d, t0)

cell_num:20813, time:35.002134948064595

Comparison between simulated and analytical developmental dynamics

cell_number = np.array(system.n)
fig, ax = plt.subplots(figsize=(3.8, 3.2))
show_tp = np.arange(0, len(system.t),1000)
ax.plot(np.array(system.t)[show_tp], cell_number[show_tp, 0], '--', lw=3, label='Stem\nsimulated', c='#9098d9')
ax.plot(np.array(system.t)[show_tp], cell_number[show_tp, 1], '--', lw=3, label='Non-stem\nsimulated', c='#ed9e44')
ax.plot(sol.t, sol.y[0], label='Stem', c='#9098d9', lw=6, alpha=0.4)
ax.plot(sol.t, sol.y[1], label='Non-stem', c='#ed9e44', lw=6, alpha=0.4)
ax.legend(loc=2,fontsize=10)
ax.ticklabel_format (style='sci', scilimits= (-1,2), axis='y')
ax.spines['right'].set_visible(False)
ax.spines['top'].set_visible(False)
ax.set_xlabel('Time')
ax.set_ylabel('Cell number')

Text(0, 0.5, 'Cell number')

Derive mutational character matrix from simulated data

seqtab = np.array([i.seq for i in system.Stemcells] + [i.seq for i in system.Diffcells])
cell_names = np.array([f'<{i.gen}_{i.cellid}>' for i in system.Stemcells] + [f'<{i.gen}_{i.cellid}>' for i in system.Diffcells])
sel_cells = np.random.choice(range(seqtab.shape[0]), 500, replace=False)
seqtab = seqtab[sel_cells]
cell_names = cell_names[sel_cells]

Visualize phylogenetic tree

tree_nwk = spx.utils.reconstruct(system.lineage_info, list(cell_names))
tree = Phylo.read(StringIO(tree_nwk), format='newick')

fig, ax = plt.subplots(figsize=(3.5, 3.5), subplot_kw={'projection': 'polar'})

spx.tree.polar_plot(tree, ax=ax, label_leaf=False, lw=0.6)

<PolarAxes: >

2. Simulation of tissue continuous growth model

d, p, r, a, b, k, t0 = (0.2, 0.6, 1.3, 0.2, 0.375, 0.8, 8)
K = 35
t = range(K)
x0 = [100, 0]

sol = solve_ivp(spx.sim_tissue.cellnumber, t_span=(0, K), y0=x0, t_eval=range(K+1),
                method='RK45', args=(a, b, k, t0, p, r, d))
fig, ax = plt.subplots(1, 2, figsize=(7,3.2))
ax[0].plot(sol.t, sol.y[0], label='Stem cells', c='#9098d9', lw=4)
ax[0].plot(sol.t, sol.y[1], label='Non-stem cells', c='#ed9e44', lw=4)
ax[0].legend(loc=2)
ax[1].plot(spx.sim_tissue.bt(np.array(t), a, b, k, t0), lw=3, c='tab:gray')
ax[0].set_xlabel('Time')
ax[1].set_xlabel('Time')
ax[0].set_ylabel('Cell number')
ax[1].set_ylabel(r'$\beta(t)$')
ax[1].set_ylim([0,1.05])

ax[0].ticklabel_format (style='sci', scilimits= (-1,2), axis='y')
plt.tight_layout()
ax[1].spines['right'].set_visible(False)
ax[1].spines['top'].set_visible(False)
ax[0].spines['right'].set_visible(False)
ax[0].spines['top'].set_visible(False)

mu = 2
system = spx.sim_tissue.simulation(x0, K, mu, a, b, p, r, k, d, t0)

cell_num:19419, time:35.0001496265336741

cell_number = np.array(system.n)
fig, ax = plt.subplots(figsize=(3.8, 3.2))
show_tp = np.arange(0, len(system.t),1000)
ax.plot(np.array(system.t)[np.arange(0, len(system.t),1000)], cell_number[np.arange(0, len(system.t),1000), 0], '--', lw=3, label='Stem\nsimulated', c='#9098d9')
ax.plot(np.array(system.t)[np.arange(0, len(system.t),1000)], cell_number[np.arange(0, len(system.t),1000), 1], '--', lw=3, label='Non-stem\nsimulated', c='#ed9e44')
ax.plot(sol.t, sol.y[0], label='Stem', c='#9098d9', lw=6, alpha=0.4)
ax.plot(sol.t, sol.y[1], label='Non-stem', c='#ed9e44', lw=6, alpha=0.4)
ax.legend(loc=2,fontsize=10)
ax.ticklabel_format (style='sci', scilimits= (-1,2), axis='y')
ax.spines['right'].set_visible(False)
ax.spines['top'].set_visible(False)
ax.set_xlabel('Time')
ax.set_ylabel('Cell number')

Text(0, 0.5, 'Cell number')

seqtab = np.array([i.seq for i in system.Stemcells] + [i.seq for i in system.Diffcells])
cell_names = np.array([f'<{i.gen}_{i.cellid}>' for i in system.Stemcells] + [f'<{i.gen}_{i.cellid}>' for i in system.Diffcells])
sel_cells = np.random.choice(range(seqtab.shape[0]), 500, replace=False)
seqtab = seqtab[sel_cells]
cell_names = cell_names[sel_cells]

tree_nwk = spx.utils.reconstruct(system.lineage_info, list(cell_names))
tree = Phylo.read(StringIO(tree_nwk), format='newick')
fig, ax = plt.subplots(figsize=(3.5, 3.5), subplot_kw={'projection': 'polar'})
spx.tree.polar_plot(tree, ax=ax, label_leaf=False, lw=0.6)

<PolarAxes: >