<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 Accuracy as a function taxon evenness & % atom 13C
Variable parameters
Step1:  Init
Step2:  Nestly
Step3:  Adding DESeq results table of 'NA' if file is not present
Step4:  aggregating confusion matrix data
Step5:  Plotting accuracy
Dominant taxa
Step6:  Rare taxa
Step7:  Sample size of dominants and rares
Method
Step8:  Plotting sensitivity ~ taxon_abundance
Method
Step9:  Notes
Step10:  
<ASSISTANT_TASK:>
Python Code:
workDir = '/home/nick/notebook/SIPSim/dev/bac_genome1210/'
sourceDir = os.path.join(workDir, 'atomIncorp_evenness')
buildDir = os.path.join(workDir, 'atomIncorp_evenness_abund')
R_dir = '/home/nick/notebook/SIPSim/lib/R/'
Explanation: Accuracy as a function taxon evenness & % atom 13C
Variable parameters:
Partition communities into 'dominant' and 'rare' subsets
dominant/rare cutoff: 1% relative abundance?
Assess sensitivity & specificity for each subset
Note: 
Much less time to perform if using data from past simulation run datasets
Analysis:
Sup Figures:
Same as above simulation runs, but split into 'dominant' vs 'rare'
User variables
End of explanation
import glob
from os.path import abspath
import nestly
from IPython.display import Image, display
%load_ext rpy2.ipython
%%R
library(ggplot2)
library(dplyr)
library(tidyr)
library(gridExtra)
if not os.path.isdir(buildDir):
    os.makedirs(buildDir)
Explanation: Init
End of explanation
# building tree structure
nest = nestly.Nest()
## varying params
#-- full run --#
#nest.add('percIncorp', [0, 25, 50, 100])
#nest.add('percTaxa', [1, 5, 10, 25, 50])
#nest.add('rep', range(1,21))
#-- full run --#
nest.add('abund_dist_p', ['mean:10,sigma:1','mean:10,sigma:2','mean:10,sigma:3'])
nest.add('percIncorp', [0,5,10,15,20,30,40,50,75,100])
nest.add('rep', range(1,21))
## set params
nest.add('percTaxa', [10], create_dir=False)
nest.add('np', [6], create_dir=False)
nest.add('abs', ['1e9'], create_dir=False)
nest.add('subsample_mean', [30000], create_dir=False) 
nest.add('subsample_scale', [5000], create_dir=False) 
nest.add('BD_min', [1.71], create_dir=False)
nest.add('BD_max', [1.75], create_dir=False)
nest.add('padj', [0.1], create_dir=False)
nest.add('log2', [0.25], create_dir=False)
nest.add('topTaxaToPlot', [100], create_dir=False)
## input/output files
nest.add('buildDir', [buildDir], create_dir=False)
nest.add('sourceDir', [sourceDir], create_dir=False)
nest.add('frag_file', ['ampFrags_kde_dif'], create_dir=False)
nest.add('comm_file', ['comm.txt'], create_dir=False)
nest.add('R_dir', [R_dir], create_dir=False)
# building directory tree
nest.build(buildDir)
bashFile = os.path.join(buildDir, 'SIPSimRun.sh')
%%writefile $bashFile
#!/bin/bash
export PATH={R_dir}:$PATH
# copying DESeq files
cp {sourceDir}/{abund_dist_p}/{percIncorp}/{rep}/{comm_file} .
cp {sourceDir}/{abund_dist_p}/{percIncorp}/{rep}/OTU_n2_abs{abs}_sub-norm_DESeq2 .
cp {sourceDir}/{abund_dist_p}/{percIncorp}/{rep}/{frag_file}_incorp_BD-shift.txt .
## DESeq2
DESeq2_rare-dominant.r \
    OTU_n2_abs{abs}_sub-norm_DESeq2 \
    {comm_file} \
    -o OTU_n2_abs{abs}_sub-norm_DESeq2
    
## Confusion matrix
### dominant
DESeq2_confuseMtx.r \
    {frag_file}_incorp_BD-shift.txt \
    OTU_n2_abs{abs}_sub-norm_DESeq2_dom \
    --padj {padj} \
    -o DESeq2-cMtx_dom
### rare    
DESeq2_confuseMtx.r \
    {frag_file}_incorp_BD-shift.txt \
    OTU_n2_abs{abs}_sub-norm_DESeq2_rare \
    --padj {padj} \
    -o DESeq2-cMtx_rare
!chmod 775 $bashFile
!cd $workDir; \
    nestrun -j 20 --template-file $bashFile -d atomIncorp_evenness_abund --log-file log.txt
Explanation: Nestly
End of explanation
p = os.path.join(workDir, 'atomIncorp_evenness_abund', 'mean*/*/*', 'DESeq2-cMtx_rare_byClass.csv')
rareFiles = glob.glob(p)
p = os.path.join(workDir, 'atomIncorp_evenness_abund', 'mean*/*/*', 'DESeq2-cMtx_dom_byClass.csv')
domFiles = glob.glob(p)
print len(rareFiles)
print len(domFiles)
byClass_tmp = [
['Sensitivity','NA','NA',0,0],
['Specificity','NA','NA',0,0],
['Pos Pred Value','NA','NA',0,0],
['Neg Pred Value','NA','NA',0,0],
['Prevalence','NA','NA',0,0],
['Detection Rate','NA','NA',0,0],
['Detection Prevalence','NA','NA',0,0],
['Balanced Accuracy','NA','NA',0,0]
]
byClass_tmp = pd.DataFrame(byClass_tmp, columns=['','byClass','abund_dist_p','percIncorp','rep'])
# finding files not present
rareFiles = set([x.rstrip('_rare_byClass.csv') for x in rareFiles])
domFiles = set([x.rstrip('_dom_byClass.csv') for x in domFiles])
# in rare but not in dominant list
missingFiles = list(rareFiles - domFiles)
cols = ['abund_dist_p','percIncorp','rep']
for f in missingFiles:    
    df_tmp = byClass_tmp.copy()
    x = f.split('/')
    df_tmp.loc[:,cols] = x[-4:-1]
    outFile = '_'.join([f, 'dom_byClass.csv'])
    df_tmp.to_csv(outFile)
Explanation: Adding DESeq results table of 'NA' if file is not present
End of explanation
## byClass
### dominant
!cd $workDir; \
    nestagg delim \
        -d atomIncorp_evenness_abund \
        -k abund_dist_p,percIncorp,rep \
        -o ./atomIncorp_evenness_abund/DESeq2-cMtx_dom_byClass.csv \
        DESeq2-cMtx_dom_byClass.csv
### rare
!cd $workDir; \
    nestagg delim \
        -d atomIncorp_evenness_abund \
        -k abund_dist_p,percIncorp,rep \
        -o ./atomIncorp_evenness_abund/DESeq2-cMtx_rare_byClass.csv \
        DESeq2-cMtx_rare_byClass.csv        
Explanation: aggregating confusion matrix data
End of explanation
%%R -i workDir
setwd(workDir)
byClass = read.csv('./atomIncorp_evenness_abund/DESeq2-cMtx_dom_byClass.csv')
byClass %>% head
%%R -w 700
to.keep = c('Sensitivity', 'Specificity', 'Balanced Accuracy')
byClass.f = byClass %>%
    filter(X %in% to.keep) %>%
    mutate(percIncorp = as.character(percIncorp))
byClass.f$percIncorp = factor(byClass.f$percIncorp, levels=sort(unique(as.numeric(byClass.f$percIncorp))))
ggplot(byClass.f, aes(percIncorp, byClass, fill=abund_dist_p)) +
    geom_boxplot() +
    labs(y='Value', x='Atom percent 13C') +
    facet_grid(X ~ .) +
        theme(
        text = element_text(size=16)
        )
Explanation: Plotting accuracy
Dominant taxa
End of explanation
%%R -i workDir
setwd(workDir)
byClass = read.csv('./atomIncorp_evenness_abund/DESeq2-cMtx_rare_byClass.csv')
byClass %>% head
%%R -w 700
to.keep = c('Sensitivity', 'Specificity', 'Balanced Accuracy')
byClass.f = byClass %>%
    filter(X %in% to.keep) %>%
    mutate(percIncorp = as.character(percIncorp))
byClass.f$percIncorp = factor(byClass.f$percIncorp, levels=sort(unique(as.numeric(byClass.f$percIncorp))))
ggplot(byClass.f, aes(percIncorp, byClass, fill=abund_dist_p)) +
    geom_boxplot() +
    labs(y='Value', x='Atom percent 13C') +
    facet_grid(X ~ .) +
        theme(
        text = element_text(size=16)
        )
Explanation: Rare taxa
End of explanation
%%R -i workDir
setwd(workDir)
dom.files = grep('DESeq2_dom$', list.files('./atomIncorp_evenness_abund/', recursive=TRUE), value=TRUE)
rare.files = grep('DESeq2_rare$', list.files('./atomIncorp_evenness_abund/', recursive=TRUE), value=TRUE)
read_taxa = function(file.name, dir='atomIncorp_evenness_abund'){
    f = paste(c(dir, file.name), collapse='/')
    deseq.res = readRDS(f)
    n.taxa = length(unique(rownames(deseq.res)))
    return(n.taxa)
    }
