### Notebook 10
### Data set 10: Pedicularis
### Authors: Eaton & Ree 2013
### Data Location: SRP021469

"""
Explanation: Notebook 10:
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_10/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. 

SRA link: http://trace.ncbi.nlm.nih.gov/Traces/study/?acc=SRP021469
End of explanation
"""


## IPython code
import pandas as pd
import numpy as np
import urllib2
import os

## open the SRA run table from github url
url = "https://raw.githubusercontent.com/"+\
      "dereneaton/RADmissing/master/empirical_10_SraRunTable.txt"
intable = urllib2.urlopen(url)
indata = pd.read_table(intable, sep="\t")

## print first few rows
print indata.head()

def wget_download(SRR, outdir, outname):
    """ Python function to get sra data from ncbi and write to
    outdir with a new name using bash call wget """
    
    ## get output name
    output = os.path.join(outdir, outname+".sra")
    
    ## create a call string 
    call = "wget -q -r -nH --cut-dirs=9 -O "+output+" "+\
           "ftp://ftp-trace.ncbi.nlm.nih.gov/"+\
           "sra/sra-instant/reads/ByRun/sra/SRR/"+\
           "{}/{}/{}.sra;".format(SRR[:6], SRR, SRR)
        
    ## call bash script
    ! $call 

"""
Explanation: For each ERS (individuals) get all of the ERR (sequence file accessions).
End of explanation
"""


for ID, SRR in zip(indata.Sample_Name_s, indata.Run_s):
    wget_download(SRR, "empirical_10/fastq/", ID)

%%bash
## convert sra files to fastq using fastq-dump tool
## output as gzipped into the fastq directory
fastq-dump --gzip -O empirical_10/fastq/ empirical_10/fastq/*.sra

## remove .sra files
rm empirical_10/fastq/*.sra

%%bash
ls -lh empirical_10/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_10/          ## 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_10_m4         ## 14. output name      ' params.txt
sed -i '/## 18./c\empirical_10/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
"""


%%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_10_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
## read in the data
sdat = pd.read_table("empirical_10/stats/s2.rawedit.txt", header=0, nrows=14)

## print summary stats
print sdat["passed.total"].describe()

## find which sample has the most raw data
maxraw = sdat["passed.total"].max()
print "\nmost raw data in sample:"
print sdat['sample '][sdat['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
sdat = pd.read_table("empirical_10/stats/s3.clusters.txt", header=0, nrows=14)

## print summary stats
print "summary of means\n=================="
print sdat['dpt.me'].describe()

## print summary stats
print "\nsummary of std\n=================="
print sdat['dpt.sd'].describe()

## print summary stats
print "\nsummary of proportion lowdepth\n=================="
print pd.Series(1-sdat['d>5.tot']/sdat["total"]).describe()

## find which sample has the greatest depth of retained loci
max_hiprop = (sdat["d>5.tot"]/sdat["total"]).max()
print "\nhighest coverage in sample:"
print sdat['taxa'][sdat['d>5.tot']/sdat["total"]==max_hiprop]

maxprop =(sdat['d>5.tot']/sdat['total']).max()
print "\nhighest prop coverage in sample:"
print sdat['taxa'][sdat['d>5.tot']/sdat['total']==maxprop]

import numpy as np
## print mean and std of coverage for the highest coverage sample
with open("empirical_10/clust.85/38362_rex.depths", 'rb') as indat:
    depths = np.array(indat.read().strip().split(","), dtype=int)
    
print "Means for sample 38362_rex"
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_10/clust.85/38362_rex.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="dataset10/sample=38362_rex")

## 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_10_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_10/stats/empirical_10_m4.stats


%%bash
head -n 10 empirical_10/stats/empirical_10_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_10/ \
                      -n empirical_10_m4 -s empirical_10/outfiles/empirical_10_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_10/ \
                      -n empirical_10_m2 -s empirical_10/outfiles/empirical_10_m2.phy
                      

%%bash 
head -n 20 empirical_10/RAxML_info.empirical_10_m4

%%bash 
head -n 20 empirical_10/RAxML_info.empirical_10_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_10/RAxML_bipartitions.empirical_10_m4")
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
"""


