<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 Notebook 9
Step1:  Download the sequence data
Sequence data for this study are archived on the NCBI sequence read archive (SRA). Below I read in SraRunTable.txt for this project which contains all of the information we need to download the data. 
Project DRA
Step3:  For each ERS (individuals) get all of the ERR (sequence file accessions).
Step4:  Here we pass the SRR number and the sample name to the wget_download function so that the files are saved with their sample names.
Step5:  Make a params file
Step6:  Note
Step7:  Trimming the barcode
In this data set the data were uploaded separated by sample, but with the barcode still attached to the sequences. The python code below will remove the 5bp barcode from each sequence.
Step8:  Assemble in pyrad
Step9:  Results
We are interested in the relationship between the amount of input (raw) data between any two samples, the average coverage they recover when clustered together, and the phylogenetic distances separating samples.  
Raw data amounts
The average number of raw reads per sample is 1.36M.
Step10:  Look at distributions of coverage
pyrad v.3.0.63 outputs depth information for each sample which I read in here and plot. First let's ask which sample has the highest depth of coverage. The std of coverages is pretty low in this data set compared to several others.
Step11:  Plot the coverage for the sample with highest mean coverage
Green shows the loci that were discarded and orange the loci that were retained. The majority of data were discarded for being too low of coverage.
Step12:  Print final stats table
Step13:  Infer ML phylogeny in raxml as an unrooted tree
Step14:  Plot the tree in R using ape
Step15:  Get phylo distances (GTRgamma dist)
<ASSISTANT_TASK:>
Python Code:
### Notebook 9
### Data set 9 (Ohomopterus)
### Authors: Takahashi et al. 2014
### Data Location: DRP001067
Explanation: Notebook 9:
This is an IPython notebook. Most of the code is composed of bash scripts, indicated by %%bash at the top of the cell, otherwise it is IPython code. This notebook includes code to download, assemble and analyze a published RADseq data set.
End of explanation
%%bash
## make a new directory for this analysis
mkdir -p empirical_9/fastq/
Explanation: Download the sequence data
Sequence data for this study are archived on the NCBI sequence read archive (SRA). Below I read in SraRunTable.txt for this project which contains all of the information we need to download the data. 
Project DRA: DRA001025
Study: DRP001067
SRA link: http://trace.ddbj.nig.ac.jp/DRASearch/study?acc=DRP001067
End of explanation
import os
def wget_download_ddbj(SRR, outdir):
     Python function to get sra data from ncbi and write to
    outdir with a new name using bash call wget 
    
    ## create a call string 
    call = "wget -q -r -nH --cut-dirs=9 -P "+outdir+" "+\
           "ftp://ftp.ddbj.nig.ac.jp/ddbj_database/dra/sra/ByExp/"+\
           "sra/DRX/DRX011/DRX011{:03d}".format(SRR)
    
    ## run wget call 
    ! $call
for ID in range(602,634):
    wget_download_ddbj(ID, "empirical_9/fastq/")
