

In [20]:

    
import scipy as sp
from scipy import stats
import pandas as pd
import matplotlib as mpl
import matplotlib.pyplot as plt
%matplotlib inline
%config InlineBackend.figure_format = 'svg'
exec(open('../settings.py').read(), globals())



In [21]:

    
store = pd.HDFStore('./lg_model_plot_data.hdf')



In [22]:

    
full_title = {'full': 'Full model',
              'basal_div': 'No acceleration',
              'nodiv': 'No division',
              'noact': 'No activation',
              'noflux': 'No influx',
              'noflux_noact': 'Only divisions'}

Plot full model



In [23]:

    
scale_growth = 1.0 / 1000.0



In [26]:

    
time = sp.array(store['time'])
gexpdata = store['gexpdata']
outgrowth_pop = store['outgrowth_pop']

condition = 'full'
  
fig, ax = plt.subplots(1, 2, sharex = True, figsize = (80.0/25.4, 52.0/25.4))
# fig.suptitle(full_title[condition], fontsize = 8)

# experimental data
ax[1].errorbar(sp.array(outgrowth_pop['time']), scale_growth * sp.array(outgrowth_pop['mean']),
               scale_growth * sp.array(outgrowth_pop['sem']),
               label = 'Exp.', fmt = 'o', color = 'black')

ax[0].errorbar(gexpdata['time'], gexpdata['G', 'mean'], gexpdata['G', 'sem'], fmt='o', color = 'black', label = 'Exp.')

# best prediction

L = store['best_fit_{}'.format(condition)]['L']
G = store['best_fit_{}'.format(condition)]['G']

if condition == 'full':
    Glabel = 'fit'
else:
    Glabel = 'prediction'

# ax[1].plot(time, scale_growth * L, lw = 2, color = 'black', label = 'prediction')
# ax[0].plot(time, G,                  lw = 2, color = 'black', label = Glabel)

# 3-sigma confidence interval
L_3sigma_max = store['L_3sigma_max_{0}'.format(condition)]
L_3sigma_min = store['L_3sigma_min_{0}'.format(condition)]
ax[1].fill_between(time, scale_growth * L_3sigma_min, scale_growth * L_3sigma_max, alpha = 0.2, color = 'green', label = 'Prediction')

G_3sigma_max = store['G_3sigma_max_{0}'.format(condition)]
G_3sigma_min = store['G_3sigma_min_{0}'.format(condition)]
ax[0].fill_between(time, G_3sigma_min, G_3sigma_max, alpha = 0.2, color = 'black', label = 'Fit')

# 2sigma confidence interval
L_2sigma_max = store['L_2sigma_max_{0}'.format(condition)]
L_2sigma_min = store['L_2sigma_min_{0}'.format(condition)]
ax[1].fill_between(time, scale_growth * L_2sigma_min, scale_growth * L_2sigma_max, alpha = 0.2, color = 'green')

G_2sigma_max = store['G_2sigma_max_{0}'.format(condition)]
G_2sigma_min = store['G_2sigma_min_{0}'.format(condition)]
ax[0].fill_between(time, G_2sigma_min, G_2sigma_max, alpha = 0.2, color = 'black')

# 1sigma confidence interval
L_1sigma_max = store['L_1sigma_max_{0}'.format(condition)]
L_1sigma_min = store['L_1sigma_min_{0}'.format(condition)]
ax[1].fill_between(time, scale_growth * L_1sigma_min, scale_growth * L_1sigma_max, alpha = 0.2, color = 'green')

G_1sigma_max = store['G_1sigma_max_{0}'.format(condition)]
G_1sigma_min = store['G_1sigma_min_{0}'.format(condition)]
ax[0].fill_between(time, G_1sigma_min, G_1sigma_max, alpha = 0.2, color = 'black')

# formatting
fig.patch.set_alpha(1.0)

ax[1].legend(loc = 'upper left', frameon = False, numpoints = 1)
#     ax[1].plot([0.6, 2.45], 2 * [3.37], color = 'black')
ax[1].set_xlim(-0.5, 8.5)
ax[1].set_ylim(0.0, 4.0)
ax[1].set_yticks([0, 1, 2, 3, 4])

