hypertraps-demos.Rmd

---
title: "Demonstration of HyperTraPS in R"
output: html_document
---

HyperTraPS (hypercubic transition path sampling) is a family of algorithms for inferring the dynamics of "accumulation" processes. These are processes where a set of binary features are acquired over time.

HyperTraPS will take a set of observed states, described by binary strings recording the presence/absence of each feature. It may optionally take initial states from which these states are reached, and information on the timings associated with each observation. It returns various summaries of which feature are acquired and when, and how these features influence each other.

### Loading the software

If we just want the HyperTraPS code without other dependencies, we only need Rcpp, and can use the following

```{r, message=FALSE}
library(Rcpp)
sourceCpp("hypertraps-r.cpp")
```

If we want various helper functions and ggplot functions in addition, use this

```{r, message=FALSE}
source("hypertraps.R")
```

### Simple demo 

Here we'll construct a simple synthetic dataset. The `m.1` matrix will store a set of initial states, and the `m.2` matrix will store a set of later observed states. The first row of `m.1`, for example, stores the state 00000, where no features have been acquired. The first row of `m.2` stores 10000, where the first of five features has been acquired.

The times correspond to each observation, so each transition has an associated time of 0.1 (in whatever units we are working with).

```{r}
m.1 = matrix(rep(c(0,0,0,0,0,
               1,0,0,0,0,
               1,1,0,0,0,
               1,1,1,0,0,
               1,1,1,1,0,
               0,0,0,0,0,
               0,0,0,0,1,
               0,0,0,1,1,
               0,0,1,1,1,
               0,1,1,1,1),5), byrow=TRUE, ncol=5)
m.2 = matrix(rep(c(1,0,0,0,0,
               1,1,0,0,0,
               1,1,1,0,0,
               1,1,1,1,0,
               1,1,1,1,1,
               0,0,0,0,1,
               0,0,0,1,1,
               0,0,1,1,1,
               0,1,1,1,1,
               1,1,1,1,1),5), byrow=TRUE, ncol=5)
times = rep(0.1, 50)
```

Let's run HyperTraPS with these "before" and "after" observations. By using `times` as both the start and end time arguments, we say that each transition takes precisely that associated time -- we could allow broader time windows to capture uncertainty in timing. Finally, we provide labels for the five individual features involved.

```{r}
my.post = HyperTraPS(m.2, initialstates = m.1, 
                     starttimes = times, endtimes = times,
                     featurenames = c("A", "B", "C", "D", "E")) 
```

That output takes us through the HyperTraPS process. First, the arguments provided to the function call are described, then the algorithm chosen (MH MCMC) and a summary of the input data is given. To estimate runtime, HyperTraPS reports the time taken for a single likelihood calculation, then scales this by the estimated number of calculations required to give a runtime estimate. Here this is only 1.17 seconds, but for chains long enough for satisfactory convergence (and more complicated systems) this may be dramatically longer.

Because we didn't turn off output, HyperTraPS periodically outputs information about the run as it goes. Every 100th step, it outputs: the likelihood of the current parameterisation; the step acceptance rate through the whole run; the step acceptance rate since the last output; the likelihood of the most recent proposal. These can help design efficient MCMC approaches: if the acceptance rate is too low and/or the recent proposal likelihood is much lower than the current, consider a smaller perturbation kernel. It also helps us see when the chain is "burned in" -- when the likelihood fluctuates, rather than consistently increasing, we're probably more stable.

After the MCMC run, the posterior analysis begins, and outputs a few details. These include the size of the posterior sample set being explored, the model being used, and a quick summary of the mean acquisition orderings of each feature. These are really just checks for debugging; not much useful information can be seen from them.

The default chain length (1000 steps) is good for a short demo, but we'll get more robust results if we run for longer. Let's do 10000 steps, suppressing output to avoid pages of updates. "length" specifies log10(number of steps).

```{r}
my.post = HyperTraPS(m.2, initialstates = m.1, 
                     starttimes = times, endtimes = times,
                     featurenames = c("A", "B", "C", "D", "E"),
                     length = 4, limited_output = 1) 
```

Now -- visualising the results.

*Transition graph plot.* This plot shows a set of sampled routes on the hypercubic space. Each edge is a transition, corresponding to the acquisition of a new feature. The edges are labelled by the feature acquired and the statistics of the time taken (mean +- s.d.). The edge width gives the probability flux through that transition.

