# ---
# jupyter:
#   jupytext:
#     text_representation:
#       extension: .py
#       format_name: light
#       format_version: '1.5'
#       jupytext_version: 1.14.4
#   kernelspec:
#     display_name: Python 2
#     language: python
#     name: python2
# ---

# + [markdown] slideshow={"slide_type": "slide"}
# # The ipyrad.analysis tool kit
# <NAME>  
#

# + [markdown] slideshow={"slide_type": "slide"}
# ### Install software
# All required software for this walkthrough is available on conda. 

# + slideshow={"slide_type": "-"}
# conda install -c ipyrad ipyrad structure clumpp bpp 
# conda install -c eaton-lab toytree toyplot
# conda install -c bioconda raxml

# + [markdown] slideshow={"slide_type": "slide"}
# ### Start an ipyparallel cluster
# In a separate terminal run the following command to start a cluster of engines. If working on a notebook running remotely, use the dashboard to open a new terminal. 

# +
# ipcluster start --n=4
# -

# You should then be able to connect to the engines in your notebook:

# +
## connect to the cluster
import ipyparallel as ipp
ipyclient = ipp.Client()

## print number of engines
print len(ipyclient), "connected engines"

# + [markdown] slideshow={"slide_type": "slide"}
# ### Assemble a RAD data set
# The code here is to assemble the example empirical data set from the [ipyrad tutotial](http://ipyrad.readthedocs.io/userguide.html#empirical-examples). 
# -

## import ipyrad
import ipyrad as ip

# Minimal workflow: scroll down for details.

# + slideshow={"slide_type": "-"}
## create an Assembly object
data = ip.Assembly("simdata")

## set I/O paths for the data
data.set_params("project_dir", "~/workshop")
data.set_params("raw_fastq_path", "ipsimdata/rad_example_R1_.fastq.gz")
data.set_params("barcodes_path", "ipsimdata/rad_example_barcodes.txt")

## run all steps of the Assembly
data.run("1234567")

# + [markdown] slideshow={"slide_type": "subslide"}
# ### Modify more parameters

# + slideshow={"slide_type": "-"}
## set params
data.set_params("filter_adapters", 2)
data.set_params("output_formats", "lpask")

## show params
data.get_params()

# + [markdown] slideshow={"slide_type": "subslide"}
# ### Assemble the data set
# You can run one or more steps just like in the CLI. 
# -

## run all steps of assembly
data.run("1234567")

# + [markdown] slideshow={"slide_type": "subslide"}
# ### Access assembly results
# You can easily access summary stats for the assembly as a data frame. 
# -

## summary stats
data.stats

# + [markdown] slideshow={"slide_type": "subslide"}
# ### Plot statistics

# + slideshow={"slide_type": "-"}
import toyplot

## plot barplot
c, a, m = toyplot.bars(
    data.stats.hetero_est, 
    height=250, width=500,
    )

## style the axes
a.x.ticks.locator = toyplot.locator.Explicit(
    locations=range(len(data.stats)),
    labels=data.stats.index)
a.y.label.text = "Heterozygosity"
a.y.ticks.show = True

# + [markdown] slideshow={"slide_type": "subslide"}
# ### Access result files  
# You can also access the stats files for each step, and the output files for downstream analyses. 

# +
## s2 stats file
print data.stats_files.s2

## the .loci file location
print data.outfiles.loci

# + [markdown] slideshow={"slide_type": "slide"}
# ## ipyrad.analysis tools
# The ipyrad.analysis module includes many *wrapper* tools that can be used to efficiently run evolutionary analysis tools in a notebook. 
# -

## import the toolkit
import ipyrad.analysis as ipa

# + [markdown] slideshow={"slide_type": "slide"}
# ### RAxML analysis
# Simply enter the location of the phylip file, which can be accessed from the `.outfiles` attribute of the Assembly object. You can also provide a name and output directory, and set many other optional parameters. 

