From dd1322b4ad41dcf59a1ec7d868a773cb42784e66 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Yoshiki=20V=C3=A1zquez=20Baeza?= Date: Tue, 26 Sep 2017 15:39:03 -0700 Subject: [PATCH] ENH: Adds new method `beta_phylogenetic_alt` (#143) This implementation is suggested for large datasets. --- ci/recipe/meta.yaml | 1 + q2_diversity/__init__.py | 13 +- q2_diversity/_beta/__init__.py | 12 +- q2_diversity/_beta/_method.py | 47 ++++++ q2_diversity/plugin_setup.py | 73 ++++++++- q2_diversity/tests/data/crawford.biom | Bin 0 -> 137185 bytes q2_diversity/tests/data/crawford.nwk | 1 + q2_diversity/tests/data/vaw.biom | Bin 0 -> 33800 bytes q2_diversity/tests/data/vaw.nwk | 1 + q2_diversity/tests/test_beta.py | 222 +++++++++++++++++++++++++- 10 files changed, 355 insertions(+), 15 deletions(-) create mode 100644 q2_diversity/tests/data/crawford.biom create mode 100644 q2_diversity/tests/data/crawford.nwk create mode 100644 q2_diversity/tests/data/vaw.biom create mode 100644 q2_diversity/tests/data/vaw.nwk diff --git a/ci/recipe/meta.yaml b/ci/recipe/meta.yaml index 20490046..9096e130 100644 --- a/ci/recipe/meta.yaml +++ b/ci/recipe/meta.yaml @@ -32,6 +32,7 @@ requirements: - ipywidgets - scikit-bio - scikit-learn + - psutil - biom-format >=2.1.5,<2.2.0 - unifrac >=0.9.2,<0.10.0 - qiime2 {{ release }}.* diff --git a/q2_diversity/__init__.py b/q2_diversity/__init__.py index 2515cc09..eaee8fc9 100644 --- a/q2_diversity/__init__.py +++ b/q2_diversity/__init__.py @@ -8,8 +8,9 @@ from ._alpha import (alpha, alpha_phylogenetic, alpha_group_significance, alpha_correlation, alpha_rarefaction) -from ._beta import (beta, beta_phylogenetic, bioenv, beta_group_significance, - beta_correlation, beta_rarefaction) +from ._beta import (beta, beta_phylogenetic, beta_phylogenetic_alt, bioenv, + beta_group_significance, beta_correlation, + beta_rarefaction) from ._ordination import pcoa from ._core_metrics import core_metrics_phylogenetic, core_metrics from ._filter import filter_distance_matrix @@ -19,8 +20,10 @@ __version__ = get_versions()['version'] del get_versions -__all__ = ['beta', 'beta_phylogenetic', 'alpha', 'alpha_phylogenetic', 'pcoa', - 'alpha_group_significance', 'bioenv', 'beta_group_significance', - 'alpha_correlation', 'core_metrics_phylogenetic', 'core_metrics', + +__all__ = ['beta', 'beta_phylogenetic', 'beta_phylogenetic_alt', 'alpha', + 'alpha_phylogenetic', 'pcoa', 'alpha_group_significance', 'bioenv', + 'beta_group_significance', 'alpha_correlation', + 'core_metrics_phylogenetic', 'core_metrics', 'filter_distance_matrix', 'beta_correlation', 'alpha_rarefaction', 'beta_rarefaction'] diff --git a/q2_diversity/_beta/__init__.py b/q2_diversity/_beta/__init__.py index 465fe317..134f2158 100644 --- a/q2_diversity/_beta/__init__.py +++ b/q2_diversity/_beta/__init__.py @@ -6,11 +6,13 @@ # The full license is in the file LICENSE, distributed with this software. # ---------------------------------------------------------------------------- -from ._method import (beta_phylogenetic, beta, phylogenetic_metrics, - non_phylogenetic_metrics, all_metrics) +from ._method import (beta_phylogenetic, beta_phylogenetic_alt, beta, + phylogenetic_metrics, non_phylogenetic_metrics, + all_metrics, phylogenetic_metrics_alt_dict) from ._visualizer import (bioenv, beta_group_significance, beta_correlation, beta_rarefaction) -__all__ = ['beta_phylogenetic', 'beta', 'bioenv', 'beta_group_significance', - 'phylogenetic_metrics', 'non_phylogenetic_metrics', 'all_metrics', - 'beta_correlation', 'beta_rarefaction'] +__all__ = ['beta_phylogenetic', 'beta_phylogenetic_alt', 'beta', 'bioenv', + 'beta_group_significance', 'phylogenetic_metrics', + 'non_phylogenetic_metrics', 'all_metrics', 'beta_correlation', + 'beta_rarefaction', 'phylogenetic_metrics_alt_dict'] diff --git a/q2_diversity/_beta/_method.py b/q2_diversity/_beta/_method.py index 7e616ce8..97293c5f 100644 --- a/q2_diversity/_beta/_method.py +++ b/q2_diversity/_beta/_method.py @@ -12,6 +12,13 @@ import skbio.diversity import skbio.tree import sklearn.metrics +import unifrac +import psutil + +from q2_types.feature_table import BIOMV210Format +from q2_types.tree import NewickFormat + +from functools import partial # We should consider moving these functions to scikit-bio. They're part of @@ -20,6 +27,13 @@ def phylogenetic_metrics(): return {'unweighted_unifrac', 'weighted_unifrac'} +def phylogenetic_metrics_alt_dict(): + return {'unweighted_unifrac': unifrac.unweighted, + 'weighted_unifrac': unifrac.weighted_unnormalized, + 'weighted_normalized_unifrac': unifrac.weighted_normalized, + 'generalized_unifrac': unifrac.generalized} + + def non_phylogenetic_metrics(): return {'cityblock', 'euclidean', 'seuclidean', 'sqeuclidean', 'cosine', 'correlation', 'hamming', 'jaccard', 'chebyshev', 'canberra', @@ -63,6 +77,39 @@ def beta_phylogenetic(table: biom.Table, phylogeny: skbio.TreeNode, return results +def beta_phylogenetic_alt(table: BIOMV210Format, phylogeny: NewickFormat, + metric: str, n_jobs: int=1, + variance_adjusted: bool=False, + alpha=None, + bypass_tips: bool=False) -> skbio.DistanceMatrix: + + metrics = phylogenetic_metrics_alt_dict() + generalized_unifrac = 'generalized_unifrac' + + if metric not in metrics: + raise ValueError("Unknown metric: %s" % metric) + + if alpha is not None and metric != generalized_unifrac: + raise ValueError('The alpha parameter is only allowed when the choice' + ' of metric is generalized_unifrac') + + # this behaviour is undefined, so let's avoid a seg fault + cpus = psutil.cpu_count(logical=False) + if n_jobs > cpus: + raise ValueError('The value of n_jobs cannot exceed the number of ' + 'processors (%d) available in this system.' % cpus) + + if metric == generalized_unifrac: + alpha = 1.0 if alpha is None else alpha + f = partial(metrics[metric], alpha=alpha) + else: + f = metrics[metric] + + # unifrac processes tables and trees should be filenames + return f(str(table), str(phylogeny), threads=n_jobs, + variance_adjusted=variance_adjusted, bypass_tips=bypass_tips) + + def beta(table: biom.Table, metric: str, n_jobs: int=1)-> skbio.DistanceMatrix: if metric not in non_phylogenetic_metrics(): raise ValueError("Unknown metric: %s" % metric) diff --git a/q2_diversity/plugin_setup.py b/q2_diversity/plugin_setup.py index 4918ce07..867f05cf 100644 --- a/q2_diversity/plugin_setup.py +++ b/q2_diversity/plugin_setup.py @@ -7,7 +7,7 @@ # ---------------------------------------------------------------------------- from qiime2.plugin import (Plugin, Str, Properties, MetadataCategory, Choices, - Metadata, Int, Bool, Range) + Metadata, Int, Bool, Range, Float) import q2_diversity from q2_diversity import _alpha as alpha @@ -38,7 +38,19 @@ 'and exploring community alpha and beta diversity through ' 'statistics and visualizations in the context of sample ' 'metadata.'), - short_description='Plugin for exploring community diversity.' + short_description='Plugin for exploring community diversity.', + citation_text=('Unweighted UniFrac: ' + 'Lozupone and Knight 2005 Appl Environ Microbiol; DOI: ' + '10.1128/AEM.71.12.8228-8235.2005.\n' + 'Weighted UniFrac: ' + 'Lozupone et al. 2007 Appl Environ Microbiol; DOI: ' + '10.1128/AEM.01996-06.\n' + 'Variance adjusted UniFrac: ' + 'Chang et al. BMC Bioinformatics 2011 ' + 'https://doi.org/10.1186/1471-2105-12-118.\n' + 'Generalized UniFrac: ' + 'Chen et al. 2012 Bioinformatics; DOI: ' + '10.1093/bioinformatics/bts342') ) plugin.methods.register_function( @@ -46,7 +58,7 @@ inputs={'table': FeatureTable[Frequency], 'phylogeny': Phylogeny[Rooted]}, parameters={'metric': Str % Choices(beta.phylogenetic_metrics()), - 'n_jobs': Int}, + 'n_jobs': Int % Range(1, None)}, outputs=[('distance_matrix', DistanceMatrix % Properties('phylogenetic'))], input_descriptions={ 'table': ('The feature table containing the samples over which beta ' @@ -68,6 +80,61 @@ " for all pairs of samples in a feature table.") ) + +plugin.methods.register_function( + function=q2_diversity.beta_phylogenetic_alt, + inputs={'table': FeatureTable[Frequency], + 'phylogeny': Phylogeny[Rooted]}, + parameters={'metric': Str % Choices(beta.phylogenetic_metrics_alt_dict()), + 'n_jobs': Int, + 'variance_adjusted': Bool, + 'alpha': Float % Range(0, 1, inclusive_end=True), + 'bypass_tips': Bool}, + outputs=[('distance_matrix', DistanceMatrix % Properties('phylogenetic'))], + input_descriptions={ + 'table': ('The feature table containing the samples over which beta ' + 'diversity should be computed.'), + 'phylogeny': ('Phylogenetic tree containing tip identifiers that ' + 'correspond to the feature identifiers in the table. ' + 'This tree can contain tip ids that are not present in ' + 'the table, but all feature ids in the table must be ' + 'present in this tree.') + }, + parameter_descriptions={ + 'metric': 'The beta diversity metric to be computed.', + 'n_jobs': 'The number of workers to use.', + 'variance_adjusted': ('Perform variance adjustment based on Chang et ' + 'al. BMC Bioinformatics 2011. Weights distances ' + 'based on the proportion of the relative ' + 'abundance represented between the samples at a' + ' given node under evaluation.'), + 'alpha': ('This parameter is only used when the choice of metric is ' + 'generalized_unifrac. The value of alpha controls importance' + ' of sample proportions. 1.0 is weighted normalized UniFrac.' + ' 0.0 is close to unweighted UniFrac, but only if the sample' + ' proportions are dichotomized.'), + 'bypass_tips': ('In a bifurcating tree, the tips make up about 50% of ' + 'the nodes in a tree. By ignoring them, specificity ' + 'can be traded for reduced compute time. This has the' + ' effect of collapsing the phylogeny, and is analogous' + ' (in concept) to moving from 99% to 97% OTUs') + }, + output_descriptions={'distance_matrix': 'The resulting distance matrix.'}, + name='Beta diversity (phylogenetic) - High Performance Computation', + description=("Computes a user-specified phylogenetic beta diversity metric" + " for all pairs of samples in a feature table. This " + "implementation is recommended for large datasets, otherwise " + "the results are identical to beta_phylogenetic.\n\n" + "This method is an implementation of the Strided State " + "UniFrac algorithm. Multiple variants of the UniFrac metric " + "are available, including Generalized UniFrac (Chen et al. " + "2012), Variance Adjusted UniFrac (Chang et al. 2011), " + "as well as Weighted normalized and unnormalized UniFrac " + "(Lozupone et al. 2007) and unweighted UniFrac " + "(Lozupone et al. 2005)") +) + + plugin.methods.register_function( function=q2_diversity.beta, inputs={'table': FeatureTable[Frequency]}, diff --git a/q2_diversity/tests/data/crawford.biom b/q2_diversity/tests/data/crawford.biom new file mode 100644 index 0000000000000000000000000000000000000000..fa5bf1bc21e438461861f78ad6b052c2491cd5de GIT binary patch literal 137185 zcmeFa2Y6KF*7iSi5D*lUUKJ?{WHPDrUK9ZVMWluhAW{Mep@b%&BBCHoEC?daj*3c= zP;4}5q9RI(bVMQ2Yw}-va<9qcPS(l$Ue|lR|M#6YbDiVBeXqUe`914tvuE#zX>D4x zDtbrRJ6z7IU_n=2S8?wj^*^b2{ZT5nSJJ<$-v@ulh4&Zm{%%^hE0?RV{NLp&4ZS$d zUpY^j&#PAyxyBvQ|>@338cZ~{G#>N)TKJ8Ga!%cdP2(kFpi8Sh)s&hS@POy9Vu z*pU$fukn{hv&n$?zG{3OCJ%hII_^q3FXy=VXw@W+5CA}u?)jp*I+ecegcGn576Hq%g z(*027ArWySqGFo__74gQ^Mt$dA+9Xm0|VT_wFBI>146n41l03{)(Z%(;|>iE3=2*R z$d5REcZY|&!`v=c0n!6Q-Ju~NGSran4habe4|ll=u^tv05)|Zer43Lnm!3b+;|X(z zrq`2C4-5|qmh%+h^Mr?ogiG>^k{%cw9^{ep-^qGNKxnX(vsh+5FhH1*uJrOa@gA4! zE{gYfLfkdrYp}}EpDJPnh-+9~~r#vMo-W?tk8YaB=ur8`c&ae5`dT@sZga-$> zT%{;JFeD%#Sn8oP>5|KU5Gj8d(gQ<6g2Uzf+79&m9(Q0+sMJq6iuayJxXQC09^&@M z=T;!ysfQ4WuSmM^hKEYLwpab!z|fEYr#)0gytc0(Pk2aJXsD|SohMA%Lx8lys-#Of zL!J7+k90Y|RFTVdeO^ft^5U2l& zl|h$rDcmjRsYkjrl8|62OnuUWq#-)}y#eW=9%&Bpp$$p*NFtqf)rj=K(7>Qzsf5O) zyTb#e;kq#U_|;E%u*VsTn^JsWh}$jQ!qtrRz;KUKZ_QZ`2y#39D{Zz(KYw6Is5mZH zON#ddxWhuFF|;CG;@yG5)tYrrc(7CcHl%xk!rju1Ty060c_b)6&hs$ofdSG15?x`xYWa=q|3==pmDj{lP=}=gh{?SupS!Z4v>25m{|{V>g_Ss!@|Nn zQqE4ShX%Vna-PnthsXy=d>81xekb!?sK?`aoOJ;Z0@o9C9*;XH%q@B8O7WhMK({+A z+_#R>e&h}dloPpJ-6-De4hZ&0!+na*?+yzP$K~qIy39I`-h*@*=R*Reo_mrmtw*ME zm#Y`)l17;zU9R4w2Zp;toq4Sf>Cz4Zobfh-bm{Nm0iyS1JyezhQqIUsy7XYldq37? z9CP|_6zLvmKcO;c^e5dJf1Lh1fb_r+Pk4ZwXCUdGFsWyWkA|L}_wdjl0eG6@Ju)o| z*C2|Qe8~vvaSf(;8UI5Ag?9+)o?v&FQ_pFu<@EmU2@D7dmwJw&_&`~Wg-g2~M!L+` zAt6$p;hA(9engLDT_$Vi{Bf)ct+VbL!FpJLOffE3Jn0^p-K9CZMw0FhlbKo0KZ^7q znJ}DoIGS|HOK^~U?ikXgeTF;dPas{^jp0uJ8cVvghY%+(<4BirILMiI#zXhzSE@5O zFvvB5;@vX;I^*+1(mjEJ0Zu(XL%Ou{U}rrxDU%)$;><6TNq0N#N8+cDF5_N^!}~1h z?vOwMcDbG-JmP92_e3@G|L=$6#ljGlO(l`~*4mHj{K&AG)1(_zLOn@KBk`U9MTI zho$9xHt9h!%{mLFIiyQ@!kq6fb6F1!4t3Ub^GJ7QUuT{CD(Qit0kVK|x#mMp&r7IV za_@SL<7HfO+W+gU%h#^}Dd!ukhs*dZ?f*^Er95FyITw&F{Y$1Em+LLkrJu@#;9 z?klHL{h>jwMHDafBNKqj^)~6UK9m{WIFOV?-J;~&z04)9O7Ea z@!@Wd^Zjud>(U;ayf0^6)=N%aR*){kl6PE!EZ-N4d;h5YL9d&oUth@Qy}mGF%rNiN zlil?N=VW?)A-DXw2v%CobvWlbMcOri+VQcNAhJ8(y?*D*li_^D{q^`-trBI8X6TFF(Rs8?Y?82`0UD8Nn(^YulzZk z_fGHmB4eYx7lUf|lau}z&R5p&eEGcRYa;)3uI1$S^2QB}@FpTV&LIFSpf}0=~8L$*P6-8hi^Kw6yE^Z)Kz{ohB= z<3MFpL2=9vS|i%$qlfm3j*Nu7&+a{K9LVdt?)h^z{X1<{dYQSpws7vBx$~r*uCS}J_YZoWTwXtyt6xDKmA{4ymQ7z{k>S_?`;uTgSqp0=e1ITFS{4J*j za%v!_26Ad3rv`Fr;6I=S($=AQd<&LOD*S6GU57f?I{(#6*O#Ni|5em~jyJm3&Up=o zeiKmSxrY7!+Eb1~0*7XP9#`SLb#o?u zR_XEmPpljN$H98_8m}0Yx_IIG&)bihP+kfW-mlS5bADTKy!G-MUFPqq_V*9_&b*cQ z@UA}&k6wCo?b>%|k9lHS+_^QY&&~Yhy|C?FFK!v%{_^_O?Jqxi-gD$k>h=Q@-tVx# z>8Qi27jCK5xyOc;>wCAjv24w@_7@tjI{V!hw?90!AiCn2)`7*UcRs&$?A-CU2laYt zz!wX<1=cySw8QNm>y0{m<;N5AwqCzkyH?{vN7m2pnbN+)x|FW_tFODZzWUko=PQ(X zC-)1Zdky)i(9?O{MOy6LweaSmu-k1L*Z85}vxh!=b@s!a{qa3xR&VH?`s4NCXHP|@ zCS1Dt%c3oZM_s+OZong--E(C1*b6rooY+yP_sB+1zqjwe@4wxB?bqY&KbyD0Gih1* z>u2_b{(f?Mr7xGB@0)bsM*R=I>-ORKzpMORXIJ^@E6&dOv;Km_Qq@~EKJfm5Z5{qT z@ksmg*M4pEQ|OcuOX?LW^83n^pW9x0bIg$eGvdD7cJ#@?E3b6gS#nhWQ+0OrtbKjX z(N&wPZ#nvPbWHuw@tdpmAOBhPfI7XU1l?z^T)pe*uXaq?z2VT$U&mG5e{M#F17TN6 zruP53^}%!B?y9_ZRfDj0KWuqv(*h>g#^zBU;O>of&;J4*|4*2{g(=j zd-_WLN{?Q8cG0-`b3Pio;$3ZsGaV&D#ns_`QCSD^K>0T~Ob( z|Kw-o`VNe}(kJDs&NsXFoDiD7?vCE;+~a!v(K+(=4+TSaRN44BMM&I7fRenpsr$_g?=cf~K#iYOQsWkHGr^|)7+jY3QZC|@n;bTfa(EhD+VZ+xo zdiP9&3i~cMZNFvMsTrO_@#`8^I=nGC_oOD3tG(BzdA-TonzRZm(X2*9ziRdWXz|LN zHW%({IpL0uuT?Es?Tw>zmzBTg-8R+U@A&c?KfP1u^1B_H_Z`@x#2Y_O4*7gpg}=|V zsa|j2b03dbQvcISZg72)4>+cj(Jm^#Cs`=G~t)80IIG-TA89SzUT`5{lm3$Imga3Z|! z=xq=7y7>I+^{4Vwy8i8e4axIj#+|CZ@!HNxlRM?j9o@D=K$)#=9{8-dm$IBtPU&&{ zf!PgO%vv@3tOaws`Z%H}<@A z^7FT=eKhdVvGw-MYCmFm;X~7w6=-;-W%Zc1v|MJ4-0(ZH^2wg-P4H|74JaZ|rKc&yR321!+Bj{d>3@Z#jq z^~V}ky0*1?ql?e29CxgK<7?aMPPjOAL3eID6u*$WKzNi+o?;qQ$keYqGBFF7mC}nD*^A=S%H>&V>8l zA+g#~%|mJJqV18icPnl){6GFMuA(^C0spJlp^AI2D>=V&rCmd+>?Qdhxc<~eel~Vq z&LteZJe9@ur#{}GfA@c^>rg9<8vLYY#uDjuC?7~hJLhjXHIP#SIW>?|135L2Qv?42 zHIUu`&aL&+7ToW}AAU_A9D||tAg9+~e7bDF;_M&eJWJ=iRF~&L z-OfG@`uFtbFN0*`i||bQ26^0R&!g*4+8XZ;k=h#f$}u8vyVkl(&Zu4AdfuHdS@nGwowo+=hhGBm3|(1 z4pr(?*H@~`bHjm7ecqK8?u(a)zk{9Uq3>orJRr!~pQZ$KEx$ZZ940&I$TN>=FW>o{ z^T_%2u`6G^JX{+f8*VtycRDYPm)$euN$kw~2Fb2XvVlhWbBG!*Pxl8o&(rJU7d~CK z^9YlAbDpPkUaHH!4zkfu=6!=?AA(Tld1!qMMdO2m<#}1DPrXe~b+@F~*>?fYAo$KN z`I0BSv)VVv?Qx#Rt?VsLdU*mpZfUwMR~6RfscFN!kjLXxX{LRHoP9W)_E4SABkj%EXQf6aUD6`?(w~!k^&{2eJP(~| z-ymmOE;+wmKGk?nuvCw<=ZEP0?hx4~Ms$5_#uqOS+B^Fd)uDKgY=hFa2|JfaFDa_jx?bRx;#cM8!lzB zZ;eOZl6T?mWfs==jd->vzGy0p9td1;q!31bX-PYDv13Gr-vwr&T6hh99Zt)}+g9 z;cT;&dEcN=*|GO_3jhXf%Kpdx3kShrhS7v!7?vMc^;#9+2_E!pG_yyJ(9(&_6?HFC!|57?;GUn zUm?NX{fZu^^GFhdoc)uYfUe_PkkAH8$LPwsTXy=8@#;xFk8IE3)PFb9W&fdYXa9_+ zNOx-A*$=EcpI@>bEE7x*igyRd20K!3Jz1C5?~L=kNSA$9f}Q;&GVL2AyX3hgz7NGa zyVE%RD1vm^15vg>OW!xhS3f}^fs*;m`v%GWF3xz^kIo~t=AFl*NSFPS!kzs&`m-*o zvp?Se(xukC`(O=ZUB(AzUz2Fko&DjQ{fBtpAgO=v{)>9~N85=rt~%p0wuSKZYp*Wl z8A9imdDgpM%TVaP`tgLw0?w0p-=NSyIUtLDgTjIX14SR6>HL8nXFteT(q-I`4dt@f zHz+(@rk>3E1_g!$hdAq*c$`Pqp|VGk(>ycn8|3VV$G2~g=2r#>86HA2?;9k$`Z(?9SvtS$3l-+%_c_w#dxUp= z_laX$0DL9)Y<)WdW-ztmfxY@_CS ziF8??J6{Vj?HeSkDDSv2BU5~cAZM{}kSwd6eW6~VcxQhaXFsi3tjqT!XPrD7y7J0? zkO9&sGwmBBgJiI?4xh{Ckx|uIx6dP8mdVo4Gw&NDTrx1F?;GSBA0*z{r|mU5kCb0_ z?n~b{NXr>0%k2R9Lih&9%QWuHo0;|v3Y7(c)XxHncjj@Y{lAq-mxYV)E+k!6QQ^+_ zO5Qg}##`_HLvK^O?29G?Q09Guc^y?%g7)e^hL2bYE|nY{B7u?WMR#Bko7UmBf1t#YGHGJL5k1 z_uyT4GkmTh{{+7{-i7zk_{eCdR-y+Fjvw+f^rCnd?zqT-5mEgH3!-ShK)3k!)A7RJ zDmr#ZbY%PpZ$W=0U$~nOj){}aiu!qL@Bn-*pYW#p!ulJ0&6jX`mETF%JRT%n^LPll z#c>$E)&p}KA)h&t;9DGjz}I$Swa=r7Q(omdM!Iqxhi-A5fUo5euC(?W6WuRLQgjmj zop={6pYJ5-6df=5@%c_$JPp4f-i51kU{LMQGx$4m{7JfUoP}<2`~~0Q_#3{(agL5> zj`Q&E!n<&!bEXd;u4MWg$&WA2Y4HERS8nCI0Nvud2;XXdm&oV#cNxCc5BH-h@Uv@w zR}rV>*YFo7Kd|xmGfeb zvZO0VIntG*Jamhr0(^_3BKgcw3BJWq8Gd&Cr3&JdSGlT^u3Yy)x473#d;6*57cbPkknj3C?bgDR=XXp$QSCyf>v_h!zP?AVDl6x9-8S#S)c0;(AGdVOy!B7@x^LyY z{mU*~-alwV)08Jael>LaqOETp-+AHQ>-)M)*swe3NUfMTjpwi4(R*wCPS+>&etJXa z!__8MoA}oD{e#{g<5|1*V6&k$SI<5_@BEzZJFf01wspwy0-dj4Zm@Jr_q9E1t(^Ge z#VZ%?i!E^dvFq-&%_j6XSEF&2QM2ps=z09(*wD3K9c+HCX5*TTt6iNvcE|YR8)6$@ zk6yd6L6eF#ip?&tqxtax#}{8%bs=EPng-2g9EojQcVp#^HD=Cjy=}<}_k~9;?6~l2 zmz1%ktIzCk;&O@O(=I&t-nf+ecUIrg`FPh8`;WhO;qxvVj*QqjVC%{qo)g_Jl)Lce zg~4NHH0Zjp`o?+NcQ@{Radpg^&(6hfSroS=F7C(!N9ydKfBu#8bI;E{|7y9dSGN{D z(JD4_?X-4_0>*!RsnEXa8UL_%#kDd?A#b!^6uYnb!6%k?c&*VNR|3a>y?4c}>rQZw zD$W7JF4dfP`*i08G48BRK_^Xiz1}F|#IC_IF&MD=uXG+XS?Ju9fpML9n=Sww5 zU8=cpp%b;F!}SG?{y5!vsB-|0T(N!cii4L5?S92kM!)JUP3jFfU!xOo-t#s(vAn~t z-?753oVa!K8z>*ojqJT^=Z2R) z=yxTm{LR>uF~4j-+i6TggsyqYdxZC7avBL3ntrkw96eJfmlAbjDgH~WV_RdG?VEe|b{|HU2tc0uUy z1)&$uhR5{X_2z3AA2_*wn*MRd#r`QNUOP2Ut;5F#j_&ovs|{!5b{A>2V(G!jlPBkX zxp=cKZOZo*#)Q8f{HpGEaj7HPpB&Y% zXX`^Z-`e z_rJP(-MJ0r%bXwma)-5LSKkv{?3bTYe(e2IgVLiW^@#eS;^nGWKMAe0+x_tPVq-Sm zf27cF2YR1uJ@L=h+tyC#J~b}laIMR?mv8v3(%%a%<_SDjruyrrZpK`gnD^$kU+Yi2 z@m0g+J)&Y3{CRf6g12t;y6{$&pI866cHPex*CxEU;;q2jJsX!egu6&^n*R#q_B7ef_}Rt@a=JrNMXWdL5r| zq5S7_660crfB49*PhYtmP%LWx{0UW$)W4cCV8r-!>u!u+@z#Yp3n$&k{X%wkpx=QZ z_uu#gU(MXPy&vLyR9+q7|NnY|Z*;H7K7INm=RR7ta-$;OKJ-$NBIla)xZ~Tg&s-~2 z{A{4}Z)LP^nOCO2f4y%+Az!Vf?O*zDyl=nT9ci`ae56jt|Htp!AKmL%M4vvtKDK+S zEBDh)UMw(e<%{2UZF=m%#sAFPdo6b0xG_N+a+mFWSDRd}FIyLR+uE>W+TCd=@0Cn@ z*S284^FF|(5n6gg#$!fAjE)%^Gh~daSYEiAsPNx>&z%3XPtJvc>vITt9s1w@{6?N< z;?o$MOa9KsdfP(gUexiwN%~(4PkU}7cc%SoTk8VZI?d*^OdQ2R)G2&cuhzRYV`;9& zY5n>depUzPXgM{IQv*3QkW&LWHIP#SIW>?|1OM05fbaPo{J^BEJnOopelM`DTk3U0 z_4&G`e!bN7B44-6hnmsz==J6Nr0aD`U4Lr4UZ>RGe=QFEeyHW(>&yDBTbIwm;k)kZ z{0^<>)#Lg4auJkQec|X7JxKmhvcZGBw?b(BXq@mg9UL7uB6?Wb{S^AS#c&?sa{dVw z2^bnVB0kRdT#%N7?*q}#)z9JI&Gfpr`kW`dZmj3wJe8n)XrA;w5zU9@N$(3$UGsD= z&SP8WY6LH#}g{ym~ z_#x3lV!>hC;T{&8jt{knQTO4iRTO4i4 zXO4&ATO5zTXO4EHE61axD@S|i7Dos87Dq?&nd33|7Dp%e%+Z;2<>*4Xay$;*;&=kS z#nF{~=6Djm#nBCZcH`Sqh_lAG?(nrfl&=SLi%;hPi?0`sSD*RPdO`XHe(pDYs9f4_ zWU1(VX?ZnIeMwi2Naz+vKlobDJWfZ!x5lylbUcq^1K{&GHV`^>agC#LX@16&t{f9cSB{C$t^7O#-{P1A zU&kfkXc;>uPQHc>j**+fBL|{^Pez=6kCFD4F%CDvPeGi1uMsYaOS^gZS;T4kRL#yqcfaNmq_HNLP+Gp<5ga;9DGT zk2DpQhK=xEb@2Z_`vu~ZSGm38-@~^!wvf*pTj5(A z+u$?DcG8t&2kFZ319XdHCwz-z7x~Pw8@|P{2fkl_Nq?@+?|G)Zh`ST-!sUyLAL4TT z2wm$#`SwA#_T*7Kf~A06+Yi}C+B+qFLXSAZ`x13^qUOV`+tS6<<@*1AYFNW zgKqKs4&TbxLGn3Yhu~{{@b{;~@crhcEU!l%K^$`>k*-{SK)1M#!ne4Nko#2rz}mp-Sc?aW#?>GO)rrOzp6}%U%s*&FZFrFVyH*SS9+ZM8l}%GGM_%DXz}UuidufYUZc+~s?U7- z{G!FD&oBD%WjSu@bBsEV3SaY>VX^o%PWGqFl?Kjy`uw8Br_V8J`I%3jXH=j0^tnch zPoHbla_P9C&oip7{Y#%?R9)MlKEG&j=yQx(PUg_(8Sg;*5Dr;~$nSg6X@Wl2sCkfi zTkrmr_hcVm!!E>X`Va310DkMpg#^f^bZ7vb}Zi;vUi9GO#}Z?rh|xkigqpKny3 z^QF%@YPooQqR%;Mxinw;e52}`FMY02bfmR{HKD^7}K+uW?S^Vq&8;%t7|P^!Z0^e@_1V;(QHN zpMxx!&YvMp8>&7Jsrl2q>T{5)YhLyFM~hdVd(`&F{ZOBORG-^tSmwO?zNa{F;gkn$ zpY=#rj{2l4M+4{9CLLfUAZ2EZgF*jul?Q{A39T{e0TzJT5poa zrp}&0k+D(nBNCz`y297`Q0^z8Tio5?Tlsv7e9mWg_*yP*mp$NXxwKvOBwgEOFVdBx zH*|}m4}2{Lj}sB_t#P6+9na%LBz$cj!XfLHw7JQzS-Ses=LnZC&Y8oa;46pl^^c7h z5H0scdi!aA^7+2U0q`x(f%N&3FW=`olUFny&zw(_@ANSm+kVIQUu)?nmR{Tm5JP9nbA+B7B~2pMkD?$}x#_<(LfJ;+O*8%EzUsFt%~MkmXMG=j0lso8-!$kJ-;3~dJd!+hjExx)r3=-Fv}xZp9dV`c zcj0U~EG~L*%uxB#>KFbJ;&dKzcq4`k8#5Ru)CKL!6esPlP1NA%NDb5V>B{>Gbc=Tud@G-`$>)5|fp6t=F8u8BIS+Bn^(yJgH6Oag^%{JO>vi&(>kar8*PHOQ zygF_yAYJqG7U{~d5W2;&2)@PfHu=o47{0~v4t(ZVLb`G+C0#j|LAN-T!`E`~__qSS zHU7Oz$Mg929(*4E-iNMy%JBi|%JCs|i{m5sRz5x^pYyR2zSaYef1kkD@@l*Klyv3z zjCAE#1>NHK9KOZz1^LYJC47rxHGJk+L%MRTC0#kzLAN-*f^TuGC!aYI;aeOV;QNhV zSzf=~h&bg{u1%yX*JkJz*VphZu5ZX^u5aO6T;IXh@@hN#o^<8dLb`Hng>G?dgKu$c zC!aZXz_&PlfX^H|Nmq_tq$|g6=oZHw_!h@r@|ojD_!h@L_}R6ypAe_K%Jnnp%JmC$ zi)%l8i|beNnd<<2i|aS|T3&5uzmu*U2T50sL(nab!|*MRBjhtj5`2r}5BSV+lyv1d zM!Iqwhi-A4fNyb}B%e7>!M8Y0!)J~&q$|gtq$|f+=oZIc@GXwN$!Csp@GXw>@U!bL z$%s>4<@$$o<+=df;<^am;<`jWb6tjSab1C*9oJRFG1oQHmFqfmi|Yn_i|Z!&%ykRC z#dRBgc3dfl)9+uBuPpDMNJX4=KY{+kfflds_ha#5XDi>oyT`nFa;JUqvv~8u*ZkZm4&bUNZMce z?@W3zL7H$m#3`5Ntvu<B@0G zbc>@Jd@DZ>kk9$44qxkm`&AA2T3#I=YLc$`s71PRJP6(5cnH45QJZ|`r~}{Ps0-h3 ze8}>5r2xb!uX4FbSFS+l7MBOU#T7(8a|OeB5I$I+D*EkHNP%I>Bd-&ZH|x7t)pEap)Gu6YwpL zuH-YvlkhE$Zt$7oDbkgrJL$^N1G>f06TZdKi+twj4d3GE17F7#83(ewt`&hebc?GWe2XiJeCFy8-{Kko-;XQH`Dh^Gtly8K;VZZDJq_LB8wB6t8%#d)4S{d* z4TZ1e(s?O{bj{N+(v@R4bc-VvzQqwoK68wKZ*j!KXO5AiE5|6(m18t?i(?FYiz9)2 z<`@g#;ur^?ImVN&91}=aj)~AMj%VOo9Fxdrj>+&Xjw$e&<5|*`<2llm<9X;7$5i+h z#|!YS`&g&J*ZwW*9KZWc_19GScQ+BI_X(wqZ@&9oecynlBTnl@`25ZnANLY`trz8d z8M?(e1HQ#MlYGwCEAX{k{QY4Td@Yxb1G7oje9a+UIp#vQIOf5(I9?^6Ip)K+I9`LV z`LyPt*AeG8Ubh?;KR9MswCrLeOU^goYker+o6s%31@N_e{C)K;___`jz6|&METrQl zPd?vijf>#(_tm$dYk4(Ki%D0Gcc5DwOW<31T1r0WX&HR22X1f6;b)hp6^LW5cS%>S z_n=!`@58saJ|LgDK7?;^eFR_2tNrL>(ltLTNmq_fpj#ZD!nZg+BcC}|!M8X*htC{e zkggnGlCB)9p<5hl;9DGP$!Csr@GXw7;4{a1(v>5TbmiCp-Qw5?U(3P$eiMAF-*2Yl zx!->cpZonc(3MX)z9n5bzJqRYd=KBs#}@KAA6wyTJ#fF@2H&sWXStrL5!~S_c!e_pH&@H~7;A{D~9sUg8 zYKOnj@!Ssg!{>JRD|9Wdw!;IYYo30CZgKn$-^$ZL@;Oh3;A=f_J3I_u%d6w`5z>_- ziFD=o1G>d=6u!lAjC|%e4&UN90iQWelCB)5NLP;2&@GNL@GXu%$!CtU@GXwN;Ahuw z{zjbDZ_dHj`cS^}&@H}X_!i$kG=?>DBqcC)fqX>M9qbT{z zaVLC>qZs_`@>CpgcylS&U8F15-Ow$r67aQLJg=67Z_TUs(D6L4-V0y*z2wLHd|N5# zc=K_TCS5toK({!`!ng8M4!(80q&yv;#v#uI<5Ecl_< z^5Y-p+Yh-C;_&8D&dLxh&MNRN&Z^{dJG>9R=9||K_ru4VuRm2IUEAFQq$@{t=oUu} z_!dV^@|mL+e2e2j_}TTRhY;sC?qyhy*M^TbU!Lkfu=wi2xA+3cXFfN4i!Ts9-h6rT zkgjI{VM&sknifF%TB=<`y2bSje2Z%m`OGyLzQr{KKHhx!d6sm|&vT?J z$MeuFj;Zi1ju*&hj%n~Mju+u)m!Ii~)BY)38TO@m2|nI@>$8_3SbQ_!TYNLgXTDe9 zTYR(NpTnt{iWWt{iVdw>TES z*K+VU{T6&}2g2c7r$B{jQbc^GC_*xFm#|QAOe0)g9b3Q(T&-wTmI^KNwSV_8a zd;;C#_!Pdy@fm!JV-+3G9G}D2`+jBq@7|%G?3mTpyNil<3zeR7KkgUwIg%g0Fqi8~ z_;~Y;W2+%pd~4ua?QJdj+}_r~*ZSf1_7!};aXS4vhTjG;>k)@HU)~Z)*Su|jZgFjd zujS(Qw+X(Lx6O1s=k06w+2!pUI*;V7S^S7G@-Iw8JNvxFeG4CNzPx=0!Q%TKzLmEv z1c95=l`T@Gdu@k<*bU#}*aIJLK90SlE60zd zE5|te_*OmKAfM~uCVZr|fN${?B%k>T!MFJCfUo6} zb)dIG3zM$xw+QLVQ53qxaVLC>qZoYa_sHV#^}cNBhh5^G)eP3Fa?7l@+;<_)?{~=z zarzU=-H7wcmw%kA1bnR*W~}>%M%#;k|Dw zZbWQ`e`YBdamuS)A*3r;D0GV}48E3&`%gH0?YG>2>e2Dsf9k{M{?hB`a4&|5*b_*%oa_}Y-qd~M-dd=JBC zzDGz`zILQ5-=okij`r|%9Fq3dBC3B>tX}xRH0tU=al(^jTt~#|d@Fg&5T{R8J%%`) z&m?~t<1jIILY&r<@ID$J867t)I#zBziXYM$zUD!=<01z}L^*peJ5y&D_!j@;bbjt{ zPr$eOTUR=s``eT7{W!DSU%wmTG!Md=<^K9lAx_(c=D9m`i?auOtta8k@Vj+S^0~eC zg0Fet{@5G7mP`9%AJVnGMv$%?eW6<%k?<{!e&jPp6nu-LKYY!DaAY|z4?vvqD%U{L zl`9&$#q~6Ni)#@1%rzLk#Wf_G&mT(X=k^gpK7anOY(9TD{OsyA7I9p!ainX#j(~2} zYdn0bUPqG8T%+Jy^*S29mRH;T7}B+UB#^EgW1(9dUl`pE*8; zZ*i=I&m5nSt{k6|t{k61w>Vb8w>Um0pEH9uK~PZ;c1*>3ALw65;cBumL)AY$RPdHbJ*IHp92_^ELULpKsu6J@9z&E&S~A z^Bv-t>wD6bYYTLXYb$&$7w2aid@DcO>3Gi14)~m(AD}bGPSTZQ7j%nbH+(BUd*EBw ztM|g!`WFuAAHID=^m4W9N5tuM?X>S7zBs*H?b?Sptyk$se&>si`w71CE9cM9EzV!y zTb%pJ=XUrjd@UEZ!vpZO9)vH;dFeOAY58@W_?>j+Itbn3It1V1I!r!u9f5CgCBfJ7 z>U{JE>B@1Gbmce(-QqY7-{LqyK69LeZ*iQ0pIv@VBhFekoPn?Pp?rTrxA@M&xA^`d zpZWfVZ}FXj@5h(rJary%T7J!2GU>|o4|I#`0(^_>BKgd93BJX38NQZR+us$^mE$Vu z%5e?4#c>_JmV?Lh8}P01{3adG@{z6_`Jr1J1>jrlsv!ByQ3$@paR>bD`b%NNW%qnm5yV-~XBCBS z@!kpF;w?r#^A?A1@!kdBk2lNt>2Abnel?FJNY^};gl=)&17FL<{pw!$R=+Am$8*0b z4WIi}8R*PWmUQJP2i@W*58uj91@bvR72#_=aKEYqU(2igsxs-Ck1C`qM^)$+$9?cE zj{C`Hj%x5NjtAiTwX-beq3VcJUgfGmx^mToZgJIuZ*e_HK65<;-{PtbU(2icsYAMQ z)FoXx0-#$QZuk~QAo5Gbma(zZgGUcw>ZMdXO4RCEspx|{n}ZU z<5>g5S?jrm@U=dauMu>MuQ7a!uL=3g*A%|R*9^XvOWRv>(lt*lNLP-Q&@GNu@GXwk zB{j4>B`X#y2bG*e2b$!eCvB`2l(2aWPJAf-l@Otx?CL*=l4Cf zGu=AjU4T7?IIS1q^NWj*>jYoxML9b|w>Z1Nw>TdspY!zud@UExLtWu(JxIIDa{PW0 zaaw-OUpLZ~>nZ3KS9kbYE?y_}fN!l6deZT{PUr<+*KxwtF*a(X{8N;&)rarugUi($ zaazB^l^!RV?E_!Al`jIi#n+e4D|{J#r;dbgwaUW+J|$7`_5K9m^WT5Ohu9x+ zS`XfKnZBQiZ=aF@i1S;|Iva_4?>%rD_(1p;Z!~aSelSaSbP*xnkj4TygNT;~IfD?GM70;knCr_{yz(BcWS-qu^V7 zqv31$gfHX$nz*lX4C3_u5XoQ0I1H%?h|BK2&asH|;cL0L zKfVB8>p{k!EZ1Yx5U1tW{`ey4+Fqwax42${Z*jd$K6A~0Z*k3ppB>jLh-0o(DK(H{ffzc%FX~ zzV>tB%5Wd`0ym5q!UYza?tGfT-9J z1Eako+}rRi&c*Pp{Jldy=WhvottakpOW|u?gwwn3Tt>R~ujQm`zE(iDINpVCalA)9 zbG#4V;`jhQb9_j;a(qO(a(oQk;#dja;`oGo=J*u8#qk+@=2%6#a(qs@a(n^Z;`kE2 z#j%=v=2!#Y;#dpcZ(PW79AAexzwg!Qak7g13cl8d@~ww%@g>67^6~hv0lqapY^39P zeAq<3GfzbiiW;h412)6ga%;Z6CS7^Hfo}193*XAucjR-vzK5^%!Q;Rd_*z~a2ey)~ z9NS1&j_uGbjveqVjvvTpj-Bu=j$QEm+Fh3O(Qd?P{)H>Ux_1wJtqB@5iy2X=3pT|6Zz_)mg((%l53_kyUa2z^woFH8}PC~ahPQka@*J<*( zeVu`?dF657PxxA1ZC__e*L?g%x^nyt-QqY0-{Lq=K650)w>bWR?>7#l&ntfG(hG=F zUgf$-x^i8DZgE|PujS%)b_KrG&aTq&+|I7S=XQ1-I&<70T{&(-w>WOWxAJqFe9lh_ ze60s=XQ}Z0@{{E_?aJ?6j=GpD7wO8C8@k1n2fmhz^OG08m7jccJm)7re9lh+=*&@& zbmb@n-Qu_dzLlTCV0}w>Zj>&m3jpTO8%!`{gOiepDWD`h8Zo zvfLlG0^cY3U0^s{` zWtlfO;;g&{!dGtP^FX)wg5XDn%vkgj=a3f|-v%NOTd=W7jL>qGh4K)3kX!ngADF!`LPN8oGy@I2lQzLrb#^eE|?r}m^P zM+fK@M@RS;$7AF(M<@6eM`!qc{VU7yqYL7cSGgW1UAdlsZgF*mZ*e_IK67=0Z*e^Z zKRd4Oh_lAU9`Kc0`FcXP_1@U?!FcMx=ocQAa5 zcL@2Mx1sR0Ts&^Xz}IqV-iDE`c^gi;a>PQnIO5=293#kQj(GSM$4K~^cj3r#+!%#8 ztACA#ul1pPW1w4n3GlUiJb#UaZ_Quh=y;yL#>40Dn-idGc{NWHNmq_%pj#Z1;9Gf` zOg`sn3Vf{xZg0=R_sdh3>!;@s$6U{ou3S^0TU;-|*K%=wrop%J^CBJ3`I!#CRN7?> zXTQ7d&c+qPGTb!O)AbVK{Mujl^c#kHx?YB_+{!lty2UpWzLmFE$mhJxg0J<%d7BMi z%cbMQ9MZM@%_Utq=0UeOUWJc0m*lC-qoMMtU%7?D`|n~chRz)Ckggm{ zpj#YE;alxz8Ts6Pmc!S2;P$fuzF)u3a(sRlaen!H;E)7PU}s0GsJCx zulI!sck@7Z@5sP_v~q8RZ}D$}Z{>9}`JC6U;cNLguiwDedw+p_-x0`(C+XLU?+Y4XI zrFr_1bj{N~(v{;U=oZJ%@GXvC$Y+lI@GXvC;WNhp(v{;k(v{ ze2e1?_ z`W)fP@b4g;hOgYpcLut}_a}Ud?<{?;bWRZ?(S*a{{ z*K*2ukl}fe0^~DCLHHI&A^6O32kFXDm~`bR0^Q;$3g63G9oIi*}3k8Z`D^R_*Q+D zCZFr841BA;%EIURDo483S9#L4zA8YsI4Z)oI4Y6P9F^f)997^mM^)05<37@r<9_HC zM>Y6b4xX1EfUoUAIDGBFnYXIb@jP$UfY0+*P3X#}9JNSSjt8Mz91p>_@==?7&PN^i zS`R#L)rFtkcprc`Yn|?fuiVNP2;JiIz_<8<$Y;J__!eIXd@Yxb`=O+3p2A31j&SG} zM?Lr!M}7GEy)&)7$xjma9i#z#&4c7g)~~*OmVEm$HAGx?&m%QLoYsr*`JFF5t}%SA z7v*dM-QsKt-{NdWKIf}Bd@UEZ#}@Fl9wc8`uH#!GPRD)W%kuu3R*3T(*E8H-(;B|T z+XlYH+m?LheHgyQ`v`o^gSO9hq-*qM zXU7$RIOgh0x^hKAx48Pjx45FnXRiM6Ev^CZwY=JY29mD%i6&h+o`!C541#ZQ3?`pB zhQPNthQiM-KQV~QZrmG&IBVP+4&UO9g|Fr1^-~;tYyC8Wj_37LJbYe1jfAfG)Vz%% zUGp{?y2UXDzLmEG@;Pr~;cGqcd^HZfmRH-~c+!<)0_n;z5xT|k416sI=VKClD<6~T zc+STZ_?(Ytp(~$qJV&~6JP+OCmyh?EnpZsgB8s__5Xg=a}{p#e;FV44N%WH`9+ZQ84oIVo#I^wk6qL)B2FS`o}4z z>utpOar(u@$1R4>obN!lIG4b;IG2*o?RFV_Ef@F4dHx)Jc6t5+aXL>4Uxw#Ezl6_xtD#$bYv5aaYsqK6b?_~|ui$IBbe>;N zy6||f7$%agyfar!&Ga8X>^5cCb= z{PLDzAD3_8YyBwichD`~@8Mg#Tgd0UZG~^;Z5w>8FKhhSjyNs9=5Yt<%JlMze!h)bEGTBdFU2LGJK2UAM%;w0(^_( zB7Ek!M7nZZCS5tMK({!q!nZiC!ME1G*WqhBlKIDPy_&xMy@5DgKc$V&zU$t;_3ur@ zX}t)a-}&O>Zc%xJ)Aou#hU`(;!P!=^X9@AXDuJEk8;7+`Vh`6znA7l zob@|P9{8Gn&1YWd7GFO27GHkynXdqRi?1MjEtmGMLZoZEyn}S*C=A`=C<0&0!Ph~G z!q;_*aAf#h?@l^i^5pZKR#^)6C17GWj+vC0Pvulr~5NEZ=((pC!n#VHGm7n>_!ngR! zk?akPie934nkj*g@&$76=x3A)AC8NS8Wg?#3F9KOZ(1bpV}O1knrNxJfN zgKlv=1z*d-<9c`aIxY!^Z++t|QF_qvJg)bIk2l}A-V1{ADMxS8m7@=Iiz5QQm5;vU zb3P*BYdvuP?+0JULFxZluG^y!$6WnMSFQojEv|v^Ev{(zS}x(ra=v;RaoNpRgAkY9 zd^H$xT3?dC?k(f{rak-EIvGl%LHhmJeizJ(*p^}BvA z(eG9cE@w{CFzV;vC$}n%&q~j%D zKHpi@*MhI(g>ZSlk39&TIUXWiIch_7?5&QAb*tp{#*ZunYW;qWSfq-#6# zkggm-&@GN&_!dV9`OG1Y)O%kRM;Ls+{*yku`Q7*tjyN3;g)77J+V$XTeJEdj=oViC z_!eJ7@|mv@e2cF!d@YyesR`+tr>3MUM>FUaM|1cVM+@?qqa}QcqZNF=JY_jwwMLxs zDpwoQm8&gui|b+d7S|)>GgmwK7T2Tjv*T)yI2{j#E6aJf1LCZCxFdXv_c8buZzuAZ zw=;Z;w+sC2cppa`=kW>BHIH4PTU<}Vx462I&sdMs% z@fKIws4LgY5y$gYAL!P66#-w5=lQBHd~J6!ex%Q5e%I?G>2sug`r@4H_5I*8pFaO( z@%6`XR^A4X&v_dN-;XbS9U=2^Gcz_)mZ!q@T& zZ4eOFhIT#?zUEKY z0nd=Gd7VVM=5;c3i*pKmE3eOz&v|_gzSa}3JD-Qo98*bGju%K*j%m;>ju+uu9Mj2X zj+fwD952Jqu6@oxoboEyOwyI>73db%Ech1JZ1S0F4t$GiE_^MowzGMpE61y(E6053 z7RPJwwckpA%JRPH*D227@w*?|ci;3Ih%142r``F*`F@aj6Mjj&3txsfoe~$o*Lo7( z4EIgH1z+n)xfepWxEH~<^7l6RoWI5JwOri4-huDeF0)+6EJ2*#_o4JS*VbAXO1nIjRt#jycC zb8IACIX0259GjtA9ACrNa_~Ch8~D~b;#)eN*Ad^rFNk;PAKvRz-{bGfryN^ISB|aF zEskyQEspK*Esh;@Jaha2pE-6yXO3NlzMFLA+XLO=+Y8^~`w_myw~vlzzMtST-_OvQ z?-$aQV?T6@<5&1rzd1lY_nY6~Yd!FM^E-UM`6kPC`a#4o*CEoC>o9bS>j->{D~WvO z`UAejbrimqSLc~yq-#GqPP%fOfNpV|gl}=2BA+=Nmq`uq$|f?&@GO? z;aeQ%$Y+l8@GXvH_E|brESB%fS zn$724gRk|%^Tl=e);M&7j^}acCVU=;Zb8?2()MwibgjP>=vMuu!nf*Ae@IrJ>n~TX zv=eAO@HmtkzLr<}YaY^-BQNR7kq^4XksrR6Q}U7Fcg_OjGe<%A7DplY%yCCH#}zh? zD?+;R7bRW!?}Tpg7lUu{7blpgJ9!*9q~2I8r6s4|e>@rXLXM_*G2b z6Xf_+E#DLDJ^y{S?>+v0+xM2gn&k(Fdz{aI!19B_+>WoW8Wr?UIqx|&kW&LWHIP#S zIW>?|135L2Qv*3QkW&LWHIP#SIW>?|135L2Qv*3QkW&LWHIP#SIW>?|135L2Qv*3Q zkW&LWHSnKZ1D)Ggs z!`D~XbM}Rm=Sv>@sp6>b6G~lLcl6r2h>3^xo}F~;Ordkn^@_Q6JKMGB*`JSHy7BWf8y{G6Z^9R!E~=Nb z?)1-7RyJCbYfY7ex2j$qm4DQ<&-Vv^eYV5WlygZHXI)A7`bgn}on>;hu%$6}dmfkE=Y@sW)#kB_} zRiCor?a+HtqLa>?-uc3%al4+Kc<<6YVqpQDb(%kzkIkva>GK0^BkUjxN~yD3TINTM2w%f^k#im zYSZK@m+t&#RNllTPfmQ`$E%**N6*e&JgR(RiCGIb&Ph78XKbS%PyLbaTJO9YpFO-I zIdSr_@%Nlv*>c^+{hy}fPl$Qv>b+MkhFq&~?eKHw8z$dODz@il@s!R<`%aI~d!$VA zmZVju$LCI6k+gfy&GH-SZt6XGW3%&vu5HXUar)s=$>WPPD$r==;e6L7x_Y-bfA_e? zlPBg*ed%!BUFZ? z?QWa8HUB$m`tW))S55tMcdlL2deylWJ2j!v#Xb#V-aX~1bnE)&e0v^z;daT~rCO~l zIIiX9?>2q6_rY8f+YYOF>EWq2F0LwaD{AwJJr5R}m~U9iOSLw4Ez=`kk9!iIPVTzp zpZx_dpNiYrv2lm}cU(@`KV$gAGta)=wNzq4(v4qx=P5t2`LO$jMP7Pi(~*(~izJ>* z{=&67@9JhNgKrJlG-U6f0=>r$D=@6}rFxr>lsT9$@#Ez9JfjL+ZL+fZt;9`2|0Qwt`yk&YfJ9lgQ5iNlj;7U<9< zrs}vko91lpyywk3QiGC5TbE#3Q^A9B~+H^Jg z-BTAmA5N_J=k;5;Q|fN+z4&VIm5W8Mb-%P>%K5R$;mKnz)!CGq=UUZ-`AIkLx>e%Z zfJ+l5r3PNRBcW^a#0eXpPq_W@;!)CQk~ZerFzxV19UIR)bm!%S!%feRx|H%vYKP>W zmr^F5k4bKvH1_nNsj0n_^ITf+L*qqFm!9&x^m(b4zfGV1)1rb4!hY;^_qExR=S-Tj zWPa|Hxk+74FPkx9Vsui6!tJN_O1Ku-;>P*`-zF??bAILz$pdl?Xj!-W(wo~Vp1HZK z{FA%OoH-v8*D_(_$yRlXC5-uFM}b|H&ZOM2scZ2b1$&fEteRXZU#aFRM=#yKvc=M? z1KM7CcyZUVDS=6!ot`&kWrH=j*HlfIm-IpTGq*eCJNx*3iM@VKzWnpv4<_~6+2(5J z!ad3*)|<6(_qx;jrwRqdnKn@_)5J-KFb^z`%9l3Oo7GBN+HPrqw?q|mruj-F~< z=+?CFDs{*`@#%lIyt!%7>sK0Wxlt;q?xW*-zg{tQe^Sh*BX3;U;JG}pSfc^I)sI=1 zT6lR%rA?(C9Xr3>_{LA&x;S@Jrl*d1YtzQpZq%Lf|Fw4|P)%Ly8WIqZAu2K` zldS_x6~U9i0kjTD ziiieC4n;`-AxD#tK!)>9gxq`UecpT5x?St7_9Y8)_St)%{qO(V|JnQ7|K8_>T|RvA z z1*gu`EUneH{TNiFMh)62{Hz8N+H?0s%or9v$ zqqhN{)(CP*RN8daPiNwtb@MS7^9XlV2vX;0{R`{$?a|rc#oe&M-m9QX(?~$I7&OwU z!Bc`41s@>^3@_ERA)p90wg4J+dzWVPCH-1AOuyw_4;zKo_0CYp;ZWiYy=ml1>fGsReCrh>j^F&#otQb`mI+PDj zIyrQ1a>Ms^+%Nldm(Q>i_~Hr7I2DXOa!}FfL}~* z?8)Kb^5Nl=!_|w`cqd7!xPwgLIOkysFcrhjm7}G3rAX4EF2QA%F&XS(0WK^fGHS-; z3MxklR@{m7dsj3+)WNP$Jqu`$p|K04R$P)e_x@uT88llZO)QRvrzpe+4qnfxRN!qR z(!_g}%iz&di)e;Sc65i(Oo>A6OvX|j~dm`FYa(oAnK zHdRU!@5PaPYxxg_yPO#g z>CvjsZ7GZ^hmpH!A!QH+;KhZZHWcnYM-{#AU{Bvs0z$vH;k$TU-2>Ld*mP_G0!@y-UO(Nsr-NP`z$Ws9dTBUNc9n~*%R;g+BlxR_co80hYm=KmySJwzS zK8DPMQ@7cUP-;ly8JFGU+TWNC6oCwn8@9j&C>%wfruU$9fDEr2-hekWXrn+W)+p?? z)s8A-XQ3L%moQ`b@X~BwitZ8<4HnA13lxFciQ0J-0S{iF^`NwLhgwJ2)13iXH^W;0 zJ_U;a#BqfE!FW-wld+vouF=e;lJlh#xP_7?h-*hjNCT*wOr1oX0hEGxw9x~&oL?=7 z1>tG>=@jNI*j5{@;uI)$K{${s(m5Lye3>v5BejJvV6i;00FHuUc^{^%+rn*day%Vv{itdaIV z`tt6NX&t7&=Kds2cB9?cN!sog&wHz>Lzt%_=YD(?PHY2S(Wgw`aN-;#x%}Bz`BU7q z*p1K6ZAz+|TJ`?xxX!?JZX@$|didG-+4y~>pEm30`20!^-LY!1X<3C`rqAGIh>Goq z)N~ZIq*byvPFj4o^l6fQ>HzzFygo@Uehg30Ke?iLw69vTwOB~7DZx(){@|(m>2SYQ zhJdG$0t|5}6cHyOi$k~ehwf-(7JW~Dx({DkK3ZB(S}|Nw`DibS`ty~n4_?ordoY)& znbypp!Bb1rOgm!WVB%7BB8oVHxM(mjNKHf&af4lf>Mk28O1f;YYmr)ND_t_!h3gX7 z9rMk&sXa|n93v?%)xq}KI+DUaPPzg@s7Fk$^yDZP%^-_){KvS7VIKNglENwi7R*ZGbdE<5mzSlSK3$GR;P(8llzg$S*o~_ z+@D}Y7#ni>iO7&6PWYlf!I@x9Ac+%F`iV}2WO2f-{scS1_<+-~eKbsbRK1WTsJ*~4 zrh$Vjz?;Q6#R3vHNz=p~$^B9t-spaz8vocKlg8Kq7j#;63Vr)+uOB1QzUCu6hiKOX z;$FrM&dh5(O7xp!0XL|@iITP`HZ(m~wOkv?*38U0$aQI0LP-P_lZe;;v~F5+*Q}5X z4-M~?>LAy>0Y{Fc(Mqj&blPFS}`9GxVNCiTbKu1ga~C-)<4IaM4@Ze&c4F-+%I#N<}2rEQK;q;9@c z{ZK{!uIn)qUDbM}?Eb#$PB4rG++cC2pdLiAfC83r6?l^cypIr+zw`-Q?!@HYS(`g2 z#xRF}HYWG%T3QMZnY_j1zF2E;<+-lSJrN_T-`&qfUEC?I|DvDmd@)H}PZHOsP=F`P zXh8d-!L}fEJGvb=7!)`dwDi=uUsD|t`%oQsoawN??i*k)e9}=dM9ihq>{cvREpan|#qtxY*jS@~a2At_?wyG~q4>#v9Y*NvCw$1COBn)>*b$ zk?c@?I4@i@aA<^*OQ8+FQjNYWCg8QBxfHx-o*EMdFlk@rM4%W)C4S|yI4D(4bCx2r z4`lL(o6KP^lHxYOZc3URtO0$6o(FR#6f!M1}zohi}4u zTCB?0tbx%!y?ut926qZK=z)h~h4No8( zQrIdY$#|+A&yHUt40VL3KpRvg{46gjiaHgDgP0STD-{c6+8=(&zMK@^L}VBrSTok82dFZWtHbDB>>G5N@r80giN%+&ttck99U+O! zAr#jTAG6t-ZS+&zV9Fe!{}ICcbdvTbUcB2D=Pf?KBq-h{eieH$bt7*hzfD;0%y3Fy ztJ=wr_uOK)#T`hBP&;Y*Mtt+t$0C4;sFBBrzqI)ra4F(cs#kUx#f3^S#HjiVDN`vX zZUiITNX5E^Lx&;W6yO;ouokQ)`%sm9#G4{Yuz|NhQ&r!M?!%Yyxn4I=zzOKW>pt}l zGX{21DZu$d4Dal4if5o;4K!D^ESa+qBHNN9$)VKHKkiksA?ZmG%J7V+S+^qfY3 zE|=BBDwsv-VMU_Q0twR8`K7nAP=@0W6UAJkaMY)oS;}f?jw;|5!?&6@#wxHl zVbc-GLM^#aQIq8kJG-gO)&hRRo^qXXLhx-^3VfZNKHIl@>@( zYkZV8#+kFLW_b;+&UozZ`Vt*lGQ33_{=F}-qokTEFEQTmz9rJ+Zad@dV0;P&%kQkl z`t%;btji;;(oPMb7^fzOAE}jL=j-qB?unxFb3Q4EJ6$$|vlP4n?vtM_!lQ#{Q7XK; z1mg}6PyI7~i15BWF(7?Ux{vDhO(v$w7eGO-HtugcaC*c8C-fhZ)4xks88=S_K9Qd| z%yn#7KzRZBV=h`<{F3lyvR(6cw3_=}tw19#86tW|F%|nW!qme0QUY2vqrj8j?qQe< zjd?t%3%vmDmD?2mC`H+Y3@XwC6%#_KwA?7oW%QHK<)ic*KUHWbg*#M0h}SmeHzjj? zN`_KnAGD4+aW%G@NH(5=_UeL7?kI#+!XBJhFB%t&i;`lsj(|Xd! z8r)RQhLUPnw8^(eK9IJ6L$7ktdv;DPjUuCZD<)PSlx+8cr$J+eiAM|**bF_KFlLx~ z#4v>k6Q%D%Zw&J$*YSoI;95e!`g_wj>{a9-$}TWrKt4h&QX|e`h-^b%MPAbT3?-R_ z%r$>3%*x9;bIrln1%}H9_wU1NLK=N_qR>3aN-ZWIQ=4VOsRfYPXy+qlX)j#}r9s}Rp_>d0)fccsDud}} z*wBwkT(s588<~Rdjh!L`)f$e5{8fbHY8-e?4&0(Dx-ok6pVLdyue#saNc=QjrcCO} z*dO6e$GnWfo~&1XXxFj=%c@u2$y6Jsdjg1I9>V(sM^Nm+=!xrW)*VG0bH@vgpDy!K zo=sQhFn?}ivQ36*L=01YHWCZe>ng03DIa28^cDQg!uwwlmT4~+?yV_vF>V1nlAlOh`#B>bi8Wx z7T!v4P4V4kEZq{$@U7L@QGDfQIjuD<8f+dneh*(aav zYY)$ZcBq;`ZA+up;5*PGVj>F`&}`M*rbb2_FHTh1+pk{S)uzOufeQJ<0`jun%0l(^ zixP70X}m9BD_?j!Y;3^_{zO}rUsikVW}8g*alE3mK2WgpJe&XO{C>58vtNFpCCroO z+1+ulN9rg+NxqUS(_Smgs=3zMX6Ss;ul@lUALWCG7{1wfiX9j_F~TO|sAP&3nQI^m z6C4-tI2b2Cbr^QkE}+P^bP1FYaclvhERkkN_bnTh^+QQ@-cNgO(@TvSxJ~klOCz#`|*TZ@Rem9}T@aZS3 zSDO|g;kQ2Kg~h+t(!kC_w@6z68@F_o>q+IGnez0 z7|ZD&9D%p(cjmah{Bip+t8cbd(w^<4WaLfB$T59i=HGr8#ccQH{Qv4Ww!zBAYTkIp z_gqo`^b+%nM$!?P|LkKP7dF3j|Gj;U+1@XX-HvZ3|8CB$t?b_P17lyN_83FT{34-+ zrvA=U&$nnKoI6&x ze~`u{|2oQy*IyO*b+YBmf`A193j!7dEC^T-_yZs?#zWD6B`it3?;r}}x|(@tAN0@j znfc}Nj6VQK|Aq?xJ`bJzhF4BT2-X;focX5o-++IPyc%l~6NNL=)BlBkZv8j3`OmYw zi-$It--!5@i20+xwkI2`OzU-+q@e?NGO2 zW>?pu$Sl(b_w>}$t_xfJ%eOkrJ&;upbii(2`Cl|i$_r&z^%;g!Bf}xOF!pTMRUdsD zne*wYwU1Y&J&tryH7|Y8?6u(bVy~l-86G9gA5`VJWGs8cw{4i{-m`eik)Q@2+%E?@ zzt%yMeeQe_mM?oXSJKilkLRa#B>v*J`=kFTe1I!ep1DY<<;(^#Lra1U>J=h0AT zD&bC(b{f7XDc(j?AE*)XN; zIjyByDRO)$LM;@zd5h5DFV05>qGl~}yZpUPbcU7l?M0b(k9eCK)@|-pXIqtalWr7V z(Edr%htt|}e<9~EQZ_ih*kIlEjb;-r-sXp@lM0@&SP25ZnHJpLOz=- zsG0ZiUGF37madA+@ktcAc0HG9%5w43YM+n9+}s{hmV0U9{fuwJD$84Pe)=NQVb8`V zN6PtoGmQ;JDC&Euu$#0e>bg4&5=Yi#W)Y^Qe3SY~kbyCGtuQPnKy-Xd)P$sMnSS5R zk>(~%JD$72dgFTUoYxm_$B6IZezra(^u877;}RH0_lw_#e~INf_;pXparT5}5{@fs zRty=L^N#({QCfX%Wy|wx|Lx0~t=Rk7>wq4djJ|8*<3ERblXDVB9`0Nq%~9PrQFYo| ze80(ifnsk})vU!I3TCxT3W?lj+>~-ku-`!4GLNfhzdJJ2GWWUybN`g|sNbuCs%twk zpFB-STe>3disU%%K$|V?lkaE~4<~=}dTq;sOidU@c_}z2z<>RRlGOR%9+N7|Nz`OhfbKg;; zkKbAcoiv`$;xBjyV3ueL0u}@;2v`uXAn*r3;9Yzn;2pZJ(yr{g2N#pzqo2dXhS1%N$dp!mH;DTBB!E} zQL2hYtrV$J6;(Z?Kap}o>Lu!lG^vyxD&^1vm2yCB#i@rxttu*2N@sT7TjIs*B?(S2 zek0-8nR)YO=9?dew=?GB$dO~~*KJt`()s=1gF5q(J|3d;f2d+8OzBa4(9oYk|GGy# zXs`}sT0O?sp?=+~A~Gn~{qE7DAS6jY`Irr&sDSiPJp1KAAae9zmuzVdHdjag4vrrT zLedzEr_ROFiR6eGH5MBl9g545rCgW&{=Kj=})GHW9h(nYRHUs z)wfBAQN85UtN|0Z&{r#->jCwhi>Jm+v~{8lsu}HO0h*>-OWk4jiT`3f*PHbmh>yfm zv2?sIV4RnV+}d3wWT>7#@w}*SgQ>4K6*oIS&?l5GL*JH1^i|~Ot2Em^HW)KIWV!8? zI;5*qp3u=`K51fQtr5gVM&38WOZ1_9E&7@uFPYJLN8@k?3z(B2gmYUr2e&CX(xXEB zSuO!!8-;Tl97#A1<3WI(KAV`kKFJRK8LC`Rc6spR-*U@?~FvrYmIe#hYx$4yze zAPPtMuD*Vf?~B*POewd4W`qhZM%0HiR5rV ztx=81)BuRtrW{Tah-cI<0-B$4-b>fq6buELS_94dx`TV$f_vIREsZU~a7$~eN~d{D zxH>w{p6Gt>Y>4#Q;`dOPe3WmY_SoW^)y9j~QUH_>5+5q(q*?utSVNS{jo$<7^TzLm z4SC~R#HyZ8Pc#0u1IX1ML}`F9DYRnvfeA1HCcp%k025#WOn?b60Vco%m;k`Z&J#y8 z@Sv>_;0pc$PIdshJ&WvAnRgUDlDmrGnUm#7)yW=sceCPgf2;D{?1A@pDLLX?o;1t1 z)pmQJjwtHTPE2IT_LmsxOZ3LaK=`BS)S~^)vLjLk7hkOdK4?I>8=dNb=Jy9fxZ8t)5l=a#Dyhsj z=~8r~e!zNP6RL*LLo}pTh*+8K@hY{|pqmf)I!8RqbPf9A{X_DN3fWctsMqW9%0i#j35|Cn?dD`;_aQtRKBECjSkGBEmEPF`YPB6?e$T(M*5|FaajO z1eieS5>We4ueDK`-d?yZv=0q86_%s@WoBPtxwVp-Inx`3$=ZLSmk!?64&7R)>DpiQ zMfj6V*-N*6xjC2~zY>DU4{ExFn308JdVi&{N$;;Pb~U_m`=g&O-fo|$-wM8O;dQA; zIG?dspZ6JO^!X|Sv@ao$IrF=Sj?6{D+o)FtV-uaB{ZohLzx(RJuaosK5cZrlYoKkc zYWi^ot{s^DtmD?-)$NH54WBkXytxBrRpEU?Nu_bD;pZPfJ06+1yltxX>)6ocZQ$Dp zI+pD*PRX)I^m)y=pjUz?3auh(0E^f3{*K@N`sTsEcNgvrPV5Nn-(Ba0&;RiJxmLQ8 z+rb2w025#WOkmjr9^VdCC#*E$T(mH>C>xS)mbC)*rwE=-lcpRpPvx;c^_b-q>%++o z^{PXGXC3TNRs$UJaWoTP0!)AjFaah|x&++zPjv5qe6Ki}U!?uhRGY#mi?X&KOW)r; zpZycr)$Bc;we6pXlhpuw1`aa;Ccp%k025#WB~QR@|3vpZ$)|HrPWQ`Kzz%iiu%ar9 zq}PUd>`>449Mo$LJ$=QteN;He4N@1wT&6jn+5&^tz9IoeNfC-c?0k{1V{XYoui+CTUepq3YMOoXA)v|xWui0GE z(f`@@PsGV;z>*vs!33B96JP>NfC-d30k{1V-A^T--h&qLK1$Cy1yvSFX2a+H9(1Qe zN5A}LZUS+#DxRCnAtt~Cm;e)C0!&~H5O6wAYCjL>cAlTlldgfm@Dxq3pC_erboaBl mO6go?S8ntZhn&W*{TwNe-zP>~78#w!?=+r6CI~A%e*XvBCayjJ literal 0 HcmV?d00001 diff --git a/q2_diversity/tests/data/vaw.nwk b/q2_diversity/tests/data/vaw.nwk new file mode 100644 index 00000000..1ba14a5b --- /dev/null +++ b/q2_diversity/tests/data/vaw.nwk @@ -0,0 +1 @@ +(GG_OTU_1:1,(GG_OTU_2:1,GG_OTU_3:1):1,(GG_OTU_5:1,GG_OTU_4:1):1); diff --git a/q2_diversity/tests/test_beta.py b/q2_diversity/tests/test_beta.py index 8e9f38a5..dd055deb 100644 --- a/q2_diversity/tests/test_beta.py +++ b/q2_diversity/tests/test_beta.py @@ -21,9 +21,11 @@ import pandas as pd import scipy.misc import qiime2 +from qiime2.plugin.testing import TestPluginBase -from q2_diversity import (beta, beta_phylogenetic, bioenv, - beta_group_significance, beta_correlation, + +from q2_diversity import (beta, beta_phylogenetic, beta_phylogenetic_alt, + bioenv, beta_group_significance, beta_correlation, beta_rarefaction) from q2_diversity._beta._visualizer import (_get_distance_boxplot_data, _metadata_distance, @@ -161,6 +163,222 @@ def test_beta_phylogenetic_weighted_unifrac_threads_error(self): metric='weighted_unifrac', n_jobs=-1) +class BetaDiversityAltTests(TestPluginBase): + # Note that some of these tests replicate the cases in biocore/unifrac + package = 'q2_diversity.tests' + + def test_beta_unweighted(self): + bt_fp = self.get_data_path('crawford.biom') + tree_fp = self.get_data_path('crawford.nwk') + + actual = beta_phylogenetic_alt(table=bt_fp, + phylogeny=tree_fp, + metric='unweighted_unifrac') + + # computed with beta-phylogenetic + data = np.array([0.71836067, 0.71317361, 0.69746044, 0.62587207, + 0.72826674, 0.72065895, 0.72640581, 0.73606053, + 0.70302967, 0.73407301, 0.6548042, 0.71547381, + 0.78397813, 0.72318399, 0.76138933, 0.61041275, + 0.62331299, 0.71848305, 0.70416337, 0.75258475, + 0.79249029, 0.64392779, 0.70052733, 0.69832716, + 0.77818938, 0.72959894, 0.75782689, 0.71005144, + 0.75065046, 0.78944369, 0.63593642, 0.71283615, + 0.58314638, 0.69200762, 0.68972056, 0.71514083]) + ids = ('10084.PC.481', '10084.PC.593', '10084.PC.356', '10084.PC.355', + '10084.PC.354', '10084.PC.636', '10084.PC.635', '10084.PC.607', + '10084.PC.634') + expected = skbio.DistanceMatrix(data, ids=ids) + + self.assertEqual(actual.ids, expected.ids) + for id1 in actual.ids: + for id2 in actual.ids: + npt.assert_almost_equal(actual[id1, id2], expected[id1, id2]) + + def test_beta_unweighted_parallel(self): + bt_fp = self.get_data_path('crawford.biom') + tree_fp = self.get_data_path('crawford.nwk') + + actual = beta_phylogenetic_alt(table=bt_fp, + phylogeny=tree_fp, + metric='unweighted_unifrac', + n_jobs=2) + + # computed with beta-phylogenetic + data = np.array([0.71836067, 0.71317361, 0.69746044, 0.62587207, + 0.72826674, 0.72065895, 0.72640581, 0.73606053, + 0.70302967, 0.73407301, 0.6548042, 0.71547381, + 0.78397813, 0.72318399, 0.76138933, 0.61041275, + 0.62331299, 0.71848305, 0.70416337, 0.75258475, + 0.79249029, 0.64392779, 0.70052733, 0.69832716, + 0.77818938, 0.72959894, 0.75782689, 0.71005144, + 0.75065046, 0.78944369, 0.63593642, 0.71283615, + 0.58314638, 0.69200762, 0.68972056, 0.71514083]) + ids = ('10084.PC.481', '10084.PC.593', '10084.PC.356', '10084.PC.355', + '10084.PC.354', '10084.PC.636', '10084.PC.635', '10084.PC.607', + '10084.PC.634') + expected = skbio.DistanceMatrix(data, ids=ids) + + self.assertEqual(actual.ids, expected.ids) + for id1 in actual.ids: + for id2 in actual.ids: + npt.assert_almost_equal(actual[id1, id2], expected[id1, id2]) + + def test_beta_weighted(self): + bt_fp = self.get_data_path('crawford.biom') + tree_fp = self.get_data_path('crawford.nwk') + + actual = beta_phylogenetic_alt(table=bt_fp, + phylogeny=tree_fp, + metric='weighted_unifrac') + + # computed with beta-phylogenetic (weighted_unifrac) + data = np.array([0.44656238, 0.23771096, 0.30489123, 0.23446002, + 0.65723575, 0.44911772, 0.381904, 0.69144829, + 0.39611776, 0.36568012, 0.53377975, 0.48908025, + 0.35155196, 0.28318669, 0.57376916, 0.23395746, + 0.24658122, 0.60271637, 0.39802552, 0.36567394, + 0.68062701, 0.36862049, 0.48350632, 0.33024631, + 0.33266697, 0.53464744, 0.74605075, 0.53951035, + 0.49680733, 0.79178838, 0.37109012, 0.52629343, + 0.22118218, 0.32400805, 0.43189708, 0.59705893]) + ids = ('10084.PC.481', '10084.PC.593', '10084.PC.356', '10084.PC.355', + '10084.PC.354', '10084.PC.636', '10084.PC.635', '10084.PC.607', + '10084.PC.634') + expected = skbio.DistanceMatrix(data, ids=ids) + + self.assertEqual(actual.ids, expected.ids) + for id1 in actual.ids: + for id2 in actual.ids: + npt.assert_almost_equal(actual[id1, id2], expected[id1, id2]) + + def test_variance_adjusted_normalized(self): + bt_fp = self.get_data_path('vaw.biom') + tree_fp = self.get_data_path('vaw.nwk') + + actual = beta_phylogenetic_alt(table=bt_fp, + phylogeny=tree_fp, + metric='weighted_normalized_unifrac', + variance_adjusted=True) + + data = np.array([[0.0000000, 0.4086040, 0.6240185, 0.4639481, + 0.2857143, 0.2766318], + [0.4086040, 0.0000000, 0.3798594, 0.6884992, + 0.6807616, 0.4735781], + [0.6240185, 0.3798594, 0.0000000, 0.7713254, + 0.8812897, 0.5047114], + [0.4639481, 0.6884992, 0.7713254, 0.0000000, + 0.6666667, 0.2709298], + [0.2857143, 0.6807616, 0.8812897, 0.6666667, + 0.0000000, 0.4735991], + [0.2766318, 0.4735781, 0.5047114, 0.2709298, + 0.4735991, 0.0000000]]) + ids = ('Sample1', 'Sample2', 'Sample3', 'Sample4', 'Sample5', + 'Sample6') + expected = skbio.DistanceMatrix(data, ids=ids) + + self.assertEqual(actual.ids, expected.ids) + for id1 in actual.ids: + for id2 in actual.ids: + npt.assert_almost_equal(actual[id1, id2], expected[id1, id2]) + + def test_generalized_unifrac(self): + bt_fp = self.get_data_path('vaw.biom') + tree_fp = self.get_data_path('vaw.nwk') + + actual = beta_phylogenetic_alt(table=bt_fp, + phylogeny=tree_fp, + metric='generalized_unifrac', + alpha=0.5) + + data = np.array([[0.0000000, 0.4040518, 0.6285560, 0.5869439, + 0.4082483, 0.2995673], + [0.4040518, 0.0000000, 0.4160597, 0.7071068, + 0.7302479, 0.4860856], + [0.6285560, 0.4160597, 0.0000000, 0.8005220, + 0.9073159, 0.5218198], + [0.5869439, 0.7071068, 0.8005220, 0.0000000, + 0.4117216, 0.3485667], + [0.4082483, 0.7302479, 0.9073159, 0.4117216, + 0.0000000, 0.6188282], + [0.2995673, 0.4860856, 0.5218198, 0.3485667, + 0.6188282, 0.0000000]]) + ids = ('Sample1', 'Sample2', 'Sample3', 'Sample4', 'Sample5', + 'Sample6') + expected = skbio.DistanceMatrix(data, ids=ids) + + self.assertEqual(actual.ids, expected.ids) + for id1 in actual.ids: + for id2 in actual.ids: + npt.assert_almost_equal(actual[id1, id2], expected[id1, id2]) + + def test_generalized_unifrac_no_alpha(self): + bt_fp = self.get_data_path('crawford.biom') + tree_fp = self.get_data_path('crawford.nwk') + + actual = beta_phylogenetic_alt(table=bt_fp, + phylogeny=tree_fp, + metric='generalized_unifrac', + alpha=None) + + # alpha=1 should be equal to weighted normalized UniFrac + data = np.array([0.2821874, 0.16148405, 0.20186143, 0.1634832, + 0.40351108, 0.29135056, 0.24790944, 0.41967404, + 0.24642185, 0.22218489, 0.34007547, 0.27722011, + 0.20963881, 0.16897221, 0.3217958, 0.15237816, + 0.16899207, 0.36445044, 0.25408941, 0.23358681, + 0.4069374, 0.24615927, 0.28573888, 0.20578184, + 0.20742006, 0.31249151, 0.46169893, 0.35294595, + 0.32522355, 0.48437103, 0.21534558, 0.30558908, + 0.12091004, 0.19817777, 0.24792853, 0.34293674]) + ids = ('10084.PC.481', '10084.PC.593', '10084.PC.356', '10084.PC.355', + '10084.PC.354', '10084.PC.636', '10084.PC.635', '10084.PC.607', + '10084.PC.634') + expected = skbio.DistanceMatrix(data, ids=ids) + + self.assertEqual(actual.ids, expected.ids) + for id1 in actual.ids: + for id2 in actual.ids: + npt.assert_almost_equal(actual[id1, id2], expected[id1, id2]) + + def test_beta_phylogenetic_alpha_on_non_generalized(self): + bt_fp = self.get_data_path('crawford.biom') + tree_fp = self.get_data_path('tree.nwk') + + with self.assertRaisesRegex(ValueError, 'The alpha parameter is only ' + 'allowed when the choice of metric is ' + 'generalized_unifrac'): + beta_phylogenetic_alt(table=bt_fp, phylogeny=tree_fp, + metric='unweighted_unifrac', + alpha=0.11) + + def test_beta_phylogenetic_non_phylo_metric(self): + bt_fp = self.get_data_path('crawford.biom') + tree_fp = self.get_data_path('tree.nwk') + + with self.assertRaises(ValueError): + beta_phylogenetic_alt(table=bt_fp, phylogeny=tree_fp, + metric='braycurtis') + + def test_beta_phylogenetic_unknown_metric(self): + bt_fp = self.get_data_path('crawford.biom') + tree_fp = self.get_data_path('tree.nwk') + + with self.assertRaises(ValueError): + beta_phylogenetic_alt(table=bt_fp, phylogeny=tree_fp, + metric='not-a-metric') + + def test_beta_phylogenetic_too_many_jobs(self): + bt_fp = self.get_data_path('crawford.biom') + tree_fp = self.get_data_path('tree.nwk') + + with self.assertRaises(ValueError): + # cannot guarantee that this will always be true, but it would be + # odd to see a machine with these many CPUs + beta_phylogenetic_alt(table=bt_fp, phylogeny=tree_fp, + metric='unweighted_unifrac', n_jobs=11117) + + class BioenvTests(unittest.TestCase): def test_bioenv(self):