Getting Started
===============

We present a model to analyze populations of cells where we have one or more phenotypic measurements from, such as cell lifetime, cell fates, cell shape, migration, etc.
The data should be in the form of a lineage binary tree, where each node represents one cell with its measurements.
We have provided a way to create synthetic data in the same format to test the performance of our model.
tHMM uses expectation-maximization (EM) algorithm to cluster the cells based solely on their measurements and their relationship with other cells in a lineage.
Each cluster, aka state, represents a sub-population of cells that belong to a distribution for their observations.
For instance, considering cell lifetime as a measurement, cells that belong to state 1 have the lifetime duration that comes from a Gamma(s1, k1), and cells that belong to state2 have the lifetime duration that comes from a Gamma(s2, k2).
This way, we can quantify phenotypic heterogeneity within a population of cells by revealing each cell's hidden state. Our model, incorporating the cell-cell relationship within each lineage is superior to other clustering methods such as K-means.


We first introduce how to synthesize these populations by working up from the basic unit; cells. We then introduce how to synthesize lineages, which are just hierarchical lineage tree groupings of cells based on their family history. Our model ultimately analyzes populations which are groups of one or more lineages that share the same states.


1. Lineage Tree Data Model
---------------------------

Lineage trees in tHMM are represented at the lineage level using a **SciPy CSR sparse array** (representing parent-to-daughter connectivity) alongside **NumPy arrays** for observations, hidden states, and cell lifetimes.

Transition matrices define the probabilities of how cells divide and switch states across generations.

The transition matrix defines the rate at which cells change from one state to another. Specifically, if an element of a transition matrix :math:`T` at row :math:`i` and column :math:`j` is defined as :math:`T_{i,j}`, then

.. math:: T_{i,j} = \mathbb{P}(z_{\text{daughter}} = j | z_{\text{parent}} = i).

Indexing for states starts at :math:`0`. Usually the number of states is represented as the capital letter :math:`K` and indexed by :math:`k`. For most examples, we deal with two states, i.e., :math:`K=2`.

.. code:: python

    import numpy as np

    # State transition probability matrix
    T = np.array([[0.75, 0.25],
                  [0.15, 0.85]], dtype="float")

Note that the rows of the transition matrix must sum to 1 by the Law of Total Probability.

--------------

2. Creating a synthetic lineage
-------------------------------

.. code:: python

    from lineage.LineageTree import LineageTree
    from lineage.states.StateDistributionGamma import StateDistribution

2.1. Defining the :math:`\pi` initial probability vector and :math:`T` stochastic transition rate matrix
--------------------------------------------------------------------------------------------------------

Before, we “hard-coded” that the first cell in our lineage should be
state :math:`0`. In a Markov model, this first state (the state of the
root cell), like the states of the daughter cells, are probabilistically
expressed. These probabilities are stored in the :math:`\pi` initial
probability vector. In particular, if an element of the initial
probability vector , :math:`\pi`, at index :math:`i`, is defined as
:math:`\pi_{i}`, then

.. math:: \pi_{i}=\mathbb{P}(z_{0}=i).

\ We require for :math:`\pi` a :math:`K\times 1` list of probabilities.
These probabilities must add up to :math:`1` and they should be either
in a :math:`1`-dimensional list or a :math:`1`-dimensional numpy array.
An example is shown below.

.. code:: python

    # pi: the initial probability vector
    pi = np.array([0.6, 0.4], dtype="float")
    # Recall that this means that the first cell in our lineage in generation 1 
    # has a 60% change of being state 0 and a 40% chance of being state 1.
    # The values of this vector have to add up to 1 because of the 
    # Law of Total Probability.
    
    # T: transition probability matrix
    T = np.array([[0.75, 0.25],
                  [0.25, 0.75]], dtype="float")

2.2. Defining the :math:`E` emissions matrix using state distributions
----------------------------------------------------------------------

The emission matrix :math:`E` is a little more complicated to define
because this is where the user has complete freedom in defining what
type of observation(s) they care about. In particular, the user has to
first begin with defining what physical observation they will want to
extract from images of their cells, or test on synthetically created
lineages. For example, if one is observing kinematics or physics, they
might want to use the Gaussian distribution parameterized by a mean and
covariance to model their observations (velocity, acceleration, etc.).
If one wanted to model lifetimes of cell, one could utilize a Gamma 
distribution with a shape and scale parameter.
These distributions can then be combined into a multivariate
distribution.

Ultimately, the user needs to provide three things based on the
phenotype they wish to observe, model, and predict:

1. a *probability distribution function*: a function that returns a
   **likelihood** when given a **single random observation** and
   **parameters** describing the distribution
2. a *random variable*: a function that returns **random observations**
   from the distribution when given **parameters** describing the
   distribution
3. an *estimator*: a function that returns **parameters** that describe a
   distribution when given **random observations**

These three things fundamentally define any probability distribution.
For more information about how to define these functions by example, please see "2.stateDistribution.rst".

