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

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

<p>Having done some research, I could only find one functional software package, the R-package <a href="https://github.com/KlausVigo/phangorn" rel="nofollow noreferrer">phangorn</a> which implements both, the SH- as well as SOWH-test and thought to share my R-workflow:</p>



<p><strong>Reading in multiple-alignment (MLN). Here: amino-acids (AA)</strong></p>

<pre class="lang-r prettyprint-override"><code> library(phangorn)
 mln &lt;- read.phyDat("my_mln.fasta", format = "fasta", type = "AA")
</code></pre>

<p><strong>Partitioning MLN into two parts (part A and B) given a predefined breakpoint (here 100)</strong></p>

<pre class="lang-r prettyprint-override"><code>breakpoint &lt;- 100
 # partA: 1:breakpoing
 ( mln_partA &lt;- 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 &lt;- subset(mln, select = (breakpoint + 1):length(attr(mln, "index")), site.pattern = FALSE) )
 # 58 sequences with 1702 character and 720 different site patterns.
</code></pre>

<p><strong>Constructing tree(s) separately from part A and B using Neighbourhood Joining (NJ) method</strong></p>

<pre class="lang-r prettyprint-override"><code> ## Note, other tree building algorithms might be used
 ( treeNJ_partA &lt;- NJ(dist.ml(mln_partA)) )
 # Phylogenetic tree with 58 tips and 56 internal nodes.
 # Unrooted; includes branch lengths.
 ( treeNJ_partB &lt;- NJ(dist.ml(mln_partB)) )
 # Phylogenetic tree with 58 tips and 56 internal nodes.
 # Unrooted; includes branch lengths.
</code></pre>

<p><strong>Using topology of previously constructed trees and compute optimal likelihood based on partA data</strong></p>

<pre class="lang-r prettyprint-override"><code> ( treeNJ_partA_fit1 &lt;- optim.pml(pml(treeNJ_partA, mln_partA)) )
 # optimize edge weights:  -1402.126 --&gt; -1379.043
 #
 #  loglikelihood: -1379.041
 #
 # unconstrained loglikelihood: -302.6292
 ( treeNJ_partB_fit1 &lt;- optim.pml(pml(treeNJ_partB, mln_partA)) )
 # optimize edge weights:  -1633.646 --&gt; -1530.745
 #
 #  loglikelihood: -1530.682
 #
 # unconstrained loglikelihood: -302.6292
</code></pre>

<p><strong>Performing SH-test for for both trees</strong></p>

