<jupyter_start><jupyter_text>Figure 8:  Sampling cell observations from 'true' process time In this notebook, we simulate a system with the bifurcating topology $AB(C/D)$, plot the expression over observed process times, and visualize the corresponding 'pseudotime' assignments for each cell.  Dependencies<jupyter_code>import time
t1 = time.time()
## uncomment this cell if run in colab ##

%%capture
!pip install --upgrade git+https://github.com/meichenfang/velocyto.py.git
!pip install scanpy umap-learn


!git clone https://github.com/pachterlab/GFCP_2022.git
%cd GFCP_2022
from vis import * 

## run locally ##
# %cd ../../GFCP_2021
# from vis import * 
# %cd notebooks
import warnings
warnings.filterwarnings("ignore")

from os.path import exists

plt.rcParams.update({'font.size': 22})<jupyter_output><empty_output><jupyter_text>Simulations Generate or load simulation data<jupyter_code>!mkdir -p data
nCells,nGenes,T,tau,topo =  2000,100,10,[0,3,7],"ab(c/d)"
topo_ = 'abcd'

meta = nCells,nGenes,T,tau,topo
filename="data/simdata_"+topo_

tsim_ = time.time()

from numpy.random import MT19937
from numpy.random import RandomState, SeedSequence
rs = RandomState(MT19937(SeedSequence(2812)))


szfig = [1,1]
N = np.prod(szfig)
for i_ in range(N):
    simdata = simulate_occup_meas(nCells,nGenes,T,tau,topo,seed=18)
    sim_export(simdata,meta,filename+str(i_))

dtime = time.time()-tsim_
print('Sim runtime: {:.1f} seconds.'.format(dtime))<jupyter_output>Sim runtime: 18.2 seconds.<jupyter_text>Plot expression over observed time<jupyter_code>import seaborn as sns
import matplotlib
%matplotlib inline
def getXtheo(vlm,geneind,meta):
  nCells,nGenes,T,tau,topo = meta
  j=geneind

  tvec = vlm.ca['time']
  gamma = vlm.ra['gamma']
  beta = vlm.ra['beta']
  dict_nk = {"aba":2,"ab(a/c)":3,"ab(c/d)":4}
  dict_Kval = {"aba":[0],"ab(a/c)":[0,2],"ab(c/d)":[2,3]}
  n_K = dict_nk[topo]
  K = np.vstack([vlm.ra['k'+str(i)] for i in range(n_K)]).T

  branches = dict_Kval[topo]
  n_branches = len(branches)

  i=0
  Xtheo = [eval_x(get_cell_spec_K(K[j],branches[i]),tau,t,beta[j],gamma[j]) for t in tvec]
  Xtheo = np.asarray(Xtheo)
  Xtheo[np.isnan(Xtheo)]=0
  Xtheo[np.isinf(Xtheo)]=0

  return Xtheo ,tvec
def addJitter(Xtheo,celltype,ind=1):
  types = np.unique(celltype)
  for i in types:
    if i < 3:
      Xtheo[celltype == i,ind] = Xtheo[celltype == i,ind] + np.random.normal(0,1,Xtheo[celltype == i,ind].shape[0])
    else:
      Xtheo[celltype == i,ind] = Xtheo[celltype == i,ind] + np.random.normal(-1,1,Xtheo[celltype == i,ind].shape[0])

  return Xtheo
#Get X theoretical (not imputed values by velocity)
vlm = vcy.VelocytoLoom(loom_filepath=filename+str(0)+'.loom')
tvec = vlm.ca['time']
tvec
geneind = 8
#for i in geneind:

from matplotlib.ticker import (MultipleLocator, FormatStrFormatter,
                               AutoMinorLocator)
import random
from matplotlib import cm

X, t = getXtheo(vlm,geneind,meta)
#Add jitter
X[:,1] = X[:,1] + np.random.normal(0,0.5,X.shape[0])
X[vlm.ca['celltype'] == 3,1] = -1*X[vlm.ca['celltype'] == 3,1]


fig, ax = plt.subplots(nrows = 1, ncols = 1, figsize=(8,3))
names = ['Cell Group '+str(i) for i in vlm.ca['celltype']]
sns.scatterplot(t,X[:,1] + 5,hue=pd.Categorical(names),alpha=0.7,palette=sns.color_palette("rocket")[::-1][0:4])

# ---- Figure aesthetics ---
handles, labels = plt.gca().get_legend_handles_labels()
by_label = dict(zip(labels, handles))
ax.legend(by_label.values(), by_label.keys(),loc='center left',bbox_to_anchor=(1, 0.5),prop={'size': 10},frameon=False,ncol=1)

# Hide the right and top spines
ax.spines['right'].set_visible(False)
ax.spines['top'].set_visible(False)

ax.set_xlabel('Sampled Process Times (t)',fontsize=12)
ax.set_ylabel('Expression (a.u.)',fontsize=12)


