uscbiostats / aphylo

Statistical inference of genetic functions in phylogenetic trees

USCbiostats

Name	Description
AnnoQR	R client wrap for AnnoQ API (https://github.com/blueOwl/AnnoQR)
aphylo	Statistical inference of genetic functions in phylogenetic trees
bayesnetworks	C++ program to fit Bayesian networks, illustrated with simulated data
BinaryDosage	Converts VCF files to a binary format
causnet	What the Package Does (One Line, Title Case)
fdrci	Permutation-Based FDR Point and Confidence Interval Estimation
fmcmc	A friendly MCMC framework
GxEScanR	R version of GxEScan
HiLDA	An R package for inferring the mutational exposures difference between groups
hJAM	hJAM is a hierarchical model which unifies the framework of Mendelian Randomization and Transcriptome-wide association studies.
iMutSig	A web application to identify the most similar mutational signature using shiny
jsPhyloSVG	htmlwidgets for the jsPhyloSVG JavaScript library
LUCIDus	Latent and Unknown Cluster Analysis using Integrated Data
MergeBinaryDosage	R package for merging binary dosage files
partition	A fast and flexible framework for agglomerative partitioning in R
pfamscanr	An R client for EMBL-EBI’s PfamScan API
polygons	Flexible functions for computing polygons coordinates in R
rphyloxml	Read and write phyloXML files in R
selectKSigs	Selection of K in finding the number of mutational signatures
slurmR	slurmR: A Lightweight Wrapper for Slurm
xrnet	R Package for Hierarchical Regularized Regression to Incorporate External Data
xtune	An R package for Lasso and Ridge Regression with differential penalization based on prior knowledge

aphylo's Issues

Before starting simulations

• write tests for aphylo_formula.
• write more tests for aphylo_mcmc using formulas (compare against raw MCMC)
• rewrite predict.aphylo_estimates so that it uses the loglikelihood function created from the model call.
• Include a warning when fixing parameters in the MLE estimation.

Fix parameters in the likelihood

Create functions of type: obj_fun(p, dat, prior), obj_fun_fix_eta(p, dat, prior), obj_fun_fix_psi(p, dat, prior)

Plot loglikelihood

For the plot.aphylo_estimates method, use the modified version of the likelihood, i.e. using priors. Right now it is plotting the solution without considering the priors.

Extending the MCMC function

While we are intending to make the MCMC function its own R package, and hence, having its own repo, we will start working on the improvements right here. The idea of this Github-issue is to list all suggestions (wishes) for the MCMC function. For now, the devel road has the following:

1. Implement Tempering MCMC.

2. Implement Adaptative MCMC (global and local).

3. Parallel MCMC (multiple chains).

Besides, we would like to include some benchmarks with the following packages: rjags, rstan, mcmc, etc. on some standard problems

Implement the MAP estimator

Phylogenetic Trees Visualization using the ggtree package

In Bioconductor, there's an R package designed for visualization of annotated phylogenetic trees

ggtree: a phylogenetic tree viewer for different types of tree annotations

It might be worth taking a look at it.

Numerical instability and computational complexity of the LogLike function

After a few tests of the phylo_mcmc function I incidentally realized something that Duncan mentioned some time ago, the complexity of the algorithm is bad. On raw, and perhaps sloppy, terms, as a function of the number of 'functions' -p-, the algorithm performs (p*2^p) computations, ergo, for p = 1, 2, 4 we get

1x2^1 = 2
2x2^2 = 8
4x2^4 = 64

That said, I need to look more carefully to the LogLike function itself, and furthermore, to the probabilities function in C++ to see if we can improve efficiency. This is specially important for the MCMC case.

Here are some tests with the current state of the function:

> benchmark(
+     `1 fun`=with(dat1, LogLike(experiment, offspring, noffspring, c(.2,.2), c(.2,.2), .5)),
+     `2 fun`=with(dat2, LogLike(experiment, offspring, noffspring, c(.2,.2), c(.2,.2), .5)),
+     `3 fun`=with(dat3, LogLike(experiment, offspring, noffspring, c(.2,.2), c(.2,.2), .5)), 
+     replications=100, relative="elapsed")[,1:4]
   test replications elapsed relative
1 1 fun          100   0.005      1.0
2 2 fun          100   0.011      2.2
3 3 fun          100   0.027      5.4

The instability part comes from the fact that, as P increases (and n as well), the size of the rootnode probabilities tend to zero, which is another issue that I encountered throughout the process. That makes the LogLike undefined. Need to take a look on that too.

Need to change the name of the package

Considering CRAN as a Task view on Phylogenetics (https://cran.r-project.org/web/views/Phylogenetics.html), the current name is too generic and should be more informative.

Visualization tools

Web-based

• Suggested by Paul Thomas, we can use iTOL. Seems to have API.

• This is a javascript SVG based visualization tool jsPhyloSVG. Seems to have API, and moreover, suitable to be used with Shiny. Elsevier uses this. So far, this is the tool that I like the most.

• SILVA Tree Viewer is a webapp.

• TreeLink is a webapp.

• PhyD3 makes interactive visualizations. Suitable for Shiny.

In R

• The ape::plot.phylo function has several options that I was not aware of. Now that I realized that the coordinates of the nodes are stored in .PlotPhyloEnv$last_plot.phylo as xx and yy, I can use this as input to create more personalized visualizations.

• The ggtree R package, currently used here, is not working properly when handling the annotations.

Re-implement multi-function trees likelihood calculation

Use `ape::as.phylo` methods

The phylo class of the ape package seems to be a very popular storage method for phylogenetic trees, so it is worthwhile either storing this package's objects as phylo objects or create conversion functions to take advantage of the ape package.

leave-one-out cross validation

Need to implement this to assess the quality of the predictions.

`rdrop_annotations` returns with error when tip.annotations has a 9/NA

phylo_mcmc has a bug when control = list(nchains) is greater than 1

the function must account for this. Also, it'll be nice if we include the gelman test.

Using reorder(x, "postorder") already gives us a peeling sequence!

Just figured out that ape::reorder(., "postorder") can be mapped to the po-tree that we are creating with aphylo. This means, we don't need to have a special class of object to store the trees, it suffices to pass the second column of x$edges after calling the reorder function and use that as the peeling sequence. This is both encouraging and a shame since: (1) it makes more compatible still with the ape package, and (2) means that all the development on PO trees that I've done is useless 😭 (but that's OK!).