```{r}
plotHypercube.sampledgraph2(my.post)
```

There are a lot of options for this plot type. Here are some of them. Respectively: (i) no time labels; (ii) small time labels; (iii) different label angles; (iv) cube truncated two steps after the root; (v) total timings from start, not timing of each transition; (vi) different edge style.

```{r fig.asp=1, fig.width=9}
ggarrange(plotHypercube.sampledgraph2(my.post, no.times=TRUE) + theme(legend.position = "none"),
          plotHypercube.sampledgraph2(my.post, no.times=TRUE, small.times=TRUE) + theme(legend.position = "none"),
          plotHypercube.sampledgraph2(my.post, edge.label.angle="none") + theme(legend.position = "none"),
          plotHypercube.sampledgraph2(my.post, truncate=2) + theme(legend.position = "none"),
          plotHypercube.sampledgraph2(my.post, use.timediffs=FALSE) + theme(legend.position = "none"),
          plotHypercube.sampledgraph2(my.post, use.arc=TRUE) + theme(legend.position = "none")
          )
```


*Bubble/motif plot.* These plots summarise the inferred dynamics of features, forgetting specific states. The size of a circle at point x,y, or the height of bar y at point x, gives the probability that feature y is acquired at ordinal time x. You can choose which is more beautiful/informative.

```{r}
ggarrange(plotHypercube.bubbles(my.post),
          plotHypercube.motifs(my.post))
```

*Motif time series plot.* Given a set of sampling times, this plot gives the set of likely states of the system at those times. As with the motifs above, the height of a bar gives its probability at that time.

```{r}
ggarrange(plotHypercube.motifseries(my.post, t.set=c(0,0.1,0.5,1,5)))
```

*Timing histograms.* Histograms of acquisition time for each feature, and probability that each feature is not acquired within a given time threshold (here, 1). By default this would be plotted on a log timescale; we'll linearise for this simple system.

```{r}
plotHypercube.timehists(my.post, t.thresh=1, log.time=FALSE)
```

*Time series plot.* Each step up the vertical axis corresponds to the acquisition of a new feature. The horizontal axis gives the time of the corresponding event, and the colour labels the feature that was acquired. Plotted with and without logged time axis.

```{r}
ggarrange(plotHypercube.timeseries(my.post),
          plotHypercube.timeseries(my.post, log.time=FALSE))
```

*Influences between features.* How each acquired trait (horizontal axis) influences the propensity for each other trait (vertical axis) to be acquired. Red is repression, blue is promotion; the transparency gives the width of the associated posterior. Here we see that 1 and 5 strongly cross-repress (being the first steps on the two competing pathways) and each feature promotes the next step on its pathway (wrapping the diagonal).

```{r}
plotHypercube.influences(my.post)
```

This can be perhaps more intuitively, and certainly more generally, be represented as a graph of influences: here we also impose a threshold on the coefficient of variation of the associated posterior, to report only the more "reliable" influences.

```{r}
plotHypercube.influencegraph(my.post, cv.thresh = 0.7)
```

*Likelihood trace.* Last in this list, but perhaps the first to view, the trace of likelihood estimates through the course of the MCMC search (or optimisation process; see below) offers a first diagnostic about the performance of the inference process. There are two black lines and one red. The red line shows the likelihood estimate that is actually used in the MCMC process. The two black lines show periodic, independent recalculations of the estimate. In the perfect case they will totally overlap and resemble uncorrelated white noise with no systematic trend.

```{r}
plotHypercube.lik.trace(my.post)
```

Some things that can go wrong: 
* The two black lines look very different. The path sampling process is not given consistent estimates for the likelihood. Use more walkers.
* Red line is systematically above black line. A rare extreme estimate of the likelihood has been fixed and has frozen the simulation. Use more walkers as above for more precise estimation; consider using auxiliary pseudo-marginal MCMC.
* Systematic trend in likelihood: gradual increase, changing the baseline level. Chain has not equilibrated and is taking a long time to locate good parameter sets; consider running longer chains, or increasing step size (but see below).
* Systematic trend in likelihood: punctuated increases. Chain has not equilibrated and is getting stuck in local optima; consider running longer chains, or decreasing step size (but see above).

### Predictions of unseen and future behaviour

We can use the inferred transition network to make predictions about what will happen next in a given state. Let's query the inferred hypercube to see what will happen next when we are part of the way down one pathway.

```{r}
prediction.step = predictNextStep(my.post, c(1,1,0,0,0))
plotHypercube.prediction(prediction.step)
```

The graphic reports the likely next steps, as a word cloud and a set of probabilities. Here, we're making the strong prediction that the next step after 11000 is 11100.

We can also make predictions about hidden values in new observations. First let's mask an observation and ask HyperTraPS to "fill in the blanks" of 1????.

```{r}
prediction.hidden = predictHiddenVals(my.post, c(1,2,2,2,2))
plotHypercube.prediction(prediction.hidden)
``` 

Here, we see a word cloud of possible states that could correspond to our observation, and a set of probabilities associated with each *feature* being 1. That is, it's highly likely that feature 2 is 1, and less likely that feature 5 is 1. But doesn't this assume something about how far we've come on the hypercube?

Yes -- by default this assumes that all "levels" of the hypercube -- all counts of "1" that are possible, given our uncertainty -- are equally likely. But we might be more interested in specifying something about how many features we're likely to have acquired. We can do this by specifying weights for each "level":

```{r}
prediction.hidden = predictHiddenVals(my.post, c(1,2,2,2,2), level.weight=c(0,0,1,0,0,0))
plotHypercube.prediction(prediction.hidden)
```

Here, we're saying that we believe there's only one more "1" in the ?s -- and we see a correspondingly strong prediction about where that "1" is. We can be more agnostic, and the predictions are shaped accordingly, like below. We might want our prediction to mirror how many "levels" we saw in our original data, for example.

```{r}
predictHiddenVals(my.post, c(1,2,2,2,2), level.weight=c(0,0,1,1,1,0))
plotHypercube.prediction(prediction.hidden)
```

### Uncertainty in observation times

HyperTraPS-CT allows uncertainty in transition timings to be captured. In the run above, we specified a precise timing for each transition, by giving a zero-width time window for each observation (start times = end times). This time window is flexible. We can allow it to have infinite width, in which case absolute timings are meaningless and we just infer the orderings of events. Or we can specify a finite time window, allowing a transition to take any time within that window, to capture uncertainty in observation timings.

```{r}
# precisely specified timings, as above
my.post.time.precise = HyperTraPS(m.2, initialstates = m.1, 
                               starttimes = times, endtimes = times,  
                               limited_output = 1,
                               featurenames = c("A", "B", "C", "D", "E")); 
# infinite width time window for transitions (equivalent to just inferring ordering)
my.post.time.inf = HyperTraPS(m.2, initialstates = m.1, 
                                starttimes = times*0, endtimes = times*Inf, 
                                limited_output = 1,
                                featurenames = c("A", "B", "C", "D", "E"));
# finite time window for each uncertain transition time
my.post.time.uncertain = HyperTraPS(m.2, initialstates = m.1, 
                     starttimes = times*0.25, endtimes = times*4, 
                     limited_output = 1,
                     featurenames = c("A", "B", "C", "D", "E")); 
ggarrange(plotHypercube.timehists(my.post.time.precise, t.thresh=3), 
          plotHypercube.timehists(my.post.time.uncertain, t.thresh=3),
          plotHypercube.timehists(my.post.time.inf, t.thresh=3),
          nrow=3)
plotHypercube.sampledgraph2(my.post.time.precise, thresh=0.1, use.arc=FALSE, edge.label.size=3) + theme(legend.position="none") + expand_limits(x = c(-0.1, 1.1))
plotHypercube.sampledgraph2(my.post.time.uncertain, thresh=0.1, use.arc=FALSE, edge.label.size=3) + theme(legend.position="none") + expand_limits(x = c(-0.1, 1.1))
plotHypercube.sampledgraph2(my.post.time.inf, thresh=0.1, use.arc=FALSE, edge.label.size=3) + theme(legend.position="none") + expand_limits(x = c(-0.1, 1.1))
```

The summary plots here reflects the different inferences about acquisition timescales that come from the approaches with different time windows. 

### Input/output between R data structure and file format

We can write the output of HyperTraPS to a set of files for storage and later retrieval

