#!/usr/bin/env python
# coding: utf-8

# # scEU-seq organoid
# 
# This tutorial uses the intestine organoid data from Battich, et al (2020). This tutorial is the second one of the two tutorials for demonstrating how dynamo can use used to analyze the scEU-seq data. Please refer the [cell cycle](https://dynamo-release.readthedocs.io/en/latest/scEU_seq_rpe1_analysis_kinetic.html) tutorial for details on how to analyze the cell cycle dataset.
# 

# In[2]:


import warnings
warnings.filterwarnings('ignore')

import dynamo as dyn
import anndata
import pandas as pd
import numpy as np
import scipy.sparse

from anndata import AnnData
from scipy.sparse import csr_matrix

# dyn.get_all_dependencies_version()


# ## Load data

# In[3]:


organoid = dyn.sample_data.scEU_seq_organoid()


# In[4]:


organoid


# In[5]:


# mapping:
cell_mapper = {
    '1': 'Enterocytes',
    '2': 'Enterocytes',
    '3': 'Enteroendocrine',
    '4': 'Enteroendocrine progenitor',
    '5': 'Tuft cells',
    '6': 'TA cells',
    '7': 'TA cells',
    '8': 'Stem cells',
    '9': 'Paneth cells',
    '10': 'Goblet cells',
    '11': 'Stem cells',
 }

organoid.obs['cell_type'] = organoid.obs.som_cluster_id.map(cell_mapper).astype('str')


# ## typical dynamo analysis workflow
# 

# In[6]:


dyn.pl.basic_stats(organoid)


# In[7]:


organoid


# In[8]:


organoid.obs


# In[9]:


organoid.obs.groupby(['exp_type', 'time']).agg('count')


# In[10]:


adata = organoid.copy()
adata.obs.time = adata.obs.time.astype('str')
adata.obs.loc[adata.obs['time'] == 'dmso', 'time'] = -1
adata.obs['time'] = adata.obs['time'].astype(float)
adata = adata[adata.obs.time != -1, :]
adata = adata[adata.obs.exp_type == 'Pulse', :]
adata.layers['new'], adata.layers['total'] = adata.layers['ul'] + adata.layers['sl'], adata.layers['su'] + adata.layers['sl'] + adata.layers['uu'] + adata.layers['ul']
del adata.layers['uu'], adata.layers['ul'], adata.layers['su'], adata.layers['sl']
dyn.pp.recipe_monocle(adata, n_top_genes=1000, total_layers=False)
# preprocessor = dyn.pp.Preprocessor(cell_cycle_score_enable=True)
# preprocessor.config_monocle_recipe(adata, n_top_genes=1000)                                  
# preprocessor.preprocess_adata_monocle(adata)
dyn.pl.basic_stats(adata)
dyn.pl.show_fraction(organoid)


# In[11]:


adata.obs.time = adata.obs.time/60


# In[12]:


adata.obs.time  = adata.obs.time.astype('float')
dyn.tl.dynamics(adata, model='deterministic', tkey='time', assumption_mRNA='ss')

dyn.tl.reduceDimension(adata)


# In[13]:


dyn.tl.cell_velocities(adata, ekey='M_t', vkey='velocity_T', enforce=True)


# In[14]:


adata.obsm['X_umap_ori'] = adata.obs.loc[:, ['rotated_umap1', 'rotated_umap2']].values.astype(float)


# ## Visualize time-resolved vector flow learned with dynamo 

# In[ ]:


dyn.tl.cell_velocities(adata, basis='umap_ori')

dyn.pl.streamline_plot(adata, color='cell_type', basis='umap_ori')


# In[ ]:


dyn.pl.streamline_plot(adata, color='cell_cycle_phase', basis='umap_ori')


# In[17]:


adata.var_names[adata.var.use_for_transition][:5]


# In[ ]:


dyn.pl.phase_portraits(adata, genes=['Brat1', 'Ccnd2', 'Ckmt1', 'Pdgfb', 'Gpa33'],
                       color='som_cluster_id', basis='umap_ori')


# ## Animate intestine organoid differentiation

# In[42]:


dyn.vf.VectorField(adata, basis='umap_ori')


# In[43]:


progenitor = adata.obs_names[adata.obs.cell_type == 'Stem cells']
len(progenitor)


# In[44]:


np.random.seed(19491001)

from matplotlib import animation
info_genes = adata.var_names[adata.var.use_for_transition]
dyn.pd.fate(adata, basis='umap_ori', init_cells=progenitor[:100], interpolation_num=100,  direction='forward',
   inverse_transform=False, average=False)


# In[ ]:


get_ipython().run_cell_magic('capture', '', "import matplotlib.pyplot as plt\n\nfig, ax = plt.subplots()\nax = dyn.pl.topography(adata, basis='umap_ori', color='cell_type', ax=ax, save_show_or_return='return',  figsize=(24, 24))\nax.set_aspect(0.8)\n")


# In[ ]:


get_ipython().run_cell_magic('capture', '', "adata.obs['time'] = adata.obs.time.astype('float')\ninstance = dyn.mv.StreamFuncAnim(adata=adata, basis='umap_ori', color='cell_type', ax=ax)\n")


# In[ ]: