<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 Quick guide to the ipyrad API
Getting Started
Welcome! This tutorial will introduce you to the basics of working with ipyrad to assemble RADseq data.
Note
Step1:  Assembly objects
Assembly objects are used by ipyrad to access data stored on disk and to manipulate it. Each biological sample in a data set is represented as a Sample object, and a set of Samples is stored inside an Assembly object.
The Assembly object has functions to assemble data, and to view or plot the resulting assembled files and statistics. Assembly objects can be copied or merged to allow branching events where different parameters can subsequently be applied to different Assembly objects. Examples of this are shown in the workflow. 
Every analysis begins by creating at least one Assembly object. Here we will name it "data1".
Step2:  The printout tells us that we created the object data1, and also that it found 4 engines on our system that can be used for computation. An engine is simply a CPU. When working on a single machine it will usually be easiest to simply let the Assembly object connect to all available local engines. However, on HPC clusters you may need to modify the controller or the number of engines, as shown below
Step3:  For more information about connecting CPUs for parallelization see ipyparallel setup.
Modifying Assembly parameters
The arguments get_params() and set_params() are used to view and modify parameter settings of an Assembly object, respectively.
Step4:  To get more detailed information about each parameter use ip.get_params_info(), or look up their funcion in the documentation (Parameters). To quickly look up the proper formatting for a parameter,  you can use ip.get_params_info(N), where N is the number of a parameter. Example
Step5:  Sample Objects
Each biological sample in a data set is represented as a Sample object. Sample objects are created during step1() of an analysis, at which time they are linked to an Assembly object. 
When getting started raw data should be in one of two forms
Step6:  Step 1
Step7:  Step 2
Step8:  Branching Assembly objects
Let's imagine at this point that we are interested in clustering our data at two different clustering thresholds. We will try 0.90 and 0.85. First we need to make a copy/branch of the Assembly object. This will inherit the locations of the data linked in the first object, but diverge in any future applications to the object. Thus, the two Assembly objects can share the same working directory, and inherit shared files, but will diverge in creating new files linked to only one or the other. You can view the directories linked to an Assembly object with the .dirs argument, shown below. The prefix_outname (param 14) of the new object is automatically set to the Assembly object name.
Step9:  Step 3
Step10:  Branched Assembly objects
And you can see below that the two Assembly objects are now working with several shared directories (working, fastq, edits) but with different clust directories (clust_0.85 and clust_0.9).
Step11:  Saving stats outputs
Example
Step12:  Example of plotting with ipyrad
There are a a few simple plotting functions in ipyrad useful for visualizing results. These are in the module ipyrad.plotting. Below is an interactive plot for visualizing the distributions of coverages across the 12 samples in the test data set.
Step13:  Step 4
Step14:  Step 5
Step15:  Quick parameter explanations are always on-hand
Step16:  Log history
A common problem at the end of an analysis, or while troubleshooting it, is that you find you've completely forgotten which parameters you used at what point, and when you changed them. Documenting or executing code inside Jupyter notebooks (like the one you're reading right now) is a great way to keep track of all of this. In addition, ipyrad also stores a log history which time stamps all modifications to Assembly objects.
Step17:  Saving Assembly objects
Assembly objects can be saved and loaded so that interactive analyses can be started, stopped, and returned to quite easily. The format of these saved files is a serialized 'dill' object used by Python. Individual Sample objects are saved within Assembly objects. These objects to not contain the actual sequence data, but only link to it, and so are not very large. The information contained includes parameters and the log of Assembly objects, and the statistics and state of Sample objects. Assembly objects are autosaved each time an assembly step function is called, but you can also create your own checkpoints with the save command.
<ASSISTANT_TASK:>
Python Code:
import ipyrad as ip
Explanation: Quick guide to the ipyrad API
Getting Started
Welcome! This tutorial will introduce you to the basics of working with ipyrad to assemble RADseq data.
Note: this tutorial was created in a Jupyter Notebook and assumes that you’re following-along in a notebook of your own. If you installed ipyrad then you will also have jupyter installed by default. For a little background on how jupyter notebooks see Notebooks. All of the code below is written in Python. Follow along in your own notebook. 
To begin, we're going to import the ipyrad module and name it ip so that we can refer to it more easily:
End of explanation
## create an Assembly object named data1. 
data1 = ip.Assembly("data1")
Explanation: Assembly objects
Assembly objects are used by ipyrad to access data stored on disk and to manipulate it. Each biological sample in a data set is represented as a Sample object, and a set of Samples is stored inside an Assembly object.
The Assembly object has functions to assemble data, and to view or plot the resulting assembled files and statistics. Assembly objects can be copied or merged to allow branching events where different parameters can subsequently be applied to different Assembly objects. Examples of this are shown in the workflow. 
Every analysis begins by creating at least one Assembly object. Here we will name it "data1".
End of explanation
## create an Assembly object linked to 8 engines using MPI
data1 = ip.Assembly("data1", N=4, controller="MPI")
Explanation: The printout tells us that we created the object data1, and also that it found 4 engines on our system that can be used for computation. An engine is simply a CPU. When working on a single machine it will usually be easiest to simply let the Assembly object connect to all available local engines. However, on HPC clusters you may need to modify the controller or the number of engines, as shown below:
End of explanation
## setting/modifying parameters for this Assembly object
data1.set_params('project_dir', "./test_rad")
data1.set_params('raw_fastq_path', "./data/sim_rad_test_R1_.fastq.gz")
data1.set_params('barcodes_path', "./data/sim_rad_test_barcodes.txt")
data1.set_params('filter_adapters', 0)
data1.set_params('datatype', 'rad')
## print the parameters for `data`
data1.get_params()
Explanation: For more information about connecting CPUs for parallelization see ipyparallel setup.
Modifying Assembly parameters
The arguments get_params() and set_params() are used to view and modify parameter settings of an Assembly object, respectively.
End of explanation
ip.get_params_info(10)
Explanation: To get more detailed information about each parameter use ip.get_params_info(), or look up their funcion in the documentation (Parameters). To quickly look up the proper formatting for a parameter,  you can use ip.get_params_info(N), where N is the number of a parameter. Example:
End of explanation
## This would link fastq files from the 'sorted_fastq_path' if present
## Here it does nothing b/c there are no files in the sorted_fastq_path
data1.link_fastqs()
Explanation: Sample Objects
Each biological sample in a data set is represented as a Sample object. Sample objects are created during step1() of an analysis, at which time they are linked to an Assembly object. 
When getting started raw data should be in one of two forms: 
+ Non-demultiplexed data files (accompanied by a barcode file).
+ Demultiplexed data files. 
Note: For additional ways to add raw data files to a data set see link_fastqs. 
If the data are already demultiplexed then fastq files can be linked directly to the Assembly object, which in turn will create new Sample objects from them, or link them to existing Sample objects based on the file names (or pair of fastq files for paired data files). The files may be gzip compressed. If the data are not demultiplexed then you will have to run the step1 function below to demultiplex the raw data.
End of explanation
## run step 1 to demultiplex the data
data1.step1()
## print the results for each Sample in data1
print data1.stats
## remove the lane control sequence
#data1.samples.pop("FGXCONTROL")
Explanation: Step 1: Demultiplexing raw data files
Step1 uses barcode information to demultiplex data files found in param 2 ['raw_fastq_path']. It will create a Sample object for each barcoded sample. Below we use the step1() function to demultiplex. The stats attribute of an Assembly object is returned as a pandas data frame.
End of explanation
## example of ways to run step 2 to filter and trim reads
#data1.step2(["1A_0"])              ## run on a single sample
#data1.step2(["1B_0", "1C_0"])      ## run on one or more samples
data1.step2(force=True)             ## run on all samples, overwrite finished
## print the results
print data1.stats
#data1.samples["veitchii"].files
Explanation: Step 2: Filter reads
If for some reason we wanted to execute on just a subsample of our data, we could do this by selecting only certain samples to call the step2 function on. Because step2 is a function of data, it will always execute with the parameters that are linked to data.
End of explanation
## create a copy of our Assembly object
data2 = data1.branch(newname="data2")
## set clustering threshold to 0.90
data2.set_params(11, 0.90)
## look at inherited parameters
data2.get_params()
Explanation: Branching Assembly objects
Let's imagine at this point that we are interested in clustering our data at two different clustering thresholds. We will try 0.90 and 0.85. First we need to make a copy/branch of the Assembly object. This will inherit the locations of the data linked in the first object, but diverge in any future applications to the object. Thus, the two Assembly objects can share the same working directory, and inherit shared files, but will diverge in creating new files linked to only one or the other. You can view the directories linked to an Assembly object with the .dirs argument, shown below. The prefix_outname (param 14) of the new object is automatically set to the Assembly object name.
End of explanation
import ipyrad as ip
data1 = ip.load_assembly("test_rad/data1")
## run step 3 to cluster reads within samples using vsearch
data1.step3(force=True)
## print the results
print data1.stats
## run step 3 to cluster reads in data2 at 0.90 sequence similarity
data2.step3(force=True) 
## print the results
print data2.stats
Explanation: Step 3: clustering within-samples
End of explanation
print "data1 directories:"
for (i,j) in data1.dirs.items():
    print "{}\t{}".format(i, j)
    
