<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 <span style="color
Step1:  Pre-filter data from other programs (e.g., FreeBayes, GATK)
You can use the program bcftools to pre-filter your data to exclude indels and low quality SNPs. If you ran the conda install commands above then you will have all of the required tools installed. To achieve the format that ipyrad expects you will need to exclude indel containing SNPs (this may change in the future). Further quality filtering is optional. 
The example below reduced the size of a VCF data file from 29Gb to 80Mb! VCF contains a lot of information that you do not need to retain through all of your analyses. We will keep only the final genotype calls.  
Note that the code below is bash script. You can run this from a terminal, or in a jupyter notebook by appending the (%%bash) header like below.
Step2:  A peek at the cleaned VCF file
Step3:  Converting clean VCF to HDF5
Here I using a VCF file from whole geome data for 20 monkey's from an unpublished study (in progress). It contains >6M SNPs all from chromosome 1. Because many SNPs are close together and thus tightly linked we will likely wish to take linkage into account in our downstream analyses.
The ipyrad analysis tools can do this by encoding linkage block information into the HDF5 file. Here we encode ld_block_size of 20K bp. This breaks the 1 scaffold (chromosome) into about 10K linkage blocks. See the example below of this information being used in an ipyrad PCA analysis.
Step4:  Downstream analyses
The data file now contains 6M SNPs across 20 samples and N linkage blocks. By default the PCA tool subsamples a single SNP per linkage block. To explore variation over multiple random subsamplings we can use the nreplicates argument.
Step5:  Run a single PCA analysis from subsampled unlinked SNPs
Step6:  Run multiple PCAs over replicates of subsampled SNPs
Here you can see the results for a different 10K SNPs that are sampled in each replicate iteration. If the signal in the data is robust then we should expect to see the points clustering at a similar place across replicates. Internally ipyrad will rotate axes to ensure the replicate plots align despite axes swapping (which is arbitrary in PCA space). You can see this provides a better view of uncertainty in our estimates than the plot above (and it looks cool!)
<ASSISTANT_TASK:>
Python Code:
# conda install ipyrad -c bioconda 
# conda install htslib -c bioconda
# conda install bcftools -c bioconda
import ipyrad.analysis as ipa
import pandas as pd
Explanation: <span style="color:gray">ipyrad-analysis toolkit:</span> vcf_to_hdf5
View as notebook
Many genome assembly tools will write variant SNP calls to the VCF format (variant call format). This is a plain text file that stores variant calls relative to a reference genome in tabular format. It includes a lot of additional information about the quality of SNP calls, etc., but is not very easy to read or efficient to parse. To make analyses run a bit faster ipyrad uses a simplified format to store this information in the form of an HDF5 database. You can easily convert any VCF file to this HDF5 format using the ipa.vcf_to_hdf5() tool. 
This tool includes an added benefit of allowing you to enter an (optional) ld_block_size argument when creating the file which will store information that can be used downstream by many other tools to subsample SNPs and perform bootstrap resampling in a way that reduces the effects of linkage among SNPs. If your data are assembled RAD data then the ld_block_size is not required, since we can simply use RAD loci as the linkage blocks. But if you want to combine reference-mapped RAD loci located nearby in the genome as being on the same linkage block then you can enter a value such as 50,000 to create 50Kb linkage block that will join many RAD loci together and sample only 1 SNP per block in each bootstrap replicate. If your data are not RAD data, e.g., whole genome data, then the ld_block_size argument will be required in order to  encode linkage information as discrete blocks into your database. 
Required software
If you are converting a VCF file assembled from some other tool (e.g., GATK, freebayes, etc.) then you will need to install the htslib and bcftools software and use them as described below.
End of explanation
%%bash
# compress the VCF file if not already done (creates .vcf.gz)
bgzip data.vcf
# tabix index the compressed VCF (creates .vcf.gz.tbi)
tabix data.vcf.gz
# remove multi-allelic SNPs and INDELs and PIPE to next command
bcftools view -m2 -M2 -i'CIGAR="1X" & QUAL>30' data.vcf.gz -Ou | 
    # remove extra annotations/formatting info and save to new .vcf
    bcftools annotate -x FORMAT,INFO  > data.cleaned.vcf
    
# recompress the final file (create .vcf.gz)
bgzip data.cleaned.vcf
Explanation: Pre-filter data from other programs (e.g., FreeBayes, GATK)
You can use the program bcftools to pre-filter your data to exclude indels and low quality SNPs. If you ran the conda install commands above then you will have all of the required tools installed. To achieve the format that ipyrad expects you will need to exclude indel containing SNPs (this may change in the future). Further quality filtering is optional. 
The example below reduced the size of a VCF data file from 29Gb to 80Mb! VCF contains a lot of information that you do not need to retain through all of your analyses. We will keep only the final genotype calls.  
Note that the code below is bash script. You can run this from a terminal, or in a jupyter notebook by appending the (%%bash) header like below.
End of explanation
# load the VCF as an datafram
dfchunks = pd.read_csv(
    "/home/deren/Documents/ipyrad/sandbox/Macaque-Chr1.clean.vcf.gz",
    sep="\t", 
    skiprows=1000, 
    chunksize=1000,
)
# show first few rows of first dataframe chunk
next(dfchunks).head()
Explanation: A peek at the cleaned VCF file
End of explanation
# init a conversion tool
converter = ipa.vcf_to_hdf5(
    name="Macaque_LD20K",
    data="/home/deren/Documents/ipyrad/sandbox/Macaque-Chr1.clean.vcf.gz",
    ld_block_size=20000,
)
# run the converter
converter.run()
Explanation: Converting clean VCF to HDF5
Here I using a VCF file from whole geome data for 20 monkey's from an unpublished study (in progress). It contains >6M SNPs all from chromosome 1. Because many SNPs are close together and thus tightly linked we will likely wish to take linkage into account in our downstream analyses.
The ipyrad analysis tools can do this by encoding linkage block information into the HDF5 file. Here we encode ld_block_size of 20K bp. This breaks the 1 scaffold (chromosome) into about 10K linkage blocks. See the example below of this information being used in an ipyrad PCA analysis.
End of explanation
# init a PCA tool and filter to allow no missing data
pca = ipa.pca(
    data="./analysis-vcf2hdf5/Macaque_LD20K.snps.hdf5",
    mincov=1.0, 
)
Explanation: Downstream analyses
The data file now contains 6M SNPs across 20 samples and N linkage blocks. By default the PCA tool subsamples a single SNP per linkage block. To explore variation over multiple random subsamplings we can use the nreplicates argument.
End of explanation
pca.run_and_plot_2D(0, 1, seed=123);
Explanation: Run a single PCA analysis from subsampled unlinked SNPs
End of explanation
pca.run_and_plot_2D(0, 1, seed=123, nreplicates=25);
Explanation: Run multiple PCAs over replicates of subsampled SNPs
Here you can see the results for a different 10K SNPs that are sampled in each replicate iteration. If the signal in the data is robust then we should expect to see the points clustering at a similar place across replicates. Internally ipyrad will rotate axes to ensure the replicate plots align despite axes swapping (which is arbitrary in PCA space). You can see this provides a better view of uncertainty in our estimates than the plot above (and it looks cool!)
End of explanation
<END_TASK>