ax[0].legend(loc = (0.4, 0.05), frameon = False, numpoints = 1)
ax[0].set_yticks(sp.arange(0.7, 1.01, 0.1))
ax[0].set_yticks(sp.arange(0.75, 1.06, 0.1), minor = True)
ax[0].set_ylim(0.79, 1.01)

#     ax[1].get_yaxis().set_label_coords(-0.3, 0.5)
#     ax[0].get_yaxis().set_label_coords(-0.3, 0.5)


ax[0].set_xlabel('Time (days)')
ax[1].set_xlabel('Time (days)')
ax[1].set_ylabel('Outgrowth (mm)')
ax[0].set_ylabel('Growth fraction')
plt.tight_layout()
plt.savefig('../../figure_plots/lg_model_{}.svg'.format(condition))

plt.show()

Full model without data



In [6]:

    
time = sp.array(store['time'])
gexpdata = store['gexpdata']
outgrowth_pop = store['outgrowth_pop']

condition = 'full'
  
fig, ax = plt.subplots(1, 2, sharex = True, figsize = (80.0/25.4, 40.0/25.4))
# fig.suptitle(full_title[condition], fontsize = 8)

# experimental data
# ax[1].errorbar(sp.array(outgrowth_pop['time']), scale_growth * sp.array(outgrowth_pop['mean']),
#                scale_growth * sp.array(outgrowth_pop['sem']),
#                label = 'data', fmt = 'o', color = 'black')

# ax[0].errorbar(gexpdata['time'], gexpdata['G', 'mean'], gexpdata['G', 'sem'], fmt='o', color = 'black', label = 'data')

# best prediction

L = store['best_fit_{}'.format(condition)]['L']
G = store['best_fit_{}'.format(condition)]['G']

if condition == 'full':
    Glabel = 'fit'
else:
    Glabel = 'prediction'

# ax[1].plot(time, scale_growth * L, lw = 2, color = 'black', label = 'prediction')
# ax[0].plot(time, G,                  lw = 2, color = 'black', label = Glabel)

# 3-sigma confidence interval
L_3sigma_max = store['L_3sigma_max_{0}'.format(condition)]
L_3sigma_min = store['L_3sigma_min_{0}'.format(condition)]
ax[1].fill_between(time, scale_growth * L_3sigma_min, scale_growth * L_3sigma_max, alpha = 0.2, color = 'green', label = 'prediction')

G_3sigma_max = store['G_3sigma_max_{0}'.format(condition)]
G_3sigma_min = store['G_3sigma_min_{0}'.format(condition)]
ax[0].fill_between(time, G_3sigma_min, G_3sigma_max, alpha = 0.2, color = 'black', label = 'fit')

# 2sigma confidence interval
L_2sigma_max = store['L_2sigma_max_{0}'.format(condition)]
L_2sigma_min = store['L_2sigma_min_{0}'.format(condition)]
ax[1].fill_between(time, scale_growth * L_2sigma_min, scale_growth * L_2sigma_max, alpha = 0.2, color = 'green')

G_2sigma_max = store['G_2sigma_max_{0}'.format(condition)]
G_2sigma_min = store['G_2sigma_min_{0}'.format(condition)]
ax[0].fill_between(time, G_2sigma_min, G_2sigma_max, alpha = 0.2, color = 'black')

# 1sigma confidence interval
L_1sigma_max = store['L_1sigma_max_{0}'.format(condition)]
L_1sigma_min = store['L_1sigma_min_{0}'.format(condition)]
ax[1].fill_between(time, scale_growth * L_1sigma_min, scale_growth * L_1sigma_max, alpha = 0.2, color = 'green')

G_1sigma_max = store['G_1sigma_max_{0}'.format(condition)]
G_1sigma_min = store['G_1sigma_min_{0}'.format(condition)]
ax[0].fill_between(time, G_1sigma_min, G_1sigma_max, alpha = 0.2, color = 'black')

# formatting
fig.patch.set_alpha(1.0)

ax[1].legend(loc = 'upper left', frameon = False, numpoints = 1)
#     ax[1].plot([0.6, 2.45], 2 * [3.37], color = 'black')
ax[1].set_xlim(-0.5, 8.5)
ax[1].set_ylim(0.0, 4.0)
ax[1].set_yticks([0, 1, 2, 3, 4])