```{r}
writeHyperinf(my.post, "simpledemo", my.post$L, postlabel = "simpledemo", fulloutput=TRUE)
```

We can retrieve output from files, which may have been previously written as above, or may have come from running HyperTraPS at the command line.

```{r}
my.post.r = readHyperinf("simpledemo", postlabel = "simpledemo", fulloutput=TRUE)
```

### Other functional examples

HyperTraPS allows substantial flexibility in how evolutionary pathways are inferred, how likelihoods are estimated, and other aspects of the approach. 

If we want to sacrifice some accuracy in estimating likelihood for computational speed, we can run an example with fewer walkers sampling pathways (`walkers`)

```{r}
my.post.sparse = HyperTraPS(m.2, initialstates = m.1, 
                            starttimes = times, endtimes = times,
                            featurenames = c("A", "B", "C", "D", "E"), 
                            walkers = 2,
                            limited_output = 1)
```

The original, discrete-time HyperTraPS approach is recovered if we don't use timing information

```{r fig.width=9}
my.post.dt = HyperTraPS(m.2, initialstates = m.1, 
                        featurenames = c("A", "B", "C", "D", "E"),
                        limited_output = 1) 
plotHypercube.summary(my.post.dt, continuous.time = FALSE)
```
          
We can ask HyperTraPS to perform stepwise regularisation after fitting a model, pruning extraneous parameters (`regularise`). The regularisation plot shows how the information criterion behaves as we prune back parameters -- the minimum gives us our optimal model.

```{r fig.width=9}
my.post.regularise = HyperTraPS(m.2, initialstates = m.1, regularise = 1,
                                walkers = 20,
                                limited_output = 1)
plotHypercube.regularisation(my.post.regularise)
plotHypercube.summary(my.post.regularise)
```

We can use simulated annealing to get a maximum likelihood estimate rather than a Bayesian picture (`sa`)

```{r fig.width=9}
my.post.sa = HyperTraPS(m.2, initialstates = m.1, sa = 1,
                        limited_output = 1)
plotHypercube.summary(my.post.sa)
```

We can also use "phenotype landscape inference" -- an unbiased sampling method, unlike the (corrected) bias sampling in HyperTraPS (`PLI`). This takes longer but is more flexible with uncertain data

```{r fig.width=9}
my.post.pli = HyperTraPS(m.2, initialstates = m.1, pli = 1,
                         limited_output = 1)
plotHypercube.summary(my.post.pli)
```

HyperTraPS supports different parameter structures (`model`). By default we use an "L2" parameterisation, where each feature can individually influence the acquisition of every other feature. "L1" has features completely independent; "L0" has all features identical. "L3" allows *pairs* of features to influence each feature's acquisition; "L4" allows *triplets* of features to influence each feature's acquisition. This approach, labelled -1, allows every edge on the hypercube to have its own independent parameters -- corresponding to the most flexible possible picture, where arbitrary sets of features influence other features. We then regularise as above (`regularise`) to remove extraneous parameters

```{r}
my.post.bigmodel.regularise = HyperTraPS(m.2, initialstates = m.1, model = -1,
                                         regularise = 1, walkers = 20,
                                         limited_output = 1)
plotHypercube.regularisation(my.post.bigmodel.regularise)
```

Here's another example of model choice. The data is now generated by a process where *pairs* of features influence the acquisition of other features. To capture these interactions, an "L^2" picture (`model = 2`, where each feature *individually* influences each other feature) is insufficient -- an "L^3" picture (`model = 3`) allowing pair influence is needed. The full parameterisation, where every edge on the hypercube transition network (`model = -1`) has an independent parameter, can also capture this behaviour.

