workDir = '/home/nick/notebook/SIPSim/dev/bac_genome1210/'
genomeDir = '/home/nick/notebook/SIPSim/dev/bac_genome1210/genomes/'
R_dir = '/home/nick/notebook/SIPSim/lib/R/'

"""
Explanation: Goal
Question: how is incorporator identification accuracy affected by the percent isotope incorporation of taxa?


Using genome dataset created in the "dataset" notebook


Simulates isotope dilution or short incubations

Method  
25% taxa incorporate
incorporation % same for all incorporators
incorporation % treatments: 10, 20, 40, 60, 80, 100%
Total treatments: 6
Title: perc_incorp



User variables
End of explanation
"""


import glob
from os.path import abspath
import nestly

%load_ext rpy2.ipython

%%R
library(ggplot2)
library(dplyr)
library(tidyr)
library(gridExtra)

"""
Explanation: Init
End of explanation
"""


# building tree structure
nest = nestly.Nest()

## varying params
nest.add('percIncorp', range(0,101,20))

## set params
nest.add('np_many', [24], create_dir=False)
nest.add('np_few', [8], create_dir=False)
nest.add('percTaxa', [25], create_dir=False)
nest.add('abs', ['1e10'], create_dir=False)
nest.add('subsample', [20000], create_dir=False) 
nest.add('BD_min', [1.71], create_dir=False)
nest.add('BD_max', [2], 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('fileName', ['ampFrags'], create_dir=False)
nest.add('genome_index', [os.path.join(genomeDir, 'genome_index.txt')], create_dir=False)
nest.add('genome_dir', [genomeDir], create_dir=False)
nest.add('primers', [os.path.join(workDir, '../', '515F-806R.fna')], create_dir=False)
nest.add('R_dir', [R_dir], create_dir=False)

# building directory tree
buildDir = os.path.join(workDir, 'percIncorpUnif')
nest.build(buildDir)

bashFile = os.path.join(workDir, 'SIPSimRun.sh')

%%writefile $bashFile
#!/bin/bash

# simulating fragments
SIPSim fragments \
    {genome_index} \
    --fp {genome_dir} \
    --fr {primers} \
    --fld skewed-normal,9000,2500,-5 \
    --flr None,None \
    --nf 10000 \
    --np {np_many} \
    2> {fileName}.log \
    > {fileName}.pkl
    
# converting to kde object
SIPSim fragment_kde \
    {fileName}.pkl \
    > {fileName}_kde.pkl
    
# adding diffusion
SIPSim diffusion \
    {fileName}_kde.pkl \
    --np {np_many} \
    > {fileName}_kde_dif.pkl
    
# creating a community file
SIPSim gradientComms \
    {genome_index} \
    --n_comm 2 \
    > comm.txt

# making incorp file
SIPSim incorpConfigExample \
  --percTaxa {percTaxa} \
  --percIncorpUnif {percIncorp} \
  > {percTaxa}_{percIncorp}.config

# adding isotope incorporation to BD distribution
SIPSim isoIncorp \
    {fileName}_kde_dif.pkl \
    {percTaxa}_{percIncorp}.config \
    --comm comm.txt \
    --np {np_many} \
    > {fileName}_kde_dif_incorp.pkl
 
# calculating BD shift from isotope incorporation
SIPSim BD_shift \
    {fileName}_kde_dif.pkl \
    {fileName}_kde_dif_incorp.pkl \
    --np {np_few} \
    > {fileName}_kde_dif_incorp_BD-shift.txt
       
# simulating gradient fractions
SIPSim fractions \
    comm.txt \
    > fracs.txt

# simulating an OTU table
SIPSim OTU_table \
    {fileName}_kde_dif_incorp.pkl \
    comm.txt \
    fracs.txt \
    --abs {abs} \
    --np {np_few} \
    > OTU_abs{abs}.txt
    
# subsampling from the OTU table (simulating sequencing of the DNA pool)
SIPSim OTU_subsample \
    --dist_params low:{subsample},high:{subsample} \
    OTU_abs{abs}.txt \
    > OTU_n2_abs{abs}_sub{subsample}.txt

# making a wide table
SIPSim OTU_wideLong -w \
    OTU_n2_abs{abs}_sub{subsample}.txt \
    > OTU_n2_abs{abs}_sub{subsample}_w.txt
    
# making metadata (phyloseq: sample_data)
SIPSim OTU_sampleData \
    OTU_n2_abs{abs}_sub{subsample}.txt \
    > OTU_n2_abs{abs}_sub{subsample}_meta.txt
  
#-- R analysis --#
export PATH={R_dir}:$PATH
# plotting taxon abundances
OTU_taxonAbund.r \
    OTU_n2_abs{abs}_sub{subsample}.txt \
    -r {topTaxaToPlot} \
    -o OTU_n2_abs{abs}_sub{subsample} 
# running DeSeq2 and making confusion matrix on predicting incorporators
## making phyloseq object from OTU table
phyloseq_make.r \
    OTU_n2_abs{abs}_sub{subsample}_w.txt \
    -s OTU_n2_abs{abs}_sub{subsample}_meta.txt \
    > OTU_n2_abs{abs}_sub{subsample}.physeq
## filtering phyloseq object to just taxa/samples of interest
phyloseq_edit.r \
    OTU_n2_abs{abs}_sub{subsample}.physeq \
    --BD_min {BD_min} --BD_max {BD_max} \
    > OTU_n2_abs{abs}_sub{subsample}_filt.physeq
## making ordination
phyloseq_ordination.r \
    OTU_n2_abs{abs}_sub{subsample}_filt.physeq \
    OTU_n2_abs{abs}_sub{subsample}_bray-NMDS.pdf
## DESeq2
phyloseq_DESeq2.r \
    OTU_n2_abs{abs}_sub{subsample}_filt.physeq \
    > OTU_n2_abs{abs}_sub{subsample}_DESeq2
## Confusion matrix
DESeq2_confuseMtx.r \
    {fileName}_kde_dif_incorp_BD-shift.txt \
    OTU_n2_abs{abs}_sub{subsample}_DESeq2 \
    --padj {padj} --log2 {log2} 

%%writefile $bashFile
#!/bin/bash

## just R analysis

#-- R analysis --#
export PATH={R_dir}:$PATH
# plotting taxon abundances
OTU_taxonAbund.r \
    OTU_n2_abs{abs}_sub{subsample}.txt \
    -r {topTaxaToPlot} \
    -o OTU_n2_abs{abs}_sub{subsample} 
# running DeSeq2 and making confusion matrix on predicting incorporators
## making phyloseq object from OTU table
phyloseq_make.r \
    OTU_n2_abs{abs}_sub{subsample}_w.txt \
    -s OTU_n2_abs{abs}_sub{subsample}_meta.txt \
    > OTU_n2_abs{abs}_sub{subsample}.physeq
## filtering phyloseq object to just taxa/samples of interest
phyloseq_edit.r \
    OTU_n2_abs{abs}_sub{subsample}.physeq \
    --BD_min {BD_min} --BD_max {BD_max} \
    > OTU_n2_abs{abs}_sub{subsample}_filt.physeq
## making ordination
phyloseq_ordination.r \
    OTU_n2_abs{abs}_sub{subsample}_filt.physeq \
    OTU_n2_abs{abs}_sub{subsample}_bray-NMDS.pdf
## DESeq2
phyloseq_DESeq2.r \
    OTU_n2_abs{abs}_sub{subsample}_filt.physeq \
    --log2 {log2} \
    > OTU_n2_abs{abs}_sub{subsample}_DESeq2
## Confusion matrix
DESeq2_confuseMtx.r \
    {fileName}_kde_dif_incorp_BD-shift.txt \
    OTU_n2_abs{abs}_sub{subsample}_DESeq2 \
    --padjBH {padj} 

!chmod 775 $bashFile

!cd $workDir; \
    nestrun -j 6 --template-file $bashFile -d percIncorpUnif --log-file logR.txt

# aggregating confusion matrix data
## table
!cd $workDir; \
    nestagg delim \
        -d percIncorpUnif \
        -k percIncorp \
        -o ./percIncorpUnif/DESeq2-cMtx_table.csv \
        DESeq2-cMtx_table.csv

## overall
!cd $workDir; \
    nestagg delim \
        -d percIncorpUnif \
        -k percIncorp \
        -o ./percIncorpUnif/DESeq2-cMtx_overall.csv \
        DESeq2-cMtx_overall.csv
    
## byClass
!cd $workDir; \
    nestagg delim \
        -d percIncorpUnif \
        -k percIncorp \
        -o ./percIncorpUnif/DESeq2-cMtx_byClass.csv \
        DESeq2-cMtx_byClass.csv

"""
Explanation: Using nestly: different incorporation percentages
End of explanation
"""


%%R -i workDir -w 600 -h 600
setwd(workDir)

byClass = read.csv('./percIncorpUnif/DESeq2-cMtx_byClass.csv')

byClass$byClass[is.na(byClass$byClass)] = 0

byClass$percIncorp = factor(byClass$percIncorp, levels=as.character(unique(sort(byClass$percIncorp))))
p = ggplot(byClass, aes(X, byClass, color=percIncorp)) +
    geom_point() +
    labs(y='value') +
    theme_bw() +
    theme(
        text=element_text(size=16),
        axis.text.x=element_text(angle=90, hjust=1, vjust=0.5),
        axis.title.x=element_blank()
        )
p

%%R -w 600 -h 400

p + scale_y_continuous(limits=c(0,0.15))

"""
Explanation: Plotting results
End of explanation
"""


%%R -i workDir

setwd(workDir)

tbl.c = read.csv('percIncorpUnif/50/DESeq2-cMtx_data.csv')

tbl.c.f = tbl.c %>% 
    filter(incorp.known==TRUE & incorp.pred==FALSE) %>%
    select(taxon)

write.table(tbl.c.f, file='percIncorpUnif/50/DESeq2-cMtx_FN', row.names=F, col.names=F, quote=F)

!cd $workDir'percIncorpUnif/50/'; \
    $R_dir/OTU_taxonAbund.r OTU_n2_abs1e8_sub30000.txt -t DESeq2-cMtx_FN -l -o DESeq2-cMtx_FN --width 14

%%R -i workDir

setwd(workDir)

tbl.c = read.csv('percIncorpUnif/50/DESeq2-cMtx_data.csv')

tbl.c.f = tbl.c %>% 
    filter(incorp.known==FALSE & incorp.pred==TRUE) %>%
    select(taxon)

write.table(tbl.c.f, file='percIncorpUnif/50/DESeq2-cMtx_FP', row.names=F, col.names=F, quote=F)

!cd $workDir'percIncorpUnif/50/'; \
    $R_dir/OTU_taxonAbund.r OTU_n2_abs1e8_sub30000.txt -t DESeq2-cMtx_FP -l -o DESeq2-cMtx_FP --width 14

%%R
tbl.c %>% filter(taxon %in% tbl.c.f$taxon)

"""
Explanation: Plotting false negative taxa (actual incorporators)
End of explanation
"""


%%R -i workDir

setwd(workDir)

tbl.c = read.csv('percIncorpUnif/100/DESeq2-cMtx_data.csv')
tbl.otu = read.delim('percIncorpUnif/100/OTU_n2_abs1e10_sub20000.txt', sep='\t')

%%R
# OTU total counts
tbl.otu.sum = tbl.otu %>%
    group_by(library, taxon) %>%
    summarize(total_count = sum(count))
tbl.otu.sum %>% head

%%R
# 
label.tp.fn = function(known, pred){
    if(is.na(known) | is.na(pred)){
        return(NA)
        } else
    if(known==TRUE & pred==TRUE){
        return('TP')
        } else
    if(known==TRUE & pred==FALSE){
        return('FN')
        } else {
        return(NA)
        }
    }

tbl.c.tp.fn = tbl.c %>% 
    mutate(tp.fn = mapply(label.tp.fn, incorp.known, incorp.pred)) %>%
    filter(! is.na(tp.fn))

tbl.tp.fn = inner_join(tbl.c.tp.fn, tbl.otu.sum, c('taxon' = 'taxon')) 
tbl.tp.fn %>% head

%%R
# how many TP & FN?
tbl.tp.fn %>% 
    group_by(library, tp.fn) %>%
    summarize(n = n())

%%R -h 350
tbl.tp.fn$library = as.character(tbl.tp.fn$library)
ggplot(tbl.tp.fn, aes(library, total_count, color=tp.fn)) +
    geom_boxplot() +
    labs(y='Total count') +
    theme(
        text = element_text(size=18)
        )  

%%R -h 700 -w 900
tbl.tp.fn$library = as.character(tbl.tp.fn$library)
ggplot(tbl.tp.fn, aes(taxon, total_count, group=taxon, color=library)) +
    geom_point(size=3) +
    geom_line() +
    scale_y_log10() +
    labs(y='Total count') +
    facet_grid(. ~ tp.fn, scales='free_x') +
    theme(
        text = element_text(size=18),
        axis.text.x = element_text(angle=90, hjust=1)
        )  

"""
Explanation: Enrichment of TP for abundant incorporators?

What is the abundance distribution of TP and FP?
Are more abundant incorporators being detected more than low abundant taxa
End of explanation
"""


%%R
# OTU total counts

heavy.cut = 1.71

tbl.otu.sum = tbl.otu %>%
    filter(! grepl('inf', fraction)) %>%
    separate(fraction, into = c('BD_min','BD_max'), sep='-', convert=TRUE) %>%
    filter(BD_min >= heavy.cut & BD_max <= 2) %>%
    group_by(library, taxon) %>%
    summarize(total_count = sum(count))
tbl.otu.sum %>% head

%%R

tbl.tp.fn = inner_join(tbl.c.tp.fn, tbl.otu.sum, c('taxon' = 'taxon')) 
tbl.tp.fn %>% head

%%R -h 350
tbl.tp.fn$library = as.character(tbl.tp.fn$library)
ggplot(tbl.tp.fn, aes(library, total_count, color=tp.fn)) +
    geom_boxplot() +
    labs(y='Total count') +
    theme(
        text = element_text(size=18)
        )  

%%R -h 700 -w 900
tbl.tp.fn$library = as.character(tbl.tp.fn$library)
ggplot(tbl.tp.fn, aes(taxon, total_count, group=taxon, color=library)) +
    geom_point(size=3) +
    geom_line() +
    scale_y_log10() +
    labs(y='Total count') +
    facet_grid(. ~ tp.fn, scales='free_x') +
    theme(
        text = element_text(size=18),
        axis.text.x = element_text(angle=90, hjust=1)
        )  

"""
Explanation: Notes:

FNs can be abundant
TP can have a lower abundance in the 'treatment' library

Just looking at the 'heavy' fractions
End of explanation
"""


# building tree structure
from os.path import abspath
nest = nestly.Nest()

##  values
vals = [str(x) for x in range(1,5)]
nest.add('vals', vals)


## input files
nest.add('--np', [1], create_dir=False)

buildDir = '/home/nick/t/nestly/' #os.path.join(workDir, 'vals')
nest.build(buildDir)

%%writefile /home/nick/t/example.sh
#!/bin/bash

export TIME='elapsed,maxmem,exitstatus\n%e,%M,%x'

echo {--np} > {--np}_test.txt

!cd /home/nick/t/; \
    chmod 777 example.sh 

!cd /home/nick/t/; \
    nestrun -j 2 --template-file example.sh -d nestly

"""
Explanation: Notes:

The TPs (for the most part) are dramatically different in abundance between control and treatment


Sandbox
End of explanation
"""


