<jupyter_start><jupyter_text>`SuchLinkedTrees`I didn't want to write this either.  Working with linked treesIf you are interested in studying how two groups of organisms interact (or, rather, have interacted over evolutionary time), you will find yourself with two trees of distinct groups of taxa that are linked by a matrix of interaction observations. This is sometimes called a 'dueling trees' problem.![dueling trees](dueling_trees.png) If the trees happen to have the same number of taxa, and the interaction matrix happens to be a unit matrix, then you can compute the distance matrixfor each of your tres and use the [Mantel test](https://en.wikipedia.org/wiki/Mantel_test) to compare them. However, this is a pretty special case. Hommola *et al.* describe a method extends the Mantel test in this paper here :* [A Permutation Test of Host–Parasite Cospeciation](http://dx.doi.org/10.1093/molbev/msp062). *Molecular Biology and Evolution*, Vol. 26, No. 7. (01 July 2009), pp. 1457-1468, by <NAME>, <NAME>, <NAME>, <NAME><jupyter_code>%load_ext Cython
%pylab inline
from SuchTree import SuchTree
import pandas as pd
import numpy as np
import seaborn
from SuchTree import SuchLinkedTrees, pearson
T1 = SuchTree( 'SuchTree/tests/test.tree' )
T2 = SuchTree( 'http://edhar.genomecenter.ucdavis.edu/~russell/fishpoo/fishpoo2_p200_c2_unique_2_clustalo_fasttree.tree' )
links = pd.read_csv( 'http://edhar.genomecenter.ucdavis.edu/~russell/fishpoo/fishpoo2_p200_c2_host_count_table.tsv', 
                    sep='\t', index_col='Host')
links.index = map( lambda x : x.replace(' ','_'), links.index )<jupyter_output><empty_output><jupyter_text>Let your trees be linked!Now, we create a `SuchLinkedTrees` object, which connects the host and guest trees to the link matrix (`True` or `False` values derived from thecount matrix). This pre-indexes the table for fast access later, so creatingobject takes a little while.<jupyter_code>%%time
SLT = SuchLinkedTrees( T1, T2, links )<jupyter_output>139761682196456 allocating columns in 3446384640
bulding default subset.
bulding default link list.
CPU times: user 9min 41s, sys: 2.54 s, total: 9min 43s
Wall time: 10min<jupyter_text>What's goodies are inside?Let's look at how `SuchLinkedTrees` slices the link data by clade.We'll pick a random clade...<jupyter_code>SLT.TreeB.get_leafs( 7027 )<jupyter_output><empty_output><jupyter_text>This one looks good.Now, we tell our `SuchLinkedTrees` object to subset itself using that clade.The default subset is the root node of the guest tree.<jupyter_code>SLT.subset( 7027 )
print 'subset size  :', SLT.subset_size
print 'subset links :', SLT.subset_n_links
print 'link pairs   :', ( SLT.subset_n_links * ( SLT.subset_n_links -1 ) ) / 2<jupyter_output>subset size  : 89
subset links : 245
link pairs   : 29890<jupyter_text>print SLT.col_idsprint SLT.subset_columnsprint SLT.subset_leafs  Paired distancesWe can test whether or not the guest organisms are likely to have co-diversified with the host organisms by computing the distances between each possible pair oflinks in the link matrix through each of the two trees, and then testing how wellthose two sets of distances correlate with one another.<jupyter_code>result = SLT.linked_distances()
seaborn.jointplot( result['TreeA'], result['TreeB'] )<jupyter_output><empty_output><jupyter_text>Yeesh. Lousy correlation, isn't it?Anyway, thanks to `SuchTree`, it only takes a tiny sliver of a second to compute the 29,890 distances through each of the twotrees. Unfortunately, for $n$ links, the number of link pairs is$$\frac{n (n-1)}{2}$$This is $\mathcal{O}(n^2)$ scaling, which is $bad$. The biggest clade is $$ \frac{ 54327 \times 54326 }{2}= 1,475,684,301 $$link pairs. And that's just one particularly big clade!  Solution : samplingWe can avoid this problem (with some care) using `SuchLinkedTrees.sample_linked_distances()`.<jupyter_code>result_sampled = SLT.sample_linked_distances(sigma=0.05, n=10000, buckets=10)
result_sampled<jupyter_output><empty_output><jupyter_text>Clade `7027` is not particularly big, so we've actually over-sampled it by a little more than $10\times$.It shoudln't much different, though.<jupyter_code>seaborn.jointplot( result_sampled['TreeA'], result_sampled['TreeB'] )<jupyter_output><empty_output><jupyter_text>But, of couse, sampled distributions don't look exactly like the distributionsfrom which they were taken, so let's have a look at 'em.<jupyter_code>seaborn.kdeplot(result['TreeB'])
seaborn.kdeplot(result_sampled['TreeB'][100000:400000])<jupyter_output><empty_output><jupyter_text>Down to businessNow that you've seen the major features of `SuchLinkedTrees`, let's use it to do something useful.How about my main dissertation question? *Are there any co-diversifying bacteria in my fish?*<jupyter_code>import pyprind