ax[0].legend(loc = (0.4, 0.05), frameon = False, numpoints = 1)
ax[0].set_yticks(sp.arange(0.7, 1.01, 0.1))
ax[0].set_yticks(sp.arange(0.75, 1.06, 0.1), minor = True)
ax[0].set_ylim(0.79, 1.01)

#     ax[1].get_yaxis().set_label_coords(-0.3, 0.5)
#     ax[0].get_yaxis().set_label_coords(-0.3, 0.5)


ax[0].set_xlabel('time (days)')
ax[1].set_xlabel('time (days)')
ax[1].set_ylabel('outgrowth (mm)')
ax[0].set_ylabel('growth fraction')
plt.tight_layout()
plt.savefig('../../figure_plots/lg_model_without_data_{}.svg'.format(condition))

plt.show()

Plot scenarios



In [7]:

    
time = sp.array(store['time'])
gexpdata = store['gexpdata']
outgrowth_pop = store['outgrowth_pop']
conditions = list(store['conditions'])
conditions.remove('full')


for condition in conditions:

  
    fig, ax = plt.subplots(2, sharex = True, figsize = (33.0/25.4, 65.0/25.4))
    fig.suptitle(full_title[condition], fontsize = 8)

    # experimental data
    ax[1].errorbar(sp.array(outgrowth_pop['time']), scale_growth * sp.array(outgrowth_pop['mean']),
                   scale_growth * sp.array(outgrowth_pop['sem']),
                   label = '_data', fmt = 'o', color = 'black')
    
    ax[0].errorbar(gexpdata['time'], gexpdata['G', 'mean'], gexpdata['G', 'sem'], fmt='o', color = 'black', label = '_data')
    
#     # best prediction
    
#     L = store['best_fit_{}'.format(condition)]['L']
#     G = store['best_fit_{}'.format(condition)]['G']
    
#     if condition == 'full':
#         Glabel = 'fit'
#     else:
#         Glabel = 'prediction'
        
#     ax[1].plot(time, scale_growth * L, lw = 2, color = 'black', label = 'prediction')
#     ax[0].plot(time, G,                  lw = 2, color = 'black', label = Glabel)
        
    # 3-sigma confidence interval
    L_3sigma_max = store['L_3sigma_max_{0}'.format(condition)]
    L_3sigma_min = store['L_3sigma_min_{0}'.format(condition)]
    ax[1].fill_between(time, scale_growth * L_3sigma_min, scale_growth * L_3sigma_max, alpha = 0.2, color = 'green')
    
    G_3sigma_max = store['G_3sigma_max_{0}'.format(condition)]
    G_3sigma_min = store['G_3sigma_min_{0}'.format(condition)]
    ax[0].fill_between(time, G_3sigma_min, G_3sigma_max, alpha = 0.2, color = 'green')
    
    # 2sigma confidence interval
    L_2sigma_max = store['L_2sigma_max_{0}'.format(condition)]
    L_2sigma_min = store['L_2sigma_min_{0}'.format(condition)]
    ax[1].fill_between(time, scale_growth * L_2sigma_min, scale_growth * L_2sigma_max, alpha = 0.2, color = 'green')
    
    G_2sigma_max = store['G_2sigma_max_{0}'.format(condition)]
    G_2sigma_min = store['G_2sigma_min_{0}'.format(condition)]
    ax[0].fill_between(time, G_2sigma_min, G_2sigma_max, alpha = 0.2, color = 'green')
    
    # 1sigma confidence interval
    L_1sigma_max = store['L_1sigma_max_{0}'.format(condition)]
    L_1sigma_min = store['L_1sigma_min_{0}'.format(condition)]
    ax[1].fill_between(time, scale_growth * L_1sigma_min, scale_growth * L_1sigma_max, alpha = 0.2, color = 'green')
    
    G_1sigma_max = store['G_1sigma_max_{0}'.format(condition)]
    G_1sigma_min = store['G_1sigma_min_{0}'.format(condition)]
    ax[0].fill_between(time, G_1sigma_min, G_1sigma_max, alpha = 0.2, color = 'green')
    
    # formatting
    fig.patch.set_alpha(1.0)
    