dom.taxa = list()
for(f in dom.files){
    dom.taxa[[f]] = read_taxa(f)
    }
rare.taxa = list()
for(f in rare.files){
    rare.taxa[[f]] = read_taxa(f)
    }
%%R 
tbl.dom = as.data.frame(do.call(rbind, dom.taxa))
tbl.dom$abund = 'Dominant'
tbl.rare = as.data.frame(do.call(rbind, rare.taxa))
tbl.rare$abund = 'Rare'
tbl.dom.rare = rbind(tbl.dom, tbl.rare)
colnames(tbl.dom.rare)[1] = 'N.taxa'
tbl.dom.rare$file = rownames(tbl.dom.rare) 
tbl.dom.rare = tbl.dom.rare %>% 
    separate(file, c('abund_dist_p','percIncorp','rep','file'), sep='/') 
tbl.dom.rare %>% head
%%R -w 700 -h 400
tbl.dom.rare$percIncorp = factor(tbl.dom.rare$percIncorp, levels=sort(unique(as.numeric(tbl.dom.rare$percIncorp))))
#tbl.dom.rare$abund_dist_p = factor(tbl.dom.rare$abund_dist_p, levels=sort(unique(as.numeric(tbl.dom.rare$abund_dist_p))))
ggplot(tbl.dom.rare, aes(percIncorp, N.taxa, fill=abund_dist_p)) +
    geom_boxplot(position='dodge') +
    labs(y='Number of taxa', x='Atom percent 13C') +
    facet_grid(abund ~ ., scales='free_y') +
    theme(
        text = element_text(size=16)
        )
Explanation: Sample size of dominants and rares
Method:
Get number of taxa from dominant/rare-parsed DESeq2 object files
End of explanation
%%bash -s "$R_dir"
find atomIncorp_evenness_abund/ -name "*DESeq2" | \
    perl -pe 's/(.+\/).+/\1/' | \
    xargs -P 30 -I % bash -c \
    "$1DESeq2_addInfo.r %OTU_n2_abs1e9_sub-norm_DESeq2 %comm.txt %ampFrags_kde_dif_incorp_BD-shift.txt > %OTU_n2_abs1e9_sub-norm_DESeq2_info"
    
%%R -i workDir
setwd(workDir)
deseq.files = grep('_DESeq2_info$', list.files('./atomIncorp_evenness_abund/', recursive=TRUE), value=TRUE)
tbl.l = list()
for(f in deseq.files){
    f = paste(c('./atomIncorp_evenness_abund/', f), collapse='/')
    tbl.l[[f]] = readRDS(f) %>% mutate(file = f)
    }
tbl = do.call(rbind, tbl.l)
rownames(tbl) = seq(1,nrow(tbl))
tbl %>% head
%%R
tbl.e = tbl %>%
    mutate(file = gsub('.+//','', file)) %>%
    separate(file, c('abund_dist_p','percIncorp','rep','file'), sep='/') #%>%
    #unite(percIncorp_percTaxa, percIncorp, percTaxa, sep='_', remove=FALSE)
tbl.e %>% head
%%R
BD.shift.cut = 0.05
padj.cut = 0.1
clsfy = function(guess,known){
    if(is.na(guess) | is.na(known)){
        return(NA)
        }
    if(guess == TRUE){
        if(guess == known){
            return('TP')
            } else {
            return('FP')
            }
        } else
    if(guess == FALSE){
        if(guess == known){
            return('TN')
            } else {
            return('FN')
            }
        } else {
        stop('Error: true or false needed')
        }
    }
tbl.e = tbl.e %>% 
    mutate(true_incorp = BD_shift > BD.shift.cut,
           DESeq_incorp = padj < padj.cut,
           cls = mapply(clsfy, DESeq_incorp, true_incorp)) 
    
tbl.e %>% head    
%%R
tbl.e$taxon = reorder(tbl.e$taxon, dense_rank(tbl.e$mean_rel_abund_perc))
tbl.e$percIncorp = factor(tbl.e$percIncorp, levels=sort(unique(as.numeric(tbl.e$percIncorp))))
#tbl.e$percTaxa = factor(tbl.e$percTaxa, levels=sort(unique(as.numeric(tbl.e$percTaxa))))
%%R -w 1000
ggplot(tbl.e, aes(taxon, mean_rel_abund_perc, color=cls)) +
    geom_point(alpha=0.5) +
    facet_grid(abund_dist_p ~ percIncorp) +
    theme(
        text = element_text(size=16),
        axis.text.x = element_blank()
        )
