diff --git a/previews/PR12/.documenter-siteinfo.json b/previews/PR12/.documenter-siteinfo.json
index ca0fe71..8de2fd2 100644
--- a/previews/PR12/.documenter-siteinfo.json
+++ b/previews/PR12/.documenter-siteinfo.json
@@ -1 +1 @@
-{"documenter":{"julia_version":"1.10.3","generation_timestamp":"2024-05-17T03:10:50","documenter_version":"1.4.1"}}
\ No newline at end of file
+{"documenter":{"julia_version":"1.10.3","generation_timestamp":"2024-05-17T12:07:52","documenter_version":"1.4.1"}}
\ No newline at end of file
diff --git a/previews/PR12/index.html b/previews/PR12/index.html
index 5194537..6a6101c 100644
--- a/previews/PR12/index.html
+++ b/previews/PR12/index.html
@@ -1,2 +1,2 @@
<!DOCTYPE html>
-<html lang="en"><head><meta charset="UTF-8"/><meta name="viewport" content="width=device-width, initial-scale=1.0"/><title>home · PhyloCoalSimulations.jl</title><meta name="title" content="home · PhyloCoalSimulations.jl"/><meta property="og:title" content="home · PhyloCoalSimulations.jl"/><meta property="twitter:title" content="home · PhyloCoalSimulations.jl"/><meta name="description" content="Documentation for PhyloCoalSimulations.jl."/><meta property="og:description" content="Documentation for PhyloCoalSimulations.jl."/><meta property="twitter:description" content="Documentation for PhyloCoalSimulations.jl."/><meta property="og:url" content="https://cecileane.github.io/PhyloCoalSimulations.jl/stable/"/><meta property="twitter:url" content="https://cecileane.github.io/PhyloCoalSimulations.jl/stable/"/><link rel="canonical" href="https://cecileane.github.io/PhyloCoalSimulations.jl/stable/"/><script data-outdated-warner src="assets/warner.js"></script><link href="https://cdnjs.cloudflare.com/ajax/libs/lato-font/3.0.0/css/lato-font.min.css" rel="stylesheet" type="text/css"/><link href="https://cdnjs.cloudflare.com/ajax/libs/juliamono/0.050/juliamono.min.css" rel="stylesheet" type="text/css"/><link href="https://cdnjs.cloudflare.com/ajax/libs/font-awesome/6.4.2/css/fontawesome.min.css" rel="stylesheet" type="text/css"/><link href="https://cdnjs.cloudflare.com/ajax/libs/font-awesome/6.4.2/css/solid.min.css" rel="stylesheet" type="text/css"/><link href="https://cdnjs.cloudflare.com/ajax/libs/font-awesome/6.4.2/css/brands.min.css" rel="stylesheet" type="text/css"/><link href="https://cdnjs.cloudflare.com/ajax/libs/KaTeX/0.16.8/katex.min.css" rel="stylesheet" type="text/css"/><script>documenterBaseURL="."</script><script src="https://cdnjs.cloudflare.com/ajax/libs/require.js/2.3.6/require.min.js" data-main="assets/documenter.js"></script><script src="search_index.js"></script><script src="siteinfo.js"></script><script src="../versions.js"></script><link class="docs-theme-link" rel="stylesheet" type="text/css" href="assets/themes/documenter-dark.css" data-theme-name="documenter-dark" data-theme-primary-dark/><link class="docs-theme-link" rel="stylesheet" type="text/css" href="assets/themes/documenter-light.css" data-theme-name="documenter-light" data-theme-primary/><script src="assets/themeswap.js"></script></head><body><div id="documenter"><nav class="docs-sidebar"><div class="docs-package-name"><span class="docs-autofit"><a href>PhyloCoalSimulations.jl</a></span></div><button class="docs-search-query input is-rounded is-small is-clickable my-2 mx-auto py-1 px-2" id="documenter-search-query">Search docs (Ctrl + /)</button><ul class="docs-menu"><li class="is-active"><a class="tocitem" href>home</a></li><li><span class="tocitem">manual</span><ul><li><a class="tocitem" href="man/getting_started/">getting started</a></li><li><a class="tocitem" href="man/mapping_genetree_to_network/">mapping gene trees into the species network</a></li><li><a class="tocitem" href="man/converting_coal2generation_units/">converting between units</a></li><li><a class="tocitem" href="man/correlated_inheritance/">correlated inheritance</a></li><li><a class="tocitem" href="man/more_examples/">more examples</a></li></ul></li><li><span class="tocitem">library</span><ul><li><a class="tocitem" href="lib/public/">public</a></li><li><a class="tocitem" href="lib/internal/">internals</a></li></ul></li></ul><div class="docs-version-selector field has-addons"><div class="control"><span class="docs-label button is-static is-size-7">Version</span></div><div class="docs-selector control is-expanded"><div class="select is-fullwidth is-size-7"><select id="documenter-version-selector"></select></div></div></div></nav><div class="docs-main"><header class="docs-navbar"><a class="docs-sidebar-button docs-navbar-link fa-solid fa-bars is-hidden-desktop" id="documenter-sidebar-button" href="#"></a><nav class="breadcrumb"><ul class="is-hidden-mobile"><li class="is-active"><a href>home</a></li></ul><ul class="is-hidden-tablet"><li class="is-active"><a href>home</a></li></ul></nav><div class="docs-right"><a class="docs-navbar-link" href="https://github.com/cecileane/PhyloCoalSimulations.jl" title="View the repository on GitHub"><span class="docs-icon fa-brands"></span><span class="docs-label is-hidden-touch">GitHub</span></a><a class="docs-navbar-link" href="https://github.com/cecileane/PhyloCoalSimulations.jl/blob/main/docs/src/index.md" title="Edit source on GitHub"><span class="docs-icon fa-solid"></span></a><a class="docs-settings-button docs-navbar-link fa-solid fa-gear" id="documenter-settings-button" href="#" title="Settings"></a><a class="docs-article-toggle-button fa-solid fa-chevron-up" id="documenter-article-toggle-button" href="javascript:;" title="Collapse all docstrings"></a></div></header><article class="content" id="documenter-page"><h1 id="PhyloCoalSimulations"><a class="docs-heading-anchor" href="#PhyloCoalSimulations">PhyloCoalSimulations</a><a id="PhyloCoalSimulations-1"></a><a class="docs-heading-anchor-permalink" href="#PhyloCoalSimulations" title="Permalink"></a></h1><p><a href="https://github.com/cecileane/PhyloCoalSimulations.jl">PhyloCoalSimulations</a> is a <a href="http://julialang.org">Julia</a> package to simulate phylogenies under the coalescent. It depends on <a href="https://github.com/crsl4/PhyloNetworks.jl">PhyloNetworks</a> for the phylogenetic data structures, and manipulation of phylogenies.</p><p>References: please see <a href="https://github.com/cecileane/PhyloCoalSimulations.jl/blob/main/CITATION.bib">bibtex entries here</a>.</p><ul><li><p>for this package and the network coalescent model with inheritance correlation:<br/>Fogg, Allman & Ané (2023). PhyloCoalSimulations: A simulator for network multispecies coalescent models, including a new extension for the inheritance of gene flow. <a href="https://academic.oup.com/sysbio/advance-article/doi/10.1093/sysbio/syad030/7187011?guestAccessKey=616778e9-ea09-4e4e-9853-229499b5f068">Systematic Biology</a>, 72(5):1171–1179. <a href="https://doi.org/10.1093/sysbio/syad030">doi:10.1093/sysbio/syad030</a>.</p></li><li><p>for <a href="https://github.com/crsl4/PhyloNetworks.jl">PhyloNetworks</a>, which this package depends on:<br/>Solís-Lemus, Bastide & Ané (2017). PhyloNetworks: a package for phylogenetic networks. <a href="https://academic.oup.com/mbe/article/doi/10.1093/molbev/msx235/4103410/PhyloNetworks-a-package-for-phylogenetic-networks?guestAccessKey=230afceb-df28-4160-832d-aa7c73f86369">Molecular Biology and Evolution</a>, 34(12):3292–3298. <a href="https://doi.org/10.1093/molbev/msx235">doi:10.1093/molbev/msx235</a></p></li></ul><p>For a tutorial, see the manual:</p><ul><li><a href="man/getting_started/#getting-started">getting started</a></li><li class="no-marker"><ul><li><a href="man/getting_started/#installation">installation</a></li><li><a href="man/getting_started/#basic-simulation-example">basic simulation example</a></li><li class="no-marker"><ul><li><a href="man/getting_started/#example-network">example network</a></li><li><a href="man/getting_started/#basic-example:-simulate,-save-to-file,-plot">basic example: simulate, save to file, plot</a></li></ul></li><li><a href="man/getting_started/#several-individuals-per-species">several individuals per species</a></li></ul></li><li><a href="man/mapping_genetree_to_network/#mapping-gene-trees-into-the-species-phylogeny">mapping gene trees into the species phylogeny</a></li><li class="no-marker"><ul><li><a href="man/mapping_genetree_to_network/#naming-internal-nodes">naming internal nodes</a></li><li><a href="man/mapping_genetree_to_network/#cleaning-gene-trees">cleaning gene trees</a></li></ul></li><li><a href="man/converting_coal2generation_units/#converting-between-units">converting between units</a></li><li class="no-marker"><ul><li><a href="man/converting_coal2generation_units/#converting-edge-lengths-in-gene-trees">converting edge lengths in gene trees</a></li></ul></li><li><a href="man/converting_coal2generation_units/#number-of-generations-in-the-network-and-gene-trees">number of generations in the network and gene trees</a></li><li><a href="man/correlated_inheritance/#correlated-inheritance">correlated inheritance</a></li><li><a href="man/more_examples/#example-uses">example uses</a></li><li class="no-marker"><ul><li><a href="man/more_examples/#counting-deep-coalescences">counting deep coalescences</a></li><li><a href="man/more_examples/#number-of-lineages-inherited-via-gene-flow">number of lineages inherited via gene flow</a></li><li><a href="man/more_examples/#rate-variation-across-species">rate variation across species</a></li></ul></li></ul><p>For help on individual functions, see the library:</p><ul><li><a href="lib/public/#public-documentation">public documentation</a></li><li class="no-marker"><ul><li><a href="lib/public/#functions-and-types">functions & types</a></li><li><a href="lib/public/#index">index</a></li></ul></li><li><a href="lib/internal/#internal-documentation">internal documentation</a></li><li class="no-marker"><ul><li><a href="lib/internal/#functions-and-types">functions & types</a></li><li><a href="lib/internal/#index">index</a></li></ul></li></ul></article><nav class="docs-footer"><a class="docs-footer-nextpage" href="man/getting_started/">getting started »</a><div class="flexbox-break"></div><p class="footer-message">Powered by <a href="https://github.com/JuliaDocs/Documenter.jl">Documenter.jl</a> and the <a href="https://julialang.org/">Julia Programming Language</a>.</p></nav></div><div class="modal" id="documenter-settings"><div class="modal-background"></div><div class="modal-card"><header class="modal-card-head"><p class="modal-card-title">Settings</p><button class="delete"></button></header><section class="modal-card-body"><p><label class="label">Theme</label><div class="select"><select id="documenter-themepicker"><option value="auto">Automatic (OS)</option><option value="documenter-light">documenter-light</option><option value="documenter-dark">documenter-dark</option></select></div></p><hr/><p>This document was generated with <a href="https://github.com/JuliaDocs/Documenter.jl">Documenter.jl</a> version 1.4.1 on <span class="colophon-date" title="Friday 17 May 2024 03:10">Friday 17 May 2024</span>. Using Julia version 1.10.3.</p></section><footer class="modal-card-foot"></footer></div></div></div></body></html>
+<html lang="en"><head><meta charset="UTF-8"/><meta name="viewport" content="width=device-width, initial-scale=1.0"/><title>home · PhyloCoalSimulations.jl</title><meta name="title" content="home · PhyloCoalSimulations.jl"/><meta property="og:title" content="home · PhyloCoalSimulations.jl"/><meta property="twitter:title" content="home · PhyloCoalSimulations.jl"/><meta name="description" content="Documentation for PhyloCoalSimulations.jl."/><meta property="og:description" content="Documentation for PhyloCoalSimulations.jl."/><meta property="twitter:description" content="Documentation for PhyloCoalSimulations.jl."/><meta property="og:url" content="https://cecileane.github.io/PhyloCoalSimulations.jl/stable/"/><meta property="twitter:url" content="https://cecileane.github.io/PhyloCoalSimulations.jl/stable/"/><link rel="canonical" href="https://cecileane.github.io/PhyloCoalSimulations.jl/stable/"/><script data-outdated-warner src="assets/warner.js"></script><link href="https://cdnjs.cloudflare.com/ajax/libs/lato-font/3.0.0/css/lato-font.min.css" rel="stylesheet" type="text/css"/><link href="https://cdnjs.cloudflare.com/ajax/libs/juliamono/0.050/juliamono.min.css" rel="stylesheet" type="text/css"/><link href="https://cdnjs.cloudflare.com/ajax/libs/font-awesome/6.4.2/css/fontawesome.min.css" rel="stylesheet" type="text/css"/><link href="https://cdnjs.cloudflare.com/ajax/libs/font-awesome/6.4.2/css/solid.min.css" rel="stylesheet" type="text/css"/><link href="https://cdnjs.cloudflare.com/ajax/libs/font-awesome/6.4.2/css/brands.min.css" rel="stylesheet" type="text/css"/><link href="https://cdnjs.cloudflare.com/ajax/libs/KaTeX/0.16.8/katex.min.css" rel="stylesheet" type="text/css"/><script>documenterBaseURL="."</script><script src="https://cdnjs.cloudflare.com/ajax/libs/require.js/2.3.6/require.min.js" data-main="assets/documenter.js"></script><script src="search_index.js"></script><script src="siteinfo.js"></script><script src="../versions.js"></script><link class="docs-theme-link" rel="stylesheet" type="text/css" href="assets/themes/documenter-dark.css" data-theme-name="documenter-dark" data-theme-primary-dark/><link class="docs-theme-link" rel="stylesheet" type="text/css" href="assets/themes/documenter-light.css" data-theme-name="documenter-light" data-theme-primary/><script src="assets/themeswap.js"></script></head><body><div id="documenter"><nav class="docs-sidebar"><div class="docs-package-name"><span class="docs-autofit"><a href>PhyloCoalSimulations.jl</a></span></div><button class="docs-search-query input is-rounded is-small is-clickable my-2 mx-auto py-1 px-2" id="documenter-search-query">Search docs (Ctrl + /)</button><ul class="docs-menu"><li class="is-active"><a class="tocitem" href>home</a></li><li><span class="tocitem">manual</span><ul><li><a class="tocitem" href="man/getting_started/">getting started</a></li><li><a class="tocitem" href="man/mapping_genetree_to_network/">mapping gene trees into the species network</a></li><li><a class="tocitem" href="man/converting_coal2generation_units/">converting between units</a></li><li><a class="tocitem" href="man/correlated_inheritance/">correlated inheritance</a></li><li><a class="tocitem" href="man/more_examples/">more examples</a></li></ul></li><li><span class="tocitem">library</span><ul><li><a class="tocitem" href="lib/public/">public</a></li><li><a class="tocitem" href="lib/internal/">internals</a></li></ul></li></ul><div class="docs-version-selector field has-addons"><div class="control"><span class="docs-label button is-static is-size-7">Version</span></div><div class="docs-selector control is-expanded"><div class="select is-fullwidth is-size-7"><select id="documenter-version-selector"></select></div></div></div></nav><div class="docs-main"><header class="docs-navbar"><a class="docs-sidebar-button docs-navbar-link fa-solid fa-bars is-hidden-desktop" id="documenter-sidebar-button" href="#"></a><nav class="breadcrumb"><ul class="is-hidden-mobile"><li class="is-active"><a href>home</a></li></ul><ul class="is-hidden-tablet"><li class="is-active"><a href>home</a></li></ul></nav><div class="docs-right"><a class="docs-navbar-link" href="https://github.com/cecileane/PhyloCoalSimulations.jl" title="View the repository on GitHub"><span class="docs-icon fa-brands"></span><span class="docs-label is-hidden-touch">GitHub</span></a><a class="docs-navbar-link" href="https://github.com/cecileane/PhyloCoalSimulations.jl/blob/main/docs/src/index.md" title="Edit source on GitHub"><span class="docs-icon fa-solid"></span></a><a class="docs-settings-button docs-navbar-link fa-solid fa-gear" id="documenter-settings-button" href="#" title="Settings"></a><a class="docs-article-toggle-button fa-solid fa-chevron-up" id="documenter-article-toggle-button" href="javascript:;" title="Collapse all docstrings"></a></div></header><article class="content" id="documenter-page"><h1 id="PhyloCoalSimulations"><a class="docs-heading-anchor" href="#PhyloCoalSimulations">PhyloCoalSimulations</a><a id="PhyloCoalSimulations-1"></a><a class="docs-heading-anchor-permalink" href="#PhyloCoalSimulations" title="Permalink"></a></h1><p><a href="https://github.com/cecileane/PhyloCoalSimulations.jl">PhyloCoalSimulations</a> is a <a href="http://julialang.org">Julia</a> package to simulate phylogenies under the coalescent. It depends on <a href="https://github.com/crsl4/PhyloNetworks.jl">PhyloNetworks</a> for the phylogenetic data structures, and manipulation of phylogenies.</p><p>References: please see <a href="https://github.com/cecileane/PhyloCoalSimulations.jl/blob/main/CITATION.bib">bibtex entries here</a>.</p><ul><li><p>for this package and the network coalescent model with inheritance correlation:<br/>Fogg, Allman & Ané (2023). PhyloCoalSimulations: A simulator for network multispecies coalescent models, including a new extension for the inheritance of gene flow. <a href="https://academic.oup.com/sysbio/advance-article/doi/10.1093/sysbio/syad030/7187011?guestAccessKey=616778e9-ea09-4e4e-9853-229499b5f068">Systematic Biology</a>, 72(5):1171–1179. <a href="https://doi.org/10.1093/sysbio/syad030">doi:10.1093/sysbio/syad030</a>.</p></li><li><p>for <a href="https://github.com/crsl4/PhyloNetworks.jl">PhyloNetworks</a>, which this package depends on:<br/>Solís-Lemus, Bastide & Ané (2017). PhyloNetworks: a package for phylogenetic networks. <a href="https://academic.oup.com/mbe/article/doi/10.1093/molbev/msx235/4103410/PhyloNetworks-a-package-for-phylogenetic-networks?guestAccessKey=230afceb-df28-4160-832d-aa7c73f86369">Molecular Biology and Evolution</a>, 34(12):3292–3298. <a href="https://doi.org/10.1093/molbev/msx235">doi:10.1093/molbev/msx235</a></p></li></ul><p>For a tutorial, see the manual:</p><ul><li><a href="man/getting_started/#getting-started">getting started</a></li><li class="no-marker"><ul><li><a href="man/getting_started/#installation">installation</a></li><li><a href="man/getting_started/#basic-simulation-example">basic simulation example</a></li><li class="no-marker"><ul><li><a href="man/getting_started/#example-network">example network</a></li><li><a href="man/getting_started/#basic-example:-simulate,-save-to-file,-plot">basic example: simulate, save to file, plot</a></li></ul></li><li><a href="man/getting_started/#several-individuals-per-species">several individuals per species</a></li></ul></li><li><a href="man/mapping_genetree_to_network/#mapping-gene-trees-into-the-species-phylogeny">mapping gene trees into the species phylogeny</a></li><li class="no-marker"><ul><li><a href="man/mapping_genetree_to_network/#naming-internal-nodes">naming internal nodes</a></li><li><a href="man/mapping_genetree_to_network/#cleaning-gene-trees">cleaning gene trees</a></li></ul></li><li><a href="man/converting_coal2generation_units/#converting-between-units">converting between units</a></li><li class="no-marker"><ul><li><a href="man/converting_coal2generation_units/#converting-edge-lengths-in-gene-trees">converting edge lengths in gene trees</a></li></ul></li><li><a href="man/converting_coal2generation_units/#number-of-generations-in-the-network-and-gene-trees">number of generations in the network and gene trees</a></li><li><a href="man/correlated_inheritance/#correlated-inheritance">correlated inheritance</a></li><li><a href="man/more_examples/#example-uses">example uses</a></li><li class="no-marker"><ul><li><a href="man/more_examples/#counting-deep-coalescences">counting deep coalescences</a></li><li><a href="man/more_examples/#number-of-lineages-inherited-via-gene-flow">number of lineages inherited via gene flow</a></li><li><a href="man/more_examples/#rate-variation-across-species">rate variation across species</a></li></ul></li></ul><p>For help on individual functions, see the library:</p><ul><li><a href="lib/public/#public-documentation">public documentation</a></li><li class="no-marker"><ul><li><a href="lib/public/#functions-and-types">functions & types</a></li><li><a href="lib/public/#index">index</a></li></ul></li><li><a href="lib/internal/#internal-documentation">internal documentation</a></li><li class="no-marker"><ul><li><a href="lib/internal/#functions-and-types">functions & types</a></li><li><a href="lib/internal/#index">index</a></li></ul></li></ul></article><nav class="docs-footer"><a class="docs-footer-nextpage" href="man/getting_started/">getting started »</a><div class="flexbox-break"></div><p class="footer-message">Powered by <a href="https://github.com/JuliaDocs/Documenter.jl">Documenter.jl</a> and the <a href="https://julialang.org/">Julia Programming Language</a>.</p></nav></div><div class="modal" id="documenter-settings"><div class="modal-background"></div><div class="modal-card"><header class="modal-card-head"><p class="modal-card-title">Settings</p><button class="delete"></button></header><section class="modal-card-body"><p><label class="label">Theme</label><div class="select"><select id="documenter-themepicker"><option value="auto">Automatic (OS)</option><option value="documenter-light">documenter-light</option><option value="documenter-dark">documenter-dark</option></select></div></p><hr/><p>This document was generated with <a href="https://github.com/JuliaDocs/Documenter.jl">Documenter.jl</a> version 1.4.1 on <span class="colophon-date" title="Friday 17 May 2024 12:07">Friday 17 May 2024</span>. Using Julia version 1.10.3.</p></section><footer class="modal-card-foot"></footer></div></div></div></body></html>
diff --git a/previews/PR12/lib/internal/index.html b/previews/PR12/lib/internal/index.html
index 6efb1fe..61a0243 100644
--- a/previews/PR12/lib/internal/index.html
+++ b/previews/PR12/lib/internal/index.html
@@ -1,7 +1,7 @@
<!DOCTYPE html>
-<html lang="en"><head><meta charset="UTF-8"/><meta name="viewport" content="width=device-width, initial-scale=1.0"/><title>internals · PhyloCoalSimulations.jl</title><meta name="title" content="internals · PhyloCoalSimulations.jl"/><meta property="og:title" content="internals · PhyloCoalSimulations.jl"/><meta property="twitter:title" content="internals · PhyloCoalSimulations.jl"/><meta name="description" content="Documentation for PhyloCoalSimulations.jl."/><meta property="og:description" content="Documentation for PhyloCoalSimulations.jl."/><meta property="twitter:description" content="Documentation for PhyloCoalSimulations.jl."/><meta property="og:url" content="https://cecileane.github.io/PhyloCoalSimulations.jl/stable/lib/internal/"/><meta property="twitter:url" content="https://cecileane.github.io/PhyloCoalSimulations.jl/stable/lib/internal/"/><link rel="canonical" href="https://cecileane.github.io/PhyloCoalSimulations.jl/stable/lib/internal/"/><script data-outdated-warner src="../../assets/warner.js"></script><link href="https://cdnjs.cloudflare.com/ajax/libs/lato-font/3.0.0/css/lato-font.min.css" rel="stylesheet" type="text/css"/><link href="https://cdnjs.cloudflare.com/ajax/libs/juliamono/0.050/juliamono.min.css" rel="stylesheet" type="text/css"/><link href="https://cdnjs.cloudflare.com/ajax/libs/font-awesome/6.4.2/css/fontawesome.min.css" rel="stylesheet" type="text/css"/><link href="https://cdnjs.cloudflare.com/ajax/libs/font-awesome/6.4.2/css/solid.min.css" rel="stylesheet" type="text/css"/><link href="https://cdnjs.cloudflare.com/ajax/libs/font-awesome/6.4.2/css/brands.min.css" rel="stylesheet" type="text/css"/><link href="https://cdnjs.cloudflare.com/ajax/libs/KaTeX/0.16.8/katex.min.css" rel="stylesheet" type="text/css"/><script>documenterBaseURL="../.."</script><script src="https://cdnjs.cloudflare.com/ajax/libs/require.js/2.3.6/require.min.js" data-main="../../assets/documenter.js"></script><script src="../../search_index.js"></script><script src="../../siteinfo.js"></script><script src="../../../versions.js"></script><link class="docs-theme-link" rel="stylesheet" type="text/css" href="../../assets/themes/documenter-dark.css" data-theme-name="documenter-dark" data-theme-primary-dark/><link class="docs-theme-link" rel="stylesheet" type="text/css" href="../../assets/themes/documenter-light.css" data-theme-name="documenter-light" data-theme-primary/><script src="../../assets/themeswap.js"></script></head><body><div id="documenter"><nav class="docs-sidebar"><div class="docs-package-name"><span class="docs-autofit"><a href="../../">PhyloCoalSimulations.jl</a></span></div><button class="docs-search-query input is-rounded is-small is-clickable my-2 mx-auto py-1 px-2" id="documenter-search-query">Search docs (Ctrl + /)</button><ul class="docs-menu"><li><a class="tocitem" href="../../">home</a></li><li><span class="tocitem">manual</span><ul><li><a class="tocitem" href="../../man/getting_started/">getting started</a></li><li><a class="tocitem" href="../../man/mapping_genetree_to_network/">mapping gene trees into the species network</a></li><li><a class="tocitem" href="../../man/converting_coal2generation_units/">converting between units</a></li><li><a class="tocitem" href="../../man/correlated_inheritance/">correlated inheritance</a></li><li><a class="tocitem" href="../../man/more_examples/">more examples</a></li></ul></li><li><span class="tocitem">library</span><ul><li><a class="tocitem" href="../public/">public</a></li><li class="is-active"><a class="tocitem" href>internals</a><ul class="internal"><li><a class="tocitem" href="#functions-and-types"><span>functions & types</span></a></li><li><a class="tocitem" href="#index"><span>index</span></a></li></ul></li></ul></li></ul><div class="docs-version-selector field has-addons"><div class="control"><span class="docs-label button is-static is-size-7">Version</span></div><div class="docs-selector control is-expanded"><div class="select is-fullwidth is-size-7"><select id="documenter-version-selector"></select></div></div></div></nav><div class="docs-main"><header class="docs-navbar"><a class="docs-sidebar-button docs-navbar-link fa-solid fa-bars is-hidden-desktop" id="documenter-sidebar-button" href="#"></a><nav class="breadcrumb"><ul class="is-hidden-mobile"><li><a class="is-disabled">library</a></li><li class="is-active"><a href>internals</a></li></ul><ul class="is-hidden-tablet"><li class="is-active"><a href>internals</a></li></ul></nav><div class="docs-right"><a class="docs-navbar-link" href="https://github.com/cecileane/PhyloCoalSimulations.jl" title="View the repository on GitHub"><span class="docs-icon fa-brands"></span><span class="docs-label is-hidden-touch">GitHub</span></a><a class="docs-navbar-link" href="https://github.com/cecileane/PhyloCoalSimulations.jl/blob/main/docs/src/lib/internal.md" title="Edit source on GitHub"><span class="docs-icon fa-solid"></span></a><a class="docs-settings-button docs-navbar-link fa-solid fa-gear" id="documenter-settings-button" href="#" title="Settings"></a><a class="docs-article-toggle-button fa-solid fa-chevron-up" id="documenter-article-toggle-button" href="javascript:;" title="Collapse all docstrings"></a></div></header><article class="content" id="documenter-page"><h1 id="internal-documentation"><a class="docs-heading-anchor" href="#internal-documentation">internal documentation</a><a id="internal-documentation-1"></a><a class="docs-heading-anchor-permalink" href="#internal-documentation" title="Permalink"></a></h1><p>Documentation for <code>PhyloCoalSimulations</code>'s internal functions. These functions are not exported and their access (API) should not be considered stable. But they can still be used, like this for example: <code>PhyloCoalSimulations.foo()</code> for a function named <code>foo()</code>.</p><h2 id="functions-and-types"><a class="docs-heading-anchor" href="#functions-and-types">functions & types</a><a id="functions-and-types-1"></a><a class="docs-heading-anchor-permalink" href="#functions-and-types" title="Permalink"></a></h2><article class="docstring"><header><a class="docstring-article-toggle-button fa-solid fa-chevron-down" href="javascript:;" title="Collapse docstring"></a><a class="docstring-binding" id="PhyloCoalSimulations.mappingnodes" href="#PhyloCoalSimulations.mappingnodes"><code>PhyloCoalSimulations.mappingnodes</code></a> — <span class="docstring-category">Type</span></header><section><div><pre><code class="language-julia hljs">mappingnodes(gene tree)</code></pre><p>Type to define an iterator over degree-2 mapping nodes in a gene tree, assuming these degree-2 nodes (other than the root) have a name to map them to nodes in a species phylogeny. See <a href="#PhyloCoalSimulations.mappingnodes"><code>ismappingnode</code></a>.</p></div><a class="docs-sourcelink" target="_blank" href="https://github.com/cecileane/PhyloCoalSimulations.jl/blob/d9bc44ccf207f472b5ee0672e997ebe5af08b3ac/src/utils.jl#L19-L25">source</a></section></article><article class="docstring"><header><a class="docstring-article-toggle-button fa-solid fa-chevron-down" href="javascript:;" title="Collapse docstring"></a><a class="docstring-binding" id="PhyloCoalSimulations.coalescence_edge-NTuple{4, Any}" href="#PhyloCoalSimulations.coalescence_edge-NTuple{4, Any}"><code>PhyloCoalSimulations.coalescence_edge</code></a> — <span class="docstring-category">Method</span></header><section><div><pre><code class="language-julia hljs">coalescence_edge(edge1, edge2, number, populationid)</code></pre><p>Create a coalescence between edges 1 and 2: with a new parent node <code>n</code> numbered <code>number</code> and a new parent edge <code>e</code> above the parent node, of length 0 and numbered <code>number</code>. Both <code>n.inCycle</code> and <code>e.inCycle</code> are set to <code>populationid</code>.</p></div><a class="docs-sourcelink" target="_blank" href="https://github.com/cecileane/PhyloCoalSimulations.jl/blob/d9bc44ccf207f472b5ee0672e997ebe5af08b3ac/src/simulatecoalescent_onepop.jl#L95-L102">source</a></section></article><article class="docstring"><header><a class="docstring-article-toggle-button fa-solid fa-chevron-down" href="javascript:;" title="Collapse docstring"></a><a class="docstring-binding" id="PhyloCoalSimulations.convert2tree!-Tuple{PhyloNetworks.Node}" href="#PhyloCoalSimulations.convert2tree!-Tuple{PhyloNetworks.Node}"><code>PhyloCoalSimulations.convert2tree!</code></a> — <span class="docstring-category">Method</span></header><section><div><pre><code class="language-julia hljs">convert2tree!(rootnode)</code></pre><p>Return a network with all nodes and edges that can be reached from <code>rootnode</code>. <strong>Warning</strong>: Assumes that edges are correctly directed (with correct <code>isChild1</code> attribute) and that the graph is a tree. This is <em>not</em> checked.</p><p>If the root node is still attached to an incomplete root edge, this edge & node are first disconnected.</p></div><a class="docs-sourcelink" target="_blank" href="https://github.com/cecileane/PhyloCoalSimulations.jl/blob/d9bc44ccf207f472b5ee0672e997ebe5af08b3ac/src/simulatecoalescent_onepop.jl#L249-L258">source</a></section></article><article class="docstring"><header><a class="docstring-article-toggle-button fa-solid fa-chevron-down" href="javascript:;" title="Collapse docstring"></a><a class="docstring-binding" id="PhyloCoalSimulations.get_rootedgenumber-Tuple{Any}" href="#PhyloCoalSimulations.get_rootedgenumber-Tuple{Any}"><code>PhyloCoalSimulations.get_rootedgenumber</code></a> — <span class="docstring-category">Method</span></header><section><div><pre><code class="language-julia hljs">get_rootedgenumber(network)</code></pre><p>1 + maximum of 0 and of all the network's edge numbers: can be used as a unique identifier of the edge above the network's root.</p></div><a class="docs-sourcelink" target="_blank" href="https://github.com/cecileane/PhyloCoalSimulations.jl/blob/d9bc44ccf207f472b5ee0672e997ebe5af08b3ac/src/simulatecoalescent_network.jl#L1-L6">source</a></section></article><article class="docstring"><header><a class="docstring-article-toggle-button fa-solid fa-chevron-down" href="javascript:;" title="Collapse docstring"></a><a class="docstring-binding" id="PhyloCoalSimulations.initializetip" href="#PhyloCoalSimulations.initializetip"><code>PhyloCoalSimulations.initializetip</code></a> — <span class="docstring-category">Function</span></header><section><div><pre><code class="language-julia hljs">initializetip(species::AbstractString, individual::AbstractString,
- number::Integer, delim=""::AbstractString)</code></pre><p>Create a leaf node and a pendant edge of length 0, incident to each other, both numbered <code>number</code>. Return the pendant edge. The leaf name is made by concatenating <code>species</code>, <code>delim</code> and <code>individual</code>.</p></div><a class="docs-sourcelink" target="_blank" href="https://github.com/cecileane/PhyloCoalSimulations.jl/blob/d9bc44ccf207f472b5ee0672e997ebe5af08b3ac/src/simulatecoalescent_onepop.jl#L119-L126">source</a></section></article><article class="docstring"><header><a class="docstring-article-toggle-button fa-solid fa-chevron-down" href="javascript:;" title="Collapse docstring"></a><a class="docstring-binding" id="PhyloCoalSimulations.initializetipforest" href="#PhyloCoalSimulations.initializetipforest"><code>PhyloCoalSimulations.initializetipforest</code></a> — <span class="docstring-category">Function</span></header><section><div><pre><code class="language-julia hljs">initializetipforest(speciesnode::Node, nindividuals::Integer,
- number::Integer, delim)</code></pre><p>Vector of pendant leaf edges, with leaves named after <code>speciesnode</code>, and numbered with consecutive number IDs starting at <code>number</code>. If nindividuals is 1, then the leaf name is simply the species name. Otherwise, then the leaf names include the individual number and the default delimiter is <code>_</code>. For example, if the species name is <code>s</code> then leaf names are: <code>s_1</code>, <code>s_2</code>, etc. by default. Pendant leaf edges have inCycle set to the number of the corresponding edge in the species network.</p></div><a class="docs-sourcelink" target="_blank" href="https://github.com/cecileane/PhyloCoalSimulations.jl/blob/d9bc44ccf207f472b5ee0672e997ebe5af08b3ac/src/simulatecoalescent_onepop.jl#L142-L154">source</a></section></article><article class="docstring"><header><a class="docstring-article-toggle-button fa-solid fa-chevron-down" href="javascript:;" title="Collapse docstring"></a><a class="docstring-binding" id="PhyloCoalSimulations.ismappingnode-Tuple{PhyloNetworks.Node}" href="#PhyloCoalSimulations.ismappingnode-Tuple{PhyloNetworks.Node}"><code>PhyloCoalSimulations.ismappingnode</code></a> — <span class="docstring-category">Method</span></header><section><div><pre><code class="language-julia hljs">ismappingnode(node)</code></pre><p>Boolean: true if <code>node</code> is of degree 2, has a single child, and has a name. (The root is of degree-2 but is not a mapping node).</p></div><a class="docs-sourcelink" target="_blank" href="https://github.com/cecileane/PhyloCoalSimulations.jl/blob/d9bc44ccf207f472b5ee0672e997ebe5af08b3ac/src/utils.jl#L11-L16">source</a></section></article><article class="docstring"><header><a class="docstring-article-toggle-button fa-solid fa-chevron-down" href="javascript:;" title="Collapse docstring"></a><a class="docstring-binding" id="PhyloCoalSimulations.map2population!-NTuple{4, Any}" href="#PhyloCoalSimulations.map2population!-NTuple{4, Any}"><code>PhyloCoalSimulations.map2population!</code></a> — <span class="docstring-category">Method</span></header><section><div><pre><code class="language-julia hljs">map2population!(forest, population_node, populationid, nextlineageID)</code></pre><p>Extend each incomplete edge in the forest with a new degree-2 node <code>n</code> and a new incomplete edge <code>e</code>, with the following information to map <code>n</code> and <code>e</code> into the species phylogeny:</p><ul><li><code>e.inCycle</code> is set to <code>populationid</code>, and</li><li><code>n.name</code> is set to <code>population_node.name</code> if this name is non-empty, or <code>string(population_node.number)</code> otherwise (with any negative sign replaced by the string "minus").</li></ul><p><code>e.number</code> and <code>n.number</code> are set to <code>nextlineageID</code>, which is incremented by 1 for each incomplete edge in the forest.</p><p>The forest is updated to contain the newly-created incomplete edges, replacing the old incomplete (and now complete) edges.</p><p>Output: nextlineageID, incremented by the number of newly created degree-2 lineages.</p><p><strong>example</strong></p><pre><code class="language-julia-repl hljs">julia> using PhyloNetworks; net = readTopology("(A:1,B:1);");
+<html lang="en"><head><meta charset="UTF-8"/><meta name="viewport" content="width=device-width, initial-scale=1.0"/><title>internals · PhyloCoalSimulations.jl</title><meta name="title" content="internals · PhyloCoalSimulations.jl"/><meta property="og:title" content="internals · PhyloCoalSimulations.jl"/><meta property="twitter:title" content="internals · PhyloCoalSimulations.jl"/><meta name="description" content="Documentation for PhyloCoalSimulations.jl."/><meta property="og:description" content="Documentation for PhyloCoalSimulations.jl."/><meta property="twitter:description" content="Documentation for PhyloCoalSimulations.jl."/><meta property="og:url" content="https://cecileane.github.io/PhyloCoalSimulations.jl/stable/lib/internal/"/><meta property="twitter:url" content="https://cecileane.github.io/PhyloCoalSimulations.jl/stable/lib/internal/"/><link rel="canonical" href="https://cecileane.github.io/PhyloCoalSimulations.jl/stable/lib/internal/"/><script data-outdated-warner src="../../assets/warner.js"></script><link href="https://cdnjs.cloudflare.com/ajax/libs/lato-font/3.0.0/css/lato-font.min.css" rel="stylesheet" type="text/css"/><link href="https://cdnjs.cloudflare.com/ajax/libs/juliamono/0.050/juliamono.min.css" rel="stylesheet" type="text/css"/><link href="https://cdnjs.cloudflare.com/ajax/libs/font-awesome/6.4.2/css/fontawesome.min.css" rel="stylesheet" type="text/css"/><link href="https://cdnjs.cloudflare.com/ajax/libs/font-awesome/6.4.2/css/solid.min.css" rel="stylesheet" type="text/css"/><link href="https://cdnjs.cloudflare.com/ajax/libs/font-awesome/6.4.2/css/brands.min.css" rel="stylesheet" type="text/css"/><link href="https://cdnjs.cloudflare.com/ajax/libs/KaTeX/0.16.8/katex.min.css" rel="stylesheet" type="text/css"/><script>documenterBaseURL="../.."</script><script src="https://cdnjs.cloudflare.com/ajax/libs/require.js/2.3.6/require.min.js" data-main="../../assets/documenter.js"></script><script src="../../search_index.js"></script><script src="../../siteinfo.js"></script><script src="../../../versions.js"></script><link class="docs-theme-link" rel="stylesheet" type="text/css" href="../../assets/themes/documenter-dark.css" data-theme-name="documenter-dark" data-theme-primary-dark/><link class="docs-theme-link" rel="stylesheet" type="text/css" href="../../assets/themes/documenter-light.css" data-theme-name="documenter-light" data-theme-primary/><script src="../../assets/themeswap.js"></script></head><body><div id="documenter"><nav class="docs-sidebar"><div class="docs-package-name"><span class="docs-autofit"><a href="../../">PhyloCoalSimulations.jl</a></span></div><button class="docs-search-query input is-rounded is-small is-clickable my-2 mx-auto py-1 px-2" id="documenter-search-query">Search docs (Ctrl + /)</button><ul class="docs-menu"><li><a class="tocitem" href="../../">home</a></li><li><span class="tocitem">manual</span><ul><li><a class="tocitem" href="../../man/getting_started/">getting started</a></li><li><a class="tocitem" href="../../man/mapping_genetree_to_network/">mapping gene trees into the species network</a></li><li><a class="tocitem" href="../../man/converting_coal2generation_units/">converting between units</a></li><li><a class="tocitem" href="../../man/correlated_inheritance/">correlated inheritance</a></li><li><a class="tocitem" href="../../man/more_examples/">more examples</a></li></ul></li><li><span class="tocitem">library</span><ul><li><a class="tocitem" href="../public/">public</a></li><li class="is-active"><a class="tocitem" href>internals</a><ul class="internal"><li><a class="tocitem" href="#functions-and-types"><span>functions & types</span></a></li><li><a class="tocitem" href="#index"><span>index</span></a></li></ul></li></ul></li></ul><div class="docs-version-selector field has-addons"><div class="control"><span class="docs-label button is-static is-size-7">Version</span></div><div class="docs-selector control is-expanded"><div class="select is-fullwidth is-size-7"><select id="documenter-version-selector"></select></div></div></div></nav><div class="docs-main"><header class="docs-navbar"><a class="docs-sidebar-button docs-navbar-link fa-solid fa-bars is-hidden-desktop" id="documenter-sidebar-button" href="#"></a><nav class="breadcrumb"><ul class="is-hidden-mobile"><li><a class="is-disabled">library</a></li><li class="is-active"><a href>internals</a></li></ul><ul class="is-hidden-tablet"><li class="is-active"><a href>internals</a></li></ul></nav><div class="docs-right"><a class="docs-navbar-link" href="https://github.com/cecileane/PhyloCoalSimulations.jl" title="View the repository on GitHub"><span class="docs-icon fa-brands"></span><span class="docs-label is-hidden-touch">GitHub</span></a><a class="docs-navbar-link" href="https://github.com/cecileane/PhyloCoalSimulations.jl/blob/main/docs/src/lib/internal.md" title="Edit source on GitHub"><span class="docs-icon fa-solid"></span></a><a class="docs-settings-button docs-navbar-link fa-solid fa-gear" id="documenter-settings-button" href="#" title="Settings"></a><a class="docs-article-toggle-button fa-solid fa-chevron-up" id="documenter-article-toggle-button" href="javascript:;" title="Collapse all docstrings"></a></div></header><article class="content" id="documenter-page"><h1 id="internal-documentation"><a class="docs-heading-anchor" href="#internal-documentation">internal documentation</a><a id="internal-documentation-1"></a><a class="docs-heading-anchor-permalink" href="#internal-documentation" title="Permalink"></a></h1><p>Documentation for <code>PhyloCoalSimulations</code>'s internal functions. These functions are not exported and their access (API) should not be considered stable. But they can still be used, like this for example: <code>PhyloCoalSimulations.foo()</code> for a function named <code>foo()</code>.</p><h2 id="functions-and-types"><a class="docs-heading-anchor" href="#functions-and-types">functions & types</a><a id="functions-and-types-1"></a><a class="docs-heading-anchor-permalink" href="#functions-and-types" title="Permalink"></a></h2><article class="docstring"><header><a class="docstring-article-toggle-button fa-solid fa-chevron-down" href="javascript:;" title="Collapse docstring"></a><a class="docstring-binding" id="PhyloCoalSimulations.mappingnodes" href="#PhyloCoalSimulations.mappingnodes"><code>PhyloCoalSimulations.mappingnodes</code></a> — <span class="docstring-category">Type</span></header><section><div><pre><code class="language-julia hljs">mappingnodes(gene tree)</code></pre><p>Type to define an iterator over degree-2 mapping nodes in a gene tree, assuming these degree-2 nodes (other than the root) have a name to map them to nodes in a species phylogeny. See <a href="#PhyloCoalSimulations.mappingnodes"><code>ismappingnode</code></a>.</p></div><a class="docs-sourcelink" target="_blank" href="https://github.com/cecileane/PhyloCoalSimulations.jl/blob/7ddfd0974ad019ce6407a497f7048aebe963bb36/src/utils.jl#L19-L25">source</a></section></article><article class="docstring"><header><a class="docstring-article-toggle-button fa-solid fa-chevron-down" href="javascript:;" title="Collapse docstring"></a><a class="docstring-binding" id="PhyloCoalSimulations.coalescence_edge-NTuple{4, Any}" href="#PhyloCoalSimulations.coalescence_edge-NTuple{4, Any}"><code>PhyloCoalSimulations.coalescence_edge</code></a> — <span class="docstring-category">Method</span></header><section><div><pre><code class="language-julia hljs">coalescence_edge(edge1, edge2, number, populationid)</code></pre><p>Create a coalescence between edges 1 and 2: with a new parent node <code>n</code> numbered <code>number</code> and a new parent edge <code>e</code> above the parent node, of length 0 and numbered <code>number</code>. Both <code>n.inCycle</code> and <code>e.inCycle</code> are set to <code>populationid</code>.</p></div><a class="docs-sourcelink" target="_blank" href="https://github.com/cecileane/PhyloCoalSimulations.jl/blob/7ddfd0974ad019ce6407a497f7048aebe963bb36/src/simulatecoalescent_onepop.jl#L95-L102">source</a></section></article><article class="docstring"><header><a class="docstring-article-toggle-button fa-solid fa-chevron-down" href="javascript:;" title="Collapse docstring"></a><a class="docstring-binding" id="PhyloCoalSimulations.convert2tree!-Tuple{PhyloNetworks.Node}" href="#PhyloCoalSimulations.convert2tree!-Tuple{PhyloNetworks.Node}"><code>PhyloCoalSimulations.convert2tree!</code></a> — <span class="docstring-category">Method</span></header><section><div><pre><code class="language-julia hljs">convert2tree!(rootnode)</code></pre><p>Return a network with all nodes and edges that can be reached from <code>rootnode</code>. <strong>Warning</strong>: Assumes that edges are correctly directed (with correct <code>isChild1</code> attribute) and that the graph is a tree. This is <em>not</em> checked.</p><p>If the root node is still attached to an incomplete root edge, this edge & node are first disconnected.</p></div><a class="docs-sourcelink" target="_blank" href="https://github.com/cecileane/PhyloCoalSimulations.jl/blob/7ddfd0974ad019ce6407a497f7048aebe963bb36/src/simulatecoalescent_onepop.jl#L249-L258">source</a></section></article><article class="docstring"><header><a class="docstring-article-toggle-button fa-solid fa-chevron-down" href="javascript:;" title="Collapse docstring"></a><a class="docstring-binding" id="PhyloCoalSimulations.get_rootedgenumber-Tuple{Any}" href="#PhyloCoalSimulations.get_rootedgenumber-Tuple{Any}"><code>PhyloCoalSimulations.get_rootedgenumber</code></a> — <span class="docstring-category">Method</span></header><section><div><pre><code class="language-julia hljs">get_rootedgenumber(network)</code></pre><p>1 + maximum of 0 and of all the network's edge numbers: can be used as a unique identifier of the edge above the network's root.</p></div><a class="docs-sourcelink" target="_blank" href="https://github.com/cecileane/PhyloCoalSimulations.jl/blob/7ddfd0974ad019ce6407a497f7048aebe963bb36/src/simulatecoalescent_network.jl#L1-L6">source</a></section></article><article class="docstring"><header><a class="docstring-article-toggle-button fa-solid fa-chevron-down" href="javascript:;" title="Collapse docstring"></a><a class="docstring-binding" id="PhyloCoalSimulations.initializetip" href="#PhyloCoalSimulations.initializetip"><code>PhyloCoalSimulations.initializetip</code></a> — <span class="docstring-category">Function</span></header><section><div><pre><code class="language-julia hljs">initializetip(species::AbstractString, individual::AbstractString,
+ number::Integer, delim=""::AbstractString)</code></pre><p>Create a leaf node and a pendant edge of length 0, incident to each other, both numbered <code>number</code>. Return the pendant edge. The leaf name is made by concatenating <code>species</code>, <code>delim</code> and <code>individual</code>.</p></div><a class="docs-sourcelink" target="_blank" href="https://github.com/cecileane/PhyloCoalSimulations.jl/blob/7ddfd0974ad019ce6407a497f7048aebe963bb36/src/simulatecoalescent_onepop.jl#L119-L126">source</a></section></article><article class="docstring"><header><a class="docstring-article-toggle-button fa-solid fa-chevron-down" href="javascript:;" title="Collapse docstring"></a><a class="docstring-binding" id="PhyloCoalSimulations.initializetipforest" href="#PhyloCoalSimulations.initializetipforest"><code>PhyloCoalSimulations.initializetipforest</code></a> — <span class="docstring-category">Function</span></header><section><div><pre><code class="language-julia hljs">initializetipforest(speciesnode::Node, nindividuals::Integer,
+ number::Integer, delim)</code></pre><p>Vector of pendant leaf edges, with leaves named after <code>speciesnode</code>, and numbered with consecutive number IDs starting at <code>number</code>. If nindividuals is 1, then the leaf name is simply the species name. Otherwise, then the leaf names include the individual number and the default delimiter is <code>_</code>. For example, if the species name is <code>s</code> then leaf names are: <code>s_1</code>, <code>s_2</code>, etc. by default. Pendant leaf edges have inCycle set to the number of the corresponding edge in the species network.</p></div><a class="docs-sourcelink" target="_blank" href="https://github.com/cecileane/PhyloCoalSimulations.jl/blob/7ddfd0974ad019ce6407a497f7048aebe963bb36/src/simulatecoalescent_onepop.jl#L142-L154">source</a></section></article><article class="docstring"><header><a class="docstring-article-toggle-button fa-solid fa-chevron-down" href="javascript:;" title="Collapse docstring"></a><a class="docstring-binding" id="PhyloCoalSimulations.ismappingnode-Tuple{PhyloNetworks.Node}" href="#PhyloCoalSimulations.ismappingnode-Tuple{PhyloNetworks.Node}"><code>PhyloCoalSimulations.ismappingnode</code></a> — <span class="docstring-category">Method</span></header><section><div><pre><code class="language-julia hljs">ismappingnode(node)</code></pre><p>Boolean: true if <code>node</code> is of degree 2, has a single child, and has a name. (The root is of degree-2 but is not a mapping node).</p></div><a class="docs-sourcelink" target="_blank" href="https://github.com/cecileane/PhyloCoalSimulations.jl/blob/7ddfd0974ad019ce6407a497f7048aebe963bb36/src/utils.jl#L11-L16">source</a></section></article><article class="docstring"><header><a class="docstring-article-toggle-button fa-solid fa-chevron-down" href="javascript:;" title="Collapse docstring"></a><a class="docstring-binding" id="PhyloCoalSimulations.map2population!-NTuple{4, Any}" href="#PhyloCoalSimulations.map2population!-NTuple{4, Any}"><code>PhyloCoalSimulations.map2population!</code></a> — <span class="docstring-category">Method</span></header><section><div><pre><code class="language-julia hljs">map2population!(forest, population_node, populationid, nextlineageID)</code></pre><p>Extend each incomplete edge in the forest with a new degree-2 node <code>n</code> and a new incomplete edge <code>e</code>, with the following information to map <code>n</code> and <code>e</code> into the species phylogeny:</p><ul><li><code>e.inCycle</code> is set to <code>populationid</code>, and</li><li><code>n.name</code> is set to <code>population_node.name</code> if this name is non-empty, or <code>string(population_node.number)</code> otherwise (with any negative sign replaced by the string "minus").</li></ul><p><code>e.number</code> and <code>n.number</code> are set to <code>nextlineageID</code>, which is incremented by 1 for each incomplete edge in the forest.</p><p>The forest is updated to contain the newly-created incomplete edges, replacing the old incomplete (and now complete) edges.</p><p>Output: nextlineageID, incremented by the number of newly created degree-2 lineages.</p><p><strong>example</strong></p><pre><code class="language-julia-repl hljs">julia> using PhyloNetworks; net = readTopology("(A:1,B:1);");
julia> leafA = net.node[1]; edge2A_number = net.edge[1].number;
@@ -22,7 +22,7 @@
julia> [e.node[1].name for e in f]
2-element Vector{String}:
"A"
- "A"</code></pre></div><a class="docs-sourcelink" target="_blank" href="https://github.com/cecileane/PhyloCoalSimulations.jl/blob/d9bc44ccf207f472b5ee0672e997ebe5af08b3ac/src/simulatecoalescent_onepop.jl#L173-L217">source</a></section></article><article class="docstring"><header><a class="docstring-article-toggle-button fa-solid fa-chevron-down" href="javascript:;" title="Collapse docstring"></a><a class="docstring-binding" id="PhyloCoalSimulations.population_mappedto-Tuple{Union{PhyloNetworks.EdgeT{PhyloNetworks.Node}, PhyloNetworks.Node}}" href="#PhyloCoalSimulations.population_mappedto-Tuple{Union{PhyloNetworks.EdgeT{PhyloNetworks.Node}, PhyloNetworks.Node}}"><code>PhyloCoalSimulations.population_mappedto</code></a> — <span class="docstring-category">Method</span></header><section><div><pre><code class="language-julia hljs">population_mappedto(edge or node)</code></pre><p>Identifier of the population (edge in the species network) that a gene tree's edge or a node is mapped onto, or <code>nothing</code> if not mapped. For example, coalescent nodes in gene trees map to a <em>node</em> in the species phylogeny, instead of mapping to an <em>edge</em>.</p></div><a class="docs-sourcelink" target="_blank" href="https://github.com/cecileane/PhyloCoalSimulations.jl/blob/d9bc44ccf207f472b5ee0672e997ebe5af08b3ac/src/utils.jl#L1-L8">source</a></section></article><article class="docstring"><header><a class="docstring-article-toggle-button fa-solid fa-chevron-down" href="javascript:;" title="Collapse docstring"></a><a class="docstring-binding" id="PhyloCoalSimulations.simulatecoal_onepopulation!-Tuple{AbstractVector, Vararg{Any}}" href="#PhyloCoalSimulations.simulatecoal_onepopulation!-Tuple{AbstractVector, Vararg{Any}}"><code>PhyloCoalSimulations.simulatecoal_onepopulation!</code></a> — <span class="docstring-category">Method</span></header><section><div><pre><code class="language-julia hljs">simulatecoal_onepopulation!([rng::AbstractRNG,]
+ "A"</code></pre></div><a class="docs-sourcelink" target="_blank" href="https://github.com/cecileane/PhyloCoalSimulations.jl/blob/7ddfd0974ad019ce6407a497f7048aebe963bb36/src/simulatecoalescent_onepop.jl#L173-L217">source</a></section></article><article class="docstring"><header><a class="docstring-article-toggle-button fa-solid fa-chevron-down" href="javascript:;" title="Collapse docstring"></a><a class="docstring-binding" id="PhyloCoalSimulations.population_mappedto-Tuple{Union{PhyloNetworks.EdgeT{PhyloNetworks.Node}, PhyloNetworks.Node}}" href="#PhyloCoalSimulations.population_mappedto-Tuple{Union{PhyloNetworks.EdgeT{PhyloNetworks.Node}, PhyloNetworks.Node}}"><code>PhyloCoalSimulations.population_mappedto</code></a> — <span class="docstring-category">Method</span></header><section><div><pre><code class="language-julia hljs">population_mappedto(edge or node)</code></pre><p>Identifier of the population (edge in the species network) that a gene tree's edge or a node is mapped onto, or <code>nothing</code> if not mapped. For example, coalescent nodes in gene trees map to a <em>node</em> in the species phylogeny, instead of mapping to an <em>edge</em>.</p></div><a class="docs-sourcelink" target="_blank" href="https://github.com/cecileane/PhyloCoalSimulations.jl/blob/7ddfd0974ad019ce6407a497f7048aebe963bb36/src/utils.jl#L1-L8">source</a></section></article><article class="docstring"><header><a class="docstring-article-toggle-button fa-solid fa-chevron-down" href="javascript:;" title="Collapse docstring"></a><a class="docstring-binding" id="PhyloCoalSimulations.simulatecoal_onepopulation!-Tuple{AbstractVector, Vararg{Any}}" href="#PhyloCoalSimulations.simulatecoal_onepopulation!-Tuple{AbstractVector, Vararg{Any}}"><code>PhyloCoalSimulations.simulatecoal_onepopulation!</code></a> — <span class="docstring-category">Method</span></header><section><div><pre><code class="language-julia hljs">simulatecoal_onepopulation!([rng::AbstractRNG,]
lineagelist,
population_length,
nextlineageID,
@@ -46,4 +46,4 @@
3 nodes: 2 tips, 0 hybrid nodes, 1 internal tree nodes.
tip labels: s2, s1
(s2:0.302,s1:0.302);
-</code></pre></div><a class="docs-sourcelink" target="_blank" href="https://github.com/cecileane/PhyloCoalSimulations.jl/blob/d9bc44ccf207f472b5ee0672e997ebe5af08b3ac/src/simulatecoalescent_onepop.jl#L1-L51">source</a></section></article><h2 id="index"><a class="docs-heading-anchor" href="#index">index</a><a id="index-1"></a><a class="docs-heading-anchor-permalink" href="#index" title="Permalink"></a></h2><ul><li><a href="#PhyloCoalSimulations.mappingnodes"><code>PhyloCoalSimulations.mappingnodes</code></a></li><li><a href="#PhyloCoalSimulations.coalescence_edge-NTuple{4, Any}"><code>PhyloCoalSimulations.coalescence_edge</code></a></li><li><a href="#PhyloCoalSimulations.convert2tree!-Tuple{PhyloNetworks.Node}"><code>PhyloCoalSimulations.convert2tree!</code></a></li><li><a href="#PhyloCoalSimulations.get_rootedgenumber-Tuple{Any}"><code>PhyloCoalSimulations.get_rootedgenumber</code></a></li><li><a href="#PhyloCoalSimulations.initializetip"><code>PhyloCoalSimulations.initializetip</code></a></li><li><a href="#PhyloCoalSimulations.initializetipforest"><code>PhyloCoalSimulations.initializetipforest</code></a></li><li><a href="#PhyloCoalSimulations.ismappingnode-Tuple{PhyloNetworks.Node}"><code>PhyloCoalSimulations.ismappingnode</code></a></li><li><a href="#PhyloCoalSimulations.map2population!-NTuple{4, Any}"><code>PhyloCoalSimulations.map2population!</code></a></li><li><a href="#PhyloCoalSimulations.population_mappedto-Tuple{Union{PhyloNetworks.EdgeT{PhyloNetworks.Node}, PhyloNetworks.Node}}"><code>PhyloCoalSimulations.population_mappedto</code></a></li><li><a href="#PhyloCoalSimulations.simulatecoal_onepopulation!-Tuple{AbstractVector, Vararg{Any}}"><code>PhyloCoalSimulations.simulatecoal_onepopulation!</code></a></li></ul></article><nav class="docs-footer"><a class="docs-footer-prevpage" href="../public/">« public</a><div class="flexbox-break"></div><p class="footer-message">Powered by <a href="https://github.com/JuliaDocs/Documenter.jl">Documenter.jl</a> and the <a href="https://julialang.org/">Julia Programming Language</a>.</p></nav></div><div class="modal" id="documenter-settings"><div class="modal-background"></div><div class="modal-card"><header class="modal-card-head"><p class="modal-card-title">Settings</p><button class="delete"></button></header><section class="modal-card-body"><p><label class="label">Theme</label><div class="select"><select id="documenter-themepicker"><option value="auto">Automatic (OS)</option><option value="documenter-light">documenter-light</option><option value="documenter-dark">documenter-dark</option></select></div></p><hr/><p>This document was generated with <a href="https://github.com/JuliaDocs/Documenter.jl">Documenter.jl</a> version 1.4.1 on <span class="colophon-date" title="Friday 17 May 2024 03:10">Friday 17 May 2024</span>. Using Julia version 1.10.3.</p></section><footer class="modal-card-foot"></footer></div></div></div></body></html>
+</code></pre></div><a class="docs-sourcelink" target="_blank" href="https://github.com/cecileane/PhyloCoalSimulations.jl/blob/7ddfd0974ad019ce6407a497f7048aebe963bb36/src/simulatecoalescent_onepop.jl#L1-L51">source</a></section></article><h2 id="index"><a class="docs-heading-anchor" href="#index">index</a><a id="index-1"></a><a class="docs-heading-anchor-permalink" href="#index" title="Permalink"></a></h2><ul><li><a href="#PhyloCoalSimulations.mappingnodes"><code>PhyloCoalSimulations.mappingnodes</code></a></li><li><a href="#PhyloCoalSimulations.coalescence_edge-NTuple{4, Any}"><code>PhyloCoalSimulations.coalescence_edge</code></a></li><li><a href="#PhyloCoalSimulations.convert2tree!-Tuple{PhyloNetworks.Node}"><code>PhyloCoalSimulations.convert2tree!</code></a></li><li><a href="#PhyloCoalSimulations.get_rootedgenumber-Tuple{Any}"><code>PhyloCoalSimulations.get_rootedgenumber</code></a></li><li><a href="#PhyloCoalSimulations.initializetip"><code>PhyloCoalSimulations.initializetip</code></a></li><li><a href="#PhyloCoalSimulations.initializetipforest"><code>PhyloCoalSimulations.initializetipforest</code></a></li><li><a href="#PhyloCoalSimulations.ismappingnode-Tuple{PhyloNetworks.Node}"><code>PhyloCoalSimulations.ismappingnode</code></a></li><li><a href="#PhyloCoalSimulations.map2population!-NTuple{4, Any}"><code>PhyloCoalSimulations.map2population!</code></a></li><li><a href="#PhyloCoalSimulations.population_mappedto-Tuple{Union{PhyloNetworks.EdgeT{PhyloNetworks.Node}, PhyloNetworks.Node}}"><code>PhyloCoalSimulations.population_mappedto</code></a></li><li><a href="#PhyloCoalSimulations.simulatecoal_onepopulation!-Tuple{AbstractVector, Vararg{Any}}"><code>PhyloCoalSimulations.simulatecoal_onepopulation!</code></a></li></ul></article><nav class="docs-footer"><a class="docs-footer-prevpage" href="../public/">« public</a><div class="flexbox-break"></div><p class="footer-message">Powered by <a href="https://github.com/JuliaDocs/Documenter.jl">Documenter.jl</a> and the <a href="https://julialang.org/">Julia Programming Language</a>.</p></nav></div><div class="modal" id="documenter-settings"><div class="modal-background"></div><div class="modal-card"><header class="modal-card-head"><p class="modal-card-title">Settings</p><button class="delete"></button></header><section class="modal-card-body"><p><label class="label">Theme</label><div class="select"><select id="documenter-themepicker"><option value="auto">Automatic (OS)</option><option value="documenter-light">documenter-light</option><option value="documenter-dark">documenter-dark</option></select></div></p><hr/><p>This document was generated with <a href="https://github.com/JuliaDocs/Documenter.jl">Documenter.jl</a> version 1.4.1 on <span class="colophon-date" title="Friday 17 May 2024 12:07">Friday 17 May 2024</span>. Using Julia version 1.10.3.</p></section><footer class="modal-card-foot"></footer></div></div></div></body></html>
diff --git a/previews/PR12/lib/public/index.html b/previews/PR12/lib/public/index.html
index caa37d5..16ed441 100644
--- a/previews/PR12/lib/public/index.html
+++ b/previews/PR12/lib/public/index.html
@@ -80,7 +80,7 @@
(tree_edge_number = 6, pop_edge_number = 3)
(tree_edge_number = 3, pop_edge_number = 1)
(tree_edge_number = 1, pop_edge_number = 1)
- (tree_edge_number = 2, pop_edge_number = 1)</code></pre></div><a class="docs-sourcelink" target="_blank" href="https://github.com/cecileane/PhyloCoalSimulations.jl/blob/d9bc44ccf207f472b5ee0672e997ebe5af08b3ac/src/simulatecoalescent_network.jl#L9-L162">source</a></section></article><article class="docstring"><header><a class="docstring-article-toggle-button fa-solid fa-chevron-down" href="javascript:;" title="Collapse docstring"></a><a class="docstring-binding" id="PhyloCoalSimulations.simulatecoalescent-Tuple{Random.AbstractRNG, PhyloNetworks.HybridNetwork, Integer, Any, Any}" href="#PhyloCoalSimulations.simulatecoalescent-Tuple{Random.AbstractRNG, PhyloNetworks.HybridNetwork, Integer, Any, Any}"><code>PhyloCoalSimulations.simulatecoalescent</code></a> — <span class="docstring-category">Method</span></header><section><div><pre><code class="language-julia hljs">simulatecoalescent([rng::AbstractRNG,] net, nloci, nindividuals, populationsize;
+ (tree_edge_number = 2, pop_edge_number = 1)</code></pre></div><a class="docs-sourcelink" target="_blank" href="https://github.com/cecileane/PhyloCoalSimulations.jl/blob/7ddfd0974ad019ce6407a497f7048aebe963bb36/src/simulatecoalescent_network.jl#L9-L162">source</a></section></article><article class="docstring"><header><a class="docstring-article-toggle-button fa-solid fa-chevron-down" href="javascript:;" title="Collapse docstring"></a><a class="docstring-binding" id="PhyloCoalSimulations.simulatecoalescent-Tuple{Random.AbstractRNG, PhyloNetworks.HybridNetwork, Integer, Any, Any}" href="#PhyloCoalSimulations.simulatecoalescent-Tuple{Random.AbstractRNG, PhyloNetworks.HybridNetwork, Integer, Any, Any}"><code>PhyloCoalSimulations.simulatecoalescent</code></a> — <span class="docstring-category">Method</span></header><section><div><pre><code class="language-julia hljs">simulatecoalescent([rng::AbstractRNG,] net, nloci, nindividuals, populationsize;
nodemapping=false, round_generationnumber=true,
inheritancecorrelation=0.0)</code></pre><p>Simulate <code>nloci</code> gene trees with <code>nindividuals</code> from each species under the multispecies network coalescent, along network <code>net</code>, whose branch lengths are assumed to be in <strong>number of generations</strong>. <code>populationsize</code> should be a single number, assumed to be the (haploid) effective population size Nₑ, constant across the species phylogeny. Alternatively, <code>populationsize</code> can be a dictionary mapping the number of each edge in <code>net</code> to its Nₑ, including an extra edge number for the population above the root of the network.</p><p>Coalescent units are then calculated as <code>u=g/Nₑ</code> where <code>g</code> is the edge length in <code>net</code> (number of generations), and the coalescent model is applied using the infinite-population-size approximation.</p><p>Output: vector of gene trees with edge lengths in number of generations, calculated as <code>g=uNₑ</code> and then rounded to be an integer, unless <code>round_generationnumber</code> is false.</p><div class="admonition is-warning"><header class="admonition-header">Warning</header><div class="admonition-body"><p>When <code>populationsize</code> Nₑ is not provided as input, all edge lengths are in coalescent units. When <code>populationsize</code> is given as an argument, all edge lengths are in number of generations. The second method (using # generation and Nₑ as input) is a wrapper around the first (using coalescent units).</p></div></div><pre><code class="language-julia-repl hljs">julia> using PhyloNetworks
@@ -104,4 +104,4 @@
julia> writeMultiTopology(genetrees, stdout) # branch lengths: number of generations
(B:546.0,A:546.0);
(B:3155.0,A:3155.0);
-</code></pre></div><a class="docs-sourcelink" target="_blank" href="https://github.com/cecileane/PhyloCoalSimulations.jl/blob/d9bc44ccf207f472b5ee0672e997ebe5af08b3ac/src/simulatecoalescent_network.jl#L295-L349">source</a></section></article><h2 id="index"><a class="docs-heading-anchor" href="#index">index</a><a id="index-1"></a><a class="docs-heading-anchor-permalink" href="#index" title="Permalink"></a></h2><ul><li><a href="#PhyloCoalSimulations.simulatecoalescent-Tuple{Random.AbstractRNG, PhyloNetworks.HybridNetwork, Integer, Any, Any}"><code>PhyloCoalSimulations.simulatecoalescent</code></a></li><li><a href="#PhyloCoalSimulations.simulatecoalescent-Tuple{PhyloNetworks.HybridNetwork, Vararg{Any}}"><code>PhyloCoalSimulations.simulatecoalescent</code></a></li></ul></article><nav class="docs-footer"><a class="docs-footer-prevpage" href="../../man/more_examples/">« more examples</a><a class="docs-footer-nextpage" href="../internal/">internals »</a><div class="flexbox-break"></div><p class="footer-message">Powered by <a href="https://github.com/JuliaDocs/Documenter.jl">Documenter.jl</a> and the <a href="https://julialang.org/">Julia Programming Language</a>.</p></nav></div><div class="modal" id="documenter-settings"><div class="modal-background"></div><div class="modal-card"><header class="modal-card-head"><p class="modal-card-title">Settings</p><button class="delete"></button></header><section class="modal-card-body"><p><label class="label">Theme</label><div class="select"><select id="documenter-themepicker"><option value="auto">Automatic (OS)</option><option value="documenter-light">documenter-light</option><option value="documenter-dark">documenter-dark</option></select></div></p><hr/><p>This document was generated with <a href="https://github.com/JuliaDocs/Documenter.jl">Documenter.jl</a> version 1.4.1 on <span class="colophon-date" title="Friday 17 May 2024 03:10">Friday 17 May 2024</span>. Using Julia version 1.10.3.</p></section><footer class="modal-card-foot"></footer></div></div></div></body></html>
+</code></pre></div><a class="docs-sourcelink" target="_blank" href="https://github.com/cecileane/PhyloCoalSimulations.jl/blob/7ddfd0974ad019ce6407a497f7048aebe963bb36/src/simulatecoalescent_network.jl#L295-L349">source</a></section></article><h2 id="index"><a class="docs-heading-anchor" href="#index">index</a><a id="index-1"></a><a class="docs-heading-anchor-permalink" href="#index" title="Permalink"></a></h2><ul><li><a href="#PhyloCoalSimulations.simulatecoalescent-Tuple{Random.AbstractRNG, PhyloNetworks.HybridNetwork, Integer, Any, Any}"><code>PhyloCoalSimulations.simulatecoalescent</code></a></li><li><a href="#PhyloCoalSimulations.simulatecoalescent-Tuple{PhyloNetworks.HybridNetwork, Vararg{Any}}"><code>PhyloCoalSimulations.simulatecoalescent</code></a></li></ul></article><nav class="docs-footer"><a class="docs-footer-prevpage" href="../../man/more_examples/">« more examples</a><a class="docs-footer-nextpage" href="../internal/">internals »</a><div class="flexbox-break"></div><p class="footer-message">Powered by <a href="https://github.com/JuliaDocs/Documenter.jl">Documenter.jl</a> and the <a href="https://julialang.org/">Julia Programming Language</a>.</p></nav></div><div class="modal" id="documenter-settings"><div class="modal-background"></div><div class="modal-card"><header class="modal-card-head"><p class="modal-card-title">Settings</p><button class="delete"></button></header><section class="modal-card-body"><p><label class="label">Theme</label><div class="select"><select id="documenter-themepicker"><option value="auto">Automatic (OS)</option><option value="documenter-light">documenter-light</option><option value="documenter-dark">documenter-dark</option></select></div></p><hr/><p>This document was generated with <a href="https://github.com/JuliaDocs/Documenter.jl">Documenter.jl</a> version 1.4.1 on <span class="colophon-date" title="Friday 17 May 2024 12:07">Friday 17 May 2024</span>. Using Julia version 1.10.3.</p></section><footer class="modal-card-foot"></footer></div></div></div></body></html>
diff --git a/previews/PR12/man/converting_coal2generation_units/index.html b/previews/PR12/man/converting_coal2generation_units/index.html
index 51a3f91..b2bb1d8 100644
--- a/previews/PR12/man/converting_coal2generation_units/index.html
+++ b/previews/PR12/man/converting_coal2generation_units/index.html
@@ -27,4 +27,4 @@
R"mtext"("red: Ne values", side=1, line=-1.5, col="red4");
R"mtext"("black: edge lengths", side=1, line=-0.5);</code></pre><p><img src="../../assets/figures/net_genNe.svg" alt="species net with Ne"/></p><p>To simulate gene trees with edge lengths in generations, we can use a convenience wrapper function that takes <strong>Nₑ as an extra input</strong> to:</p><ul><li>convert edge lengths to coalescent units in the species phylogeny,</li><li>simulate gene trees with lengths in coalescent units, then</li><li>convert gene trees to have lengths in number of generations:</li></ul><pre><code class="language-julia-repl hljs" style="display:block;">julia> genetree_gen = simulatecoalescent(net_gen,3,1, Ne; nodemapping=true);</code><code class="nohighlight hljs ansi" style="display:block;"></code><br/><code class="language-julia-repl hljs" style="display:block;">julia> writeMultiTopology(genetree_gen, stdout) # 3 gene trees, lengths in #generations</code><code class="nohighlight hljs ansi" style="display:block;">((((B:200.0)H1:700.0)i1:600.0)i3:758.0,(((C:900.0)i1:600.0)i3:132.0,((A:1000.0)i2:500.0)i3:132.0):625.0);
(((((B:200.0)H1:700.0)i1:600.0)i3:2119.0,((A:1000.0)i2:500.0)i3:2119.0):971.0,((C:900.0)i1:600.0)i3:3089.0);
-(((C:900.0)i1:600.0)i3:1203.0,(((A:1000.0)i2:500.0)i3:724.0,(((B:200.0)H1:700.0)i1:600.0)i3:724.0):479.0);</code></pre><div class="admonition is-warning"><header class="admonition-header">Warning</header><div class="admonition-body"><p>When Nₑ is given as an extra input to <code>simulatecoalescent</code>, edge lengths in the network are assumed to be in number of generations. If Nₑ is <em>not</em> given as input, then edge lengths are assumed to be in coalescent units.</p></div></div></article><nav class="docs-footer"><a class="docs-footer-prevpage" href="../mapping_genetree_to_network/">« mapping gene trees into the species network</a><a class="docs-footer-nextpage" href="../correlated_inheritance/">correlated inheritance »</a><div class="flexbox-break"></div><p class="footer-message">Powered by <a href="https://github.com/JuliaDocs/Documenter.jl">Documenter.jl</a> and the <a href="https://julialang.org/">Julia Programming Language</a>.</p></nav></div><div class="modal" id="documenter-settings"><div class="modal-background"></div><div class="modal-card"><header class="modal-card-head"><p class="modal-card-title">Settings</p><button class="delete"></button></header><section class="modal-card-body"><p><label class="label">Theme</label><div class="select"><select id="documenter-themepicker"><option value="auto">Automatic (OS)</option><option value="documenter-light">documenter-light</option><option value="documenter-dark">documenter-dark</option></select></div></p><hr/><p>This document was generated with <a href="https://github.com/JuliaDocs/Documenter.jl">Documenter.jl</a> version 1.4.1 on <span class="colophon-date" title="Friday 17 May 2024 03:10">Friday 17 May 2024</span>. Using Julia version 1.10.3.</p></section><footer class="modal-card-foot"></footer></div></div></div></body></html>
+(((C:900.0)i1:600.0)i3:1203.0,(((A:1000.0)i2:500.0)i3:724.0,(((B:200.0)H1:700.0)i1:600.0)i3:724.0):479.0);</code></pre><div class="admonition is-warning"><header class="admonition-header">Warning</header><div class="admonition-body"><p>When Nₑ is given as an extra input to <code>simulatecoalescent</code>, edge lengths in the network are assumed to be in number of generations. If Nₑ is <em>not</em> given as input, then edge lengths are assumed to be in coalescent units.</p></div></div></article><nav class="docs-footer"><a class="docs-footer-prevpage" href="../mapping_genetree_to_network/">« mapping gene trees into the species network</a><a class="docs-footer-nextpage" href="../correlated_inheritance/">correlated inheritance »</a><div class="flexbox-break"></div><p class="footer-message">Powered by <a href="https://github.com/JuliaDocs/Documenter.jl">Documenter.jl</a> and the <a href="https://julialang.org/">Julia Programming Language</a>.</p></nav></div><div class="modal" id="documenter-settings"><div class="modal-background"></div><div class="modal-card"><header class="modal-card-head"><p class="modal-card-title">Settings</p><button class="delete"></button></header><section class="modal-card-body"><p><label class="label">Theme</label><div class="select"><select id="documenter-themepicker"><option value="auto">Automatic (OS)</option><option value="documenter-light">documenter-light</option><option value="documenter-dark">documenter-dark</option></select></div></p><hr/><p>This document was generated with <a href="https://github.com/JuliaDocs/Documenter.jl">Documenter.jl</a> version 1.4.1 on <span class="colophon-date" title="Friday 17 May 2024 12:07">Friday 17 May 2024</span>. Using Julia version 1.10.3.</p></section><footer class="modal-card-foot"></footer></div></div></div></body></html>
diff --git a/previews/PR12/man/correlated_inheritance/index.html b/previews/PR12/man/correlated_inheritance/index.html
index 46c3756..7395324 100644
--- a/previews/PR12/man/correlated_inheritance/index.html
+++ b/previews/PR12/man/correlated_inheritance/index.html
@@ -15,4 +15,4 @@
plot(gt1, shownodelabel=true, edgelabel=el1, edgelabelcolor=el1.label, tipoffset=0.1);
plot(gt2, shownodelabel=true, edgelabel=el2, edgelabelcolor=el2.label, tipoffset=0.1);
-plot(gt3, shownodelabel=true, edgelabel=el3, edgelabelcolor=el3.label, tipoffset=0.1);</code></pre><p><img src="../../assets/figures/net1taxon_3genetrees.svg" alt="3 gene trees on 1-taxon network with inheritance correlation"/></p><p>In all cases, any lineage has a probability γ=0.6 to come from species edge 2 (labeled in red), and probability γ=0.4 to come from species edge 3 (in green). When the inheritance correlation r increases, lineages have an increased preference to come from the same parent as other lineages (at the same locus).</p></article><nav class="docs-footer"><a class="docs-footer-prevpage" href="../converting_coal2generation_units/">« converting between units</a><a class="docs-footer-nextpage" href="../more_examples/">more examples »</a><div class="flexbox-break"></div><p class="footer-message">Powered by <a href="https://github.com/JuliaDocs/Documenter.jl">Documenter.jl</a> and the <a href="https://julialang.org/">Julia Programming Language</a>.</p></nav></div><div class="modal" id="documenter-settings"><div class="modal-background"></div><div class="modal-card"><header class="modal-card-head"><p class="modal-card-title">Settings</p><button class="delete"></button></header><section class="modal-card-body"><p><label class="label">Theme</label><div class="select"><select id="documenter-themepicker"><option value="auto">Automatic (OS)</option><option value="documenter-light">documenter-light</option><option value="documenter-dark">documenter-dark</option></select></div></p><hr/><p>This document was generated with <a href="https://github.com/JuliaDocs/Documenter.jl">Documenter.jl</a> version 1.4.1 on <span class="colophon-date" title="Friday 17 May 2024 03:10">Friday 17 May 2024</span>. Using Julia version 1.10.3.</p></section><footer class="modal-card-foot"></footer></div></div></div></body></html>
+plot(gt3, shownodelabel=true, edgelabel=el3, edgelabelcolor=el3.label, tipoffset=0.1);</code></pre><p><img src="../../assets/figures/net1taxon_3genetrees.svg" alt="3 gene trees on 1-taxon network with inheritance correlation"/></p><p>In all cases, any lineage has a probability γ=0.6 to come from species edge 2 (labeled in red), and probability γ=0.4 to come from species edge 3 (in green). When the inheritance correlation r increases, lineages have an increased preference to come from the same parent as other lineages (at the same locus).</p></article><nav class="docs-footer"><a class="docs-footer-prevpage" href="../converting_coal2generation_units/">« converting between units</a><a class="docs-footer-nextpage" href="../more_examples/">more examples »</a><div class="flexbox-break"></div><p class="footer-message">Powered by <a href="https://github.com/JuliaDocs/Documenter.jl">Documenter.jl</a> and the <a href="https://julialang.org/">Julia Programming Language</a>.</p></nav></div><div class="modal" id="documenter-settings"><div class="modal-background"></div><div class="modal-card"><header class="modal-card-head"><p class="modal-card-title">Settings</p><button class="delete"></button></header><section class="modal-card-body"><p><label class="label">Theme</label><div class="select"><select id="documenter-themepicker"><option value="auto">Automatic (OS)</option><option value="documenter-light">documenter-light</option><option value="documenter-dark">documenter-dark</option></select></div></p><hr/><p>This document was generated with <a href="https://github.com/JuliaDocs/Documenter.jl">Documenter.jl</a> version 1.4.1 on <span class="colophon-date" title="Friday 17 May 2024 12:07">Friday 17 May 2024</span>. Using Julia version 1.10.3.</p></section><footer class="modal-card-foot"></footer></div></div></div></body></html>
diff --git a/previews/PR12/man/getting_started/index.html b/previews/PR12/man/getting_started/index.html
index ba069ee..4f0783c 100644
--- a/previews/PR12/man/getting_started/index.html
+++ b/previews/PR12/man/getting_started/index.html
@@ -31,4 +31,4 @@
tip labels: B_3, C_3, C_2, C_1, ...
(((B_3:0.964,C_3:0.964):0.214,(C_2:0.756,C_1:0.756):0.422):1.125,((B_1:0.595,B_2:0.595):0.849,(A_1:0.079,(A_2:0.054,A_3:0.054):0.025):1.565):0.659);</code></pre><p>We can also ask for varying numbers of individuals. For example, we simulate below 2 individuals in A and 1 individual in each of B and C, using a dictionary to map species to their number of individuals:</p><pre><code class="language-julia-repl hljs" style="display:block;">julia> genetrees = simulatecoalescent(net, 1, Dict("A"=>2, "B"=>1, "C"=>1));</code><code class="nohighlight hljs ansi" style="display:block;"></code><br/><code class="language-julia-repl hljs" style="display:block;">julia> writeTopology(genetrees[1])</code><code class="nohighlight hljs ansi" style="display:block;">"(C:1.7342183562262905,(B:1.1125494690619213,(A_2:0.39810435449985226,A_1:0.39810435449985226):0.9144451145620689):0.42166888716436934);"</code></pre><p>We can set 0 individuals within a species to simulate missing data.</p><pre><code class="language-julia-repl hljs" style="display:block;">julia> genetrees = simulatecoalescent(net, 3, Dict("A"=>2, "B"=>1, "C"=>0));</code><code class="nohighlight hljs ansi" style="display:block;"></code><br/><code class="language-julia-repl hljs" style="display:block;">julia> writeMultiTopology(genetrees, stdout)</code><code class="nohighlight hljs ansi" style="display:block;">((A_1:0.6217550029603106,A_2:0.6217550029603106):1.8934561254937003,B:2.515211128454011);
((A_1:0.020216077041873223,A_2:0.020216077041873223):1.6564835196459733,B:1.6766995966878466);
-(B:1.3598531402719405,(A_2:0.5258264787632165,A_1:0.5258264787632165):1.034026661508724);</code></pre></article><nav class="docs-footer"><a class="docs-footer-prevpage" href="../../">« home</a><a class="docs-footer-nextpage" href="../mapping_genetree_to_network/">mapping gene trees into the species network »</a><div class="flexbox-break"></div><p class="footer-message">Powered by <a href="https://github.com/JuliaDocs/Documenter.jl">Documenter.jl</a> and the <a href="https://julialang.org/">Julia Programming Language</a>.</p></nav></div><div class="modal" id="documenter-settings"><div class="modal-background"></div><div class="modal-card"><header class="modal-card-head"><p class="modal-card-title">Settings</p><button class="delete"></button></header><section class="modal-card-body"><p><label class="label">Theme</label><div class="select"><select id="documenter-themepicker"><option value="auto">Automatic (OS)</option><option value="documenter-light">documenter-light</option><option value="documenter-dark">documenter-dark</option></select></div></p><hr/><p>This document was generated with <a href="https://github.com/JuliaDocs/Documenter.jl">Documenter.jl</a> version 1.4.1 on <span class="colophon-date" title="Friday 17 May 2024 03:10">Friday 17 May 2024</span>. Using Julia version 1.10.3.</p></section><footer class="modal-card-foot"></footer></div></div></div></body></html>
+(B:1.3598531402719405,(A_2:0.5258264787632165,A_1:0.5258264787632165):1.034026661508724);</code></pre></article><nav class="docs-footer"><a class="docs-footer-prevpage" href="../../">« home</a><a class="docs-footer-nextpage" href="../mapping_genetree_to_network/">mapping gene trees into the species network »</a><div class="flexbox-break"></div><p class="footer-message">Powered by <a href="https://github.com/JuliaDocs/Documenter.jl">Documenter.jl</a> and the <a href="https://julialang.org/">Julia Programming Language</a>.</p></nav></div><div class="modal" id="documenter-settings"><div class="modal-background"></div><div class="modal-card"><header class="modal-card-head"><p class="modal-card-title">Settings</p><button class="delete"></button></header><section class="modal-card-body"><p><label class="label">Theme</label><div class="select"><select id="documenter-themepicker"><option value="auto">Automatic (OS)</option><option value="documenter-light">documenter-light</option><option value="documenter-dark">documenter-dark</option></select></div></p><hr/><p>This document was generated with <a href="https://github.com/JuliaDocs/Documenter.jl">Documenter.jl</a> version 1.4.1 on <span class="colophon-date" title="Friday 17 May 2024 12:07">Friday 17 May 2024</span>. Using Julia version 1.10.3.</p></section><footer class="modal-card-foot"></footer></div></div></div></body></html>
diff --git a/previews/PR12/man/mapping_genetree_to_network/index.html b/previews/PR12/man/mapping_genetree_to_network/index.html
index afdd2ac..a476158 100644
--- a/previews/PR12/man/mapping_genetree_to_network/index.html
+++ b/previews/PR12/man/mapping_genetree_to_network/index.html
@@ -6,4 +6,4 @@
4 edges
5 nodes: 3 tips, 0 hybrid nodes, 2 internal tree nodes.
tip labels: B, C, A
-((C:1.695,A:1.695):0.979,B:2.474);</code></pre><p>The option <code>true</code> is to keep the root, even if it's of degree 2.</p></article><nav class="docs-footer"><a class="docs-footer-prevpage" href="../getting_started/">« getting started</a><a class="docs-footer-nextpage" href="../converting_coal2generation_units/">converting between units »</a><div class="flexbox-break"></div><p class="footer-message">Powered by <a href="https://github.com/JuliaDocs/Documenter.jl">Documenter.jl</a> and the <a href="https://julialang.org/">Julia Programming Language</a>.</p></nav></div><div class="modal" id="documenter-settings"><div class="modal-background"></div><div class="modal-card"><header class="modal-card-head"><p class="modal-card-title">Settings</p><button class="delete"></button></header><section class="modal-card-body"><p><label class="label">Theme</label><div class="select"><select id="documenter-themepicker"><option value="auto">Automatic (OS)</option><option value="documenter-light">documenter-light</option><option value="documenter-dark">documenter-dark</option></select></div></p><hr/><p>This document was generated with <a href="https://github.com/JuliaDocs/Documenter.jl">Documenter.jl</a> version 1.4.1 on <span class="colophon-date" title="Friday 17 May 2024 03:10">Friday 17 May 2024</span>. Using Julia version 1.10.3.</p></section><footer class="modal-card-foot"></footer></div></div></div></body></html>
+((C:1.695,A:1.695):0.979,B:2.474);</code></pre><p>The option <code>true</code> is to keep the root, even if it's of degree 2.</p></article><nav class="docs-footer"><a class="docs-footer-prevpage" href="../getting_started/">« getting started</a><a class="docs-footer-nextpage" href="../converting_coal2generation_units/">converting between units »</a><div class="flexbox-break"></div><p class="footer-message">Powered by <a href="https://github.com/JuliaDocs/Documenter.jl">Documenter.jl</a> and the <a href="https://julialang.org/">Julia Programming Language</a>.</p></nav></div><div class="modal" id="documenter-settings"><div class="modal-background"></div><div class="modal-card"><header class="modal-card-head"><p class="modal-card-title">Settings</p><button class="delete"></button></header><section class="modal-card-body"><p><label class="label">Theme</label><div class="select"><select id="documenter-themepicker"><option value="auto">Automatic (OS)</option><option value="documenter-light">documenter-light</option><option value="documenter-dark">documenter-dark</option></select></div></p><hr/><p>This document was generated with <a href="https://github.com/JuliaDocs/Documenter.jl">Documenter.jl</a> version 1.4.1 on <span class="colophon-date" title="Friday 17 May 2024 12:07">Friday 17 May 2024</span>. Using Julia version 1.10.3.</p></section><footer class="modal-card-foot"></footer></div></div></div></body></html>
diff --git a/previews/PR12/man/more_examples/index.html b/previews/PR12/man/more_examples/index.html
index 7bffb3d..1625bf3 100644
--- a/previews/PR12/man/more_examples/index.html
+++ b/previews/PR12/man/more_examples/index.html
@@ -1,5 +1,5 @@
<!DOCTYPE html>
-<html lang="en"><head><meta charset="UTF-8"/><meta name="viewport" content="width=device-width, initial-scale=1.0"/><title>more examples · PhyloCoalSimulations.jl</title><meta name="title" content="more examples · PhyloCoalSimulations.jl"/><meta property="og:title" content="more examples · PhyloCoalSimulations.jl"/><meta property="twitter:title" content="more examples · PhyloCoalSimulations.jl"/><meta name="description" content="Documentation for PhyloCoalSimulations.jl."/><meta property="og:description" content="Documentation for PhyloCoalSimulations.jl."/><meta property="twitter:description" content="Documentation for PhyloCoalSimulations.jl."/><meta property="og:url" content="https://cecileane.github.io/PhyloCoalSimulations.jl/stable/man/more_examples/"/><meta property="twitter:url" content="https://cecileane.github.io/PhyloCoalSimulations.jl/stable/man/more_examples/"/><link rel="canonical" href="https://cecileane.github.io/PhyloCoalSimulations.jl/stable/man/more_examples/"/><script data-outdated-warner src="../../assets/warner.js"></script><link href="https://cdnjs.cloudflare.com/ajax/libs/lato-font/3.0.0/css/lato-font.min.css" rel="stylesheet" type="text/css"/><link href="https://cdnjs.cloudflare.com/ajax/libs/juliamono/0.050/juliamono.min.css" rel="stylesheet" type="text/css"/><link href="https://cdnjs.cloudflare.com/ajax/libs/font-awesome/6.4.2/css/fontawesome.min.css" rel="stylesheet" type="text/css"/><link href="https://cdnjs.cloudflare.com/ajax/libs/font-awesome/6.4.2/css/solid.min.css" rel="stylesheet" type="text/css"/><link href="https://cdnjs.cloudflare.com/ajax/libs/font-awesome/6.4.2/css/brands.min.css" rel="stylesheet" type="text/css"/><link href="https://cdnjs.cloudflare.com/ajax/libs/KaTeX/0.16.8/katex.min.css" rel="stylesheet" type="text/css"/><script>documenterBaseURL="../.."</script><script src="https://cdnjs.cloudflare.com/ajax/libs/require.js/2.3.6/require.min.js" data-main="../../assets/documenter.js"></script><script src="../../search_index.js"></script><script src="../../siteinfo.js"></script><script src="../../../versions.js"></script><link class="docs-theme-link" rel="stylesheet" type="text/css" href="../../assets/themes/documenter-dark.css" data-theme-name="documenter-dark" data-theme-primary-dark/><link class="docs-theme-link" rel="stylesheet" type="text/css" href="../../assets/themes/documenter-light.css" data-theme-name="documenter-light" data-theme-primary/><script src="../../assets/themeswap.js"></script></head><body><div id="documenter"><nav class="docs-sidebar"><div class="docs-package-name"><span class="docs-autofit"><a href="../../">PhyloCoalSimulations.jl</a></span></div><button class="docs-search-query input is-rounded is-small is-clickable my-2 mx-auto py-1 px-2" id="documenter-search-query">Search docs (Ctrl + /)</button><ul class="docs-menu"><li><a class="tocitem" href="../../">home</a></li><li><span class="tocitem">manual</span><ul><li><a class="tocitem" href="../getting_started/">getting started</a></li><li><a class="tocitem" href="../mapping_genetree_to_network/">mapping gene trees into the species network</a></li><li><a class="tocitem" href="../converting_coal2generation_units/">converting between units</a></li><li><a class="tocitem" href="../correlated_inheritance/">correlated inheritance</a></li><li class="is-active"><a class="tocitem" href>more examples</a><ul class="internal"><li><a class="tocitem" href="#counting-deep-coalescences"><span>counting deep coalescences</span></a></li><li><a class="tocitem" href="#number-of-lineages-inherited-via-gene-flow"><span>number of lineages inherited via gene flow</span></a></li><li><a class="tocitem" href="#rate-variation-across-species"><span>rate variation across species</span></a></li></ul></li></ul></li><li><span class="tocitem">library</span><ul><li><a class="tocitem" href="../../lib/public/">public</a></li><li><a class="tocitem" href="../../lib/internal/">internals</a></li></ul></li></ul><div class="docs-version-selector field has-addons"><div class="control"><span class="docs-label button is-static is-size-7">Version</span></div><div class="docs-selector control is-expanded"><div class="select is-fullwidth is-size-7"><select id="documenter-version-selector"></select></div></div></div></nav><div class="docs-main"><header class="docs-navbar"><a class="docs-sidebar-button docs-navbar-link fa-solid fa-bars is-hidden-desktop" id="documenter-sidebar-button" href="#"></a><nav class="breadcrumb"><ul class="is-hidden-mobile"><li><a class="is-disabled">manual</a></li><li class="is-active"><a href>more examples</a></li></ul><ul class="is-hidden-tablet"><li class="is-active"><a href>more examples</a></li></ul></nav><div class="docs-right"><a class="docs-navbar-link" href="https://github.com/cecileane/PhyloCoalSimulations.jl" title="View the repository on GitHub"><span class="docs-icon fa-brands"></span><span class="docs-label is-hidden-touch">GitHub</span></a><a class="docs-navbar-link" href="https://github.com/cecileane/PhyloCoalSimulations.jl/blob/main/docs/src/man/more_examples.md" title="Edit source on GitHub"><span class="docs-icon fa-solid"></span></a><a class="docs-settings-button docs-navbar-link fa-solid fa-gear" id="documenter-settings-button" href="#" title="Settings"></a><a class="docs-article-toggle-button fa-solid fa-chevron-up" id="documenter-article-toggle-button" href="javascript:;" title="Collapse all docstrings"></a></div></header><article class="content" id="documenter-page"><h1 id="example-uses"><a class="docs-heading-anchor" href="#example-uses">example uses</a><a id="example-uses-1"></a><a class="docs-heading-anchor-permalink" href="#example-uses" title="Permalink"></a></h1><p>We are re-using the same network and simulated tree as before:</p><pre><code class="language-julia-repl hljs" style="display:block;">julia> writeTopology(net)</code><code class="nohighlight hljs ansi" style="display:block;">"((C:0.9,(B:0.2)#H1:0.7::0.6)I1:0.6,(#H1:0.6::0.4,A:1.0)I2:0.5)I3;"</code><br/><code class="language-julia-repl hljs" style="display:block;">julia> writeTopology(tree, round=true)</code><code class="nohighlight hljs ansi" style="display:block;">"((((C:0.9)I1:0.6)I3:0.64,((A:1.0)I2:0.5)I3:0.64):1.424,(((B:0.2)H1:0.6)I2:0.5)I3:2.064);"</code></pre><p><img src="../../assets/figures/genetree_example1.svg" alt="example 1, same as in mapping section"/></p><h2 id="counting-deep-coalescences"><a class="docs-heading-anchor" href="#counting-deep-coalescences">counting deep coalescences</a><a id="counting-deep-coalescences-1"></a><a class="docs-heading-anchor-permalink" href="#counting-deep-coalescences" title="Permalink"></a></h2><p>The number of deep coalescences can be quantified as the number of "extra" lineages due to incomplete lineage sorting, that can be calculated from embedding the gene tree into the species phylogeny (see <a href="https://doi.org/10.1093/sysbio/46.3.523">Maddison 1997</a> for species trees). For an edge in the network, say edge 7 going from I2 to I3 going in back in time, lineage sorting is complete if all the gene lineages entering the edge (at I2) coalesce into a single gene lineage by the time they exit the edge (at I3). If they don't, the number of extra lineages is <code>k-1</code> where <code>k</code> is the number of lineages "exiting" the edge, for that particular edge in the species network and that particular gene tree. The total number of extra lineages, for a given gene tree, is the sum across all edges in the species phylogeny.</p><p>In our gene tree above, we can count the number of lineages that exit each species edge using the degree-2 mapping nodes, then count how many lineages are "extra". We do so below using utilities <a href="../../lib/internal/#PhyloCoalSimulations.mappingnodes"><code>mappingnodes</code></a> to iterate over degree-2 mapping nodes and <a href="../../lib/internal/#PhyloCoalSimulations.population_mappedto-Tuple{Union{PhyloNetworks.EdgeT{PhyloNetworks.Node}, PhyloNetworks.Node}}"><code>population_mappedto</code></a> to extract the mapping information.</p><pre><code class="language-julia-repl hljs" style="display:block;">julia> # dictionary to store the count of extra lineages exiting each network edge. initialized at 0
+<html lang="en"><head><meta charset="UTF-8"/><meta name="viewport" content="width=device-width, initial-scale=1.0"/><title>more examples · PhyloCoalSimulations.jl</title><meta name="title" content="more examples · PhyloCoalSimulations.jl"/><meta property="og:title" content="more examples · PhyloCoalSimulations.jl"/><meta property="twitter:title" content="more examples · PhyloCoalSimulations.jl"/><meta name="description" content="Documentation for PhyloCoalSimulations.jl."/><meta property="og:description" content="Documentation for PhyloCoalSimulations.jl."/><meta property="twitter:description" content="Documentation for PhyloCoalSimulations.jl."/><meta property="og:url" content="https://cecileane.github.io/PhyloCoalSimulations.jl/stable/man/more_examples/"/><meta property="twitter:url" content="https://cecileane.github.io/PhyloCoalSimulations.jl/stable/man/more_examples/"/><link rel="canonical" href="https://cecileane.github.io/PhyloCoalSimulations.jl/stable/man/more_examples/"/><script data-outdated-warner src="../../assets/warner.js"></script><link href="https://cdnjs.cloudflare.com/ajax/libs/lato-font/3.0.0/css/lato-font.min.css" rel="stylesheet" type="text/css"/><link href="https://cdnjs.cloudflare.com/ajax/libs/juliamono/0.050/juliamono.min.css" rel="stylesheet" type="text/css"/><link href="https://cdnjs.cloudflare.com/ajax/libs/font-awesome/6.4.2/css/fontawesome.min.css" rel="stylesheet" type="text/css"/><link href="https://cdnjs.cloudflare.com/ajax/libs/font-awesome/6.4.2/css/solid.min.css" rel="stylesheet" type="text/css"/><link href="https://cdnjs.cloudflare.com/ajax/libs/font-awesome/6.4.2/css/brands.min.css" rel="stylesheet" type="text/css"/><link href="https://cdnjs.cloudflare.com/ajax/libs/KaTeX/0.16.8/katex.min.css" rel="stylesheet" type="text/css"/><script>documenterBaseURL="../.."</script><script src="https://cdnjs.cloudflare.com/ajax/libs/require.js/2.3.6/require.min.js" data-main="../../assets/documenter.js"></script><script src="../../search_index.js"></script><script src="../../siteinfo.js"></script><script src="../../../versions.js"></script><link class="docs-theme-link" rel="stylesheet" type="text/css" href="../../assets/themes/documenter-dark.css" data-theme-name="documenter-dark" data-theme-primary-dark/><link class="docs-theme-link" rel="stylesheet" type="text/css" href="../../assets/themes/documenter-light.css" data-theme-name="documenter-light" data-theme-primary/><script src="../../assets/themeswap.js"></script></head><body><div id="documenter"><nav class="docs-sidebar"><div class="docs-package-name"><span class="docs-autofit"><a href="../../">PhyloCoalSimulations.jl</a></span></div><button class="docs-search-query input is-rounded is-small is-clickable my-2 mx-auto py-1 px-2" id="documenter-search-query">Search docs (Ctrl + /)</button><ul class="docs-menu"><li><a class="tocitem" href="../../">home</a></li><li><span class="tocitem">manual</span><ul><li><a class="tocitem" href="../getting_started/">getting started</a></li><li><a class="tocitem" href="../mapping_genetree_to_network/">mapping gene trees into the species network</a></li><li><a class="tocitem" href="../converting_coal2generation_units/">converting between units</a></li><li><a class="tocitem" href="../correlated_inheritance/">correlated inheritance</a></li><li class="is-active"><a class="tocitem" href>more examples</a><ul class="internal"><li><a class="tocitem" href="#counting-deep-coalescences"><span>counting deep coalescences</span></a></li><li><a class="tocitem" href="#number-of-lineages-inherited-via-gene-flow"><span>number of lineages inherited via gene flow</span></a></li><li><a class="tocitem" href="#rate-variation-across-species"><span>rate variation across species</span></a></li></ul></li></ul></li><li><span class="tocitem">library</span><ul><li><a class="tocitem" href="../../lib/public/">public</a></li><li><a class="tocitem" href="../../lib/internal/">internals</a></li></ul></li></ul><div class="docs-version-selector field has-addons"><div class="control"><span class="docs-label button is-static is-size-7">Version</span></div><div class="docs-selector control is-expanded"><div class="select is-fullwidth is-size-7"><select id="documenter-version-selector"></select></div></div></div></nav><div class="docs-main"><header class="docs-navbar"><a class="docs-sidebar-button docs-navbar-link fa-solid fa-bars is-hidden-desktop" id="documenter-sidebar-button" href="#"></a><nav class="breadcrumb"><ul class="is-hidden-mobile"><li><a class="is-disabled">manual</a></li><li class="is-active"><a href>more examples</a></li></ul><ul class="is-hidden-tablet"><li class="is-active"><a href>more examples</a></li></ul></nav><div class="docs-right"><a class="docs-navbar-link" href="https://github.com/cecileane/PhyloCoalSimulations.jl" title="View the repository on GitHub"><span class="docs-icon fa-brands"></span><span class="docs-label is-hidden-touch">GitHub</span></a><a class="docs-navbar-link" href="https://github.com/cecileane/PhyloCoalSimulations.jl/blob/main/docs/src/man/more_examples.md" title="Edit source on GitHub"><span class="docs-icon fa-solid"></span></a><a class="docs-settings-button docs-navbar-link fa-solid fa-gear" id="documenter-settings-button" href="#" title="Settings"></a><a class="docs-article-toggle-button fa-solid fa-chevron-up" id="documenter-article-toggle-button" href="javascript:;" title="Collapse all docstrings"></a></div></header><article class="content" id="documenter-page"><h1 id="example-uses"><a class="docs-heading-anchor" href="#example-uses">example uses</a><a id="example-uses-1"></a><a class="docs-heading-anchor-permalink" href="#example-uses" title="Permalink"></a></h1><p>We are re-using the same network and simulated tree as before:</p><pre><code class="language-julia-repl hljs" style="display:block;">julia> writeTopology(net)</code><code class="nohighlight hljs ansi" style="display:block;">"((C:0.9,(B:0.2)#H1:0.7::0.6)i1:0.6,(#H1:0.6::0.4,A:1.0)i2:0.5)i3;"</code><br/><code class="language-julia-repl hljs" style="display:block;">julia> writeTopology(tree, round=true)</code><code class="nohighlight hljs ansi" style="display:block;">"((((B:0.2)H1:0.6)i2:0.5)i3:1.174,(((C:0.9)i1:0.6)i3:0.195,((A:1.0)i2:0.5)i3:0.195):0.979);"</code></pre><p><img src="../../assets/figures/genetree_example1.svg" alt="example 1, same as in mapping section"/></p><h2 id="counting-deep-coalescences"><a class="docs-heading-anchor" href="#counting-deep-coalescences">counting deep coalescences</a><a id="counting-deep-coalescences-1"></a><a class="docs-heading-anchor-permalink" href="#counting-deep-coalescences" title="Permalink"></a></h2><p>The number of deep coalescences can be quantified as the number of "extra" lineages due to incomplete lineage sorting, that can be calculated from embedding the gene tree into the species phylogeny (see <a href="https://doi.org/10.1093/sysbio/46.3.523">Maddison 1997</a> for species trees). For an edge in the network, say edge 7 going from i2 to i3 going in back in time, lineage sorting is complete if all the gene lineages entering the edge (at i2) coalesce into a single gene lineage by the time they exit the edge (at i3). If they don't, the number of extra lineages is <code>k-1</code> where <code>k</code> is the number of lineages "exiting" the edge, for that particular edge in the species network and that particular gene tree. The total number of extra lineages, for a given gene tree, is the sum across all edges in the species phylogeny.</p><p>In our gene tree above, we can count the number of lineages that exit each species edge using the degree-2 mapping nodes, then count how many lineages are "extra". We do so below using utilities <a href="../../lib/internal/#PhyloCoalSimulations.mappingnodes"><code>mappingnodes</code></a> to iterate over degree-2 mapping nodes and <a href="../../lib/internal/#PhyloCoalSimulations.population_mappedto-Tuple{Union{PhyloNetworks.EdgeT{PhyloNetworks.Node}, PhyloNetworks.Node}}"><code>population_mappedto</code></a> to extract the mapping information.</p><pre><code class="language-julia-repl hljs" style="display:block;">julia> # dictionary to store the count of extra lineages exiting each network edge. initialized at 0
edge_count = Dict(e.number => 0 for e in net.edge)</code><code class="nohighlight hljs ansi" style="display:block;">Dict{Int64, Int64} with 7 entries:
5 => 0
4 => 0
@@ -21,7 +21,7 @@
7 => 2
2 => 1
3 => 0
- 1 => 1</code></pre><p>From this, we see two interesting things.</p><ul><li>0 lineages exited edge number 3 in the species network: it's the hybrid edge from H1 to I3 (going back in time). That's because the only lineage at H1 was interited from I2, so there weren't any lineage evolving through edge 3.</li><li>2 lineages exited edge number 7 (going from I2 to I3 back in time), so that's 1 extra lineage. All other edges look as expected, with a single gene lineage exiting from them.</li></ul><p>We can now calculate the total number of extra lineages:</p><pre><code class="language-julia-repl hljs" style="display:block;">julia> filter!(p -> p.second > 0, edge_count) # filter out edges without any gene lineage</code><code class="nohighlight hljs ansi" style="display:block;">Dict{Int64, Int64} with 6 entries:
+ 1 => 1</code></pre><p>From this, we see two interesting things.</p><ul><li>0 lineages exited edge number 3 in the species network: it's the hybrid edge from H1 to i3 (going back in time). That's because the only lineage at H1 was interited from i2, so there weren't any lineage evolving through edge 3.</li><li>2 lineages exited edge number 7 (going from i2 to i3 back in time), so that's 1 extra lineage. All other edges look as expected, with a single gene lineage exiting from them.</li></ul><p>We can now calculate the total number of extra lineages:</p><pre><code class="language-julia-repl hljs" style="display:block;">julia> filter!(p -> p.second > 0, edge_count) # filter out edges without any gene lineage</code><code class="nohighlight hljs ansi" style="display:block;">Dict{Int64, Int64} with 6 entries:
5 => 1
4 => 1
6 => 1
@@ -33,25 +33,22 @@
0
1
0
- 0</code><br/><code class="language-julia-repl hljs" style="display:block;">julia> deepcoalescence = sum(values(edge_count)) # sum 'extras' over all edges in the network</code><code class="nohighlight hljs ansi" style="display:block;">1</code></pre><p>On the particular gene tree we simulated, we counted 1 deep coalescence.</p><h2 id="number-of-lineages-inherited-via-gene-flow"><a class="docs-heading-anchor" href="#number-of-lineages-inherited-via-gene-flow">number of lineages inherited via gene flow</a><a id="number-of-lineages-inherited-via-gene-flow-1"></a><a class="docs-heading-anchor-permalink" href="#number-of-lineages-inherited-via-gene-flow" title="Permalink"></a></h2><p>Our network has inheritance γ=0.4 on the minor edge, which we'll call the "gene flow" edge, and γ=0.6 on the major hybrid edge, parent to H1 on the major tree. But we may be interested in the realized proportion of lineages inherited from each parent at H1, realized in the gene trees we actually simulated. To do so, we can count the number of gene lineages that are mapped to each hybrid edge in the network.</p><p>This mapping is stored in the edge attribute <code>.inCycle</code> internally, but it's best to access it via the function <a href="../../lib/internal/#PhyloCoalSimulations.population_mappedto-Tuple{Union{PhyloNetworks.EdgeT{PhyloNetworks.Node}, PhyloNetworks.Node}}"><code>population_mappedto</code></a> (as the internal representation may change). From the plot above, the minor "gene flow" edge is edge number 5 and the major hybrid edge has number 3. So we can count the gene lineages inherited via gene flow as the number of gene tree edges with <code>inCycle</code> equal to 5.</p><p>If the gene trees have been saved to a file and later read from this file, then the <code>.inCycle</code> attributes are no longer stored in memory. In this case, we can retrieve the mapping information by the internal node names. The edges going through gene flow are those whose child node is named "H1" and parent node is named "I2".</p><p>We use the first option with the <code>.inCycle</code> attribute below. We get that our one simulated gene tree was indeed inherited via gene flow:</p><pre><code class="language-julia-repl hljs" style="display:block;">julia> sum(e.inCycle == 5 for e in tree.edge) # or:</code><code class="nohighlight hljs ansi" style="display:block;">1</code><br/><code class="language-julia-repl hljs" style="display:block;">julia> sum(PCS.population_mappedto(e) == 5 for e in tree.edge)
- # or define a function to do this for any edge, so we can re-use later:</code><code class="nohighlight hljs ansi" style="display:block;">1</code><br/><code class="language-julia-repl hljs" style="display:block;">julia> nlineages_through(edgeID, gt) = sum(PCS.population_mappedto(e) == edgeID for e in gt.edge);</code><code class="nohighlight hljs ansi" style="display:block;"></code><br/><code class="language-julia-repl hljs" style="display:block;">julia> nlineages_through(5, tree) # same as before!</code><code class="nohighlight hljs ansi" style="display:block;">1</code><br/><code class="language-julia-repl hljs" style="display:block;">julia> nlineages_through(3, tree) # lineages that went through edge 3, the major edge.</code><code class="nohighlight hljs ansi" style="display:block;">0</code></pre><p>To make this more interesting, we can simulate many gene trees then count how many of their lineages were inherited via gene flow. If we ask for 2 individuals in species B, then each gene may have 2 lineages that enter the hybrid node H1, if the two B individuals fail to coalesce. In that case, it's possible that one individual lineage was inherited via gene flow, and the other not. We'll calculate the gene flow proportion among all these lineages. This proportion should be close (but not exactly equal) to the theoretical γ=0.4 from the network model.</p><pre><code class="language-julia-repl hljs" style="display:block;">julia> ngenes = 100;</code><code class="nohighlight hljs ansi" style="display:block;"></code><br/><code class="language-julia-repl hljs" style="display:block;">julia> genetrees = simulatecoalescent(net, ngenes, Dict("B"=>2, "A"=>1, "C"=>1); nodemapping=true);</code><code class="nohighlight hljs ansi" style="display:block;"></code><br/><code class="language-julia-repl hljs" style="display:block;">julia> length(genetrees)</code><code class="nohighlight hljs ansi" style="display:block;">100</code><br/><code class="language-julia-repl hljs" style="display:block;">julia> nlineages_geneflow = sum(nlineages_through(5,gt) for gt in genetrees)</code><code class="nohighlight hljs ansi" style="display:block;">70</code><br/><code class="language-julia-repl hljs" style="display:block;">julia> nlineages_major = sum(nlineages_through(3,gt) for gt in genetrees)
- # realized γ, close to 0.4:</code><code class="nohighlight hljs ansi" style="display:block;">129</code><br/><code class="language-julia-repl hljs" style="display:block;">julia> proportion_geneflow = nlineages_geneflow / (nlineages_geneflow + nlineages_major)</code><code class="nohighlight hljs ansi" style="display:block;">0.35175879396984927</code></pre><h2 id="rate-variation-across-species"><a class="docs-heading-anchor" href="#rate-variation-across-species">rate variation across species</a><a id="rate-variation-across-species-1"></a><a class="docs-heading-anchor-permalink" href="#rate-variation-across-species" title="Permalink"></a></h2><p>The gene trees resulting from <code>simulatecoalescent</code> have their edge lengths in coalescent units. One may want to convert them to substitutions per site, so as to simulate molecular sequences along these gene trees. The mapping information is important to allow for different rates of molecular evolution across different species. Here is an example to do this.</p><p>We will use the <code>Distributions</code> package to simulate rates from a log-normal distribution across species, that is, across edges in the species network.</p><pre><code class="language-julia-repl hljs" style="display:block;">julia> using Distributions</code><code class="nohighlight hljs ansi" style="display:block;"></code><br/><code class="language-julia-repl hljs" style="display:block;">julia> lognormal_rate_dist = LogNormal(-0.125, 0.5) # μ = -σ²/2 to get a mean of 1.</code><code class="nohighlight hljs ansi" style="display:block;">Distributions.LogNormal{Float64}(μ=-0.125, σ=0.5)</code><br/><code class="language-julia-repl hljs" style="display:block;">julia> networkedge_rate = Dict(e.number => rand(lognormal_rate_dist) for e in net.edge)
- # add entry for the edge above the network's root. Find its number first.</code><code class="nohighlight hljs ansi" style="display:block;">Dict{Int64, Float64} with 7 entries:
- 5 => 0.68631
- 4 => 1.23458
- 6 => 1.84251
- 7 => 0.438
- 2 => 1.43724
- 3 => 1.09321
- 1 => 0.580741</code><br/><code class="language-julia-repl hljs" style="display:block;">julia> rootedgenumber = PhyloCoalSimulations.get_rootedgenumber(net)</code><code class="nohighlight hljs ansi" style="display:block;">8</code><br/><code class="language-julia-repl hljs" style="display:block;">julia> push!(networkedge_rate, rootedgenumber => rand(lognormal_rate_dist))</code><code class="nohighlight hljs ansi" style="display:block;">Dict{Int64, Float64} with 8 entries:
- 5 => 0.68631
- 4 => 1.23458
- 6 => 1.84251
- 7 => 0.438
- 2 => 1.43724
- 8 => 1.39978
- 3 => 1.09321
- 1 => 0.580741</code><br/><code class="language-julia-repl hljs" style="display:block;">julia> writeTopology(tree, round=true, digits=4) # before rate variation
- # multiply the length of each gene lineage by the rate of the species edge it maps into</code><code class="nohighlight hljs ansi" style="display:block;">"((((C:0.9)I1:0.6)I3:0.64,((A:1.0)I2:0.5)I3:0.64):1.4235,(((B:0.2)H1:0.6)I2:0.5)I3:2.0636);"</code><br/><code class="language-julia-repl hljs" style="display:block;">julia> for e in tree.edge
+ 0</code><br/><code class="language-julia-repl hljs" style="display:block;">julia> deepcoalescence = sum(values(edge_count)) # sum 'extras' over all edges in the network</code><code class="nohighlight hljs ansi" style="display:block;">1</code></pre><p>On the particular gene tree we simulated, we counted 1 deep coalescence.</p><h2 id="number-of-lineages-inherited-via-gene-flow"><a class="docs-heading-anchor" href="#number-of-lineages-inherited-via-gene-flow">number of lineages inherited via gene flow</a><a id="number-of-lineages-inherited-via-gene-flow-1"></a><a class="docs-heading-anchor-permalink" href="#number-of-lineages-inherited-via-gene-flow" title="Permalink"></a></h2><p>Our network has inheritance γ=0.4 on the minor edge, which we'll call the "gene flow" edge, and γ=0.6 on the major hybrid edge, parent to H1 on the major tree. But we may be interested in the realized proportion of lineages inherited from each parent at H1, realized in the gene trees we actually simulated. To do so, we can count the number of gene lineages that are mapped to each hybrid edge in the network.</p><p>This mapping is stored in the edge attribute <code>.inCycle</code> internally, but it's best to access it via the function <a href="../../lib/internal/#PhyloCoalSimulations.population_mappedto-Tuple{Union{PhyloNetworks.EdgeT{PhyloNetworks.Node}, PhyloNetworks.Node}}"><code>population_mappedto</code></a> (as the internal representation may change). From the plot above, the minor "gene flow" edge is edge number 5 and the major hybrid edge has number 3. So we can count the gene lineages inherited via gene flow as the number of gene tree edges with <code>inCycle</code> equal to 5.</p><p>If the gene trees have been saved to a file and later read from this file, then the <code>.inCycle</code> attributes are no longer stored in memory. In this case, we can retrieve the mapping information by the internal node names. The edges going through gene flow are those whose child node is named "H1" and parent node is named "i2".</p><p>We use the first option with the <code>.inCycle</code> attribute below. We get that our one simulated gene tree was indeed inherited via gene flow:</p><pre><code class="language-julia-repl hljs" style="display:block;">julia> sum(e.inCycle == 5 for e in tree.edge) # or:</code><code class="nohighlight hljs ansi" style="display:block;">1</code><br/><code class="language-julia-repl hljs" style="display:block;">julia> sum(PCS.population_mappedto(e) == 5 for e in tree.edge)</code><code class="nohighlight hljs ansi" style="display:block;">1</code></pre><p>Next we define a function to do this for any edge, so we can re-use later:</p><pre><code class="language-julia-repl hljs" style="display:block;">julia> nlineages_through(edgeID, gt) = sum(PCS.population_mappedto(e) == edgeID for e in gt.edge);</code><code class="nohighlight hljs ansi" style="display:block;"></code><br/><code class="language-julia-repl hljs" style="display:block;">julia> nlineages_through(5, tree) # same as before: now done via our new function</code><code class="nohighlight hljs ansi" style="display:block;">1</code><br/><code class="language-julia-repl hljs" style="display:block;">julia> nlineages_through(3, tree) # lineages that went through edge 3, the major edge.</code><code class="nohighlight hljs ansi" style="display:block;">0</code></pre><p>To make this more interesting, we can simulate many gene trees then count how many of their lineages were inherited via gene flow. If we ask for 2 individuals in species B, then each gene may have 2 lineages that enter the hybrid node H1, if the two B individuals fail to coalesce. In that case, it's possible that one individual lineage was inherited via gene flow, and the other not. We'll calculate the gene flow proportion among all these lineages. This proportion should be close (but not exactly equal) to the theoretical γ=0.4 from the network model.</p><pre><code class="language-julia-repl hljs" style="display:block;">julia> ngenes = 100;</code><code class="nohighlight hljs ansi" style="display:block;"></code><br/><code class="language-julia-repl hljs" style="display:block;">julia> genetrees = simulatecoalescent(net, ngenes, Dict("B"=>2, "A"=>1, "C"=>1); nodemapping=true);</code><code class="nohighlight hljs ansi" style="display:block;"></code><br/><code class="language-julia-repl hljs" style="display:block;">julia> length(genetrees)</code><code class="nohighlight hljs ansi" style="display:block;">100</code><br/><code class="language-julia-repl hljs" style="display:block;">julia> nlineages_geneflow = sum(nlineages_through(5,gt) for gt in genetrees)</code><code class="nohighlight hljs ansi" style="display:block;">72</code><br/><code class="language-julia-repl hljs" style="display:block;">julia> nlineages_major = sum(nlineages_through(3,gt) for gt in genetrees)</code><code class="nohighlight hljs ansi" style="display:block;">128</code><br/><code class="language-julia-repl hljs" style="display:block;">julia> proportion_geneflow = nlineages_geneflow / (nlineages_geneflow + nlineages_major)
+ # realized γ, close to 0.4</code><code class="nohighlight hljs ansi" style="display:block;">0.36</code></pre><h2 id="rate-variation-across-species"><a class="docs-heading-anchor" href="#rate-variation-across-species">rate variation across species</a><a id="rate-variation-across-species-1"></a><a class="docs-heading-anchor-permalink" href="#rate-variation-across-species" title="Permalink"></a></h2><p>The gene trees resulting from <code>simulatecoalescent</code> have their edge lengths in coalescent units. One may want to convert them to substitutions per site, so as to simulate molecular sequences along these gene trees. The mapping information is important to allow for different rates of molecular evolution across different species. Here is an example to do this.</p><p>We will use the <code>Distributions</code> package to simulate rates from a log-normal distribution across species, that is, across edges in the species network.</p><pre><code class="language-julia-repl hljs" style="display:block;">julia> using Distributions</code><code class="nohighlight hljs ansi" style="display:block;"></code><br/><code class="language-julia-repl hljs" style="display:block;">julia> lognormal_rate_dist = LogNormal(-0.125, 0.5) # μ = -σ²/2 to get a mean of 1.</code><code class="nohighlight hljs ansi" style="display:block;">Distributions.LogNormal{Float64}(μ=-0.125, σ=0.5)</code><br/><code class="language-julia-repl hljs" style="display:block;">julia> networkedge_rate = Dict(e.number => rand(lognormal_rate_dist) for e in net.edge)</code><code class="nohighlight hljs ansi" style="display:block;">Dict{Int64, Float64} with 7 entries:
+ 5 => 1.64446
+ 4 => 1.55236
+ 6 => 0.351044
+ 7 => 0.551339
+ 2 => 0.672204
+ 3 => 1.75424
+ 1 => 0.548284</code></pre><p>Next we want to simulate a rate for the ancestral edge above the network's root. We find its number first, then add a entry to our dictionary of rates</p><pre><code class="language-julia-repl hljs" style="display:block;">julia> rootedgenumber = PhyloCoalSimulations.get_rootedgenumber(net)</code><code class="nohighlight hljs ansi" style="display:block;">8</code><br/><code class="language-julia-repl hljs" style="display:block;">julia> push!(networkedge_rate, rootedgenumber => rand(lognormal_rate_dist))</code><code class="nohighlight hljs ansi" style="display:block;">Dict{Int64, Float64} with 8 entries:
+ 5 => 1.64446
+ 4 => 1.55236
+ 6 => 0.351044
+ 7 => 0.551339
+ 2 => 0.672204
+ 8 => 1.11207
+ 3 => 1.75424
+ 1 => 0.548284</code><br/><code class="language-julia-repl hljs" style="display:block;">julia> writeTopology(tree, round=true, digits=4) # before rate variation</code><code class="nohighlight hljs ansi" style="display:block;">"((((B:0.2)H1:0.6)i2:0.5)i3:1.1744,(((C:0.9)i1:0.6)i3:0.195,((A:1.0)i2:0.5)i3:0.195):0.9793);"</code></pre><p>Finally, we multiply the length of each gene lineage by the rate of the species edge it maps into:</p><pre><code class="language-julia-repl hljs" style="display:block;">julia> for e in tree.edge
e.length *= networkedge_rate[e.inCycle]
- end</code><code class="nohighlight hljs ansi" style="display:block;"></code><br/><code class="language-julia-repl hljs" style="display:block;">julia> writeTopology(tree, round=true, digits=4) # after rate variation</code><code class="nohighlight hljs ansi" style="display:block;">"((((C:0.5227)I1:0.7407)I3:0.8959,((A:1.8425)I2:0.219)I3:0.8959):1.9926,(((B:0.2874)H1:0.4118)I2:0.219)I3:2.8885);"</code></pre></article><nav class="docs-footer"><a class="docs-footer-prevpage" href="../correlated_inheritance/">« correlated inheritance</a><a class="docs-footer-nextpage" href="../../lib/public/">public »</a><div class="flexbox-break"></div><p class="footer-message">Powered by <a href="https://github.com/JuliaDocs/Documenter.jl">Documenter.jl</a> and the <a href="https://julialang.org/">Julia Programming Language</a>.</p></nav></div><div class="modal" id="documenter-settings"><div class="modal-background"></div><div class="modal-card"><header class="modal-card-head"><p class="modal-card-title">Settings</p><button class="delete"></button></header><section class="modal-card-body"><p><label class="label">Theme</label><div class="select"><select id="documenter-themepicker"><option value="auto">Automatic (OS)</option><option value="documenter-light">documenter-light</option><option value="documenter-dark">documenter-dark</option></select></div></p><hr/><p>This document was generated with <a href="https://github.com/JuliaDocs/Documenter.jl">Documenter.jl</a> version 1.4.1 on <span class="colophon-date" title="Friday 17 May 2024 03:10">Friday 17 May 2024</span>. Using Julia version 1.10.3.</p></section><footer class="modal-card-foot"></footer></div></div></div></body></html>
+ end</code><code class="nohighlight hljs ansi" style="display:block;"></code><br/><code class="language-julia-repl hljs" style="display:block;">julia> writeTopology(tree, round=true, digits=4) # after rate variation</code><code class="nohighlight hljs ansi" style="display:block;">"((((B:0.1344)H1:0.9867)i2:0.2757)i3:1.306,(((C:0.4935)i1:0.9314)i3:0.2169,((A:0.351)i2:0.2757)i3:0.2169):1.0891);"</code></pre></article><nav class="docs-footer"><a class="docs-footer-prevpage" href="../correlated_inheritance/">« correlated inheritance</a><a class="docs-footer-nextpage" href="../../lib/public/">public »</a><div class="flexbox-break"></div><p class="footer-message">Powered by <a href="https://github.com/JuliaDocs/Documenter.jl">Documenter.jl</a> and the <a href="https://julialang.org/">Julia Programming Language</a>.</p></nav></div><div class="modal" id="documenter-settings"><div class="modal-background"></div><div class="modal-card"><header class="modal-card-head"><p class="modal-card-title">Settings</p><button class="delete"></button></header><section class="modal-card-body"><p><label class="label">Theme</label><div class="select"><select id="documenter-themepicker"><option value="auto">Automatic (OS)</option><option value="documenter-light">documenter-light</option><option value="documenter-dark">documenter-dark</option></select></div></p><hr/><p>This document was generated with <a href="https://github.com/JuliaDocs/Documenter.jl">Documenter.jl</a> version 1.4.1 on <span class="colophon-date" title="Friday 17 May 2024 12:07">Friday 17 May 2024</span>. Using Julia version 1.10.3.</p></section><footer class="modal-card-foot"></footer></div></div></div></body></html>
diff --git a/previews/PR12/search_index.js b/previews/PR12/search_index.js
index f0f5877..3a83b3a 100644
--- a/previews/PR12/search_index.js
+++ b/previews/PR12/search_index.js
@@ -1,3 +1,3 @@
var documenterSearchIndex = {"docs":
-[{"location":"man/mapping_genetree_to_network/","page":"mapping gene trees into the species network","title":"mapping gene trees into the species network","text":"using PhyloNetworks, PhyloPlots, PhyloCoalSimulations, RCall, DataFrames, StableRNGs\nfigpath = joinpath(\"..\", \"assets\", \"figures\"); mkpath(figpath)\nfigname(x) = joinpath(figpath, x)","category":"page"},{"location":"man/mapping_genetree_to_network/#mapping-gene-trees-into-the-species-phylogeny","page":"mapping gene trees into the species network","title":"mapping gene trees into the species phylogeny","text":"","category":"section"},{"location":"man/mapping_genetree_to_network/","page":"mapping gene trees into the species network","title":"mapping gene trees into the species network","text":"Nodes in gene trees can be mapped to a node or an edge in the species phylogeny. Edges in gene trees can be mapped to an edge in the species phylogeny. Attributes of nodes and edges are used to carry this mapping information, as detailed in the documentation of function simulatecoalescent.","category":"page"},{"location":"man/mapping_genetree_to_network/","page":"mapping gene trees into the species network","title":"mapping gene trees into the species network","text":"We give examples below of how we may use this mapping information.","category":"page"},{"location":"man/mapping_genetree_to_network/#naming-internal-nodes","page":"mapping gene trees into the species network","title":"naming internal nodes","text":"","category":"section"},{"location":"man/mapping_genetree_to_network/","page":"mapping gene trees into the species network","title":"mapping gene trees into the species network","text":"First, it's useful to name internal nodes in the network, to which we can later map nodes in the gene tree.","category":"page"},{"location":"man/mapping_genetree_to_network/","page":"mapping gene trees into the species network","title":"mapping gene trees into the species network","text":"net = readTopology(\"((C:0.9,(B:0.2)#H1:0.7::0.6):0.6,(#H1:0.6,A:1):0.5);\");\nPhyloNetworks.nameinternalnodes!(net, \"i\"); # \"i\" is a prefix to name internal nodes\nwriteTopology(net)","category":"page"},{"location":"man/mapping_genetree_to_network/","page":"mapping gene trees into the species network","title":"mapping gene trees into the species network","text":"Notice the extra node names in the species phylogeny: i1, i2, and i3 at the root.","category":"page"},{"location":"man/mapping_genetree_to_network/","page":"mapping gene trees into the species network","title":"mapping gene trees into the species network","text":"Next, we use the option nodemapping=true when simulating gene trees, to ask for extra degree-2 nodes in gene trees. These nodes are created each time that a gene tree lineage crosses a node in the species phylogeny, that is, each time that a gene tree lineage crosses a speciation node or a hybridization node. We'll simulate a single tree here each time.","category":"page"},{"location":"man/mapping_genetree_to_network/","page":"mapping gene trees into the species network","title":"mapping gene trees into the species network","text":"using StableRNGs; rng = StableRNG(7); # to replicate randomness, but default RNG is better\ntree_regular = simulatecoalescent(rng, net,1,1)[1];\nwriteTopology(tree_regular, round=true) # regular nodes only\nrng = StableRNG(7); # to replicate the same coalescent simulation\ntree = simulatecoalescent(rng, net,1,1; nodemapping=true)[1];\nwriteTopology(tree, round=true) # extra degree-2 nodes for mapping","category":"page"},{"location":"man/mapping_genetree_to_network/","page":"mapping gene trees into the species network","title":"mapping gene trees into the species network","text":"Notice that regular nodes in the gene tree (nodes that we get without the nodemapping option) don't have names. With the nodemapping option, there are many new nodes, all of degree-2, named after internal nodes in the network (i1, i2, i3). These extra degree-2 nodes and their names are sufficient to map the gene tree into the species network.","category":"page"},{"location":"man/mapping_genetree_to_network/","page":"mapping gene trees into the species network","title":"mapping gene trees into the species network","text":"The network is shown on the left below, with edges annotated by their numbers.","category":"page"},{"location":"man/mapping_genetree_to_network/","page":"mapping gene trees into the species network","title":"mapping gene trees into the species network","text":"R\"svg\"(figname(\"genetree_example1.svg\"), width=7.5, height=3); # hide\nR\"par\"(mar=[.1,.2,.1,.2], oma=[0,0,0,0.8]); R\"layout\"([1 2]); # hide\nplot(net, showedgenumber=true, shownodelabel=true, tipoffset=0.1);\nR\"mtext\"(\"species network\", side=3, line=-1); # hide\nR\"mtext\"(\"grey: population edge number\", side=1, line=-1, cex=0.9); # hide\nplot(tree, edgelabel=DataFrame(number=[e.number for e in tree.edge],\n label=[e.inCycle for e in tree.edge]),\n edgelabelcolor=\"red4\", shownodelabel=true, tipoffset=0.1);\nR\"mtext\"(\"gene tree\", side=3, line=-1); # hide\nR\"mtext\"(\"red (edge inCycle value): population a gene edge maps into\", side=1, line=-2, cex=0.9); # hide\nR\"mtext\"(\"black (node names): speciation/reticulation a node maps to\", side=1, line=-1, cex=0.9); # hide\nR\"dev.off()\" # hide\nnothing # hide","category":"page"},{"location":"man/mapping_genetree_to_network/","page":"mapping gene trees into the species network","title":"mapping gene trees into the species network","text":"(Image: example 1: degree-2 node names in gene tree)","category":"page"},{"location":"man/mapping_genetree_to_network/","page":"mapping gene trees into the species network","title":"mapping gene trees into the species network","text":"In the gene tree (right), each lineage is annotated by the network edge it maps into. Degree-2 nodes appear via their names, such that each horizontal line represents a series of gene lineages, separated from each other by degree-2 nodes. For example, the horizontal line tracing B's ancestry back in time maps into the network like this:","category":"page"},{"location":"man/mapping_genetree_to_network/","page":"mapping gene trees into the species network","title":"mapping gene trees into the species network","text":"from B, go back along edge 2\nmeet hybrid node H1, was inherited from the minor hybrid edge 5,\nfrom speciation node i2, trace back along edge 7,\nmeet the network's root node i3, trace back along the network's root edge 8 before coalescing with the ancestor of the other lineages (which have already coalesced by then).","category":"page"},{"location":"man/mapping_genetree_to_network/#cleaning-gene-trees","page":"mapping gene trees into the species network","title":"cleaning gene trees","text":"","category":"section"},{"location":"man/mapping_genetree_to_network/","page":"mapping gene trees into the species network","title":"mapping gene trees into the species network","text":"Almost all examples below use this mapping information via the extra degree-2 nodes and the extra edges between these nodes.","category":"page"},{"location":"man/mapping_genetree_to_network/","page":"mapping gene trees into the species network","title":"mapping gene trees into the species network","text":"But we may want to \"clean\" gene trees of their degree-2 nodes at some point. This can be done with the PhyloNetworks utility removedegree2nodes!, like this:","category":"page"},{"location":"man/mapping_genetree_to_network/","page":"mapping gene trees into the species network","title":"mapping gene trees into the species network","text":"PhyloNetworks.removedegree2nodes!(tree, true)","category":"page"},{"location":"man/mapping_genetree_to_network/","page":"mapping gene trees into the species network","title":"mapping gene trees into the species network","text":"The option true is to keep the root, even if it's of degree 2.","category":"page"},{"location":"man/converting_coal2generation_units/","page":"converting between units","title":"converting between units","text":"using PhyloNetworks, PhyloPlots, PhyloCoalSimulations, StableRNGs\nfigpath = joinpath(\"..\", \"assets\", \"figures\"); mkpath(figpath)\nfigname(x) = joinpath(figpath, x)\nnet = readTopology(\"((C:0.9,(B:0.2)#H1:0.7::0.6)i1:0.6,(#H1:0.6::0.4,A:1.0)i2:0.5)i3;\");","category":"page"},{"location":"man/converting_coal2generation_units/#converting-between-units","page":"converting between units","title":"converting between units","text":"","category":"section"},{"location":"man/converting_coal2generation_units/#converting-edge-lengths-in-gene-trees","page":"converting between units","title":"converting edge lengths in gene trees","text":"","category":"section"},{"location":"man/converting_coal2generation_units/","page":"converting between units","title":"converting between units","text":"In earlier examples, edge lengths in gene trees are simulated in coalescent units. These lengths can be converted into numbers of generations by multiplying by the effective population size Nₑ, since coalescent units are u = g/Nₑ. This can be done with different Nₑ's across different edges in the network, including the root edge.","category":"page"},{"location":"man/converting_coal2generation_units/","page":"converting between units","title":"converting between units","text":"info: diploid versus haploid Nₑ\nThe formula u = g/Nₑ uses the haploid effective population size Nₑ. For diploid taxa and autosomes, the haploid population size should be twice the diploid population size, for example.","category":"page"},{"location":"man/converting_coal2generation_units/","page":"converting between units","title":"converting between units","text":"Here is an example using the same network and simulated gene tree as earlier.","category":"page"},{"location":"man/converting_coal2generation_units/","page":"converting between units","title":"converting between units","text":"writeTopology(net)\nrng = StableRNG(7); # as in mapping section\ntree = simulatecoalescent(rng, net,1,1; nodemapping=true)[1];\nwriteTopology(tree, round=true)","category":"page"},{"location":"man/converting_coal2generation_units/","page":"converting between units","title":"converting between units","text":"(Image: example 1, same as in mapping section)","category":"page"},{"location":"man/converting_coal2generation_units/","page":"converting between units","title":"converting between units","text":"Let's set Nₑ to 1,000 in all populations, including the population above the root, except in edge 6. For this edge 6 (population leading to species A), let's set its Nₑ to 10,000.","category":"page"},{"location":"man/converting_coal2generation_units/","page":"converting between units","title":"converting between units","text":"Ne = Dict(e.number => 1_000 for e in net.edge);\nrootedgenumber = PhyloCoalSimulations.get_rootedgenumber(net)\npush!(Ne, rootedgenumber => 1_000); # add Ne for the edge above the network's root\nNe[6] = 10_000; # higher population size for the edge to species A\nNe\nwriteTopology(tree, round=true) # lengths in coalescent units: before unit conversion\n# convert edge lengths in gene tree from coalescent units to # generations\nfor e in tree.edge\n e.length = round(e.length * Ne[e.inCycle]) # round: to get integers\nend\nwriteTopology(tree, round=true) # lengths in # of generations","category":"page"},{"location":"man/converting_coal2generation_units/","page":"converting between units","title":"converting between units","text":"Note that the simulation model assumes an infinite-Nₑ approximation, so the rescaling of edge lengths from coalescent units to number of generations will be imperfect for very small populations size. With the extreme Nₑ=1, coalescences should be immediate in a single generation back in time: g=1. Using the approximation, the simulated number of generations will typically be between 0-3 generations. But this is an extreme case, and the approximation should be very good even for moderate Nₑ's.","category":"page"},{"location":"man/converting_coal2generation_units/#number-of-generations-in-the-network-and-gene-trees","page":"converting between units","title":"number of generations in the network and gene trees","text":"","category":"section"},{"location":"man/converting_coal2generation_units/","page":"converting between units","title":"converting between units","text":"If our input network has edge lengths in number of generations, then we need extra information to simulate under the coalescent: we need the effective size of each population. If the population size is constant, then g generations correspond to u = g/Nₑ coalescent units. If the population size varies along a single population edge, then the coalescence rate on that edge is determined by the harmonic mean overlineN_e of the population size: u = int_0^g 1N_e(t) dt = goverlineN_e.","category":"page"},{"location":"man/converting_coal2generation_units/","page":"converting between units","title":"converting between units","text":"Let's assume we have a network with number of generations as edge lengths:","category":"page"},{"location":"man/converting_coal2generation_units/","page":"converting between units","title":"converting between units","text":"net_gen = readTopology(\"((C:900,(B:200)#H1:700::0.6)i1:600,(#H1:600::0.4,A:1000)i2:500)i3;\");","category":"page"},{"location":"man/converting_coal2generation_units/","page":"converting between units","title":"converting between units","text":"and that we have a dictionary listing the (harmonic mean) population size along each edge of the species network, and also along the root edge above the network. Below, we simulate a population size for each population, from a uniform distribution between 1,000 and 1,500.","category":"page"},{"location":"man/converting_coal2generation_units/","page":"converting between units","title":"converting between units","text":"using Distributions\n# uniform distribution between 1000 and 1500, that we can draw from later\nNe_distribution() = rand(rng, DiscreteUniform(1000, 1500));\nNe = Dict(e.number => Ne_distribution() for e in net_gen.edge);\nrootedgenumber = PhyloCoalSimulations.get_rootedgenumber(net_gen)\npush!(Ne, rootedgenumber => Ne_distribution()); # Nₑ above the root\nNe","category":"page"},{"location":"man/converting_coal2generation_units/","page":"converting between units","title":"converting between units","text":"To visually check that we built our Nₑ dictionary correctly (and that our edge lengths are in generations, as we think), we can do this:","category":"page"},{"location":"man/converting_coal2generation_units/","page":"converting between units","title":"converting between units","text":"using DataFrames, RCall\nR\"svg\"(figname(\"net_genNe.svg\"), width=3.5, height=3); # hide\nR\"par\"(mar=[.5,1,.0,.1]); # hide\nplot(net_gen, tipoffset=0.1, showedgelength=true, edgelabelcolor=\"red4\",\n edgelabel=DataFrame(n=[e.number for e in net_gen.edge],\n l=[Ne[e.number] for e in net_gen.edge]));\nR\"text\"(x=1, y=2.5, Ne[rootedgenumber], adj=1, col=\"red4\");\nR\"mtext\"(\"red: Ne values\", side=1, line=-1.5, col=\"red4\");\nR\"mtext\"(\"black: edge lengths\", side=1, line=-0.5);\nR\"dev.off()\" # hide\nnothing # hide","category":"page"},{"location":"man/converting_coal2generation_units/","page":"converting between units","title":"converting between units","text":"(Image: species net with Ne)","category":"page"},{"location":"man/converting_coal2generation_units/","page":"converting between units","title":"converting between units","text":"To simulate gene trees with edge lengths in generations, we can use a convenience wrapper function that takes Nₑ as an extra input to:","category":"page"},{"location":"man/converting_coal2generation_units/","page":"converting between units","title":"converting between units","text":"convert edge lengths to coalescent units in the species phylogeny,\nsimulate gene trees with lengths in coalescent units, then\nconvert gene trees to have lengths in number of generations:","category":"page"},{"location":"man/converting_coal2generation_units/","page":"converting between units","title":"converting between units","text":"genetree_gen = simulatecoalescent(net_gen,3,1, Ne; nodemapping=true);\nwriteMultiTopology(genetree_gen, stdout) # 3 gene trees, lengths in #generations","category":"page"},{"location":"man/converting_coal2generation_units/","page":"converting between units","title":"converting between units","text":"warning: Warning\nWhen Nₑ is given as an extra input to simulatecoalescent, edge lengths in the network are assumed to be in number of generations. If Nₑ is not given as input, then edge lengths are assumed to be in coalescent units.","category":"page"},{"location":"man/correlated_inheritance/","page":"correlated inheritance","title":"correlated inheritance","text":"using PhyloNetworks, PhyloPlots, PhyloCoalSimulations\nfigpath = joinpath(\"..\", \"assets\", \"figures\"); mkpath(figpath)\nfigname(x) = joinpath(figpath, x)","category":"page"},{"location":"man/correlated_inheritance/#correlated-inheritance","page":"correlated inheritance","title":"correlated inheritance","text":"","category":"section"},{"location":"man/correlated_inheritance/","page":"correlated inheritance","title":"correlated inheritance","text":"The typical assumption made by most network inference methods is that of independent inheritance at a hybrid node. It means that if multiple lineages (of a given locus) at present at a hybrid node, then each one is inherited from each parent according to the γ inheritance probabilities independently of the other lineages.","category":"page"},{"location":"man/correlated_inheritance/","page":"correlated inheritance","title":"correlated inheritance","text":"Another model is that the lineages at a hybrid node are all inherited from the same parent, still choosing a (common) parent according to the γ inheritance probabilities. This model has full (and positive) correlation between lineages, and was used by Gerard, Gibbs & Kubatko (2011) for example.","category":"page"},{"location":"man/correlated_inheritance/","page":"correlated inheritance","title":"correlated inheritance","text":"The other extreme might be interesting for modelling allopolyploid events: of 2 lineages of the same locus at a hybrid node, exactly 1 of them comes from one parent, and the other comes from the other parent. This model would have negative correlation between lineages.","category":"page"},{"location":"man/correlated_inheritance/","page":"correlated inheritance","title":"correlated inheritance","text":"By default, simulatecoalescent uses the traditional model, with independent lineages.","category":"page"},{"location":"man/correlated_inheritance/","page":"correlated inheritance","title":"correlated inheritance","text":"It also has an option to simulate lineage inheritance with positive correlation, under a coalescent model described in Fogg, Allman & Ané (2023). For this, lineages' parents are drawn according to a Dirichlet process with base distribution determined by the γ values, and with concentration parameter α = (1-r)/r, that is, r = 1/(1+α), where r is the input inheritance correlation. More specifically, consider 2 individuals (alleles) at a given locus, that have not coalesced yet and are present at a give hybrid node. According to the coalescent model with correlated inheritance, the second individual is inherited from the same parent as the first individual with probability r. And with probability 1-r, the second individual is inherited from any parent (including the parent chosen by the first individual) based on their γ values.","category":"page"},{"location":"man/correlated_inheritance/","page":"correlated inheritance","title":"correlated inheritance","text":"This model can be the result of different loci evolving according to different γ inheritance values. For example, loci that are under selection in the environment where gene flow occurred may be more likely to be passed through gene flow, whereas loci that are involved in reproduction barriers might be less likely to be passed through gene flow. This would result in different sets of γ values across different loci, and correlated inheritance overall.","category":"page"},{"location":"man/correlated_inheritance/","page":"correlated inheritance","title":"correlated inheritance","text":"At one extreme, with correlation r=1 we have α = 0 and all lineages inherit from the same (randomly sampled) parent. This is the common inheritance model. The independence model corresponds to the other extreme r=0 and α infinite.","category":"page"},{"location":"man/correlated_inheritance/","page":"correlated inheritance","title":"correlated inheritance","text":"The same correlation r (or concentration α) parameter is used at all hybrid nodes, but the Dirichlet process is applied independently across hybrid nodes.","category":"page"},{"location":"man/correlated_inheritance/","page":"correlated inheritance","title":"correlated inheritance","text":"The effect of inheritance correlation can be seen on this species network with a single species \"A\" from which we will sample 6 individuals:","category":"page"},{"location":"man/correlated_inheritance/","page":"correlated inheritance","title":"correlated inheritance","text":"net = readTopology(\"((A:0.0)#H1:3::0.6,#H1:3)root;\")\nusing RCall # hide\nR\"svg\"(figname(\"net1taxon.svg\"), width=3, height=2); # hide\nR\"par\"(mar=[0,0,0,0]); # hide\nplot(net, shownodelabel=true, showgamma=true, showedgenumber=true, tipoffset=0.1);\nR\"dev.off()\" # hide\nnothing # hide","category":"page"},{"location":"man/correlated_inheritance/","page":"correlated inheritance","title":"correlated inheritance","text":"(Image: 1-taxon network with 2-cycle)","category":"page"},{"location":"man/correlated_inheritance/","page":"correlated inheritance","title":"correlated inheritance","text":"We will now simulate 3 gene trees, each with 6 individuals sampled from A: with the independence model (equivalent to r=0 correlation), and then with r=0.5 and r=0.99 inheritance correlations.","category":"page"},{"location":"man/correlated_inheritance/","page":"correlated inheritance","title":"correlated inheritance","text":"using Random; Random.seed!(12);\ngt1 = simulatecoalescent(net, 1, 6; inheritancecorrelation=0, nodemapping=true)[1];\ngt2 = simulatecoalescent(net, 1, 6; inheritancecorrelation=0.5, nodemapping=true)[1];\ngt3 = simulatecoalescent(net, 1, 6; inheritancecorrelation=0.99, nodemapping=true)[1];\n\n# build data frames to label each edge with the species lineage it evolved in\nusing DataFrames\nspeciespath(phy) = DataFrame(\n number = [e.number for e in phy.edge],\n label = [e.inCycle for e in phy.edge]\n)\nel1 = speciespath(gt1); el2 = speciespath(gt2); el3 = speciespath(gt3);\n\nR\"svg\"(figname(\"net1taxon_3genetrees.svg\"), width=8, height=3); # hide\nR\"par\"(mar=[0,0,0,0], oma=[0,0,2,0]); R\"layout\"([1 2 3])# hide\nplot(gt1, shownodelabel=true, edgelabel=el1, edgelabelcolor=el1.label, tipoffset=0.1);\nplot(gt2, shownodelabel=true, edgelabel=el2, edgelabelcolor=el2.label, tipoffset=0.1);\nplot(gt3, shownodelabel=true, edgelabel=el3, edgelabelcolor=el3.label, tipoffset=0.1);\nR\"mtext\"(\"gene tree \" .* [\"1\",\"2\",\"3\"] .* \"\\n(r=\" .* [\"0\",\"0.5\",\"0.99\"] .* \")\", side=3, line=-1.1, outer=true, at=[0.2,0.5,0.8]); # hide\nR\"dev.off()\" # hide\nnothing # hide","category":"page"},{"location":"man/correlated_inheritance/","page":"correlated inheritance","title":"correlated inheritance","text":"(Image: 3 gene trees on 1-taxon network with inheritance correlation)","category":"page"},{"location":"man/correlated_inheritance/","page":"correlated inheritance","title":"correlated inheritance","text":"In all cases, any lineage has a probability γ=0.6 to come from species edge 2 (labeled in red), and probability γ=0.4 to come from species edge 3 (in green). When the inheritance correlation r increases, lineages have an increased preference to come from the same parent as other lineages (at the same locus).","category":"page"},{"location":"man/more_examples/","page":"more examples","title":"more examples","text":"using PhyloNetworks, PhyloCoalSimulations\nnet = readTopology(\"((C:0.9,(B:0.2)#H1:0.7::0.6)I1:0.6,(#H1:0.6::0.4,A:1.0)I2:0.5)I3;\");\nusing Random; Random.seed!(261); # as in mapping block\ntree = simulatecoalescent(net,1,1; nodemapping=true)[1];","category":"page"},{"location":"man/more_examples/#example-uses","page":"more examples","title":"example uses","text":"","category":"section"},{"location":"man/more_examples/","page":"more examples","title":"more examples","text":"We are re-using the same network and simulated tree as before:","category":"page"},{"location":"man/more_examples/","page":"more examples","title":"more examples","text":"writeTopology(net)\nwriteTopology(tree, round=true)","category":"page"},{"location":"man/more_examples/","page":"more examples","title":"more examples","text":"(Image: example 1, same as in mapping section)","category":"page"},{"location":"man/more_examples/#counting-deep-coalescences","page":"more examples","title":"counting deep coalescences","text":"","category":"section"},{"location":"man/more_examples/","page":"more examples","title":"more examples","text":"The number of deep coalescences can be quantified as the number of \"extra\" lineages due to incomplete lineage sorting, that can be calculated from embedding the gene tree into the species phylogeny (see Maddison 1997 for species trees). For an edge in the network, say edge 7 going from I2 to I3 going in back in time, lineage sorting is complete if all the gene lineages entering the edge (at I2) coalesce into a single gene lineage by the time they exit the edge (at I3). If they don't, the number of extra lineages is k-1 where k is the number of lineages \"exiting\" the edge, for that particular edge in the species network and that particular gene tree. The total number of extra lineages, for a given gene tree, is the sum across all edges in the species phylogeny.","category":"page"},{"location":"man/more_examples/","page":"more examples","title":"more examples","text":"In our gene tree above, we can count the number of lineages that exit each species edge using the degree-2 mapping nodes, then count how many lineages are \"extra\". We do so below using utilities mappingnodes to iterate over degree-2 mapping nodes and population_mappedto to extract the mapping information.","category":"page"},{"location":"man/more_examples/","page":"more examples","title":"more examples","text":"# dictionary to store the count of extra lineages exiting each network edge. initialized at 0\nedge_count = Dict(e.number => 0 for e in net.edge)\nconst PCS = PhyloCoalSimulations; # for lazy typing!\nfor n in PCS.mappingnodes(tree) # iterate over degree-2 mapping nodes in the gene tree\n child = getchildedge(n)\n popid = PCS.population_mappedto(child) # number of species edge that 'n' came from\n # sanity check below\n isnothing(popid) && error(\"\"\"population ID not found for the child edge of\n node number $(n.number) mapping to species node $(n.name).\"\"\")\n edge_count[popid] += 1 # increment by 1 the number of lineages exiting population edge 'popid'\nend\nedge_count","category":"page"},{"location":"man/more_examples/","page":"more examples","title":"more examples","text":"From this, we see two interesting things.","category":"page"},{"location":"man/more_examples/","page":"more examples","title":"more examples","text":"0 lineages exited edge number 3 in the species network: it's the hybrid edge from H1 to I3 (going back in time). That's because the only lineage at H1 was interited from I2, so there weren't any lineage evolving through edge 3.\n2 lineages exited edge number 7 (going from I2 to I3 back in time), so that's 1 extra lineage. All other edges look as expected, with a single gene lineage exiting from them.","category":"page"},{"location":"man/more_examples/","page":"more examples","title":"more examples","text":"We can now calculate the total number of extra lineages:","category":"page"},{"location":"man/more_examples/","page":"more examples","title":"more examples","text":"filter!(p -> p.second > 0, edge_count) # filter out edges without any gene lineage\nmap!(k -> k-1, values(edge_count)) # calculate number of \"extras\": k-1\ndeepcoalescence = sum(values(edge_count)) # sum 'extras' over all edges in the network","category":"page"},{"location":"man/more_examples/","page":"more examples","title":"more examples","text":"On the particular gene tree we simulated, we counted 1 deep coalescence.","category":"page"},{"location":"man/more_examples/#number-of-lineages-inherited-via-gene-flow","page":"more examples","title":"number of lineages inherited via gene flow","text":"","category":"section"},{"location":"man/more_examples/","page":"more examples","title":"more examples","text":"Our network has inheritance γ=0.4 on the minor edge, which we'll call the \"gene flow\" edge, and γ=0.6 on the major hybrid edge, parent to H1 on the major tree. But we may be interested in the realized proportion of lineages inherited from each parent at H1, realized in the gene trees we actually simulated. To do so, we can count the number of gene lineages that are mapped to each hybrid edge in the network.","category":"page"},{"location":"man/more_examples/","page":"more examples","title":"more examples","text":"This mapping is stored in the edge attribute .inCycle internally, but it's best to access it via the function population_mappedto (as the internal representation may change). From the plot above, the minor \"gene flow\" edge is edge number 5 and the major hybrid edge has number 3. So we can count the gene lineages inherited via gene flow as the number of gene tree edges with inCycle equal to 5.","category":"page"},{"location":"man/more_examples/","page":"more examples","title":"more examples","text":"If the gene trees have been saved to a file and later read from this file, then the .inCycle attributes are no longer stored in memory. In this case, we can retrieve the mapping information by the internal node names. The edges going through gene flow are those whose child node is named \"H1\" and parent node is named \"I2\".","category":"page"},{"location":"man/more_examples/","page":"more examples","title":"more examples","text":"We use the first option with the .inCycle attribute below. We get that our one simulated gene tree was indeed inherited via gene flow:","category":"page"},{"location":"man/more_examples/","page":"more examples","title":"more examples","text":"sum(e.inCycle == 5 for e in tree.edge) # or:\nsum(PCS.population_mappedto(e) == 5 for e in tree.edge)\n# or define a function to do this for any edge, so we can re-use later:\nnlineages_through(edgeID, gt) = sum(PCS.population_mappedto(e) == edgeID for e in gt.edge);\nnlineages_through(5, tree) # same as before!\nnlineages_through(3, tree) # lineages that went through edge 3, the major edge.","category":"page"},{"location":"man/more_examples/","page":"more examples","title":"more examples","text":"To make this more interesting, we can simulate many gene trees then count how many of their lineages were inherited via gene flow. If we ask for 2 individuals in species B, then each gene may have 2 lineages that enter the hybrid node H1, if the two B individuals fail to coalesce. In that case, it's possible that one individual lineage was inherited via gene flow, and the other not. We'll calculate the gene flow proportion among all these lineages. This proportion should be close (but not exactly equal) to the theoretical γ=0.4 from the network model.","category":"page"},{"location":"man/more_examples/","page":"more examples","title":"more examples","text":"ngenes = 100;\ngenetrees = simulatecoalescent(net, ngenes, Dict(\"B\"=>2, \"A\"=>1, \"C\"=>1); nodemapping=true);\nlength(genetrees)\nnlineages_geneflow = sum(nlineages_through(5,gt) for gt in genetrees)\nnlineages_major = sum(nlineages_through(3,gt) for gt in genetrees)\n# realized γ, close to 0.4:\nproportion_geneflow = nlineages_geneflow / (nlineages_geneflow + nlineages_major)","category":"page"},{"location":"man/more_examples/#rate-variation-across-species","page":"more examples","title":"rate variation across species","text":"","category":"section"},{"location":"man/more_examples/","page":"more examples","title":"more examples","text":"The gene trees resulting from simulatecoalescent have their edge lengths in coalescent units. One may want to convert them to substitutions per site, so as to simulate molecular sequences along these gene trees. The mapping information is important to allow for different rates of molecular evolution across different species. Here is an example to do this.","category":"page"},{"location":"man/more_examples/","page":"more examples","title":"more examples","text":"We will use the Distributions package to simulate rates from a log-normal distribution across species, that is, across edges in the species network.","category":"page"},{"location":"man/more_examples/","page":"more examples","title":"more examples","text":"using Distributions\nlognormal_rate_dist = LogNormal(-0.125, 0.5) # μ = -σ²/2 to get a mean of 1.\nnetworkedge_rate = Dict(e.number => rand(lognormal_rate_dist) for e in net.edge)\n# add entry for the edge above the network's root. Find its number first.\nrootedgenumber = PhyloCoalSimulations.get_rootedgenumber(net)\npush!(networkedge_rate, rootedgenumber => rand(lognormal_rate_dist))\nwriteTopology(tree, round=true, digits=4) # before rate variation\n# multiply the length of each gene lineage by the rate of the species edge it maps into\nfor e in tree.edge\n e.length *= networkedge_rate[e.inCycle]\nend\nwriteTopology(tree, round=true, digits=4) # after rate variation","category":"page"},{"location":"man/getting_started/","page":"getting started","title":"getting started","text":"using PhyloNetworks, PhyloCoalSimulations, RCall\nfigpath = joinpath(\"..\", \"assets\", \"figures\"); mkpath(figpath)\nfigname(x) = joinpath(figpath, x)\nusing Random; Random.seed!(432)","category":"page"},{"location":"man/getting_started/#getting-started","page":"getting started","title":"getting started","text":"","category":"section"},{"location":"man/getting_started/#installation","page":"getting started","title":"installation","text":"","category":"section"},{"location":"man/getting_started/","page":"getting started","title":"getting started","text":"To install Julia see here and to install Julia packages, see here. To install PhyloCoalSimulations in the Julia REPL (as well as PhyloNetworks for many utilities), enter package mode with ] and do this:","category":"page"},{"location":"man/getting_started/","page":"getting started","title":"getting started","text":"add PhyloCoalSimulations\nadd PhyloNetworks","category":"page"},{"location":"man/getting_started/","page":"getting started","title":"getting started","text":"or do this in julian mode:","category":"page"},{"location":"man/getting_started/","page":"getting started","title":"getting started","text":"using Pkg\nPkg.add(\"PhyloCoalSimulations\")\nPkg.add(\"PhyloNetworks\")","category":"page"},{"location":"man/getting_started/#basic-simulation-example","page":"getting started","title":"basic simulation example","text":"","category":"section"},{"location":"man/getting_started/#example-network","page":"getting started","title":"example network","text":"","category":"section"},{"location":"man/getting_started/","page":"getting started","title":"getting started","text":"For a basic example, we use a simple 3-species network plotted below. On the left, the plot shows the edge numbers (black) and the γ inheritance values (blue). On the right, the length of horizontal lines are proportional to edge lengths, and the plot shows the edge length values.","category":"page"},{"location":"man/getting_started/","page":"getting started","title":"getting started","text":"using PhyloNetworks\nnet = readTopology(\"((C:0.9,(B:0.2)#H1:0.7::0.6):0.6,(#H1:0.6,A:1):0.5);\");\nusing PhyloPlots\nR\"svg\"(figname(\"net3taxa.svg\"), width=6, height=3); # hide\nR\"par\"(mar=[.1,.1,.1,.1]); R\"layout\"([1 2]); # hide\nplot(net, showedgenumber=true, shownodelabel=true, showgamma=true, tipoffset=0.1);\nR\"mtext\"(\"in grey: edge numbers\", side=1, line=-1); # hide\nplot(net, showedgelength=true, useedgelength=true, tipoffset=0.1);\nR\"mtext\"(\"in black: edge lengths\", side=1, line=-1); # hide\nR\"dev.off()\" # hide\nnothing # hide","category":"page"},{"location":"man/getting_started/","page":"getting started","title":"getting started","text":"(Image: 3-taxon network)","category":"page"},{"location":"man/getting_started/","page":"getting started","title":"getting started","text":"Note that this example network is not time consistent: the length of the path from the root to the hybridization node H1 is different depending if we go through the major edge (0.6+0.7=1.3) or the minor edge (0.5+0.6=1.1).","category":"page"},{"location":"man/getting_started/","page":"getting started","title":"getting started","text":"Coalescent simulations can be performed along such networks, also along non-ultrametric networks. If the network is ultrametric (time-consistent, and with all tips at the same distance from the root), then gene trees will also be ultrametric.","category":"page"},{"location":"man/getting_started/","page":"getting started","title":"getting started","text":"note: rooting\nThe input network needs to be rooted correctly, as some properties of gene trees may depend on the root. If a network was estimated using a method that estimates a semidirected network (where the root placement is arbitrary) then make sure to root this network appropriately, e.g. using an outgroup, before simulating genes along this network.","category":"page"},{"location":"man/getting_started/#basic-example:-simulate,-save-to-file,-plot","page":"getting started","title":"basic example: simulate, save to file, plot","text":"","category":"section"},{"location":"man/getting_started/","page":"getting started","title":"getting started","text":"We use simulatecoalescent to simulate gene trees along this network. Below, we simulate 2 gene trees, with 1 individual per species.","category":"page"},{"location":"man/getting_started/","page":"getting started","title":"getting started","text":"trees = simulatecoalescent(net, 2, 1)","category":"page"},{"location":"man/getting_started/","page":"getting started","title":"getting started","text":"Branch lengths are assumed to be in coalescent units in the species network (number of generations / effective population size), and edge lengths in gene trees are also in coalescent units.","category":"page"},{"location":"man/getting_started/","page":"getting started","title":"getting started","text":"We can work with these gene trees within Julia with downstream code, and/or we can save them to a file:","category":"page"},{"location":"man/getting_started/","page":"getting started","title":"getting started","text":"writeMultiTopology(trees, stdout) # write them to standout output (screen here)\nwriteMultiTopology(trees, \"genetrees.phy\") # warning: will overwrite \"genetrees.phy\" if this file existed\nrm(\"genetrees.phy\") # hide","category":"page"},{"location":"man/getting_started/","page":"getting started","title":"getting started","text":"Let's plot these 2 gene trees. In the plots below, we annotate each edge with its attribute that tracked the network edge on which the coalescent event occured (where the gene tree lineage originated, going back in time). For example, the gene lineage that ends in A is always mapped to network edge 6, which is the number of the external edge to A in the network (see plot of network above on the left).","category":"page"},{"location":"man/getting_started/","page":"getting started","title":"getting started","text":"R\"svg\"(figname(\"genetrees_gettingstarted_1.svg\"), width=6, height=3); # hide\nR\"par\"(mar=[.1,.1,.1,.1]); R\"layout\"([1 2]); # hide\nusing DataFrames\nfor i in 1:2\n gt = trees[i]\n plot(gt, tipoffset=0.1,\n edgelabel=DataFrame(number = [e.number for e in gt.edge],\n label = [e.inCycle for e in gt.edge]));\n R\"mtext\"(\"gene $i\", line=-1) # hide\nend\nR\"mtext\"(\"numbers: network edge each gene lineage maps to, at time of coalescence.\\n8 = number of edge above the network root\", side=1, line=-1, outer=true); # hide\nR\"dev.off()\" # hide\nnothing # hide","category":"page"},{"location":"man/getting_started/","page":"getting started","title":"getting started","text":"(Image: 3-taxon network)","category":"page"},{"location":"man/getting_started/#several-individuals-per-species","page":"getting started","title":"several individuals per species","text":"","category":"section"},{"location":"man/getting_started/","page":"getting started","title":"getting started","text":"We can ask for more individuals. To simulate 3 individuals / species for example:","category":"page"},{"location":"man/getting_started/","page":"getting started","title":"getting started","text":"simulatecoalescent(net, 1, 3) # 1 gene tree only. 3 individuals in each species","category":"page"},{"location":"man/getting_started/","page":"getting started","title":"getting started","text":"We can also ask for varying numbers of individuals. For example, we simulate below 2 individuals in A and 1 individual in each of B and C, using a dictionary to map species to their number of individuals:","category":"page"},{"location":"man/getting_started/","page":"getting started","title":"getting started","text":"genetrees = simulatecoalescent(net, 1, Dict(\"A\"=>2, \"B\"=>1, \"C\"=>1));\nwriteTopology(genetrees[1])","category":"page"},{"location":"man/getting_started/","page":"getting started","title":"getting started","text":"We can set 0 individuals within a species to simulate missing data.","category":"page"},{"location":"man/getting_started/","page":"getting started","title":"getting started","text":"genetrees = simulatecoalescent(net, 3, Dict(\"A\"=>2, \"B\"=>1, \"C\"=>0));\nwriteMultiTopology(genetrees, stdout)","category":"page"},{"location":"lib/public/#public-documentation","page":"public","title":"public documentation","text":"","category":"section"},{"location":"lib/public/","page":"public","title":"public","text":"Documentation for PhyloCoalSimulations's public (exported) functions. Most functions are internal (not exported).","category":"page"},{"location":"lib/public/#functions-and-types","page":"public","title":"functions & types","text":"","category":"section"},{"location":"lib/public/","page":"public","title":"public","text":"Modules = [PhyloCoalSimulations]\nPrivate = false\nOrder = [:function,:type]","category":"page"},{"location":"lib/public/#PhyloCoalSimulations.simulatecoalescent-Tuple{PhyloNetworks.HybridNetwork, Vararg{Any}}","page":"public","title":"PhyloCoalSimulations.simulatecoalescent","text":"simulatecoalescent([rng::AbstractRNG,] net, nloci, nindividuals;\n nodemapping=false, inheritancecorrelation=0.0)\n\nSimulate nloci gene trees with nindividuals from each species under the multispecies network coalescent, along network net whose branch lengths are assumed to be in coalescent units (ratio: number of generations / effective population size). The coalescent model uses the infinite-population-size approximation.\n\nThe random number generator rng is optional.\n\nOutput: vector of gene trees, of length nloci.\n\nnindividuals can be a single integer, or a dictionary listing the number of individuals to be simulated for each species.\n\nIf nodemapping is true, each simulated gene tree is augmented with degree-2 nodes that can be mapped to speciation or hybridization events. The mapping of gene tree nodes & edges to network edges is carried by the .inCycle attribute. The mapping of gene tree nodes to network nodes is carried by the .name attribute. Namely:\n\nA degree-3 node (1 parent + 2 children) represents a coalescent event that occurred along a population edge in net. Its .inCycle attribute is set to the number of that network population edge. Its parent edge has its .inCycle attribute also set to the number of the population edge that it originated from.\nThe gene tree's root node (of degree 2) represents a coalescent event along the network's root edge. Its .inCycle attribute is the number assigned to the network's root edge, which is set by get_rootedgenumber as the maximum edge number + 1.\nA leaf (or degree-1 node) represents an individual. It maps to a species in net. The individual leaf name is set to the species name if nindividuals is 1. Otherwise, its name is set to speciesname_i where i is the individual number in that species. Its inCycle attribute is the default -1.\nA non-root degree-2 node represents a speciation or hybridization and maps to a population node in net. Its inCycle attribute is the default -1. Its name is set to network node name, if it exists. If the network node has no name, the gene tree node is given a name built from the network node number.\n\nBy default, lineages at a hybrid node come from a parent (chosen according to inheritance probabilities γ) independently across lineages. Positive dependence can be simulated with option inheritancecorrelation. For example, if this correlation is set to 1, then all lineages inherit from the same (randomly sampled) parent. More generally, the lineages' parents are simulated according to a Dirichlet process with base distribution determined by the γ values, and with concentration parameter α = (1-r)/r, that is, r = 1/(1+α), where r is the input inheritance correlation. For more details about this model, please read the package manual or refer to Fogg, Allman & Ané (2023).\n\nAssumptions:\n\nnet must have non-missing edge lengths and γ values.\nIf nindividuals is a dictionary, it must have a key for all species, with the same spelling of species names in its keys as in the tip labels of net.\n\nexamples\n\njulia> using PhyloNetworks\n\njulia> net = readTopology(\"(A:1,B:1);\"); # branch lengths of 1 coalescent unit\n\njulia> using Random; Random.seed!(54321); # for replicability of examples below\n\njulia> simulatecoalescent(net, 2, 1) # 2 gene trees, 1 individual/species\n2-element Vector{HybridNetwork}:\n PhyloNetworks.HybridNetwork, Rooted Network\n2 edges\n3 nodes: 2 tips, 0 hybrid nodes, 1 internal tree nodes.\ntip labels: B, A\n(B:1.023,A:1.023);\n\n PhyloNetworks.HybridNetwork, Rooted Network\n2 edges\n3 nodes: 2 tips, 0 hybrid nodes, 1 internal tree nodes.\ntip labels: B, A\n(B:2.328,A:2.328);\n\n\njulia> simulatecoalescent(net, 1, 3)[1] # 1 gene tree, 3 individuals/species\nPhyloNetworks.HybridNetwork, Rooted Network\n10 edges\n11 nodes: 6 tips, 0 hybrid nodes, 5 internal tree nodes.\ntip labels: B_2, B_1, B_3, A_3, ...\n(((B_2:0.12,B_1:0.12):2.39,B_3:2.51):0.692,(A_3:0.518,(A_2:0.461,A_1:0.461):0.057):2.684);\n\n\njulia> simulatecoalescent(net, 1, Dict(\"A\"=>2, \"B\"=>1))[1] # 2 individuals in A, 1 in B\nPhyloNetworks.HybridNetwork, Rooted Network\n4 edges\n5 nodes: 3 tips, 0 hybrid nodes, 2 internal tree nodes.\ntip labels: B, A_2, A_1\n(B:2.801,(A_2:0.344,A_1:0.344):2.457);\n\nIn the next example, we use a custom random number generator (RNG), to show how to do so. In this example we use a stable RNG to make the example reproducible across julia versions. However, this particular RNG is not recommended. The RNG used by default is better (e.g. much more efficient).\n\njulia> # using Pkg; Pkg.add(\"StableRNGs\") # to install StableRNGs if not done earlier\n\njulia> using StableRNGs\n\njulia> rng = StableRNG(791);\n\njulia> tree1 = simulatecoalescent(rng, net,2,2; nodemapping=true)[1]; # first gene tree only\n\njulia> writeTopology(tree1, round=true)\n\"(((B_2:0.82,B_1:0.82):0.18)minus2:0.992,((A_1:1.0)minus2:0.427,(A_2:1.0)minus2:0.427):0.565);\"\n\njulia> PhyloNetworks.nameinternalnodes!(net, \"i\"); writeTopology(net)\n\"(A:1.0,B:1.0)i1;\"\n\njulia> tree1 = simulatecoalescent(rng, net,2,2; nodemapping=true)[1]; writeTopology(tree1, round=true)\n\"((B_2:1.0)i1:0.621,((B_1:1.0)i1:0.018,((A_1:0.61,A_2:0.61):0.39)i1:0.018):0.604);\"\n\njulia> printNodes(net)\nnode leaf hybrid hasHybEdge name inCycle edges'numbers\n1 true false false A -1 1 \n2 true false false B -1 2 \n-2 false false false i1 -1 1 2 \n\njulia> printNodes(tree1)\nnode leaf hybrid hasHybEdge name inCycle edges'numbers\n10 false false false 3 8 9 \n8 false false false i1 -1 5 8 \n5 true false false B_2 -1 5 \n9 false false false 3 7 6 9 \n7 false false false i1 -1 4 7 \n4 true false false B_1 -1 4 \n6 false false false i1 -1 3 6 \n3 false false false 1 1 2 3 \n1 true false false A_1 -1 1 \n2 true false false A_2 -1 2 \n\njulia> [(tree_edge_number = e.number, pop_edge_number = e.inCycle) for e in tree1.edge]\n9-element Vector{@NamedTuple{tree_edge_number::Int64, pop_edge_number::Int64}}:\n (tree_edge_number = 8, pop_edge_number = 3)\n (tree_edge_number = 5, pop_edge_number = 2)\n (tree_edge_number = 9, pop_edge_number = 3)\n (tree_edge_number = 7, pop_edge_number = 3)\n (tree_edge_number = 4, pop_edge_number = 2)\n (tree_edge_number = 6, pop_edge_number = 3)\n (tree_edge_number = 3, pop_edge_number = 1)\n (tree_edge_number = 1, pop_edge_number = 1)\n (tree_edge_number = 2, pop_edge_number = 1)\n\n\n\n\n\n","category":"method"},{"location":"lib/public/#PhyloCoalSimulations.simulatecoalescent-Tuple{Random.AbstractRNG, PhyloNetworks.HybridNetwork, Integer, Any, Any}","page":"public","title":"PhyloCoalSimulations.simulatecoalescent","text":"simulatecoalescent([rng::AbstractRNG,] net, nloci, nindividuals, populationsize;\n nodemapping=false, round_generationnumber=true,\n inheritancecorrelation=0.0)\n\nSimulate nloci gene trees with nindividuals from each species under the multispecies network coalescent, along network net, whose branch lengths are assumed to be in number of generations. populationsize should be a single number, assumed to be the (haploid) effective population size Nₑ, constant across the species phylogeny. Alternatively, populationsize can be a dictionary mapping the number of each edge in net to its Nₑ, including an extra edge number for the population above the root of the network.\n\nCoalescent units are then calculated as u=g/Nₑ where g is the edge length in net (number of generations), and the coalescent model is applied using the infinite-population-size approximation.\n\nOutput: vector of gene trees with edge lengths in number of generations, calculated as g=uNₑ and then rounded to be an integer, unless round_generationnumber is false.\n\nwarning: Warning\nWhen populationsize Nₑ is not provided as input, all edge lengths are in coalescent units. When populationsize is given as an argument, all edge lengths are in number of generations. The second method (using # generation and Nₑ as input) is a wrapper around the first (using coalescent units).\n\njulia> using PhyloNetworks\n\njulia> net = readTopology(\"(A:500,B:500);\"); # branch lengths of 100 generations\n\njulia> Ne = Dict(e.number => 1_000 for e in net.edge);\n\njulia> rootedgenumber = PhyloCoalSimulations.get_rootedgenumber(net)\n3\n\njulia> push!(Ne, rootedgenumber => 2_000) # Ne for population above the root\nDict{Int64, Int64} with 3 entries:\n 2 => 1000\n 3 => 2000\n 1 => 1000\n\njulia> using Random; Random.seed!(54321); # for replicability of example below\n\njulia> genetrees = simulatecoalescent(net, 2, 1, Ne);\n\njulia> writeMultiTopology(genetrees, stdout) # branch lengths: number of generations\n(B:546.0,A:546.0);\n(B:3155.0,A:3155.0);\n\n\n\n\n\n\n","category":"method"},{"location":"lib/public/#index","page":"public","title":"index","text":"","category":"section"},{"location":"lib/public/","page":"public","title":"public","text":"Pages = [\"public.md\"]","category":"page"},{"location":"","page":"home","title":"home","text":"CurrentModule = PhyloCoalSimulations","category":"page"},{"location":"#PhyloCoalSimulations","page":"home","title":"PhyloCoalSimulations","text":"","category":"section"},{"location":"","page":"home","title":"home","text":"PhyloCoalSimulations is a Julia package to simulate phylogenies under the coalescent. It depends on PhyloNetworks for the phylogenetic data structures, and manipulation of phylogenies.","category":"page"},{"location":"","page":"home","title":"home","text":"References: please see bibtex entries here.","category":"page"},{"location":"","page":"home","title":"home","text":"for this package and the network coalescent model with inheritance correlation:\nFogg, Allman & Ané (2023). PhyloCoalSimulations: A simulator for network multispecies coalescent models, including a new extension for the inheritance of gene flow. Systematic Biology, 72(5):1171–1179. doi:10.1093/sysbio/syad030.\nfor PhyloNetworks, which this package depends on:\nSolís-Lemus, Bastide & Ané (2017). PhyloNetworks: a package for phylogenetic networks. Molecular Biology and Evolution, 34(12):3292–3298. doi:10.1093/molbev/msx235","category":"page"},{"location":"","page":"home","title":"home","text":"For a tutorial, see the manual:","category":"page"},{"location":"","page":"home","title":"home","text":"Pages = [\n \"man/getting_started.md\",\n \"man/mapping_genetree_to_network.md\",\n \"man/converting_coal2generation_units.md\",\n \"man/correlated_inheritance.md\",\n \"man/more_examples.md\",\n]\nDepth = 3","category":"page"},{"location":"","page":"home","title":"home","text":"For help on individual functions, see the library:","category":"page"},{"location":"","page":"home","title":"home","text":"Pages = [\n \"lib/public.md\",\n \"lib/internal.md\",\n]\nDepth = 3","category":"page"},{"location":"lib/internal/#internal-documentation","page":"internals","title":"internal documentation","text":"","category":"section"},{"location":"lib/internal/","page":"internals","title":"internals","text":"Documentation for PhyloCoalSimulations's internal functions. These functions are not exported and their access (API) should not be considered stable. But they can still be used, like this for example: PhyloCoalSimulations.foo() for a function named foo().","category":"page"},{"location":"lib/internal/#functions-and-types","page":"internals","title":"functions & types","text":"","category":"section"},{"location":"lib/internal/","page":"internals","title":"internals","text":"Modules = [PhyloCoalSimulations]\nPublic = false\nOrder = [:type,:function]","category":"page"},{"location":"lib/internal/#PhyloCoalSimulations.mappingnodes","page":"internals","title":"PhyloCoalSimulations.mappingnodes","text":"mappingnodes(gene tree)\n\nType to define an iterator over degree-2 mapping nodes in a gene tree, assuming these degree-2 nodes (other than the root) have a name to map them to nodes in a species phylogeny. See ismappingnode.\n\n\n\n\n\n","category":"type"},{"location":"lib/internal/#PhyloCoalSimulations.coalescence_edge-NTuple{4, Any}","page":"internals","title":"PhyloCoalSimulations.coalescence_edge","text":"coalescence_edge(edge1, edge2, number, populationid)\n\nCreate a coalescence between edges 1 and 2: with a new parent node n numbered number and a new parent edge e above the parent node, of length 0 and numbered number. Both n.inCycle and e.inCycle are set to populationid.\n\n\n\n\n\n","category":"method"},{"location":"lib/internal/#PhyloCoalSimulations.convert2tree!-Tuple{PhyloNetworks.Node}","page":"internals","title":"PhyloCoalSimulations.convert2tree!","text":"convert2tree!(rootnode)\n\nReturn a network with all nodes and edges that can be reached from rootnode. Warning: Assumes that edges are correctly directed (with correct isChild1 attribute) and that the graph is a tree. This is not checked.\n\nIf the root node is still attached to an incomplete root edge, this edge & node are first disconnected.\n\n\n\n\n\n","category":"method"},{"location":"lib/internal/#PhyloCoalSimulations.get_rootedgenumber-Tuple{Any}","page":"internals","title":"PhyloCoalSimulations.get_rootedgenumber","text":"get_rootedgenumber(network)\n\n1 + maximum of 0 and of all the network's edge numbers: can be used as a unique identifier of the edge above the network's root.\n\n\n\n\n\n","category":"method"},{"location":"lib/internal/#PhyloCoalSimulations.initializetip","page":"internals","title":"PhyloCoalSimulations.initializetip","text":"initializetip(species::AbstractString, individual::AbstractString,\n number::Integer, delim=\"\"::AbstractString)\n\nCreate a leaf node and a pendant edge of length 0, incident to each other, both numbered number. Return the pendant edge. The leaf name is made by concatenating species, delim and individual.\n\n\n\n\n\n","category":"function"},{"location":"lib/internal/#PhyloCoalSimulations.initializetipforest","page":"internals","title":"PhyloCoalSimulations.initializetipforest","text":"initializetipforest(speciesnode::Node, nindividuals::Integer,\n number::Integer, delim)\n\nVector of pendant leaf edges, with leaves named after speciesnode, and numbered with consecutive number IDs starting at number. If nindividuals is 1, then the leaf name is simply the species name. Otherwise, then the leaf names include the individual number and the default delimiter is _. For example, if the species name is s then leaf names are: s_1, s_2, etc. by default. Pendant leaf edges have inCycle set to the number of the corresponding edge in the species network.\n\n\n\n\n\n","category":"function"},{"location":"lib/internal/#PhyloCoalSimulations.ismappingnode-Tuple{PhyloNetworks.Node}","page":"internals","title":"PhyloCoalSimulations.ismappingnode","text":"ismappingnode(node)\n\nBoolean: true if node is of degree 2, has a single child, and has a name. (The root is of degree-2 but is not a mapping node).\n\n\n\n\n\n","category":"method"},{"location":"lib/internal/#PhyloCoalSimulations.map2population!-NTuple{4, Any}","page":"internals","title":"PhyloCoalSimulations.map2population!","text":"map2population!(forest, population_node, populationid, nextlineageID)\n\nExtend each incomplete edge in the forest with a new degree-2 node n and a new incomplete edge e, with the following information to map n and e into the species phylogeny:\n\ne.inCycle is set to populationid, and\nn.name is set to population_node.name if this name is non-empty, or string(population_node.number) otherwise (with any negative sign replaced by the string \"minus\").\n\ne.number and n.number are set to nextlineageID, which is incremented by 1 for each incomplete edge in the forest.\n\nThe forest is updated to contain the newly-created incomplete edges, replacing the old incomplete (and now complete) edges.\n\nOutput: nextlineageID, incremented by the number of newly created degree-2 lineages.\n\nexample\n\njulia> using PhyloNetworks; net = readTopology(\"(A:1,B:1);\");\n\njulia> leafA = net.node[1]; edge2A_number = net.edge[1].number;\n\njulia> f = PhyloCoalSimulations.initializetipforest(leafA, 2, 4); # 2 edges, numbered 4 & 5\n\njulia> PhyloCoalSimulations.map2population!(f, leafA, edge2A_number, 6)\n8\n\njulia> length(f)\n2\n\njulia> f[2]\nPhyloNetworks.EdgeT{PhyloNetworks.Node}:\n number:7\n length:0.0\n attached to 1 node(s) (parent first): 7\n\njulia> [e.node[1].name for e in f]\n2-element Vector{String}:\n \"A\"\n \"A\"\n\n\n\n\n\n","category":"method"},{"location":"lib/internal/#PhyloCoalSimulations.population_mappedto-Tuple{Union{PhyloNetworks.EdgeT{PhyloNetworks.Node}, PhyloNetworks.Node}}","page":"internals","title":"PhyloCoalSimulations.population_mappedto","text":"population_mappedto(edge or node)\n\nIdentifier of the population (edge in the species network) that a gene tree's edge or a node is mapped onto, or nothing if not mapped. For example, coalescent nodes in gene trees map to a node in the species phylogeny, instead of mapping to an edge.\n\n\n\n\n\n","category":"method"},{"location":"lib/internal/#PhyloCoalSimulations.simulatecoal_onepopulation!-Tuple{AbstractVector, Vararg{Any}}","page":"internals","title":"PhyloCoalSimulations.simulatecoal_onepopulation!","text":"simulatecoal_onepopulation!([rng::AbstractRNG,]\n lineagelist,\n population_length,\n nextlineageID,\n populationID=-1\n)\n\nSimulate the coalescent process within a single population of length given in coalescent units, starting from lineages in lineagelist. This list should be a vector of incomplete edges, that is, edges incident to a single node only.\n\nVector of incomplete edges, whose lengths have been increased, is modified in place. New nodes and their parent edges are created by coalescent events, numbered with consecutive integers starting at nextlineageID.\n\nThe random number generator rng is optional.\n\nOutput: nextlineageID, incremented by number of new lineages.\n\nIn lineages, edge lengths are also considered in coalescent units.\n\nThe newly created nodes and edges have their .inCycle attribute set to populationID, so as to track the mapping of gene lineages to populations in the species phylogeny.\n\nexamples\n\njulia> PhyloCoalSimulations.simulatecoal_onepopulation!([], 2.0, 1)\n1\n\njulia> e1 = PhyloCoalSimulations.initializetip(\"s\",\"1\",1,\"\");\n\njulia> e2 = PhyloCoalSimulations.initializetip(\"s\",\"2\",2,\"\");\n\njulia> forest = [e1,e2];\n\njulia> using Random; Random.seed!(7690);\n\njulia> PhyloCoalSimulations.simulatecoal_onepopulation!(forest, Inf, 3);\n\njulia> PhyloCoalSimulations.convert2tree!(forest[1].node[1])\nPhyloNetworks.HybridNetwork, Rooted Network\n2 edges\n3 nodes: 2 tips, 0 hybrid nodes, 1 internal tree nodes.\ntip labels: s2, s1\n(s2:0.302,s1:0.302);\n\n\n\n\n\n\n","category":"method"},{"location":"lib/internal/#index","page":"internals","title":"index","text":"","category":"section"},{"location":"lib/internal/","page":"internals","title":"internals","text":"Pages = [\"internal.md\"]","category":"page"}]
+[{"location":"man/mapping_genetree_to_network/","page":"mapping gene trees into the species network","title":"mapping gene trees into the species network","text":"using PhyloNetworks, PhyloPlots, PhyloCoalSimulations, RCall, DataFrames, StableRNGs\nfigpath = joinpath(\"..\", \"assets\", \"figures\"); mkpath(figpath)\nfigname(x) = joinpath(figpath, x)","category":"page"},{"location":"man/mapping_genetree_to_network/#mapping-gene-trees-into-the-species-phylogeny","page":"mapping gene trees into the species network","title":"mapping gene trees into the species phylogeny","text":"","category":"section"},{"location":"man/mapping_genetree_to_network/","page":"mapping gene trees into the species network","title":"mapping gene trees into the species network","text":"Nodes in gene trees can be mapped to a node or an edge in the species phylogeny. Edges in gene trees can be mapped to an edge in the species phylogeny. Attributes of nodes and edges are used to carry this mapping information, as detailed in the documentation of function simulatecoalescent.","category":"page"},{"location":"man/mapping_genetree_to_network/","page":"mapping gene trees into the species network","title":"mapping gene trees into the species network","text":"We give examples below of how we may use this mapping information.","category":"page"},{"location":"man/mapping_genetree_to_network/#naming-internal-nodes","page":"mapping gene trees into the species network","title":"naming internal nodes","text":"","category":"section"},{"location":"man/mapping_genetree_to_network/","page":"mapping gene trees into the species network","title":"mapping gene trees into the species network","text":"First, it's useful to name internal nodes in the network, to which we can later map nodes in the gene tree.","category":"page"},{"location":"man/mapping_genetree_to_network/","page":"mapping gene trees into the species network","title":"mapping gene trees into the species network","text":"net = readTopology(\"((C:0.9,(B:0.2)#H1:0.7::0.6):0.6,(#H1:0.6,A:1):0.5);\");\nPhyloNetworks.nameinternalnodes!(net, \"i\"); # \"i\" is a prefix to name internal nodes\nwriteTopology(net)","category":"page"},{"location":"man/mapping_genetree_to_network/","page":"mapping gene trees into the species network","title":"mapping gene trees into the species network","text":"Notice the extra node names in the species phylogeny: i1, i2, and i3 at the root.","category":"page"},{"location":"man/mapping_genetree_to_network/","page":"mapping gene trees into the species network","title":"mapping gene trees into the species network","text":"Next, we use the option nodemapping=true when simulating gene trees, to ask for extra degree-2 nodes in gene trees. These nodes are created each time that a gene tree lineage crosses a node in the species phylogeny, that is, each time that a gene tree lineage crosses a speciation node or a hybridization node. We'll simulate a single tree here each time.","category":"page"},{"location":"man/mapping_genetree_to_network/","page":"mapping gene trees into the species network","title":"mapping gene trees into the species network","text":"using StableRNGs; rng = StableRNG(7); # to replicate randomness, but default RNG is better\ntree_regular = simulatecoalescent(rng, net,1,1)[1];\nwriteTopology(tree_regular, round=true) # regular nodes only\nrng = StableRNG(7); # to replicate the same coalescent simulation\ntree = simulatecoalescent(rng, net,1,1; nodemapping=true)[1];\nwriteTopology(tree, round=true) # extra degree-2 nodes for mapping","category":"page"},{"location":"man/mapping_genetree_to_network/","page":"mapping gene trees into the species network","title":"mapping gene trees into the species network","text":"Notice that regular nodes in the gene tree (nodes that we get without the nodemapping option) don't have names. With the nodemapping option, there are many new nodes, all of degree-2, named after internal nodes in the network (i1, i2, i3). These extra degree-2 nodes and their names are sufficient to map the gene tree into the species network.","category":"page"},{"location":"man/mapping_genetree_to_network/","page":"mapping gene trees into the species network","title":"mapping gene trees into the species network","text":"The network is shown on the left below, with edges annotated by their numbers.","category":"page"},{"location":"man/mapping_genetree_to_network/","page":"mapping gene trees into the species network","title":"mapping gene trees into the species network","text":"R\"svg\"(figname(\"genetree_example1.svg\"), width=7.5, height=3); # hide\nR\"par\"(mar=[.1,.2,.1,.2], oma=[0,0,0,0.8]); R\"layout\"([1 2]); # hide\nplot(net, showedgenumber=true, shownodelabel=true, tipoffset=0.1);\nR\"mtext\"(\"species network\", side=3, line=-1); # hide\nR\"mtext\"(\"grey: population edge number\", side=1, line=-1, cex=0.9); # hide\nplot(tree, edgelabel=DataFrame(number=[e.number for e in tree.edge],\n label=[e.inCycle for e in tree.edge]),\n edgelabelcolor=\"red4\", shownodelabel=true, tipoffset=0.1);\nR\"mtext\"(\"gene tree\", side=3, line=-1); # hide\nR\"mtext\"(\"red (edge inCycle value): population a gene edge maps into\", side=1, line=-2, cex=0.9); # hide\nR\"mtext\"(\"black (node names): speciation/reticulation a node maps to\", side=1, line=-1, cex=0.9); # hide\nR\"dev.off()\" # hide\nnothing # hide","category":"page"},{"location":"man/mapping_genetree_to_network/","page":"mapping gene trees into the species network","title":"mapping gene trees into the species network","text":"(Image: example 1: degree-2 node names in gene tree)","category":"page"},{"location":"man/mapping_genetree_to_network/","page":"mapping gene trees into the species network","title":"mapping gene trees into the species network","text":"In the gene tree (right), each lineage is annotated by the network edge it maps into. Degree-2 nodes appear via their names, such that each horizontal line represents a series of gene lineages, separated from each other by degree-2 nodes. For example, the horizontal line tracing B's ancestry back in time maps into the network like this:","category":"page"},{"location":"man/mapping_genetree_to_network/","page":"mapping gene trees into the species network","title":"mapping gene trees into the species network","text":"from B, go back along edge 2\nmeet hybrid node H1, was inherited from the minor hybrid edge 5,\nfrom speciation node i2, trace back along edge 7,\nmeet the network's root node i3, trace back along the network's root edge 8 before coalescing with the ancestor of the other lineages (which have already coalesced by then).","category":"page"},{"location":"man/mapping_genetree_to_network/#cleaning-gene-trees","page":"mapping gene trees into the species network","title":"cleaning gene trees","text":"","category":"section"},{"location":"man/mapping_genetree_to_network/","page":"mapping gene trees into the species network","title":"mapping gene trees into the species network","text":"Almost all examples below use this mapping information via the extra degree-2 nodes and the extra edges between these nodes.","category":"page"},{"location":"man/mapping_genetree_to_network/","page":"mapping gene trees into the species network","title":"mapping gene trees into the species network","text":"But we may want to \"clean\" gene trees of their degree-2 nodes at some point. This can be done with the PhyloNetworks utility removedegree2nodes!, like this:","category":"page"},{"location":"man/mapping_genetree_to_network/","page":"mapping gene trees into the species network","title":"mapping gene trees into the species network","text":"PhyloNetworks.removedegree2nodes!(tree, true)","category":"page"},{"location":"man/mapping_genetree_to_network/","page":"mapping gene trees into the species network","title":"mapping gene trees into the species network","text":"The option true is to keep the root, even if it's of degree 2.","category":"page"},{"location":"man/converting_coal2generation_units/","page":"converting between units","title":"converting between units","text":"using PhyloNetworks, PhyloPlots, PhyloCoalSimulations, StableRNGs\nfigpath = joinpath(\"..\", \"assets\", \"figures\"); mkpath(figpath)\nfigname(x) = joinpath(figpath, x)\nnet = readTopology(\"((C:0.9,(B:0.2)#H1:0.7::0.6)i1:0.6,(#H1:0.6::0.4,A:1.0)i2:0.5)i3;\");","category":"page"},{"location":"man/converting_coal2generation_units/#converting-between-units","page":"converting between units","title":"converting between units","text":"","category":"section"},{"location":"man/converting_coal2generation_units/#converting-edge-lengths-in-gene-trees","page":"converting between units","title":"converting edge lengths in gene trees","text":"","category":"section"},{"location":"man/converting_coal2generation_units/","page":"converting between units","title":"converting between units","text":"In earlier examples, edge lengths in gene trees are simulated in coalescent units. These lengths can be converted into numbers of generations by multiplying by the effective population size Nₑ, since coalescent units are u = g/Nₑ. This can be done with different Nₑ's across different edges in the network, including the root edge.","category":"page"},{"location":"man/converting_coal2generation_units/","page":"converting between units","title":"converting between units","text":"info: diploid versus haploid Nₑ\nThe formula u = g/Nₑ uses the haploid effective population size Nₑ. For diploid taxa and autosomes, the haploid population size should be twice the diploid population size, for example.","category":"page"},{"location":"man/converting_coal2generation_units/","page":"converting between units","title":"converting between units","text":"Here is an example using the same network and simulated gene tree as earlier.","category":"page"},{"location":"man/converting_coal2generation_units/","page":"converting between units","title":"converting between units","text":"writeTopology(net)\nrng = StableRNG(7); # as in mapping section\ntree = simulatecoalescent(rng, net,1,1; nodemapping=true)[1];\nwriteTopology(tree, round=true)","category":"page"},{"location":"man/converting_coal2generation_units/","page":"converting between units","title":"converting between units","text":"(Image: example 1, same as in mapping section)","category":"page"},{"location":"man/converting_coal2generation_units/","page":"converting between units","title":"converting between units","text":"Let's set Nₑ to 1,000 in all populations, including the population above the root, except in edge 6. For this edge 6 (population leading to species A), let's set its Nₑ to 10,000.","category":"page"},{"location":"man/converting_coal2generation_units/","page":"converting between units","title":"converting between units","text":"Ne = Dict(e.number => 1_000 for e in net.edge);\nrootedgenumber = PhyloCoalSimulations.get_rootedgenumber(net)\npush!(Ne, rootedgenumber => 1_000); # add Ne for the edge above the network's root\nNe[6] = 10_000; # higher population size for the edge to species A\nNe\nwriteTopology(tree, round=true) # lengths in coalescent units: before unit conversion\n# convert edge lengths in gene tree from coalescent units to # generations\nfor e in tree.edge\n e.length = round(e.length * Ne[e.inCycle]) # round: to get integers\nend\nwriteTopology(tree, round=true) # lengths in # of generations","category":"page"},{"location":"man/converting_coal2generation_units/","page":"converting between units","title":"converting between units","text":"Note that the simulation model assumes an infinite-Nₑ approximation, so the rescaling of edge lengths from coalescent units to number of generations will be imperfect for very small populations size. With the extreme Nₑ=1, coalescences should be immediate in a single generation back in time: g=1. Using the approximation, the simulated number of generations will typically be between 0-3 generations. But this is an extreme case, and the approximation should be very good even for moderate Nₑ's.","category":"page"},{"location":"man/converting_coal2generation_units/#number-of-generations-in-the-network-and-gene-trees","page":"converting between units","title":"number of generations in the network and gene trees","text":"","category":"section"},{"location":"man/converting_coal2generation_units/","page":"converting between units","title":"converting between units","text":"If our input network has edge lengths in number of generations, then we need extra information to simulate under the coalescent: we need the effective size of each population. If the population size is constant, then g generations correspond to u = g/Nₑ coalescent units. If the population size varies along a single population edge, then the coalescence rate on that edge is determined by the harmonic mean overlineN_e of the population size: u = int_0^g 1N_e(t) dt = goverlineN_e.","category":"page"},{"location":"man/converting_coal2generation_units/","page":"converting between units","title":"converting between units","text":"Let's assume we have a network with number of generations as edge lengths:","category":"page"},{"location":"man/converting_coal2generation_units/","page":"converting between units","title":"converting between units","text":"net_gen = readTopology(\"((C:900,(B:200)#H1:700::0.6)i1:600,(#H1:600::0.4,A:1000)i2:500)i3;\");","category":"page"},{"location":"man/converting_coal2generation_units/","page":"converting between units","title":"converting between units","text":"and that we have a dictionary listing the (harmonic mean) population size along each edge of the species network, and also along the root edge above the network. Below, we simulate a population size for each population, from a uniform distribution between 1,000 and 1,500.","category":"page"},{"location":"man/converting_coal2generation_units/","page":"converting between units","title":"converting between units","text":"using Distributions\n# uniform distribution between 1000 and 1500, that we can draw from later\nNe_distribution() = rand(rng, DiscreteUniform(1000, 1500));\nNe = Dict(e.number => Ne_distribution() for e in net_gen.edge);\nrootedgenumber = PhyloCoalSimulations.get_rootedgenumber(net_gen)\npush!(Ne, rootedgenumber => Ne_distribution()); # Nₑ above the root\nNe","category":"page"},{"location":"man/converting_coal2generation_units/","page":"converting between units","title":"converting between units","text":"To visually check that we built our Nₑ dictionary correctly (and that our edge lengths are in generations, as we think), we can do this:","category":"page"},{"location":"man/converting_coal2generation_units/","page":"converting between units","title":"converting between units","text":"using DataFrames, RCall\nR\"svg\"(figname(\"net_genNe.svg\"), width=3.5, height=3); # hide\nR\"par\"(mar=[.5,1,.0,.1]); # hide\nplot(net_gen, tipoffset=0.1, showedgelength=true, edgelabelcolor=\"red4\",\n edgelabel=DataFrame(n=[e.number for e in net_gen.edge],\n l=[Ne[e.number] for e in net_gen.edge]));\nR\"text\"(x=1, y=2.5, Ne[rootedgenumber], adj=1, col=\"red4\");\nR\"mtext\"(\"red: Ne values\", side=1, line=-1.5, col=\"red4\");\nR\"mtext\"(\"black: edge lengths\", side=1, line=-0.5);\nR\"dev.off()\" # hide\nnothing # hide","category":"page"},{"location":"man/converting_coal2generation_units/","page":"converting between units","title":"converting between units","text":"(Image: species net with Ne)","category":"page"},{"location":"man/converting_coal2generation_units/","page":"converting between units","title":"converting between units","text":"To simulate gene trees with edge lengths in generations, we can use a convenience wrapper function that takes Nₑ as an extra input to:","category":"page"},{"location":"man/converting_coal2generation_units/","page":"converting between units","title":"converting between units","text":"convert edge lengths to coalescent units in the species phylogeny,\nsimulate gene trees with lengths in coalescent units, then\nconvert gene trees to have lengths in number of generations:","category":"page"},{"location":"man/converting_coal2generation_units/","page":"converting between units","title":"converting between units","text":"genetree_gen = simulatecoalescent(net_gen,3,1, Ne; nodemapping=true);\nwriteMultiTopology(genetree_gen, stdout) # 3 gene trees, lengths in #generations","category":"page"},{"location":"man/converting_coal2generation_units/","page":"converting between units","title":"converting between units","text":"warning: Warning\nWhen Nₑ is given as an extra input to simulatecoalescent, edge lengths in the network are assumed to be in number of generations. If Nₑ is not given as input, then edge lengths are assumed to be in coalescent units.","category":"page"},{"location":"man/correlated_inheritance/","page":"correlated inheritance","title":"correlated inheritance","text":"using PhyloNetworks, PhyloPlots, PhyloCoalSimulations\nfigpath = joinpath(\"..\", \"assets\", \"figures\"); mkpath(figpath)\nfigname(x) = joinpath(figpath, x)","category":"page"},{"location":"man/correlated_inheritance/#correlated-inheritance","page":"correlated inheritance","title":"correlated inheritance","text":"","category":"section"},{"location":"man/correlated_inheritance/","page":"correlated inheritance","title":"correlated inheritance","text":"The typical assumption made by most network inference methods is that of independent inheritance at a hybrid node. It means that if multiple lineages (of a given locus) at present at a hybrid node, then each one is inherited from each parent according to the γ inheritance probabilities independently of the other lineages.","category":"page"},{"location":"man/correlated_inheritance/","page":"correlated inheritance","title":"correlated inheritance","text":"Another model is that the lineages at a hybrid node are all inherited from the same parent, still choosing a (common) parent according to the γ inheritance probabilities. This model has full (and positive) correlation between lineages, and was used by Gerard, Gibbs & Kubatko (2011) for example.","category":"page"},{"location":"man/correlated_inheritance/","page":"correlated inheritance","title":"correlated inheritance","text":"The other extreme might be interesting for modelling allopolyploid events: of 2 lineages of the same locus at a hybrid node, exactly 1 of them comes from one parent, and the other comes from the other parent. This model would have negative correlation between lineages.","category":"page"},{"location":"man/correlated_inheritance/","page":"correlated inheritance","title":"correlated inheritance","text":"By default, simulatecoalescent uses the traditional model, with independent lineages.","category":"page"},{"location":"man/correlated_inheritance/","page":"correlated inheritance","title":"correlated inheritance","text":"It also has an option to simulate lineage inheritance with positive correlation, under a coalescent model described in Fogg, Allman & Ané (2023). For this, lineages' parents are drawn according to a Dirichlet process with base distribution determined by the γ values, and with concentration parameter α = (1-r)/r, that is, r = 1/(1+α), where r is the input inheritance correlation. More specifically, consider 2 individuals (alleles) at a given locus, that have not coalesced yet and are present at a give hybrid node. According to the coalescent model with correlated inheritance, the second individual is inherited from the same parent as the first individual with probability r. And with probability 1-r, the second individual is inherited from any parent (including the parent chosen by the first individual) based on their γ values.","category":"page"},{"location":"man/correlated_inheritance/","page":"correlated inheritance","title":"correlated inheritance","text":"This model can be the result of different loci evolving according to different γ inheritance values. For example, loci that are under selection in the environment where gene flow occurred may be more likely to be passed through gene flow, whereas loci that are involved in reproduction barriers might be less likely to be passed through gene flow. This would result in different sets of γ values across different loci, and correlated inheritance overall.","category":"page"},{"location":"man/correlated_inheritance/","page":"correlated inheritance","title":"correlated inheritance","text":"At one extreme, with correlation r=1 we have α = 0 and all lineages inherit from the same (randomly sampled) parent. This is the common inheritance model. The independence model corresponds to the other extreme r=0 and α infinite.","category":"page"},{"location":"man/correlated_inheritance/","page":"correlated inheritance","title":"correlated inheritance","text":"The same correlation r (or concentration α) parameter is used at all hybrid nodes, but the Dirichlet process is applied independently across hybrid nodes.","category":"page"},{"location":"man/correlated_inheritance/","page":"correlated inheritance","title":"correlated inheritance","text":"The effect of inheritance correlation can be seen on this species network with a single species \"A\" from which we will sample 6 individuals:","category":"page"},{"location":"man/correlated_inheritance/","page":"correlated inheritance","title":"correlated inheritance","text":"net = readTopology(\"((A:0.0)#H1:3::0.6,#H1:3)root;\")\nusing RCall # hide\nR\"svg\"(figname(\"net1taxon.svg\"), width=3, height=2); # hide\nR\"par\"(mar=[0,0,0,0]); # hide\nplot(net, shownodelabel=true, showgamma=true, showedgenumber=true, tipoffset=0.1);\nR\"dev.off()\" # hide\nnothing # hide","category":"page"},{"location":"man/correlated_inheritance/","page":"correlated inheritance","title":"correlated inheritance","text":"(Image: 1-taxon network with 2-cycle)","category":"page"},{"location":"man/correlated_inheritance/","page":"correlated inheritance","title":"correlated inheritance","text":"We will now simulate 3 gene trees, each with 6 individuals sampled from A: with the independence model (equivalent to r=0 correlation), and then with r=0.5 and r=0.99 inheritance correlations.","category":"page"},{"location":"man/correlated_inheritance/","page":"correlated inheritance","title":"correlated inheritance","text":"using Random; Random.seed!(12);\ngt1 = simulatecoalescent(net, 1, 6; inheritancecorrelation=0, nodemapping=true)[1];\ngt2 = simulatecoalescent(net, 1, 6; inheritancecorrelation=0.5, nodemapping=true)[1];\ngt3 = simulatecoalescent(net, 1, 6; inheritancecorrelation=0.99, nodemapping=true)[1];\n\n# build data frames to label each edge with the species lineage it evolved in\nusing DataFrames\nspeciespath(phy) = DataFrame(\n number = [e.number for e in phy.edge],\n label = [e.inCycle for e in phy.edge]\n)\nel1 = speciespath(gt1); el2 = speciespath(gt2); el3 = speciespath(gt3);\n\nR\"svg\"(figname(\"net1taxon_3genetrees.svg\"), width=8, height=3); # hide\nR\"par\"(mar=[0,0,0,0], oma=[0,0,2,0]); R\"layout\"([1 2 3])# hide\nplot(gt1, shownodelabel=true, edgelabel=el1, edgelabelcolor=el1.label, tipoffset=0.1);\nplot(gt2, shownodelabel=true, edgelabel=el2, edgelabelcolor=el2.label, tipoffset=0.1);\nplot(gt3, shownodelabel=true, edgelabel=el3, edgelabelcolor=el3.label, tipoffset=0.1);\nR\"mtext\"(\"gene tree \" .* [\"1\",\"2\",\"3\"] .* \"\\n(r=\" .* [\"0\",\"0.5\",\"0.99\"] .* \")\", side=3, line=-1.1, outer=true, at=[0.2,0.5,0.8]); # hide\nR\"dev.off()\" # hide\nnothing # hide","category":"page"},{"location":"man/correlated_inheritance/","page":"correlated inheritance","title":"correlated inheritance","text":"(Image: 3 gene trees on 1-taxon network with inheritance correlation)","category":"page"},{"location":"man/correlated_inheritance/","page":"correlated inheritance","title":"correlated inheritance","text":"In all cases, any lineage has a probability γ=0.6 to come from species edge 2 (labeled in red), and probability γ=0.4 to come from species edge 3 (in green). When the inheritance correlation r increases, lineages have an increased preference to come from the same parent as other lineages (at the same locus).","category":"page"},{"location":"man/more_examples/","page":"more examples","title":"more examples","text":"using PhyloNetworks, PhyloCoalSimulations, StableRNGs\nnet = readTopology(\"((C:0.9,(B:0.2)#H1:0.7::0.6)i1:0.6,(#H1:0.6::0.4,A:1.0)i2:0.5)i3;\");\nusing Random; Random.seed!(261) # for examples below that use default RNG\nrng = StableRNG(7); # as in mapping block\ntree = simulatecoalescent(rng, net,1,1; nodemapping=true)[1];","category":"page"},{"location":"man/more_examples/#example-uses","page":"more examples","title":"example uses","text":"","category":"section"},{"location":"man/more_examples/","page":"more examples","title":"more examples","text":"We are re-using the same network and simulated tree as before:","category":"page"},{"location":"man/more_examples/","page":"more examples","title":"more examples","text":"writeTopology(net)\nwriteTopology(tree, round=true)","category":"page"},{"location":"man/more_examples/","page":"more examples","title":"more examples","text":"(Image: example 1, same as in mapping section)","category":"page"},{"location":"man/more_examples/#counting-deep-coalescences","page":"more examples","title":"counting deep coalescences","text":"","category":"section"},{"location":"man/more_examples/","page":"more examples","title":"more examples","text":"The number of deep coalescences can be quantified as the number of \"extra\" lineages due to incomplete lineage sorting, that can be calculated from embedding the gene tree into the species phylogeny (see Maddison 1997 for species trees). For an edge in the network, say edge 7 going from i2 to i3 going in back in time, lineage sorting is complete if all the gene lineages entering the edge (at i2) coalesce into a single gene lineage by the time they exit the edge (at i3). If they don't, the number of extra lineages is k-1 where k is the number of lineages \"exiting\" the edge, for that particular edge in the species network and that particular gene tree. The total number of extra lineages, for a given gene tree, is the sum across all edges in the species phylogeny.","category":"page"},{"location":"man/more_examples/","page":"more examples","title":"more examples","text":"In our gene tree above, we can count the number of lineages that exit each species edge using the degree-2 mapping nodes, then count how many lineages are \"extra\". We do so below using utilities mappingnodes to iterate over degree-2 mapping nodes and population_mappedto to extract the mapping information.","category":"page"},{"location":"man/more_examples/","page":"more examples","title":"more examples","text":"# dictionary to store the count of extra lineages exiting each network edge. initialized at 0\nedge_count = Dict(e.number => 0 for e in net.edge)\nconst PCS = PhyloCoalSimulations; # for lazy typing!\nfor n in PCS.mappingnodes(tree) # iterate over degree-2 mapping nodes in the gene tree\n child = getchildedge(n)\n popid = PCS.population_mappedto(child) # number of species edge that 'n' came from\n # sanity check below\n isnothing(popid) && error(\"\"\"population ID not found for the child edge of\n node number $(n.number) mapping to species node $(n.name).\"\"\")\n edge_count[popid] += 1 # increment by 1 the number of lineages exiting population edge 'popid'\nend\nedge_count","category":"page"},{"location":"man/more_examples/","page":"more examples","title":"more examples","text":"From this, we see two interesting things.","category":"page"},{"location":"man/more_examples/","page":"more examples","title":"more examples","text":"0 lineages exited edge number 3 in the species network: it's the hybrid edge from H1 to i3 (going back in time). That's because the only lineage at H1 was interited from i2, so there weren't any lineage evolving through edge 3.\n2 lineages exited edge number 7 (going from i2 to i3 back in time), so that's 1 extra lineage. All other edges look as expected, with a single gene lineage exiting from them.","category":"page"},{"location":"man/more_examples/","page":"more examples","title":"more examples","text":"We can now calculate the total number of extra lineages:","category":"page"},{"location":"man/more_examples/","page":"more examples","title":"more examples","text":"filter!(p -> p.second > 0, edge_count) # filter out edges without any gene lineage\nmap!(k -> k-1, values(edge_count)) # calculate number of \"extras\": k-1\ndeepcoalescence = sum(values(edge_count)) # sum 'extras' over all edges in the network","category":"page"},{"location":"man/more_examples/","page":"more examples","title":"more examples","text":"On the particular gene tree we simulated, we counted 1 deep coalescence.","category":"page"},{"location":"man/more_examples/#number-of-lineages-inherited-via-gene-flow","page":"more examples","title":"number of lineages inherited via gene flow","text":"","category":"section"},{"location":"man/more_examples/","page":"more examples","title":"more examples","text":"Our network has inheritance γ=0.4 on the minor edge, which we'll call the \"gene flow\" edge, and γ=0.6 on the major hybrid edge, parent to H1 on the major tree. But we may be interested in the realized proportion of lineages inherited from each parent at H1, realized in the gene trees we actually simulated. To do so, we can count the number of gene lineages that are mapped to each hybrid edge in the network.","category":"page"},{"location":"man/more_examples/","page":"more examples","title":"more examples","text":"This mapping is stored in the edge attribute .inCycle internally, but it's best to access it via the function population_mappedto (as the internal representation may change). From the plot above, the minor \"gene flow\" edge is edge number 5 and the major hybrid edge has number 3. So we can count the gene lineages inherited via gene flow as the number of gene tree edges with inCycle equal to 5.","category":"page"},{"location":"man/more_examples/","page":"more examples","title":"more examples","text":"If the gene trees have been saved to a file and later read from this file, then the .inCycle attributes are no longer stored in memory. In this case, we can retrieve the mapping information by the internal node names. The edges going through gene flow are those whose child node is named \"H1\" and parent node is named \"i2\".","category":"page"},{"location":"man/more_examples/","page":"more examples","title":"more examples","text":"We use the first option with the .inCycle attribute below. We get that our one simulated gene tree was indeed inherited via gene flow:","category":"page"},{"location":"man/more_examples/","page":"more examples","title":"more examples","text":"sum(e.inCycle == 5 for e in tree.edge) # or:\nsum(PCS.population_mappedto(e) == 5 for e in tree.edge)","category":"page"},{"location":"man/more_examples/","page":"more examples","title":"more examples","text":"Next we define a function to do this for any edge, so we can re-use later:","category":"page"},{"location":"man/more_examples/","page":"more examples","title":"more examples","text":"nlineages_through(edgeID, gt) = sum(PCS.population_mappedto(e) == edgeID for e in gt.edge);\nnlineages_through(5, tree) # same as before: now done via our new function\nnlineages_through(3, tree) # lineages that went through edge 3, the major edge.","category":"page"},{"location":"man/more_examples/","page":"more examples","title":"more examples","text":"To make this more interesting, we can simulate many gene trees then count how many of their lineages were inherited via gene flow. If we ask for 2 individuals in species B, then each gene may have 2 lineages that enter the hybrid node H1, if the two B individuals fail to coalesce. In that case, it's possible that one individual lineage was inherited via gene flow, and the other not. We'll calculate the gene flow proportion among all these lineages. This proportion should be close (but not exactly equal) to the theoretical γ=0.4 from the network model.","category":"page"},{"location":"man/more_examples/","page":"more examples","title":"more examples","text":"ngenes = 100;\ngenetrees = simulatecoalescent(net, ngenes, Dict(\"B\"=>2, \"A\"=>1, \"C\"=>1); nodemapping=true);\nlength(genetrees)\nnlineages_geneflow = sum(nlineages_through(5,gt) for gt in genetrees)\nnlineages_major = sum(nlineages_through(3,gt) for gt in genetrees)\nproportion_geneflow = nlineages_geneflow / (nlineages_geneflow + nlineages_major)\n# realized γ, close to 0.4","category":"page"},{"location":"man/more_examples/#rate-variation-across-species","page":"more examples","title":"rate variation across species","text":"","category":"section"},{"location":"man/more_examples/","page":"more examples","title":"more examples","text":"The gene trees resulting from simulatecoalescent have their edge lengths in coalescent units. One may want to convert them to substitutions per site, so as to simulate molecular sequences along these gene trees. The mapping information is important to allow for different rates of molecular evolution across different species. Here is an example to do this.","category":"page"},{"location":"man/more_examples/","page":"more examples","title":"more examples","text":"We will use the Distributions package to simulate rates from a log-normal distribution across species, that is, across edges in the species network.","category":"page"},{"location":"man/more_examples/","page":"more examples","title":"more examples","text":"using Distributions\nlognormal_rate_dist = LogNormal(-0.125, 0.5) # μ = -σ²/2 to get a mean of 1.\nnetworkedge_rate = Dict(e.number => rand(lognormal_rate_dist) for e in net.edge)","category":"page"},{"location":"man/more_examples/","page":"more examples","title":"more examples","text":"Next we want to simulate a rate for the ancestral edge above the network's root. We find its number first, then add a entry to our dictionary of rates","category":"page"},{"location":"man/more_examples/","page":"more examples","title":"more examples","text":"rootedgenumber = PhyloCoalSimulations.get_rootedgenumber(net)\npush!(networkedge_rate, rootedgenumber => rand(lognormal_rate_dist))\nwriteTopology(tree, round=true, digits=4) # before rate variation","category":"page"},{"location":"man/more_examples/","page":"more examples","title":"more examples","text":"Finally, we multiply the length of each gene lineage by the rate of the species edge it maps into:","category":"page"},{"location":"man/more_examples/","page":"more examples","title":"more examples","text":"for e in tree.edge\n e.length *= networkedge_rate[e.inCycle]\nend\nwriteTopology(tree, round=true, digits=4) # after rate variation","category":"page"},{"location":"man/getting_started/","page":"getting started","title":"getting started","text":"using PhyloNetworks, PhyloCoalSimulations, RCall\nfigpath = joinpath(\"..\", \"assets\", \"figures\"); mkpath(figpath)\nfigname(x) = joinpath(figpath, x)\nusing Random; Random.seed!(432)","category":"page"},{"location":"man/getting_started/#getting-started","page":"getting started","title":"getting started","text":"","category":"section"},{"location":"man/getting_started/#installation","page":"getting started","title":"installation","text":"","category":"section"},{"location":"man/getting_started/","page":"getting started","title":"getting started","text":"To install Julia see here and to install Julia packages, see here. To install PhyloCoalSimulations in the Julia REPL (as well as PhyloNetworks for many utilities), enter package mode with ] and do this:","category":"page"},{"location":"man/getting_started/","page":"getting started","title":"getting started","text":"add PhyloCoalSimulations\nadd PhyloNetworks","category":"page"},{"location":"man/getting_started/","page":"getting started","title":"getting started","text":"or do this in julian mode:","category":"page"},{"location":"man/getting_started/","page":"getting started","title":"getting started","text":"using Pkg\nPkg.add(\"PhyloCoalSimulations\")\nPkg.add(\"PhyloNetworks\")","category":"page"},{"location":"man/getting_started/#basic-simulation-example","page":"getting started","title":"basic simulation example","text":"","category":"section"},{"location":"man/getting_started/#example-network","page":"getting started","title":"example network","text":"","category":"section"},{"location":"man/getting_started/","page":"getting started","title":"getting started","text":"For a basic example, we use a simple 3-species network plotted below. On the left, the plot shows the edge numbers (black) and the γ inheritance values (blue). On the right, the length of horizontal lines are proportional to edge lengths, and the plot shows the edge length values.","category":"page"},{"location":"man/getting_started/","page":"getting started","title":"getting started","text":"using PhyloNetworks\nnet = readTopology(\"((C:0.9,(B:0.2)#H1:0.7::0.6):0.6,(#H1:0.6,A:1):0.5);\");\nusing PhyloPlots\nR\"svg\"(figname(\"net3taxa.svg\"), width=6, height=3); # hide\nR\"par\"(mar=[.1,.1,.1,.1]); R\"layout\"([1 2]); # hide\nplot(net, showedgenumber=true, shownodelabel=true, showgamma=true, tipoffset=0.1);\nR\"mtext\"(\"in grey: edge numbers\", side=1, line=-1); # hide\nplot(net, showedgelength=true, useedgelength=true, tipoffset=0.1);\nR\"mtext\"(\"in black: edge lengths\", side=1, line=-1); # hide\nR\"dev.off()\" # hide\nnothing # hide","category":"page"},{"location":"man/getting_started/","page":"getting started","title":"getting started","text":"(Image: 3-taxon network)","category":"page"},{"location":"man/getting_started/","page":"getting started","title":"getting started","text":"Note that this example network is not time consistent: the length of the path from the root to the hybridization node H1 is different depending if we go through the major edge (0.6+0.7=1.3) or the minor edge (0.5+0.6=1.1).","category":"page"},{"location":"man/getting_started/","page":"getting started","title":"getting started","text":"Coalescent simulations can be performed along such networks, also along non-ultrametric networks. If the network is ultrametric (time-consistent, and with all tips at the same distance from the root), then gene trees will also be ultrametric.","category":"page"},{"location":"man/getting_started/","page":"getting started","title":"getting started","text":"note: rooting\nThe input network needs to be rooted correctly, as some properties of gene trees may depend on the root. If a network was estimated using a method that estimates a semidirected network (where the root placement is arbitrary) then make sure to root this network appropriately, e.g. using an outgroup, before simulating genes along this network.","category":"page"},{"location":"man/getting_started/#basic-example:-simulate,-save-to-file,-plot","page":"getting started","title":"basic example: simulate, save to file, plot","text":"","category":"section"},{"location":"man/getting_started/","page":"getting started","title":"getting started","text":"We use simulatecoalescent to simulate gene trees along this network. Below, we simulate 2 gene trees, with 1 individual per species.","category":"page"},{"location":"man/getting_started/","page":"getting started","title":"getting started","text":"trees = simulatecoalescent(net, 2, 1)","category":"page"},{"location":"man/getting_started/","page":"getting started","title":"getting started","text":"Branch lengths are assumed to be in coalescent units in the species network (number of generations / effective population size), and edge lengths in gene trees are also in coalescent units.","category":"page"},{"location":"man/getting_started/","page":"getting started","title":"getting started","text":"We can work with these gene trees within Julia with downstream code, and/or we can save them to a file:","category":"page"},{"location":"man/getting_started/","page":"getting started","title":"getting started","text":"writeMultiTopology(trees, stdout) # write them to standout output (screen here)\nwriteMultiTopology(trees, \"genetrees.phy\") # warning: will overwrite \"genetrees.phy\" if this file existed\nrm(\"genetrees.phy\") # hide","category":"page"},{"location":"man/getting_started/","page":"getting started","title":"getting started","text":"Let's plot these 2 gene trees. In the plots below, we annotate each edge with its attribute that tracked the network edge on which the coalescent event occured (where the gene tree lineage originated, going back in time). For example, the gene lineage that ends in A is always mapped to network edge 6, which is the number of the external edge to A in the network (see plot of network above on the left).","category":"page"},{"location":"man/getting_started/","page":"getting started","title":"getting started","text":"R\"svg\"(figname(\"genetrees_gettingstarted_1.svg\"), width=6, height=3); # hide\nR\"par\"(mar=[.1,.1,.1,.1]); R\"layout\"([1 2]); # hide\nusing DataFrames\nfor i in 1:2\n gt = trees[i]\n plot(gt, tipoffset=0.1,\n edgelabel=DataFrame(number = [e.number for e in gt.edge],\n label = [e.inCycle for e in gt.edge]));\n R\"mtext\"(\"gene $i\", line=-1) # hide\nend\nR\"mtext\"(\"numbers: network edge each gene lineage maps to, at time of coalescence.\\n8 = number of edge above the network root\", side=1, line=-1, outer=true); # hide\nR\"dev.off()\" # hide\nnothing # hide","category":"page"},{"location":"man/getting_started/","page":"getting started","title":"getting started","text":"(Image: 3-taxon network)","category":"page"},{"location":"man/getting_started/#several-individuals-per-species","page":"getting started","title":"several individuals per species","text":"","category":"section"},{"location":"man/getting_started/","page":"getting started","title":"getting started","text":"We can ask for more individuals. To simulate 3 individuals / species for example:","category":"page"},{"location":"man/getting_started/","page":"getting started","title":"getting started","text":"simulatecoalescent(net, 1, 3) # 1 gene tree only. 3 individuals in each species","category":"page"},{"location":"man/getting_started/","page":"getting started","title":"getting started","text":"We can also ask for varying numbers of individuals. For example, we simulate below 2 individuals in A and 1 individual in each of B and C, using a dictionary to map species to their number of individuals:","category":"page"},{"location":"man/getting_started/","page":"getting started","title":"getting started","text":"genetrees = simulatecoalescent(net, 1, Dict(\"A\"=>2, \"B\"=>1, \"C\"=>1));\nwriteTopology(genetrees[1])","category":"page"},{"location":"man/getting_started/","page":"getting started","title":"getting started","text":"We can set 0 individuals within a species to simulate missing data.","category":"page"},{"location":"man/getting_started/","page":"getting started","title":"getting started","text":"genetrees = simulatecoalescent(net, 3, Dict(\"A\"=>2, \"B\"=>1, \"C\"=>0));\nwriteMultiTopology(genetrees, stdout)","category":"page"},{"location":"lib/public/#public-documentation","page":"public","title":"public documentation","text":"","category":"section"},{"location":"lib/public/","page":"public","title":"public","text":"Documentation for PhyloCoalSimulations's public (exported) functions. Most functions are internal (not exported).","category":"page"},{"location":"lib/public/#functions-and-types","page":"public","title":"functions & types","text":"","category":"section"},{"location":"lib/public/","page":"public","title":"public","text":"Modules = [PhyloCoalSimulations]\nPrivate = false\nOrder = [:function,:type]","category":"page"},{"location":"lib/public/#PhyloCoalSimulations.simulatecoalescent-Tuple{PhyloNetworks.HybridNetwork, Vararg{Any}}","page":"public","title":"PhyloCoalSimulations.simulatecoalescent","text":"simulatecoalescent([rng::AbstractRNG,] net, nloci, nindividuals;\n nodemapping=false, inheritancecorrelation=0.0)\n\nSimulate nloci gene trees with nindividuals from each species under the multispecies network coalescent, along network net whose branch lengths are assumed to be in coalescent units (ratio: number of generations / effective population size). The coalescent model uses the infinite-population-size approximation.\n\nThe random number generator rng is optional.\n\nOutput: vector of gene trees, of length nloci.\n\nnindividuals can be a single integer, or a dictionary listing the number of individuals to be simulated for each species.\n\nIf nodemapping is true, each simulated gene tree is augmented with degree-2 nodes that can be mapped to speciation or hybridization events. The mapping of gene tree nodes & edges to network edges is carried by the .inCycle attribute. The mapping of gene tree nodes to network nodes is carried by the .name attribute. Namely:\n\nA degree-3 node (1 parent + 2 children) represents a coalescent event that occurred along a population edge in net. Its .inCycle attribute is set to the number of that network population edge. Its parent edge has its .inCycle attribute also set to the number of the population edge that it originated from.\nThe gene tree's root node (of degree 2) represents a coalescent event along the network's root edge. Its .inCycle attribute is the number assigned to the network's root edge, which is set by get_rootedgenumber as the maximum edge number + 1.\nA leaf (or degree-1 node) represents an individual. It maps to a species in net. The individual leaf name is set to the species name if nindividuals is 1. Otherwise, its name is set to speciesname_i where i is the individual number in that species. Its inCycle attribute is the default -1.\nA non-root degree-2 node represents a speciation or hybridization and maps to a population node in net. Its inCycle attribute is the default -1. Its name is set to network node name, if it exists. If the network node has no name, the gene tree node is given a name built from the network node number.\n\nBy default, lineages at a hybrid node come from a parent (chosen according to inheritance probabilities γ) independently across lineages. Positive dependence can be simulated with option inheritancecorrelation. For example, if this correlation is set to 1, then all lineages inherit from the same (randomly sampled) parent. More generally, the lineages' parents are simulated according to a Dirichlet process with base distribution determined by the γ values, and with concentration parameter α = (1-r)/r, that is, r = 1/(1+α), where r is the input inheritance correlation. For more details about this model, please read the package manual or refer to Fogg, Allman & Ané (2023).\n\nAssumptions:\n\nnet must have non-missing edge lengths and γ values.\nIf nindividuals is a dictionary, it must have a key for all species, with the same spelling of species names in its keys as in the tip labels of net.\n\nexamples\n\njulia> using PhyloNetworks\n\njulia> net = readTopology(\"(A:1,B:1);\"); # branch lengths of 1 coalescent unit\n\njulia> using Random; Random.seed!(54321); # for replicability of examples below\n\njulia> simulatecoalescent(net, 2, 1) # 2 gene trees, 1 individual/species\n2-element Vector{HybridNetwork}:\n PhyloNetworks.HybridNetwork, Rooted Network\n2 edges\n3 nodes: 2 tips, 0 hybrid nodes, 1 internal tree nodes.\ntip labels: B, A\n(B:1.023,A:1.023);\n\n PhyloNetworks.HybridNetwork, Rooted Network\n2 edges\n3 nodes: 2 tips, 0 hybrid nodes, 1 internal tree nodes.\ntip labels: B, A\n(B:2.328,A:2.328);\n\n\njulia> simulatecoalescent(net, 1, 3)[1] # 1 gene tree, 3 individuals/species\nPhyloNetworks.HybridNetwork, Rooted Network\n10 edges\n11 nodes: 6 tips, 0 hybrid nodes, 5 internal tree nodes.\ntip labels: B_2, B_1, B_3, A_3, ...\n(((B_2:0.12,B_1:0.12):2.39,B_3:2.51):0.692,(A_3:0.518,(A_2:0.461,A_1:0.461):0.057):2.684);\n\n\njulia> simulatecoalescent(net, 1, Dict(\"A\"=>2, \"B\"=>1))[1] # 2 individuals in A, 1 in B\nPhyloNetworks.HybridNetwork, Rooted Network\n4 edges\n5 nodes: 3 tips, 0 hybrid nodes, 2 internal tree nodes.\ntip labels: B, A_2, A_1\n(B:2.801,(A_2:0.344,A_1:0.344):2.457);\n\nIn the next example, we use a custom random number generator (RNG), to show how to do so. In this example we use a stable RNG to make the example reproducible across julia versions. However, this particular RNG is not recommended. The RNG used by default is better (e.g. much more efficient).\n\njulia> # using Pkg; Pkg.add(\"StableRNGs\") # to install StableRNGs if not done earlier\n\njulia> using StableRNGs\n\njulia> rng = StableRNG(791);\n\njulia> tree1 = simulatecoalescent(rng, net,2,2; nodemapping=true)[1]; # first gene tree only\n\njulia> writeTopology(tree1, round=true)\n\"(((B_2:0.82,B_1:0.82):0.18)minus2:0.992,((A_1:1.0)minus2:0.427,(A_2:1.0)minus2:0.427):0.565);\"\n\njulia> PhyloNetworks.nameinternalnodes!(net, \"i\"); writeTopology(net)\n\"(A:1.0,B:1.0)i1;\"\n\njulia> tree1 = simulatecoalescent(rng, net,2,2; nodemapping=true)[1]; writeTopology(tree1, round=true)\n\"((B_2:1.0)i1:0.621,((B_1:1.0)i1:0.018,((A_1:0.61,A_2:0.61):0.39)i1:0.018):0.604);\"\n\njulia> printNodes(net)\nnode leaf hybrid hasHybEdge name inCycle edges'numbers\n1 true false false A -1 1 \n2 true false false B -1 2 \n-2 false false false i1 -1 1 2 \n\njulia> printNodes(tree1)\nnode leaf hybrid hasHybEdge name inCycle edges'numbers\n10 false false false 3 8 9 \n8 false false false i1 -1 5 8 \n5 true false false B_2 -1 5 \n9 false false false 3 7 6 9 \n7 false false false i1 -1 4 7 \n4 true false false B_1 -1 4 \n6 false false false i1 -1 3 6 \n3 false false false 1 1 2 3 \n1 true false false A_1 -1 1 \n2 true false false A_2 -1 2 \n\njulia> [(tree_edge_number = e.number, pop_edge_number = e.inCycle) for e in tree1.edge]\n9-element Vector{@NamedTuple{tree_edge_number::Int64, pop_edge_number::Int64}}:\n (tree_edge_number = 8, pop_edge_number = 3)\n (tree_edge_number = 5, pop_edge_number = 2)\n (tree_edge_number = 9, pop_edge_number = 3)\n (tree_edge_number = 7, pop_edge_number = 3)\n (tree_edge_number = 4, pop_edge_number = 2)\n (tree_edge_number = 6, pop_edge_number = 3)\n (tree_edge_number = 3, pop_edge_number = 1)\n (tree_edge_number = 1, pop_edge_number = 1)\n (tree_edge_number = 2, pop_edge_number = 1)\n\n\n\n\n\n","category":"method"},{"location":"lib/public/#PhyloCoalSimulations.simulatecoalescent-Tuple{Random.AbstractRNG, PhyloNetworks.HybridNetwork, Integer, Any, Any}","page":"public","title":"PhyloCoalSimulations.simulatecoalescent","text":"simulatecoalescent([rng::AbstractRNG,] net, nloci, nindividuals, populationsize;\n nodemapping=false, round_generationnumber=true,\n inheritancecorrelation=0.0)\n\nSimulate nloci gene trees with nindividuals from each species under the multispecies network coalescent, along network net, whose branch lengths are assumed to be in number of generations. populationsize should be a single number, assumed to be the (haploid) effective population size Nₑ, constant across the species phylogeny. Alternatively, populationsize can be a dictionary mapping the number of each edge in net to its Nₑ, including an extra edge number for the population above the root of the network.\n\nCoalescent units are then calculated as u=g/Nₑ where g is the edge length in net (number of generations), and the coalescent model is applied using the infinite-population-size approximation.\n\nOutput: vector of gene trees with edge lengths in number of generations, calculated as g=uNₑ and then rounded to be an integer, unless round_generationnumber is false.\n\nwarning: Warning\nWhen populationsize Nₑ is not provided as input, all edge lengths are in coalescent units. When populationsize is given as an argument, all edge lengths are in number of generations. The second method (using # generation and Nₑ as input) is a wrapper around the first (using coalescent units).\n\njulia> using PhyloNetworks\n\njulia> net = readTopology(\"(A:500,B:500);\"); # branch lengths of 100 generations\n\njulia> Ne = Dict(e.number => 1_000 for e in net.edge);\n\njulia> rootedgenumber = PhyloCoalSimulations.get_rootedgenumber(net)\n3\n\njulia> push!(Ne, rootedgenumber => 2_000) # Ne for population above the root\nDict{Int64, Int64} with 3 entries:\n 2 => 1000\n 3 => 2000\n 1 => 1000\n\njulia> using Random; Random.seed!(54321); # for replicability of example below\n\njulia> genetrees = simulatecoalescent(net, 2, 1, Ne);\n\njulia> writeMultiTopology(genetrees, stdout) # branch lengths: number of generations\n(B:546.0,A:546.0);\n(B:3155.0,A:3155.0);\n\n\n\n\n\n\n","category":"method"},{"location":"lib/public/#index","page":"public","title":"index","text":"","category":"section"},{"location":"lib/public/","page":"public","title":"public","text":"Pages = [\"public.md\"]","category":"page"},{"location":"","page":"home","title":"home","text":"CurrentModule = PhyloCoalSimulations","category":"page"},{"location":"#PhyloCoalSimulations","page":"home","title":"PhyloCoalSimulations","text":"","category":"section"},{"location":"","page":"home","title":"home","text":"PhyloCoalSimulations is a Julia package to simulate phylogenies under the coalescent. It depends on PhyloNetworks for the phylogenetic data structures, and manipulation of phylogenies.","category":"page"},{"location":"","page":"home","title":"home","text":"References: please see bibtex entries here.","category":"page"},{"location":"","page":"home","title":"home","text":"for this package and the network coalescent model with inheritance correlation:\nFogg, Allman & Ané (2023). PhyloCoalSimulations: A simulator for network multispecies coalescent models, including a new extension for the inheritance of gene flow. Systematic Biology, 72(5):1171–1179. doi:10.1093/sysbio/syad030.\nfor PhyloNetworks, which this package depends on:\nSolís-Lemus, Bastide & Ané (2017). PhyloNetworks: a package for phylogenetic networks. Molecular Biology and Evolution, 34(12):3292–3298. doi:10.1093/molbev/msx235","category":"page"},{"location":"","page":"home","title":"home","text":"For a tutorial, see the manual:","category":"page"},{"location":"","page":"home","title":"home","text":"Pages = [\n \"man/getting_started.md\",\n \"man/mapping_genetree_to_network.md\",\n \"man/converting_coal2generation_units.md\",\n \"man/correlated_inheritance.md\",\n \"man/more_examples.md\",\n]\nDepth = 3","category":"page"},{"location":"","page":"home","title":"home","text":"For help on individual functions, see the library:","category":"page"},{"location":"","page":"home","title":"home","text":"Pages = [\n \"lib/public.md\",\n \"lib/internal.md\",\n]\nDepth = 3","category":"page"},{"location":"lib/internal/#internal-documentation","page":"internals","title":"internal documentation","text":"","category":"section"},{"location":"lib/internal/","page":"internals","title":"internals","text":"Documentation for PhyloCoalSimulations's internal functions. These functions are not exported and their access (API) should not be considered stable. But they can still be used, like this for example: PhyloCoalSimulations.foo() for a function named foo().","category":"page"},{"location":"lib/internal/#functions-and-types","page":"internals","title":"functions & types","text":"","category":"section"},{"location":"lib/internal/","page":"internals","title":"internals","text":"Modules = [PhyloCoalSimulations]\nPublic = false\nOrder = [:type,:function]","category":"page"},{"location":"lib/internal/#PhyloCoalSimulations.mappingnodes","page":"internals","title":"PhyloCoalSimulations.mappingnodes","text":"mappingnodes(gene tree)\n\nType to define an iterator over degree-2 mapping nodes in a gene tree, assuming these degree-2 nodes (other than the root) have a name to map them to nodes in a species phylogeny. See ismappingnode.\n\n\n\n\n\n","category":"type"},{"location":"lib/internal/#PhyloCoalSimulations.coalescence_edge-NTuple{4, Any}","page":"internals","title":"PhyloCoalSimulations.coalescence_edge","text":"coalescence_edge(edge1, edge2, number, populationid)\n\nCreate a coalescence between edges 1 and 2: with a new parent node n numbered number and a new parent edge e above the parent node, of length 0 and numbered number. Both n.inCycle and e.inCycle are set to populationid.\n\n\n\n\n\n","category":"method"},{"location":"lib/internal/#PhyloCoalSimulations.convert2tree!-Tuple{PhyloNetworks.Node}","page":"internals","title":"PhyloCoalSimulations.convert2tree!","text":"convert2tree!(rootnode)\n\nReturn a network with all nodes and edges that can be reached from rootnode. Warning: Assumes that edges are correctly directed (with correct isChild1 attribute) and that the graph is a tree. This is not checked.\n\nIf the root node is still attached to an incomplete root edge, this edge & node are first disconnected.\n\n\n\n\n\n","category":"method"},{"location":"lib/internal/#PhyloCoalSimulations.get_rootedgenumber-Tuple{Any}","page":"internals","title":"PhyloCoalSimulations.get_rootedgenumber","text":"get_rootedgenumber(network)\n\n1 + maximum of 0 and of all the network's edge numbers: can be used as a unique identifier of the edge above the network's root.\n\n\n\n\n\n","category":"method"},{"location":"lib/internal/#PhyloCoalSimulations.initializetip","page":"internals","title":"PhyloCoalSimulations.initializetip","text":"initializetip(species::AbstractString, individual::AbstractString,\n number::Integer, delim=\"\"::AbstractString)\n\nCreate a leaf node and a pendant edge of length 0, incident to each other, both numbered number. Return the pendant edge. The leaf name is made by concatenating species, delim and individual.\n\n\n\n\n\n","category":"function"},{"location":"lib/internal/#PhyloCoalSimulations.initializetipforest","page":"internals","title":"PhyloCoalSimulations.initializetipforest","text":"initializetipforest(speciesnode::Node, nindividuals::Integer,\n number::Integer, delim)\n\nVector of pendant leaf edges, with leaves named after speciesnode, and numbered with consecutive number IDs starting at number. If nindividuals is 1, then the leaf name is simply the species name. Otherwise, then the leaf names include the individual number and the default delimiter is _. For example, if the species name is s then leaf names are: s_1, s_2, etc. by default. Pendant leaf edges have inCycle set to the number of the corresponding edge in the species network.\n\n\n\n\n\n","category":"function"},{"location":"lib/internal/#PhyloCoalSimulations.ismappingnode-Tuple{PhyloNetworks.Node}","page":"internals","title":"PhyloCoalSimulations.ismappingnode","text":"ismappingnode(node)\n\nBoolean: true if node is of degree 2, has a single child, and has a name. (The root is of degree-2 but is not a mapping node).\n\n\n\n\n\n","category":"method"},{"location":"lib/internal/#PhyloCoalSimulations.map2population!-NTuple{4, Any}","page":"internals","title":"PhyloCoalSimulations.map2population!","text":"map2population!(forest, population_node, populationid, nextlineageID)\n\nExtend each incomplete edge in the forest with a new degree-2 node n and a new incomplete edge e, with the following information to map n and e into the species phylogeny:\n\ne.inCycle is set to populationid, and\nn.name is set to population_node.name if this name is non-empty, or string(population_node.number) otherwise (with any negative sign replaced by the string \"minus\").\n\ne.number and n.number are set to nextlineageID, which is incremented by 1 for each incomplete edge in the forest.\n\nThe forest is updated to contain the newly-created incomplete edges, replacing the old incomplete (and now complete) edges.\n\nOutput: nextlineageID, incremented by the number of newly created degree-2 lineages.\n\nexample\n\njulia> using PhyloNetworks; net = readTopology(\"(A:1,B:1);\");\n\njulia> leafA = net.node[1]; edge2A_number = net.edge[1].number;\n\njulia> f = PhyloCoalSimulations.initializetipforest(leafA, 2, 4); # 2 edges, numbered 4 & 5\n\njulia> PhyloCoalSimulations.map2population!(f, leafA, edge2A_number, 6)\n8\n\njulia> length(f)\n2\n\njulia> f[2]\nPhyloNetworks.EdgeT{PhyloNetworks.Node}:\n number:7\n length:0.0\n attached to 1 node(s) (parent first): 7\n\njulia> [e.node[1].name for e in f]\n2-element Vector{String}:\n \"A\"\n \"A\"\n\n\n\n\n\n","category":"method"},{"location":"lib/internal/#PhyloCoalSimulations.population_mappedto-Tuple{Union{PhyloNetworks.EdgeT{PhyloNetworks.Node}, PhyloNetworks.Node}}","page":"internals","title":"PhyloCoalSimulations.population_mappedto","text":"population_mappedto(edge or node)\n\nIdentifier of the population (edge in the species network) that a gene tree's edge or a node is mapped onto, or nothing if not mapped. For example, coalescent nodes in gene trees map to a node in the species phylogeny, instead of mapping to an edge.\n\n\n\n\n\n","category":"method"},{"location":"lib/internal/#PhyloCoalSimulations.simulatecoal_onepopulation!-Tuple{AbstractVector, Vararg{Any}}","page":"internals","title":"PhyloCoalSimulations.simulatecoal_onepopulation!","text":"simulatecoal_onepopulation!([rng::AbstractRNG,]\n lineagelist,\n population_length,\n nextlineageID,\n populationID=-1\n)\n\nSimulate the coalescent process within a single population of length given in coalescent units, starting from lineages in lineagelist. This list should be a vector of incomplete edges, that is, edges incident to a single node only.\n\nVector of incomplete edges, whose lengths have been increased, is modified in place. New nodes and their parent edges are created by coalescent events, numbered with consecutive integers starting at nextlineageID.\n\nThe random number generator rng is optional.\n\nOutput: nextlineageID, incremented by number of new lineages.\n\nIn lineages, edge lengths are also considered in coalescent units.\n\nThe newly created nodes and edges have their .inCycle attribute set to populationID, so as to track the mapping of gene lineages to populations in the species phylogeny.\n\nexamples\n\njulia> PhyloCoalSimulations.simulatecoal_onepopulation!([], 2.0, 1)\n1\n\njulia> e1 = PhyloCoalSimulations.initializetip(\"s\",\"1\",1,\"\");\n\njulia> e2 = PhyloCoalSimulations.initializetip(\"s\",\"2\",2,\"\");\n\njulia> forest = [e1,e2];\n\njulia> using Random; Random.seed!(7690);\n\njulia> PhyloCoalSimulations.simulatecoal_onepopulation!(forest, Inf, 3);\n\njulia> PhyloCoalSimulations.convert2tree!(forest[1].node[1])\nPhyloNetworks.HybridNetwork, Rooted Network\n2 edges\n3 nodes: 2 tips, 0 hybrid nodes, 1 internal tree nodes.\ntip labels: s2, s1\n(s2:0.302,s1:0.302);\n\n\n\n\n\n\n","category":"method"},{"location":"lib/internal/#index","page":"internals","title":"index","text":"","category":"section"},{"location":"lib/internal/","page":"internals","title":"internals","text":"Pages = [\"internal.md\"]","category":"page"}]
}