# + slideshow={"slide_type": "skip"}
import ipyrad as ip
import ipyparallel as ipp

data = ip.load_json("/home/deren/workshop/simdata.json")
ipyclient = ipp.Client()

# -

# Minimal workflow: scroll down for details.

# + slideshow={"slide_type": "-"}
## create a raxml object
s = ipa.raxml(
        name=data.name,
        phyfile=data.outfiles.phy, 
        workdir="~/workshop/analysis-raxml");

## run the analysis
s.run()

# + [markdown] slideshow={"slide_type": "subslide"}
# ### Modify parameters and other functions
# -

## modify params
s.params.T = 4
s.params.N = 100

# + slideshow={"slide_type": "fragment"}
## print the raxml command as a string
print s.command

# + slideshow={"slide_type": "fragment"}
## overwrite existing result with this 'name'
s.run(force=True)

# + [markdown] slideshow={"slide_type": "subslide"}
# ### Access the tree files and plot
# -

print s.trees

# + slideshow={"slide_type": "fragment"}
import toytree
tre = toytree.tree(s.trees.bipartitions)
tre.root(wildcard='3')
tre.draw(
    width=300,
    node_labels=tre.get_node_values("support"),
    node_size=20,
    );

# + [markdown] slideshow={"slide_type": "slide"}
# ### introgression (abba-baba) analysis
# The baba object can be used to set up abba-baba tests, to calculate results, and to generate plots to visualize them. 
# -

# Minimal example, scroll down for details. 

# +
## create a baba object
b = ipa.baba(data=data.outfiles.loci)

## generate tests given the rooted tree
b.tests = [
    {"p4":["3L_0"], 
     "p3":["2F_0"], 
     "p2":["1D_0"],
     "p1":["1A_0"]}]

## run jobs distributed across the cluster
b.run(ipyclient)
b.results_table

# + [markdown] slideshow={"slide_type": "subslide"}
# ### Auto-generate tests 
# Instead of writing out many tests explicitly, you can instead enter a rooted tree to the baba object and use this function to auto-generate four-taxon test fitting the tree and constraints. 

# + slideshow={"slide_type": "-"}
## init baba object
b = ipa.baba(data=data.outfiles.loci, newick=tre)

## generate all possible tests on this tree
b.generate_tests_from_tree()

## set constraints on tests
cdict = {"p4": ["3L_0"], 
         "p3": ["2E_0", "2F_0"], 
         "p2": ["1D_0"]}

## generate constrainted number of tests
b.generate_tests_from_tree(
    constraint_dict=cdict,
    constraint_exact=False,
    )

# + [markdown] slideshow={"slide_type": "subslide"}
# ### Run all tests linked to a baba object

# + slideshow={"slide_type": "-"}
## run the tests (in this case 4) linked to the baba object
b.run(ipyclient)

## show results table
b.results_table

# + [markdown] slideshow={"slide_type": "subslide"}
# ### Plot results 
# -

b.plot(
    height=350,
    pct_tree_x = 0.4,
    pct_tree_y = 0.2,
    );

# + slideshow={"slide_type": "subslide"}
### Save the plot 
import toyplot.pdf
canvas, axes, mark = b.plot(height=350, pct_tree_x=0.4, pct_tree_y=0.2)
toyplot.pdf.render(canvas, "/home/deren/workshop/abba-baba.pdf")

## save the results table
b.results_table.to_csv("~/workshop/abba-baba.csv", sep="\t")

# + [markdown] slideshow={"slide_type": "slide"}
# ### Species tree inference by phylogenetic invariants
# The program tetrad follows the algorithm of SVDquartets by inferring all possible quartet trees from a large SNP alignment and uses the program quartet maxcut (Snir et al. 2012) to infer a species tree by quartet joining. 

# +
## create a tetrad class object
tet = ipa.tetrad(
    name=data.name,
    seqfile=data.outfiles.snpsphy,
    mapfile=data.outfiles.snpsmap,
    workdir="~/workshop/analysis-tetrad",
    nboots=100
    )