An optional boolean function can be provided to “censor”
cells based on the observation. In our example, cells with a Bernoulli
observation of :math:`0`, which implies that the cell died, are excluded
from the tree. Another censoring rule we have implemented is removing cells
that were born after an experimental end time.

We have already built, as a starting example, a model that resembles
lineage trees of cancer cells. In our synthetic model, our emissions are
multivariate. This first emission is a Bernoulli observation, :math:`0`
implying death and :math:`1` implying division. The second emission is
continuous RVs and are gamma distributed. Though these can be thought of
cell lifetimes or periods in a certain cell phase, we want the user to
know that these values can really mean anything and they are completely
free in choosing what the emissions and their values mean.

Ultimately, :math:`E` is defined as a :math:`K\times 1` size list of
``stateDistribution`` objects, explained in detail in "2.stateDistribution.rst"

The following code block is a standard way to define state distrbutions
and store them in an emissions list. State distributions are
instantiated via their parameters.

.. code:: python

    # E: states are defined as StateDistribution objects
    
    # State 0 parameters corresponding to the "Resistant" cells
    bern_p0 = 0.99   # bernoulli distribution parameter
    gamma_a0 = 7     # gamma distribution shape parameter
    gamma_scale0 = 7 # gamma distribution scale parameter
    
    # State 1 parameters corresponding to the "Susceptible" cells
    bern_p1 = 0.88
    gamma_a1 = 7
    gamma_scale1 = 1
    
    state_obj0 = StateDistribution(bern_p0, gamma_a0, gamma_scale0)
    state_obj1 = StateDistribution(bern_p1, gamma_a1, gamma_scale1)
    
    E = [state_obj0, state_obj1]

The final required parameters are more obvious. The first is the
number of cells one would like in their full uncensored lineage tree. This
can be any number. Since one of our observations is time-based, we can
also add a censoring condition based on time as well. Ultimately, these
design choices are left up to the user to customize based on their state
distribution type. Without loss of generality, we provide the following
example of a full lineage tree.

.. code:: python

    lineage1 = LineageTree.rand_init(pi, T, E, desired_num_cells=2**5 - 1)
    # These are the minimal arguments required to instantiate lineages
    print(lineage1)
    print("\n")

In the lineage above, note that the cells now have observations. Also
note that you did not have to “hard-code” the first cell and its state.
The first observation in the observation list for each cell is a
Bernoulli observation which can either be 1 or 0. An observatioon of 1
implies that the cell lived. An observation of 0 implies that the cell
died. The second observation in the observation is the gamma observation
and represents the lifetime of the cell. Note that some cells live for
far longer than others. This is because one of the states has a
probability distribution with a gamma distribution that draws longer
times.

3. Analyzing a full lineage
-----------------------------------

Our project’s goal is to analyze heterogeneity. We packaged the main
capability of our codebase into one function ``Analyze_list``, which runs the
tree-hidden Markov Model on an appropriately formatted dataset. In the
following example, we analyze the full lineage from above.

.. code:: python

    from lineage.Analyze import Analyze_list

    X = [lineage1] # population just contains one lineage
    tHMMobj_list, LL, gammas = Analyze_list([X], 2) # find two states
    tHMMobj = tHMMobj_list[0]

Estimated Markov parameters (:math:`\pi`, :math:`T`, :math:`E`)


Our model is blind to the true states of the cells
(unlike the code blocks above where we knew the identity of the cells, 
in terms of their state). This model primarily has to segment or
partition the tree and its cells into the number of states we think is
present in our data, and then identify the parameters that describe each
state’s distributions. We can not only check how well it estimated the
state parameters, but also the initial probability vector :math:`\pi`
and transition matrix :math:`T` vector. Note that estimating these also
get better as more lineages are added (for the :math:`\pi` vector in
particular) and in general as more cells and more lineages are added.

.. code:: python

    print(tHMMobj.estimate.pi)

.. code:: python

    print(tHMMobj.estimate.T)

.. code:: python

    for state in range(tHMMobj.num_states):
        print("State {}:".format(state))
        print("                    estimated state:", tHMMobj.estimate.E[state])
        print("original parameters given for state:", E[state])
        print("\n")


4. Creating a population with multiple lineages:
------------------------------------------------

The following is an analysis run on a larger set of lineages. We first create 10 lineages and append them to a list to form our cell populations.
In this case, we are choosing that all lineages should have 35 cells.
``Analyze_list()`` takes in the list of populations and the number of states,
and returns the list of ``tHMMobject`` s (one per population, ``tHMMobj_list``), the likelihood (``LL``), and the gammas used internally by the EM algorithm.
The instances of ``tHMMobj`` include the information about the distributions corresponding to each state and phenotype, and,
when ``write_states=True`` is passed, the Viterbi-predicted state for each cell is stored on each lineage's ``states`` attribute.
In this case, we are running ``Analyze_list`` with 2 states, and we know it is the true number of states,
because we used ``E`` as the Emissions which we defined as a list with two ``StateDistribution`` objects.