```{r fig.asp = 1.2}
logic.mat = readLines("RawData/old-hi-order.txt")
logic.mat = do.call(rbind, lapply(strsplit(logic.mat, ""), as.numeric))
logic.mat = rbind(logic.mat, logic.mat)
logic.mat.i = readLines("RawData/old-hi-order-init.txt")
logic.mat.i = do.call(rbind, lapply(strsplit(logic.mat.i, ""), as.numeric))
logic.mat.i = rbind(logic.mat.i, logic.mat.i)
logic.starts = logic.mat.i
logic.ends = logic.mat

logic.post.m1 = HyperTraPS(logic.ends, initialstates = logic.starts, length = 4, model = -1, walkers = 20, limited_output = 1)
logic.post.1 = HyperTraPS(logic.ends, initialstates = logic.starts, length = 4, model = 1, walkers = 20, limited_output = 1)
logic.post.2 = HyperTraPS(logic.ends, initialstates = logic.starts, length = 4, model = 2, walkers = 20, limited_output = 1)
logic.post.3 = HyperTraPS(logic.ends, initialstates = logic.starts, length = 4, model = 3, walkers = 20, limited_output = 1)

ggarrange(plotHypercube.graph(logic.post.m1) + ggtitle("All edges") + theme(legend.position="none"),
          plotHypercube.graph(logic.post.1) + ggtitle("L") + theme(legend.position="none"),
          plotHypercube.graph(logic.post.2)+ ggtitle("L^2") + theme(legend.position="none"),
          plotHypercube.graph(logic.post.3)+ ggtitle("L^3") + theme(legend.position="none"))
```    

We can see that the structures inferred by the "full" and "L^3" models are the same, while the "L^2" (and inappropriate "L^1") will either introduce extraneous transitions or fail to capture the true ones.

### Prior information

The Bayesian implementations of HyperTraPS-CT allow prior information to be included in the inference process. We can do this by placing bounds on the values that the parameters in our model are allowed to take. Here, we first assign wide priors (-10 to 10 in log space) to all the 25 parameters in our model. Then we go through the parameters corresponding to the base rates for each feature: these are on the diagonal of the parameter matrix. We restrict these values to be -10, except for the base rates for feature 1 and feature 5. We thus enforce feature 1 > feature 5 >> other features.

```{r}
priors = matrix(0, ncol=2, nrow=5*5)
priors[,1] = -10
priors[,2] = 10
for(i in 0:4) {
  priors[i*5+i+1,1] = -10
  priors[i*5+i+1,2] = -10
}
priors[0*5+0+1,1] = 1
priors[0*5+0+1,2] = 1
priors[4*5+4+1,1] = 0
priors[4*5+4+1,2] = 0
```

Running and visualising the posteriors shows this effect:

```{r fig.width=9}
my.post.priors = HyperTraPS(m.2, initialstates = m.1, 
                     starttimes = times, endtimes = times, 
                     priors = priors,
                     featurenames = c("A", "B", "C", "D", "E"),
                     limited_output = 1)
plotHypercube.summary(my.post.priors)
```

### Scientific examples

A set of short-form examples from past studies -- these should run in a few minutes and give approximations to the original results. In each case we read the observed states from a file (these are present in different formats, so different curation steps are involved in each case), and feature labels from another file. Then we run HyperTraPS and plot some summaries of the output.

1. Ovarian cancer case study: traits are chromosomal aberrations, observations are independent patient samples.

```{r}
cgh.mat = readLines("RawData/ovarian.txt")
cgh.mat = do.call(rbind, lapply(strsplit(cgh.mat, ""), as.numeric))
cgh.names = as.vector(read.table("RawData/ovarian-names.txt", sep=","))[[1]]

my.post.cgh = HyperTraPS(cgh.mat, 
                        length = 3, 
                        featurenames = cgh.names,
                        limited_output = 1) 
ggarrange(plotHypercube.lik.trace(my.post.cgh), 
          plotHypercube.bubbles(my.post.cgh, reorder=TRUE), 
          plotHypercube.sampledgraph2(my.post.cgh, no.times=TRUE), nrow=3)
plotHypercube.sampledgraph2(my.post.cgh, no.times=TRUE)
```

2. C4 photosynthesis case study: traits are physical/genetic features associated with C4, observations are (incomplete) phylogenetically independent intermediate species.

This case study involves some uncertain data, which are present as "2"s in the source data -- labelling features which may be 0 or 1.

```{r}
c4.mat = as.matrix(read.table("RawData/c4-curated.csv", sep=","))
c4.names = as.vector(read.table("RawData/c4-trait-names.txt", sep=","))[[1]]

my.post.c4 = HyperTraPS(c4.mat, 
                        length = 4, 
                        featurenames = c4.names,
                        limited_output = 1) 
plotHypercube.bubbles(my.post.c4, reorder=TRUE)
```