samp = random.sample(list(t),int(len(t)/5))

ax.xaxis.set_minor_locator(matplotlib.ticker.FixedLocator(samp))

# We change the fontsize of minor ticks label 
ax.tick_params(axis='both', which='major', labelsize=12)
ax.tick_params(axis='both', which='minor', size=3, direction='in')

fig.tight_layout()
plt.savefig('observedTimes.pdf')<jupyter_output><empty_output><jupyter_text>Make pseudotime image for cells<jupyter_code>#Find min time and 'root' cell

m = X[t==np.min(t),1]

#Get euclidean distance of all cells from that cell
dists = (X[:,1]-m)**2
dists.shape
dists = (dists-np.min(dists))/np.max(dists) #0 to 1 like pseudotime
fig, ax = plt.subplots(nrows = 1, ncols = 1, figsize=(8,1.5))
names = ['Cell Group '+str(i) for i in vlm.ca['celltype']]
#ys = [0]*len(names) + np.random.normal(0,0.05,len(names))
sns.scatterplot(dists ,
                vlm.ca['celltype']*-1,hue=pd.Categorical(names),
                alpha=0.7,palette=sns.color_palette("rocket")[::-1][0:4],
                s=70) #np.random.normal(0,0.02,len(names))**2


# ---- Figure aesthetics ---

handles, labels = plt.gca().get_legend_handles_labels()
by_label = dict(zip(labels, handles))
ax.legend(by_label.values(), by_label.keys(),loc='center left',bbox_to_anchor=(1, 0.5),prop={'size': 10},frameon=False,ncol=1)

# Hide the right and top spines
ax.spines['right'].set_visible(False)
ax.spines['top'].set_visible(False)

ax.set_xlabel('Pseudotime',fontsize=12)


#fig.tight_layout()
plt.savefig('pseudotime.pdf')
dists.shape<jupyter_output><empty_output><jupyter_text>Color cells by pseudotime<jupyter_code>#for i in geneind:


fig, ax = plt.subplots(nrows = 1, ncols = 1, figsize=(6.5,3))

sns.scatterplot(t,X[:,1] + 5,hue=dists,alpha=0.7,palette=sns.color_palette("viridis",as_cmap=True),legend=False)

# ---- Figure aesthetics ---
#handles, labels = plt.gca().get_legend_handles_labels()
#by_label = dict(zip(labels, handles))
#ax.legend(by_label.values(), by_label.keys(),loc='center left',bbox_to_anchor=(1, 0.5),prop={'size': 10},frameon=False,ncol=1)

# Hide the right and top spines
ax.spines['right'].set_visible(False)
ax.spines['top'].set_visible(False)

ax.set_xlabel('Sampled Process Times (t)',fontsize=12)
ax.set_ylabel('Expression (a.u.)',fontsize=12)


samp = random.sample(list(t),int(len(t)/5))

ax.xaxis.set_minor_locator(matplotlib.ticker.FixedLocator(samp))

# We change the fontsize of minor ticks label 
ax.tick_params(axis='both', which='major', labelsize=12)
ax.tick_params(axis='both', which='minor', size=3, direction='in')

fig.tight_layout()
plt.savefig('observedTimeswithPseudo.pdf')
def plot_color_gradients(cmap_category, cmap_list):
    # Create figure and adjust figure height to number of colormaps
    nrows = len(cmap_list)
    figh = 0.35 + 0.15 + (nrows + (nrows-1)*0.1)*0.22
    fig, axes = plt.subplots(nrows=nrows, figsize=(6.4, figh))
    fig.subplots_adjust(top=1-.35/figh, bottom=.15/figh, left=0.2, right=0.99)

    axes.set_title(cmap_category + ' colormaps', fontsize=14)

    #for ax, name in zip(axes, cmap_list):
    axes.imshow(gradient, aspect='auto', cmap=plt.get_cmap(cmap_list[0]))
    axes.text(-.01, .5, cmap_list[0], va='center', ha='right', fontsize=10, transform=axes.transAxes)

    # Turn off *all* ticks & spines, not just the ones with colormaps.
    #for ax in axes:
    axes.set_axis_off()

gradient = np.linspace(0, 1, 256)
gradient = np.vstack((gradient, gradient))

plot_color_gradients('Perceptually Uniform Sequential', ['viridis'])
plt.savefig('viridis.pdf')<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.17258396161754627, 'fraction_numerical': 0.02796435915010281, 'mean_word_length': 5.520107238605898, 'pattern_counts': {'":': 6, '<': 21, '<?xml version=': 0, '>': 20, 'https://': 2, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '5', 'dir': 'jupyter-structured-clean-dedup', 'id': '382255', 'license': 'BSD-2-Clause', 'n_tokens_mistral': 3077, 'n_tokens_neox': 2894, 'n_words': 531, 'path': 'notebooks/abcd_sim.ipynb', 'repo_name': 'pachterlab/GFCP_2022', 'size': '7295'}