Explanation: For each ERS (individuals) get all of the ERR (sequence file accessions).
End of explanation
%%bash
## convert sra files to fastq using fastq-dump tool
## output as gzipped into the fastq directory
fastq-dump --gzip -O empirical_9/fastq/ empirical_9/fastq/*.sra
## remove .sra files
rm empirical_9/fastq/*.sra
%%bash
ls -lh empirical_9/fastq/
Explanation: Here we pass the SRR number and the sample name to the wget_download function so that the files are saved with their sample names.
End of explanation
%%bash
pyrad --version
%%bash
## remove old params file if it exists
rm params.txt 
## create a new default params file
pyrad -n 
Explanation: Make a params file
End of explanation
%%bash
## substitute new parameters into file
sed -i '/## 1. /c\empirical_9/           ## 1. working directory ' params.txt
sed -i '/## 6. /c\TGCAG                  ## 6. cutters ' params.txt
sed -i '/## 7. /c\20                     ## 7. N processors      ' params.txt
sed -i '/## 9. /c\6                      ## 9. NQual             ' params.txt
sed -i '/## 10./c\.85                    ## 10. clust threshold  ' params.txt
sed -i '/## 12./c\4                      ## 12. MinCov           ' params.txt
sed -i '/## 13./c\10                     ## 13. maxSH            ' params.txt
sed -i '/## 14./c\empirical_9_m4          ## 14. output name      ' params.txt
sed -i '/## 18./c\empirical_9/fastq/*.gz ## 18. data location    ' params.txt
sed -i '/## 29./c\2,2                    ## 29. trim overhang    ' params.txt
sed -i '/## 30./c\p,n,s                  ## 30. output formats   ' params.txt
cat params.txt
Explanation: Note:
The data here are from Illumina Casava <1.8, so the phred scores are offset by 64 instead of 33, so we use that in the params file below.
End of explanation
os.path.splitext(fil)
import glob
import itertools
import gzip
import os
for infile in glob.glob("empirical_9/fastq/DRR*"):
    iofile = gzip.open(infile, 'rb')
    dire, fil = os.path.split(infile)
    fastq = os.path.splitext(fil)[0]
    outhandle = os.path.join(dire, "T_"+fastq)
    outfile = open(outhandle, 'wb')
    
    data = iter(iofile)
    store = []
    while 1:
        try:
            line = data.next()
        except StopIteration:
            break
        if len(line) < 80:
            store.append(line)
        else:
            store.append(line[5:])
            
        if len(store) == 10000:
            outfile.write("".join(store))
            store = []
            
    iofile.close()
    outfile.close()        
    ! gzip $outhandle
Explanation: Trimming the barcode
In this data set the data were uploaded separated by sample, but with the barcode still attached to the sequences. The python code below will remove the 5bp barcode from each sequence.
End of explanation
%%bash
pyrad -p params.txt -s 234567 >> log.txt 2>&1  
%%bash
sed -i '/## 12./c\2                    ## 12. MinCov           ' params.txt
sed -i '/## 14./c\empirical_9_m2       ## 14. output name      ' params.txt
%%bash
pyrad -p params.txt -s 7 >> log.txt 2>&1  
Explanation: Assemble in pyrad
End of explanation
import pandas as pd
import numpy as np
## read in the data
s2dat = pd.read_table("empirical_9/stats/s2.rawedit.txt", header=0, nrows=32)
## print summary stats
print s2dat["passed.total"].describe()
## find which sample has the most raw data
maxraw = s2dat["passed.total"].max()
print "\nmost raw data in sample:"
print s2dat['sample '][s2dat['passed.total']==maxraw]
Explanation: Results
We are interested in the relationship between the amount of input (raw) data between any two samples, the average coverage they recover when clustered together, and the phylogenetic distances separating samples.  
Raw data amounts
The average number of raw reads per sample is 1.36M.
End of explanation
## read in the s3 results
s9dat = pd.read_table("empirical_9/stats/s3.clusters.txt", header=0, nrows=33)
## print summary stats
print "summary of means\n=================="
print s9dat['dpt.me'].describe()
## print summary stats
print "\nsummary of std\n=================="
print s9dat['dpt.sd'].describe()
## print summary stats
print "\nsummary of proportion lowdepth\n=================="
print pd.Series(1-s9dat['d>5.tot']/s9dat["total"]).describe()
## find which sample has the greatest depth of retained loci
max_hiprop = (s9dat["d>5.tot"]/s9dat["total"]).max()
print "\nhighest coverage in sample:"
print s9dat['taxa'][s9dat['d>5.tot']/s9dat["total"]==max_hiprop]
maxprop =(s9dat['d>5.tot']/s9dat['total']).max()
print "\nhighest prop coverage in sample:"
print s9dat['taxa'][s9dat['d>5.tot']/s9dat['total']==maxprop]
## print mean and std of coverage for the highest coverage sample
with open("empirical_9/clust.85/T_DRR021966.depths", 'rb') as indat:
    depths = np.array(indat.read().strip().split(","), dtype=int)
    
print "Means for sample T_DRR021966"
print depths.mean(), depths.std()
print depths[depths>5].mean(), depths[depths>5].std()
Explanation: Look at distributions of coverage
pyrad v.3.0.63 outputs depth information for each sample which I read in here and plot. First let's ask which sample has the highest depth of coverage. The std of coverages is pretty low in this data set compared to several others.
End of explanation
import toyplot
import toyplot.svg
import numpy as np
## read in the depth information for this sample
with open("empirical_9/clust.85/T_DRR021966.depths", 'rb') as indat:
    depths = np.array(indat.read().strip().split(","), dtype=int)
    
## make a barplot in Toyplot
canvas = toyplot.Canvas(width=350, height=300)
axes = canvas.axes(xlabel="Depth of coverage (N reads)", 
                   ylabel="N loci", 
                   label="dataset9/sample=T_DRR021966")
## select the loci with depth > 5 (kept)
keeps = depths[depths>5]
## plot kept and discarded loci
edat = np.histogram(depths, range(30)) # density=True)
kdat = np.histogram(keeps, range(30)) #, density=True)
axes.bars(edat)
axes.bars(kdat)
#toyplot.svg.render(canvas, "empirical_9_depthplot.svg")
Explanation: Plot the coverage for the sample with highest mean coverage
Green shows the loci that were discarded and orange the loci that were retained. The majority of data were discarded for being too low of coverage.
End of explanation
cat empirical_9/stats/empirical_9_m4.stats
%%bash
head -n 20 empirical_9/stats/empirical_9_m2.stats
Explanation: Print final stats table
End of explanation
%%bash
## raxml argumement w/ ...
raxmlHPC-PTHREADS-AVX -f a -m GTRGAMMA -N 100 -x 12345 -p 12345 -T 20 \
                      -w /home/deren/Documents/RADmissing/empirical_9/ \
                      -n empirical_9_m4 -s empirical_9/outfiles/empirical_9_m4.phy
                      
%%bash
## raxml argumement w/ ...
raxmlHPC-PTHREADS-AVX -f a -m GTRGAMMA -N 100 -x 12345 -p 12345 -T 20 \
                      -w /home/deren/Documents/RADmissing/empirical_9/ \
                      -n empirical_9_m2 -s empirical_9/outfiles/empirical_9_m2.phy
                      
%%bash 
head -n 20 empirical_9/RAxML_info.empirical_9_m4
%%bash 
head -n 20 empirical_9/RAxML_info.empirical_9_m2
Explanation: Infer ML phylogeny in raxml as an unrooted tree
End of explanation
%load_ext rpy2.ipython
%%R -h 800 -w 800
library(ape)
tre <- read.tree("empirical_9/RAxML_bipartitions.empirical_9")
ltre <- ladderize(tre)
par(mfrow=c(1,2))
plot(ltre, use.edge.length=F)
nodelabels(ltre$node.label)
plot(ltre, type='u')
Explanation: Plot the tree in R using ape
End of explanation
%%R
mean(cophenetic.phylo(ltre))
Explanation: Get phylo distances (GTRgamma dist)
End of explanation
<END_TASK>