print "\ndata2 directories:"
for (i,j) in data2.dirs.items():
    print "{}\t{}".format(i, j)
## TODO, just make a [name]_stats directory in [work] for each data obj
data1.statsfiles
Explanation: Branched Assembly objects
And you can see below that the two Assembly objects are now working with several shared directories (working, fastq, edits) but with different clust directories (clust_0.85 and clust_0.9).
End of explanation
data1.stats.to_csv("data1_results.csv", sep="\t")
data1.stats.to_latex("data1_results.tex")
Explanation: Saving stats outputs
Example: two simple ways to save the stats data frame to a file.
End of explanation
import ipyrad.plotting as iplot
## plot for one or more selected samples
#iplot.depthplot(data1, ["1A_0", "1B_0"])
## plot for all samples in data1
iplot.depthplot(data1)
## save plot as pdf and html
#iplot.depthplot(data1, outprefix="testfig")
Explanation: Example of plotting with ipyrad
There are a a few simple plotting functions in ipyrad useful for visualizing results. These are in the module ipyrad.plotting. Below is an interactive plot for visualizing the distributions of coverages across the 12 samples in the test data set.
End of explanation
## run step 4
data1.step4() 
## print the results
print data1.stats
Explanation: Step 4: Joint estimation of heterozygosity and error rate
End of explanation
#import ipyrad as ip
## reload autosaved data. In case you quit and came back 
#data1 = ip.load_dataobj("test_rad/data1.assembly")
## run step 5
#data1.step5()
## print the results
#print data1.stats
Explanation: Step 5: Consensus base calls
End of explanation
ip.get_params_info(10)
Explanation: Quick parameter explanations are always on-hand
End of explanation
for i in data1.log:
    print i
    
print "\ndata 2 log includes its pre-branching history with data1"
for i in data2.log:
    print i
Explanation: Log history
A common problem at the end of an analysis, or while troubleshooting it, is that you find you've completely forgotten which parameters you used at what point, and when you changed them. Documenting or executing code inside Jupyter notebooks (like the one you're reading right now) is a great way to keep track of all of this. In addition, ipyrad also stores a log history which time stamps all modifications to Assembly objects.
End of explanation
## save assembly object
#ip.save_assembly("data1.p")
## load assembly object
#data = ip.load_assembly("data1.p")
#print data.name
Explanation: Saving Assembly objects
Assembly objects can be saved and loaded so that interactive analyses can be started, stopped, and returned to quite easily. The format of these saved files is a serialized 'dill' object used by Python. Individual Sample objects are saved within Assembly objects. These objects to not contain the actual sequence data, but only link to it, and so are not very large. The information contained includes parameters and the log of Assembly objects, and the statistics and state of Sample objects. Assembly objects are autosaved each time an assembly step function is called, but you can also create your own checkpoints with the save command.
End of explanation
<END_TASK>