#     ax[1].legend(loc = 'best', frameon = False, numpoints = 1)
#     ax[1].plot([0.6, 2.45], 2 * [3.37], color = 'black')
    ax[1].set_xlim(-0.5, 8.5)
    ax[1].set_ylim(0.0, 4.0)
    ax[1].set_yticks([0, 1, 2, 3, 4])
    
#     ax[0].legend(loc = (0.15,0.05), frameon = False, numpoints = 1)
    ax[0].set_yticks(sp.arange(0.7, 1.01, 0.1))
    ax[0].set_yticks(sp.arange(0.75, 1.06, 0.1), minor = True)
    ax[0].set_ylim(0.79, 1.01)
    
#     ax[1].get_yaxis().set_label_coords(-0.3, 0.5)
#     ax[0].get_yaxis().set_label_coords(-0.3, 0.5)
    
   
    ax[1].set_xlabel('Time (days)')
    ax[1].set_ylabel('Outgrowth (mm)')
    ax[0].set_ylabel('Growth fraction')
    
    plt.savefig('../../figure_plots/lg_model_{}.svg'.format(condition))

    plt.show()

Maximum outgrowth of the no acceleration scenario (number used in the results)



In [8]:

    
store['L_3sigma_max_basal_div'][8] / 1000.0









    Out[8]:





1.6916068061672069

Compare with Fei et al., 2014

Fei, Ji-Feng, Maritta Schuez, Akira Tazaki, Yuka Taniguchi, Kathleen Roensch, and Elly M. Tanaka. “CRISPR-Mediated Genomic Deletion of Sox2 in the Axolotl Shows a Requirement in Spinal Cord Neural Stem Cell Amplification during Tail Regeneration.” Stem Cell Reports 3, no. 3 (September 9, 2014): 444–59. doi:10.1016/j.stemcr.2014.06.018.



In [9]:

    
outgrowth_fei = pd.read_excel('../../data/Fei_2014_data.xlsx').iloc[:,0:3]
outgrowth_fei.columns = ['condition', 'ID', 'length']
outgrowth_fei['time'] = 6.0
outgrowth_fei.head()









    Out[9]:






  
    
      
      condition
      ID
      length
      time
    
  
  
    
      0
      GFP-gRNA#3
      #1
      616.25
      6.0
    
    
      1
      GFP-gRNA#3
      #2
      706.17
      6.0
    
    
      2
      GFP-gRNA#3
      #3
      648.37
      6.0
    
    
      3
      GFP-gRNA#3
      #4
      617.37
      6.0
    
    
      4
      GFP-gRNA#3
      #5
      637.19
      6.0



In [10]:

    
mean_outgrowth_fei = outgrowth_fei.groupby('condition').agg({'time': ['mean'], 'length': ['mean', 'sem', 'count']})



In [11]:

    
mean_outgrowth_fei









    Out[11]:






  
    
      
      length
      time
    
    
      
      mean
      sem
      count
      mean
    
    
      condition
      
      
      
      
    
  
  
    
      GFP-gRNA#3
      643.176667
      13.882569
      12
      6.0
    
    
      Sox2-gRNA#2
      504.412174
      22.598582
      23
      6.0
    
    
      Sox2-gRNA#4
      438.917917
      23.417426
      24
      6.0

Correct the SOX2 outgrowth data for the density

todo: calculate the correction factor



In [12]:

    
density_fei = pd.read_excel('../../data/Sox2_Crispr2_with Control_cell_type_quanti_6dpa-osvaldo-2.xlsx')
density_fei



In [13]:

    
density_ctr = density_fei.loc[:'Ctr_4']
density_ctr



In [14]:

    
density_SOX2ko = density_fei['SOX2-1':]
density_SOX2ko



In [15]:

    
mean_density_ctr = density_ctr.mean(axis=1).mean()
mean_density_SOX2ko = density_SOX2ko.mean(axis=1).mean()



In [16]:

    
density_correction = mean_density_SOX2ko / mean_density_ctr
density_correction









    Out[16]:





0.7058338118632893

apply correction



