4
$\begingroup$

In classic phylogenetic inference one is usually given various orthologue sequences of a given gene across various species. Those sequences are then multiple aligned and used to construct a phylogenetic tree.

Say, I split the multiple alignment of this gene into two parts (part A and B). My question now is: What is the most appropriate test to determine if part A and B share the same phylogenetic history?

Having looked through literature, it appears this can be done by constructing phylogenetic trees separately for part A and B and test for incongruence using

  • Kishino–Hasegawa (KH test);
  • Shimodaira–Hasegawa (SH-test);
  • SOWH-test etc.

as nicely summarized in Plantet et all, and Goldman et al. etc.

However the above test rely on constructed trees beforehand which might be incorrect due to inappropriate algorithms/parameters.

I was therefore wondering if there are more appropriate test maybe without having to construct trees at all?

Improve this question
$\endgroup$
9
  • 1
    $\begingroup$ I don't really see how that would be possible, off the top of my head. After all, when you construct trees, you are "operating directly on the sequences" and trees are the tool for determining evolutionary history. Could you edit your question and explain more about why you don't want to use trees? Why don't you construct the tree from your alignment? $\endgroup$
    – terdon ♦
    Jun 7, 2018 at 11:58
  • $\begingroup$ Well I didn't mean to say that I don't want to use trees, just if there are methods that don't need to compute them. I was just concerned that the way the tree is constructed might introduce a bias. I'd be happy to learn that SH-test and the like are appropriate approaches to answer the stated question. I just wanted to make sure I didn't overlook more appropriate approaches. $\endgroup$
    – Sebastian Müller
    Jun 7, 2018 at 16:02
  • $\begingroup$ I'm not sure that you can use any of the tests you've listed on fundamentally different datasets - the likelihoods will always differ in incomparable ways. One way you could do this is, if you have a priori two different hypotheses, calculate the site-specific likelihood of the whole alignment on each of the two topologies and compare the differences in likelihood across sites. If one segment of the alignment is biased towards one resolution, the difference in likelihoods between the two topologies should show that. $\endgroup$
    – NatWH
    Jun 8, 2018 at 16:00
  • $\begingroup$ Although I see in the Plantet paper the two topological-based tests, which I guess would work. The only question I have with those is how they accommodate uncertainty. It may be, for example, that in partitioning your dataset you arrive at block B containing very little phylogenetic signal for any hypothesis, and treating its point estimate as concrete is lending unjustified support to it. $\endgroup$
    – NatWH
    Jun 8, 2018 at 16:13
  • 1
    $\begingroup$ Agreed! I just read you are going to start at the Plant Science Dept in Cambridge, where I also work. Small world! We could meet to discuss in person at some point, I'm in the Baulcombe lab. Thanks again for your input anyway! $\endgroup$
    – Sebastian Müller
    Jun 11, 2018 at 10:24

1 Answer 1

Reset to default
2
$\begingroup$

Since there hasn't been any other answers, I'd thought to share my approach in case someone is having the same issue.

Gauging from the comments of tendon and NatWH, it seems using the SH-test valid stategy to tackle this problem.

Having done some research, I could only find one functional software package, the R-package phangorn which implements both, the SH- as well as SOWH-test and thought to share my R-workflow:

Reading in multiple-alignment (MLN). Here: amino-acids (AA)

 library(phangorn)
 mln <- read.phyDat("my_mln.fasta", format = "fasta", type = "AA")

Partitioning MLN into two parts (part A and B) given a predefined breakpoint (here 100)

breakpoint <- 100
 # partA: 1:breakpoing
 ( mln_partA <- subset(mln, select = 1:breakpoint, site.pattern = FALSE) )
 # 58 sequences with 100 character and 45 different site patterns.
 # The states are a r n d c q e g h i l k m f p s t w y v

 ## partB: breakpoint:end
 ( mln_partB <- subset(mln, select = (breakpoint + 1):length(attr(mln, "index")), site.pattern = FALSE) )
 # 58 sequences with 1702 character and 720 different site patterns.

Constructing tree(s) separately from part A and B using Neighbourhood Joining (NJ) method

 ## Note, other tree building algorithms might be used
 ( treeNJ_partA <- NJ(dist.ml(mln_partA)) )
 # Phylogenetic tree with 58 tips and 56 internal nodes.
 # Unrooted; includes branch lengths.
 ( treeNJ_partB <- NJ(dist.ml(mln_partB)) )
 # Phylogenetic tree with 58 tips and 56 internal nodes.
 # Unrooted; includes branch lengths.

Using topology of previously constructed trees and compute optimal likelihood based on partA data

 ( treeNJ_partA_fit1 <- optim.pml(pml(treeNJ_partA, mln_partA)) )
 # optimize edge weights:  -1402.126 --> -1379.043
 #
 #  loglikelihood: -1379.041
 #
 # unconstrained loglikelihood: -302.6292
 ( treeNJ_partB_fit1 <- optim.pml(pml(treeNJ_partB, mln_partA)) )
 # optimize edge weights:  -1633.646 --> -1530.745
 #
 #  loglikelihood: -1530.682
 #
 # unconstrained loglikelihood: -302.6292

Performing SH-test for for both trees

 ( shPart1_t12 <- SH.test(treeNJ_partA_fit1, treeNJ_partB_fit1, data = mln_partA, B = 1e5) )
 #      Trees      ln L Diff ln L p-value
 # [1,]     1 -1379.041     0.000 0.48461
 # [2,]     2 -1530.682   151.641 0.00000

As can be seen in this example, the tree constructed for partB doesn't fit partA data as much as the tree constructed from partA.

Improve this answer
$\endgroup$
7
  • $\begingroup$ Thanks for sharing! But perhaps it would be clear if the output is removed, or if the code is explained before the block $\endgroup$
    – llrs
    Jun 12, 2018 at 13:35
  • $\begingroup$ Hmm, in fact I included the output to make it clearer since I like to know intermediate results. I'll exclude a few lines though and maybe someone else can share their opinion. $\endgroup$
    – Sebastian Müller
    Jun 12, 2018 at 14:48
  • $\begingroup$ Well, it is only my opinion! But for me it would be clearer if you explained the code alongside breaking the block into smaller ones. As is I find hard to know what is going on $\endgroup$
    – llrs
    Jun 12, 2018 at 14:58
  • 1
    $\begingroup$ Agreed. I broke it into chunks, hope this improved it! $\endgroup$
    – Sebastian Müller
    Jun 12, 2018 at 15:16
  • $\begingroup$ Nitpicking: whatever model dist.ml uses for amino acids by default may not be the optimal model for the dataset (given a tree), unless you selected that elsewhere and did not show it in the code. Likewise I don't think the NJ tree of ML distances is guaranteed to give the ML tree. It appears that optim.pml does not perform topological searches by default. However, the assumptions of the SH-test call for the ML tree of the data. I would instead use an external program (RAxML, IQ-TREE), to find the overall ML topology and the ML topology of each subset, and run the SH-test on those. $\endgroup$
    – NatWH
    Jun 14, 2018 at 14:02

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

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