.. code:: python

    from lineage.Analyze import Analyze_list

    Y = []
    for _ in range(10):
        Y.append(LineageTree.rand_init(pi, T, E, desired_num_cells=35))
    tHMMobj_list, LL, gammas = Analyze_list([Y], 2, write_states=True) # find two states
    tHMMobj = tHMMobj_list[0]

.. code:: python

    print(tHMMobj.estimate.pi)

.. code:: python

    print(tHMMobj.estimate.T)

.. code:: python

    for state in range(tHMMobj.num_states):
        print("State {}:".format(state))
        print("                    estimated state:", tHMMobj.estimate.E[state])
        print("original parameters given for state:", E[state])
        print("\n")


The function ``Results()`` provides calculated features when analyzing a synthetic data.

.. code:: python

    from lineage.Analyze import Results

    results_dict = Results(tHMMobj, LL)
    print("total number of cells: ", results_dict["total_number_of_cells"])
    print("\n total number of lineages: ", results_dict["total_number_of_lineages"])
    print("\n transition matrix norm: ", results_dict["transition_matrix_similarity"])
    print("\n parameter estimtes: ", results_dict["param_estimates"])
    print("\n accuracy of state assignemnts: ", results_dict["state_similarity"])
    print("\n the distance between state 0 and state 1: ", results_dict["wasserstein"])


5. Applications - A guide to use the `tHMM` for imported experimental data.
---------------------------------------------------------------------------

As an application, we fit experimental data of cell cycle phase durations (G1 and S/G2) in response to lapatinib and gemcitabine treatments to analyze the phenotypic heterogeneity.
The data is in the form of binary tree in excel, shown beloow. We have written "lineage/LineageInputOutput.py" to properly import this data and convert it into the format usable for `tHMM`.
The following shows one lineage in the excel sheets. In each row, the difference between the two values corresponding to each cell shows the duration of G1 cell cycle phase, and 
the difference between the second value and the first value of the daughter cell shows the duration of S/G2 cell cycle phase.


.. image:: treeExcel.png
    :width: 700
    :alt: Example of a lineage tree from experimental data of G1 and S/G2 cell cycle phase durations.

.. code:: python

    from lineage.LineageInputOutput import import_exp_data
    from lineage.states.StateDistributionGaPhs import StateDistribution
    from lineage.LineageTree import LineageTree

    desired_num_states = 2 # dummy value just to initialize
    E = [StateDistribution() for _ in range(desired_num_states)]

    # control condition:
    c1 = [LineageTree(list_of_cells, E) for list_of_cells in import_exp_data(path=r"lineage/data/LineageData/AU00601_A5_1_V5.xlsx")]
    c2 = [LineageTree(list_of_cells, E) for list_of_cells in import_exp_data(path=r"lineage/data/LineageData/AU00601_A5_2_V4.xlsx")]
    c3 = [LineageTree(list_of_cells, E) for list_of_cells in import_exp_data(path=r"lineage/data/LineageData/AU00701_A5_1_V4.xlsx")]
    c4 = [LineageTree(list_of_cells, E) for list_of_cells in import_exp_data(path=r"lineage/data/LineageData/AU00801_A5_1_V4.xlsx")]
    Control = c1 + c2 + c3 + c4

    from lineage.Analyze import Analyze_list
    tHMMobj_list, LL, gammas = Analyze_list([Control], num_states=3)
    # in this example, we ran the model with 3 states.

    # finding the number of cells in the lineages:
    total_number_cells = sum([len(lineage.output_lineage) for lineage in tHMMobj_list[0].X])

    bic, dof = tHMMobj_list[0].get_BIC(LL, total_number_cells)
    print("the likelihood of having 3 states: ", LL)
    print("BIC value for this population: ", bic)
    print("The degree of freedom: ", dof)

To find out the likelihood of having different number of states we can use ``run_Analyze_over()`` 
with which we can run the model in parallel (by setting `atonce=True`) for different state numbers to minimize the run time.

To do that, we append the population for the number of states we want to analyze.
The following shows running the model for 1, 2, 3, and 4 states, in parallel, and printing the BIC value for each scenario:


.. code:: python

    from lineage.Analyze import run_Analyze_over
    import numpy as np
    desired_num_states = np.arange(1, 5)
    dataFull = []
    for _ in desired_num_states:
        dataFull.append([Control])

    # Run fitting
    output = run_Analyze_over(dataFull, desired_num_states, atonce=True)
    # output entries are (tHMMobj_list, LL, gammas) tuples, one per number of states tried
    BICs = np.array([oo[0][0].get_BIC(oo[1], total_number_cells, atonce=True)[0] for oo in output])

    print("Normalized BIC value based on the minimum: ", BICs - np.min(BICs, axis=0))