<pre class="lang-r prettyprint-override"><code> ( shPart1_t12 &lt;- 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
</code></pre>

<p>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.</p>

CC BY-SA 4.0
2018-06-12T13:35:56.373Z
Thanks for sharing! But perhaps it would be clear if the output is removed, or if the code is explained before the block
llrs
CC BY-SA 4.0
2018-06-12T14:48:24.937Z
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.
Sebastian Müller
CC BY-SA 4.0
2018-06-12T14:58:54.580Z
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
llrs
CC BY-SA 4.0
2018-06-12T15:16:21.027Z
Agreed. I broke it into chunks, hope this improved it!
Sebastian Müller
CC BY-SA 4.0
2018-06-14T14:02:27.207Z
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.
NatWH
CC BY-SA 4.0
2018-06-14T14:03:17.470Z
You could also use one of the implementations of model finder (AUTO in RAxML or MFP in IQ-TREE) to find the best-fitting model of substitution.
NatWH
CC BY-SA 4.0
2018-06-15T13:11:02.847Z
That's true. In fact I've used RAxML for the original data but used non-optimized NJ deliberately to make a more reproducible example and keeping the code simple. In fact, this underlines my original point: Testing incongruence between partitions shouldn't depend (too much) on the tree building algorithm or/and optimizations. If it does, I'd love to read about guidelines, best practices or examples on how to do this properly.
Sebastian Müller
CC BY-SA 4.0
2018-06-12T13:25:45.527Z
2018-06-12T15:15:18.260Z
2018-06-12T15:15:18.260Z
Sebastian Müller
<p>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.</p>

<p>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?</p>

<p>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 </p>

<ul>
<li>Kishino–Hasegawa (<strong>KH test</strong>); </li>
<li>Shimodaira–Hasegawa (<strong>SH-test</strong>); </li>
<li><strong>SOWH-test</strong> etc. </li>
</ul>

<p>as nicely summarized in <a href="https://www.sciencedirect.com/science/article/pii/S153204640500081X?via%3Dihub" rel="nofollow noreferrer">Plantet et all,</a> and <a href="https://doi.org/10.1080/106351500750049752" rel="nofollow noreferrer">Goldman et al.</a> etc.</p>

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

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

CC BY-SA 4.0
2018-06-07T11:58:14.560Z
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?
terdon
CC BY-SA 4.0
2018-06-07T16:02:29.640Z
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.
Sebastian Müller
CC BY-SA 4.0
2018-06-08T16:00:18.973Z
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.
NatWH
CC BY-SA 4.0
2018-06-08T16:13:28.530Z
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.
NatWH
CC BY-SA 4.0
2018-06-11T09:18:40.877Z
Good points @NatWH. Accommodating uncertainty is indeed an issue and I suppose one has to show that there is a sufficient amount of information bearing sites for each partitions, however I couldn't find suitable numbers anywhere. I'd be curious on any good practice guidelines or maybe papers that have done it nicely since this seems of general interests.
Sebastian Müller
CC BY-SA 4.0
2018-06-11T09:27:58.850Z
As for your first point @NatWH, I had the following in mind: Calculating the maximum likelihood trees (tree A and B) for both partitions (part A and B) separately. I'd then calculate the likelihoods and perform a incongruence test only for part A data under both tree topologies of tree A and B. Same procedure for part B as well for completeness. So the test is in using the same dataset(s) and the likelihood should be comparable, do you agree?
Sebastian Müller
CC BY-SA 4.0
2018-06-11T10:13:00.297Z
Yes I agree, that shouldn't violate the assumptions of the SH test, although I am just not familiar enough with the SOWH to say anything about that. I would still advise determining the site-specific differences in log-likelihood, though. Its helpful to know the distribution of the signal.
NatWH
CC BY-SA 4.0
2018-06-11T10:24:50.093Z
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!
Sebastian Müller
CC BY-SA 4.0
2018-06-11T13:06:10.723Z
Absolutely! I look forward to that in September.
NatWH
<phylogenetics><phylogeny><r>
Test to determine if two genes/exons share the same evolutionary histories?
# Test to determine if two genes/exons share the same evolutionary histories?

(Asked by: username_0 on 2018-06-06)

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,](https://www.sciencedirect.com/science/article/pii/S153204640500081X?via%3Dihub) and [Goldman et al.](https://doi.org/10.1080/106351500750049752) 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?

**Tags:** phylogenetics, phylogeny, r
> **Comments:**
> **username_1**: 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?  
> **username_0**: 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.  
> **username_2**: 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.  
> **username_2**: 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.  
> **username_0**: Good points @NatWH. Accommodating uncertainty is indeed an issue and I suppose one has to show that there is a sufficient amount of information bearing sites for each partitions, however I couldn't find suitable numbers anywhere. I'd be curious on any good practice guidelines or maybe papers that have done it nicely since this seems of general interests.  
> **username_0**: As for your first point @NatWH, I had the following in mind: Calculating the maximum likelihood trees (tree A and B) for both partitions (part A and B) separately. I'd then calculate the likelihoods and perform a incongruence test only for part A data under both tree topologies of tree A and B. Same procedure for part B as well for completeness. So the test is in using the same dataset(s) and the likelihood should be comparable, do you agree?  
> **username_2**: Yes I agree, that shouldn't violate the assumptions of the SH test, although I am just not familiar enough with the SOWH to say anything about that. I would still advise determining the site-specific differences in log-likelihood, though. Its helpful to know the distribution of the signal.  
> **username_0**: 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!  
> **username_2**: Absolutely! I look forward to that in September.  

---

## 1 answers

### Answer by username_0 on 2018-06-12. Score: 2

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.

> **Comments:**
> **username_3**: Thanks for sharing! But perhaps it would be clear if the output is removed, or if the code is explained before the block  
> **username_0**: 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.  
> **username_3**: 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  
> **username_0**: Agreed. I broke it into chunks, hope this improved it!  
> **username_2**: 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.  
> **username_2**: You could also use one of the implementations of model finder (AUTO in RAxML or MFP in IQ-TREE) to find the best-fitting model of substitution.  
> **username_0**: That's true. In fact I've used RAxML for the original data but used non-optimized NJ deliberately to make a more reproducible example and keeping the code simple. In fact, this underlines my original point: Testing incongruence between partitions shouldn't depend (too much) on the tree building algorithm or/and optimizations. If it does, I'd love to read about guidelines, best practices or examples on how to do this properly.  

---