p = pyprind.ProgBar( len( list( SLT.TreeB.get_internal_nodes() ) ), monitor=True, title='sampling trees...' )

big_nodes = []
table = {}
for n,node in enumerate( SLT.TreeB.get_internal_nodes() ) :
    p.update()
    SLT.subset( node )
    if SLT.subset_n_links > 4000 :
        big_nodes.append( node )
        result = SLT.sample_linked_distances( sigma=0.05, n=1000, buckets=100)
    else :
        result = SLT.linked_distances()
    table[node] = { 'n_leafs'    : SLT.subset_size, 
                    'n_links'    : SLT.subset_n_links,
                    'n_pairs'    : result['n_pairs'],
                    'n_samples'  : result['n_samples'],
                    'deviatnon_a': result['deviation_a'],
                    'deviation_b': result['deviation_b'],
                    'r'          : pearson( result['TreeA'], result['TreeB'] ) }<jupyter_output>sampling trees...
0%                          100%
[##############################] | ETA: 00:00:00
Total time elapsed: 06:00:03<jupyter_text>If exhaustively measure correlations for clades with fewer than 4000 links,and used sampled distances for clades with more than 4000 links, it takesabout six hours on one core.It can be threaded, but I haven't gotten around to trying that yet.<jupyter_code>C = pd.DataFrame( table ).T
seaborn.jointplot( 'n_links', 'r', data=C )<jupyter_output><empty_output><jupyter_text>Most of this is garbage. Let's focus on the good bits.<jupyter_code>seaborn.jointplot( 'n_links', 'r', data=C.query('n_leafs > 5 and r > 0.05')  )
CC = C.query('n_leafs > 5 and r > 0.05').sort_values('r', ascending=False)
print CC.shape
CC.head()<jupyter_output>(398, 7)<jupyter_text>Correlation smorrilationWell, the only reason I'm using Pearson's $r$ is because it happens to be an $\mathcal{O}(n)$ algorithm. It also assumes that the distributions of the things you are testing are normal.We already saw that they are *not* normal.We really should use a rank-order correlation test like Kendall's $\tau$ because they don't make any assumptions about the distributions, butthese are all $\mathcal{O}(n^2)$ algorithms.Once again, the massive diversity of microbes forces us to really sweat thedetails. So, let's just use Pearson's $r$ to out what the correlations *probably* are, and then find Kendall's $\tau$ for the ones that look interesting.<jupyter_code>from scipy.stats import kendalltau, pearsonr

pearson_p   = {}
kendall_tau = {}
kendall_p   = {}

for n,node in enumerate( CC.index ) :
    SLT.subset(node)
    result = SLT.linked_distances()
    p_r,p_p = pearsonr(   result['TreeA'], result['TreeB'] )
    k_t,k_p = kendalltau( result['TreeA'], result['TreeB'] )
    pearson_p[node]  = p_p
    kendall_tau[node] = k_t
    kendall_p[node]   = k_p

CC['pearson_p'] = pd.Series(pearson_p)
CC['kendall_tau'] = pd.Series(kendall_tau)
CC['kendall_p'] = pd.Series(kendall_p)
CC.head()
seaborn.jointplot( 'n_links', 'kendall_tau', data=CC )<jupyter_output><empty_output><jupyter_text>The Moment of TruthWell, there's one clade up there that has a $\tau$ of 0.394 with a $p$ of 0.040588.Could this be... it?<jupyter_code>SLT.subset( 79047 )
result = SLT.linked_distances()
seaborn.jointplot( result['TreeA'], result['TreeB'] )<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.11039283938339135, 'fraction_numerical': 0.027473893585280955, 'mean_word_length': 4.8679795769511305, 'pattern_counts': {'":': 0, '<': 53, '<?xml version=': 0, '>': 58, 'https://': 1, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '13', 'dir': 'jupyter-structured-clean-dedup', 'id': '458806', 'license': 'BSD-3-Clause', 'n_tokens_mistral': 2825, 'n_tokens_neox': 2585, 'n_words': 837, 'path': 'docs/slides2.ipynb', 'repo_name': 'ryneches/SuchTree', 'size': '8073'}