2
$\begingroup$

I am initiating a project reproducing MEGA style trees using R to find the best DNA protein models under AIC/BIC. I have read the R AICs and BIC documentation but I'm struggling to implement this.

Could I implement AICs and BICs without having to specify the number of sequences in the fasta file?

My initial code that works the same as MEGA for aligning and building a phylogenetic tree where I want to view AICs and BIC (as per reference picture):

library(seqinr)
library(adegenet)
library(ape)
library(ggtree)
library(DECIPHER)
library(Biostrings)
library(viridis)
library(ggplot2)
library(msa)

mySequenceFile <- system.file ("sequences", "allSequences.fasta", package ="msa")
mySequence <- readDNAStringSet(mySequenceFile)

mySequence

myAlignment <- msa(mySequence, "Muscle")

myAlignment

alignment <- msaConvert(myAlignment, "seqinr::alignment")
distMatrix <- dist.alignment(alignnment, "similarity")
clustering <- hclust(distMatrix)

plot(clustering)

dendrogram <- as.dendrogram(clustering)
phylotree <- as.phylo(clustering)

plot(phylotree, type="radial")

Please ignore some libraries as I previously ran into some "-" deletion issues when I was writing the code

$\endgroup$
0

1 Answer 1

1
$\begingroup$

Having looked carefully at the question, I understand the 'gap in the market' but all this code will do is generate a provisional tree. I would suggest looking at protein models first because you need models to perform an AIC/BIC. So far you haven't generated used or described a protein model and without that there's nothing to compare via AIC.

To answer the direct question you should seek to avoid removing the sample size from the AIC and base the work on the equation below. If you have large samples sizes n can be removed (which is what you are wanting).

The work you should read is jModeltest which is here, https://github.com/ddarriba/jmodeltest2. This originated from Modeltest by David Pasado.

Their code on GitHub will tell you exactly how to implement AIC and BIC. I recommend ditching R in this instance. If done via subprocessing (below) its not hard.

Their work concerns nucleotides and never really hit the mainstream, but was used by some.

The key criticism is like your approach that the tree generated was based on neighbor-joining. The tree will influence the matrix so there is a danger that the ML or Bayesian phylogeny resulting from the model will be biased by neighbour-joining. Secondly, it was not computationally expensive to do this entirely by ML or Bayesian approaches because it's a 4*4 matrix.

Protein phylogeny is very different, obviously the matrix is bigger 20*20 the other issue is that it's a 'static matrix'. There are approaches where they "combine" matrices but I don't know the algorithm. Thus we know exactly how the protein matrices are generated, "combined" needs a clearer explanation and how this is obtained could be a good application for AIC. Thus could provide useful clarity. There is in addition.the potential to generate new matrices and assess them, but the concept of introducing meaningful dynamism into the matrix does have potential and this is absent (except for cryptic combined options). AIC would then assess whether this is helpful. It might not be hugely helpful in traditional static matrices, but would compare them against further hypotheses of mutation and innovative approaches. If used deployed correctly could be cool because it allows you to experiment and ask the question what precisely does this protein matrix mean for my data?

More generally stretching the utility of AIC type metrics is very interesting in 'bleeding edge' machine learning (admittedly its a bit futuristic).


I suspect (could be wrong) jModelTest is coded in Java, an alternative would be to subprocess via Python using relevant compiled code. This was how the original Modeltest was written, essentially driven through the leading maximum likelihood package of the time known as PAUP. Python does implement AIC/BIC via the statsmodel library, scikit-learn (machine learning) or use the equation directly.

enter image description here

k= number of independent variables to build model L= maximum likelihood estimate of model .. for small sample sizes (n)

R is not known for subprocessing and I am not aware of any ML/Bayes libraries for phylogeny.

In summary, this is a good project and it is worth investing the time because it has a clear, meaningful output and potential to extend into more innovative experimental approaches. I for one would use established code (i.e. written by someone else) that implemented AIC/BIC for proteins matrices and would be fairly straight forward (for you) to write. How you maximise its use in deploying it to help understand protein matrices needs creative thinking (remember its parameters verse likelihood).

Good luck.

$\endgroup$
4
  • 1
    $\begingroup$ Maybe instead of using ML I can try to play around DL as well with some good datasets. Thanks for the the guidance, some things make more sense now and I will try to do everything using python instead. $\endgroup$ Feb 13 at 7:26
  • $\begingroup$ Cool, DL ... deep learning? Scikit-learn does machine learning, not quite the same. ML/ AIC looks a good staging post (a reviewer will ask for BIC I'm afraid). Please do remember to upvote. $\endgroup$
    – M__
    Feb 13 at 15:21
  • 1
    $\begingroup$ Yes I meant deep learning. I played around some python scripts last night and R as well. Will be creating a GitHub repository later on this week with my progress. Unfortunately I can't upvote yet due to my reputation but I choose it as the answer. Thanks again. $\endgroup$ Feb 13 at 21:08
  • $\begingroup$ Deep learning isn't trivial and you'll likely encounter skepticism, ML sci-kit lean is widely accepted and much easier. $\endgroup$
    – M__
    Feb 14 at 3:01

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.