## run the analysis
tet.run(ipyclient)
# -

# ### Access result tetrad trees and draw

# + slideshow={"slide_type": "subslide"}
tet.trees

# + slideshow={"slide_type": "subslide"}
## load unrooted result tree with toytree and draw
tre = toytree.tree(tet.trees.cons)
tre.draw(
    node_labels=tre.get_node_values("support"),
    node_size=20,
    );
# -

# ### Infer a species tree with BPP

# +
# conda install bpp -c ipyrad

# +
## setup: define how samples group into 'species'
IMAP = {
    "1": ["1A_0", "1B_0", "1C_0"],
    "D": ["1D_0"],
    "2": ["2F_0", "2E_0", "2G_0"],
    "H": ["2H_0"],
    "3": ["3J_0", "3I_0", "3K_0"],
    "L": ["3L_0"],    
}

## setup: define a guidetree
GUIDE = "(((1,D),(2,H)),(3,L));"

## init a bpp object 
bpp = ipa.bpp(
    locifile=data.outfiles.loci, 
    imap=IMAP, 
    guidetree=GUIDE,
    workdir="~/workshop/analysis-bpp"
    );

## submit jobs to run on the cluster
bpp.submit_bpp_jobs("A00", nreps=2, ipyclient=ipyclient)

# + [markdown] slideshow={"slide_type": "subslide"}
# ### Set parameters and filters
# You can define all of the parameter settings that will be used in the BPP .ctl file by modifying the .params attributes. Similarly, you can modify which loci will be included in the analysis using the .filters attributes. 

# +
## set some parameters 
bpp.params.burnin = 1000
bpp.params.nsample = 5000
bpp.params.infer_sptree = 1
bpp.params.infer_delimit = 0 

## set some filters
bpp.filters.maxloci = 200
bpp.filters.minsnps = 2

## submit jobs to run on the cluster
bpp.submit_bpp_jobs("A00", nreps=2, ipyclient=ipyclient)

# + [markdown] slideshow={"slide_type": "subslide"}
# ### Track running jobs
# Unlike some of the other ipyrad.analysis tools, the bpp object does not "block" while the jobs are running. Meaning that after it sends jobs to run on the cluster you can continue to interact with the notebook. This is useful since BPP is not multi-threaded, so you will likely want to submit many different types of jobs. You can check on running jobs like below. 

# +
## a list of submitted jobs
print bpp.asyncs

## a list of result files produced by jobs
print bpp.files

# + [markdown] slideshow={"slide_type": "slide"}
# ### Structure analyses

# + slideshow={"slide_type": "skip"}
import ipyrad as ip
import ipyrad.analysis as ipa
import ipyparallel as ipp

data = ip.load_json("/home/deren/workshop/simdata.json")
ipyclient = ipp.Client()

# +
# conda install structure -c ipyrad
# conda install clumpp -c ipyrad

# +
## create a structure class object
s = ipa.structure(
    name=data.name,
    strfile=data.outfiles.str,
    mapfile=data.outfiles.snpsmap,
    workdir="~/workshop/analysis-structure",    
    );

s.mainparams.burnin = 100
s.mainparams.numreps = 1000

## submit jobs to run on the cluster
for kpop in [2, 3, 4, 5]:
    s.submit_structure_jobs(kpop=kpop, nreps=5, ipyclient=ipyclient)

# + [markdown] slideshow={"slide_type": "subslide"}
# ### Modify parameters settings
# -

s.mainparams.burnin = 10000
s.mainparams.numreps = 100000
s.extraparams.usepopinfo = 0

# ### Summarize results with CLUMPP

# +
## get results for a single K value
s.get_clumpp_table(3)

## make a dict for all results
tables = {}
for kpop in [2, 3, 4, 5]:
    tables[kpop] = s.get_clumpp_table(kpop)
# -



StarcoderdataJupyterScriptsDedupFiltered