%%R
clsfy = function(guess,known){
    if(is.na(guess) | is.na(known)){
        return(NA)
        }
    if(guess == TRUE){
        if(guess == known){
            return(1)
            } else {
            return(NA)
            }
        } else
    if(guess == FALSE){
        if(guess == known){
            return(NA)
            } else {
            return(0)
            }
        } else {
        stop('Error: true or false needed')
        }
    }
tbl.s.c = tbl.e %>%
    mutate(TP_FN = mapply(clsfy, DESeq_incorp, true_incorp)) %>%
    filter(! is.na(TP_FN)) %>%
    group_by(percIncorp, abund_dist_p) %>%
    summarize(cor.pearson = cor(mean_rel_abund_perc, TP_FN))
tbl.s.c %>% head
%%R -w 700 -h 350
ggplot(tbl.s.c, aes(percIncorp, cor.pearson, fill=abund_dist_p)) +
    geom_bar(stat='identity', position='dodge', width=0.5) +
    labs(x='Atom percent 13C', y='r_pb', title='Point-biserial correlation coefficients') +
    theme(
        text = element_text(size=16)
        )
Explanation: Plotting sensitivity ~ taxon_abundance
Method: 
table join: DESeq res file & BD-shift
join on taxon
use cutoffs to select which taxa 'incorporated'
padj < 0.1
compare DESeq result to true result
TP,FP,TN,FN
Figure:
axes: taxon ~ abundance (pointrange: replicates)
ordered by abundance
color by classification (TP,FP,TN,FN)
facet: percIncorp ~ percTaxa
End of explanation
%%R -i workDir
setwd(workDir)
x.deseq = './validation/OTU_n2_abs1e10_sub20000_DESeq2'
deseq.res = readRDS(x.deseq)
deseq.res = as.data.frame(deseq.res)
deseq.res$taxon = rownames(deseq.res)
deseq.res %>% head
%%R -i workDir
setwd(workDir) 
x.comm = 'validation/comm.txt'
comm = read.delim(x.comm, sep='\t')
comm.s = comm %>%
    group_by(taxon_name) %>%
    summarize(mean_rel_abund_perc = mean(rel_abund_perc)) %>%
    mutate(rank = dense_rank(mean_rel_abund_perc))
comm.s %>% head 
%%R -i workDir
setwd(workDir) 
x.BD = 'validation/ampFrags_kde_dif_incorp_BD-shift.txt'
BD.shift = read.delim(x.BD, sep='\t')
BD.shift = BD.shift %>% 
    filter(lib2 == 2)
BD.shift %>% head 
%%R
tbl.j = inner_join(deseq.res, BD.shift, c('taxon' = 'taxon'))
tbl.j = inner_join(tbl.j, comm.s, c('taxon' = 'taxon_name'))
tbl.j %>% head
%%R
BD.shift.cut = 0.05
padj.cut = 0.1
clsfy = function(guess,known){
    if(guess == TRUE){
        if(guess == known){
            return('TP')
            } else {
            return('FP')
            }
        } else
    if(guess == FALSE){
        if(guess == known){
            return('TN')
            } else {
            return('FN')
            }
        } else {
        stop('Error: true or false needed')
        }
    }
tbl.j = tbl.j %>% 
    mutate(true_incorp = BD_shift > BD.shift.cut,
           DESeq_incorp = padj < padj.cut,
           cls = mapply(clsfy, DESeq_incorp, true_incorp)) 
    
tbl.j %>% head    
%%R -w 900
tbl.j$taxon = reorder(tbl.j$taxon, tbl.j$mean_rel_abund_perc)
ggplot(tbl.j, aes(taxon, mean_rel_abund_perc, color=cls)) +
    geom_point()
Explanation: Notes:
SANDBOX
End of explanation
%%R -i workDir
setwd(workDir)
deseq.res = readRDS('./validation/OTU_n2_abs1e10_sub20000_DESeq2')
%%R
comm = read.delim('./validation/comm.txt', sep='\t')
comm.s = comm %>%
    group_by(taxon_name) %>%
    summarize(mean_rel_abund_perc = mean(rel_abund_perc))
comm.s %>% head
%%R
tmp = deseq.res %>%
    as.data.frame() 
tmp$taxon_name = rownames(tmp)
tmp = inner_join(tmp, comm.s, c('taxon_name' = 'taxon_name'))
tmp = tmp %>%
    mutate(dom_rare = mean_rel_abund_perc > 1)
tmp %>% head
%%R
tmp %>% filter(dom_rare == FALSE) %>% nrow %>% print
tmp %>% filter(dom_rare == TRUE) %>% nrow %>% print
Explanation: 
End of explanation
<END_TASK>