In [17]:

    
mean_outgrowth_fei.loc[['Sox2-gRNA#2', 'Sox2-gRNA#4'], 'length'] = density_correction * mean_outgrowth_fei.loc[['Sox2-gRNA#2', 'Sox2-gRNA#4']]
mean_outgrowth_fei









    Out[17]:






  
    
      
      length
      time
    
    
      
      mean
      sem
      count
      mean
    
    
      condition
      
      
      
      
    
  
  
    
      GFP-gRNA#3
      643.176667
      13.882569
      12.000000
      6.0
    
    
      Sox2-gRNA#2
      356.031167
      15.950843
      16.234178
      6.0
    
    
      Sox2-gRNA#4
      309.803106
      16.528811
      16.940011
      6.0

Rescale outgrowth from Fei such that their control matches our case



In [18]:

    
rescale_factor = outgrowth_pop.query('time == 6')['mean'].iloc[0] / mean_outgrowth_fei.loc['GFP-gRNA#3', ('length', 'mean')]
rescale_factor









    Out[18]:





1.9877897726389329

Plot the data



In [19]:

    
time = sp.array(store['time'])
gexpdata = store['gexpdata']
outgrowth_pop = store['outgrowth_pop']

condition = 'full'
  
fig, ax = plt.subplots(1, 2, sharex = True, sharey = True, figsize = (80.0/25.4, 52.0/25.4))

# Full model
# Fei
ax[0].errorbar(mean_outgrowth_fei.loc['GFP-gRNA#3', ('time', 'mean')], 
               scale_growth * rescale_factor * mean_outgrowth_fei.loc['GFP-gRNA#3', ('length', 'mean')], 
               scale_growth * mean_outgrowth_fei.loc['GFP-gRNA#3', ('length', 'sem')],
               label = 'Fei rescaled', fmt = 'o', color = 'black')


# best prediction

L = store['best_fit_{}'.format(condition)]['L']

# 3-sigma confidence interval
L_3sigma_max = store['L_3sigma_max_{0}'.format(condition)]
L_3sigma_min = store['L_3sigma_min_{0}'.format(condition)]
ax[0].fill_between(time, scale_growth * L_3sigma_min, scale_growth * L_3sigma_max, alpha = 0.2, color = 'green', label = 'Prediction')

# 2sigma confidence interval
L_2sigma_max = store['L_2sigma_max_{0}'.format(condition)]
L_2sigma_min = store['L_2sigma_min_{0}'.format(condition)]
ax[0].fill_between(time, scale_growth * L_2sigma_min, scale_growth * L_2sigma_max, alpha = 0.2, color = 'green')

# 1sigma confidence interval
L_1sigma_max = store['L_1sigma_max_{0}'.format(condition)]
L_1sigma_min = store['L_1sigma_min_{0}'.format(condition)]
ax[0].fill_between(time, scale_growth * L_1sigma_min, scale_growth * L_1sigma_max, alpha = 0.2, color = 'green')

# No acceleration
condition = 'basal_div'
# Fei
ax[1].errorbar(mean_outgrowth_fei.loc['Sox2-gRNA#4', ('time', 'mean')], 
               scale_growth * rescale_factor * mean_outgrowth_fei.loc['Sox2-gRNA#4', ('length', 'mean')], 
               scale_growth * mean_outgrowth_fei.loc['Sox2-gRNA#4', ('length', 'sem')],
               label = 'Fei SOX2-', fmt = 'o', color = 'black')



# best prediction

L = store['best_fit_{}'.format(condition)]['L']

# 3-sigma confidence interval
L_3sigma_max = store['L_3sigma_max_{0}'.format(condition)]
L_3sigma_min = store['L_3sigma_min_{0}'.format(condition)]
ax[1].fill_between(time, scale_growth * L_3sigma_min, scale_growth * L_3sigma_max, alpha = 0.2, color = 'green', label = 'Prediction')

# 2sigma confidence interval
L_2sigma_max = store['L_2sigma_max_{0}'.format(condition)]
L_2sigma_min = store['L_2sigma_min_{0}'.format(condition)]
ax[1].fill_between(time, scale_growth * L_2sigma_min, scale_growth * L_2sigma_max, alpha = 0.2, color = 'green')

# 1sigma confidence interval
L_1sigma_max = store['L_1sigma_max_{0}'.format(condition)]
L_1sigma_min = store['L_1sigma_min_{0}'.format(condition)]
ax[1].fill_between(time, scale_growth * L_1sigma_min, scale_growth * L_1sigma_max, alpha = 0.2, color = 'green')

# formatting
fig.patch.set_alpha(1.0)

ax[0].legend(loc = 'upper left', frameon = False, numpoints = 1)
ax[1].legend(loc = 'upper left', frameon = False, numpoints = 1)
#     ax[0].plot([0.6, 2.45], 2 * [3.37], color = 'black')
ax[0].set_xlim(-0.5, 8.5)
ax[0].set_ylim(0.0, 4.0)
ax[0].set_yticks([0, 1, 2, 3, 4])

ax[0].set_xlabel('Time (days)')
ax[1].set_xlabel('Time (days)')
ax[0].set_ylabel('Outgrowth (mm)')

ax[0].set_title('Full model', fontsize = 8)
ax[1].set_title('No acceleration', fontsize = 8)

plt.tight_layout()
plt.savefig('../../figure_plots/Fei_comparisson.svg')

plt.show()

	M2	M1	R1	R2	R3	R4	R5
Ctr_1	17.000000	14.000000	11.500000	10.500000	11.000000	13.500000	9.500000
Ctr_2	18.666667	15.000000	14.666667	13.666667	18.666667	14.000000	10.000000
Ctr_3	20.000000	14.666667	13.000000	14.666667	15.000000	11.333333	8.000000
Ctr_4	17.000000	17.333333	11.000000	11.666667	12.333333	6.500000	9.000000
SOX2-1	15.333333	9.666667	9.333333	7.666667	1.333333	NaN	NaN
SOX2-2	15.000000	15.000000	9.000000	10.000000	10.000000	2.000000	NaN
SOX2-3	14.500000	15.000000	10.000000	8.000000	8.000000	2.000000	NaN
SOX2-4	16.333333	10.333333	8.333333	7.666667	6.000000	0.666667	NaN
SOX2-5	16.000000	14.666667	12.000000	7.666667	2.000000	NaN	NaN
SOX2-6	15.000000	16.666667	12.000000	8.333333	8.666667	6.333333	1.666667
SOX2-7	13.000000	11.666667	6.666667	10.333333	3.000000	NaN	NaN

	M2	M1	R1	R2	R3	R4	R5
Ctr_1	17.000000	14.000000	11.500000	10.500000	11.000000	13.500000	9.5
Ctr_2	18.666667	15.000000	14.666667	13.666667	18.666667	14.000000	10.0
Ctr_3	20.000000	14.666667	13.000000	14.666667	15.000000	11.333333	8.0
Ctr_4	17.000000	17.333333	11.000000	11.666667	12.333333	6.500000	9.0

	M2	M1	R1	R2	R3	R4	R5
SOX2-1	15.333333	9.666667	9.333333	7.666667	1.333333	NaN	NaN
SOX2-2	15.000000	15.000000	9.000000	10.000000	10.000000	2.000000	NaN
SOX2-3	14.500000	15.000000	10.000000	8.000000	8.000000	2.000000	NaN
SOX2-4	16.333333	10.333333	8.333333	7.666667	6.000000	0.666667	NaN
SOX2-5	16.000000	14.666667	12.000000	7.666667	2.000000	NaN	NaN
SOX2-6	15.000000	16.666667	12.000000	8.333333	8.666667	6.333333	1.666667
SOX2-7	13.000000	11.666667	6.666667	10.333333	3.000000	NaN	NaN

	condition	ID	length	time
0	GFP-gRNA#3	#1	616.25	6.0
1	GFP-gRNA#3	#2	706.17	6.0
2	GFP-gRNA#3	#3	648.37	6.0
3	GFP-gRNA#3	#4	617.37	6.0
4	GFP-gRNA#3	#5	637.19	6.0

	length			time
	mean	sem	count	mean
condition
GFP-gRNA#3	643.176667	13.882569	12	6.0
Sox2-gRNA#2	504.412174	22.598582	23	6.0
Sox2-gRNA#4	438.917917	23.417426	24	6.0