<jupyter_start><jupyter_text>Accuracy as a function taxon abundance 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<jupyter_code>workDir = '/home/nick/notebook/SIPSim/dev/bac_genome1210/'
sourceDir = os.path.join(workDir, 'atomIncorp_taxaIncorp')
buildDir = os.path.join(workDir, 'atomIncorp_taxaIncorp_abund')
R_dir = '/home/nick/notebook/SIPSim/lib/R/'<jupyter_output><empty_output><jupyter_text>Init<jupyter_code>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)<jupyter_output><empty_output><jupyter_text>Nestly<jupyter_code># 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))

## set params
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}/{percIncorp}/{percTaxa}/{rep}/{comm_file} .
cp {sourceDir}/{percIncorp}/{percTaxa}/{rep}/OTU_n2_abs{abs}_sub-norm_DESeq2 .
cp {sourceDir}/{percIncorp}/{percTaxa}/{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 10 --template-file $bashFile -d atomIncorp_taxaIncorp_abund --log-file log.txt
# aggregating confusion matrix data
## byClass
### dominant
!cd $workDir; \
    nestagg delim \
        -d atomIncorp_taxaIncorp_abund \
        -k percIncorp,percTaxa,rep \
        -o ./atomIncorp_taxaIncorp_abund/DESeq2-cMtx_dom_byClass.csv \
        DESeq2-cMtx_dom_byClass.csv

### rare
!cd $workDir; \
    nestagg delim \
        -d atomIncorp_taxaIncorp_abund \
        -k percIncorp,percTaxa,rep \
        -o ./atomIncorp_taxaIncorp_abund/DESeq2-cMtx_rare_byClass.csv \
        DESeq2-cMtx_rare_byClass.csv<jupyter_output><empty_output><jupyter_text>Plotting accuracy  Dominant taxa<jupyter_code>%%R -i workDir
setwd(workDir)

byClass = read.csv('./atomIncorp_taxaIncorp_abund/DESeq2-cMtx_dom_byClass.csv')
byClass %>% head
%%R -w 700 

col2keep = c('Balanced Accuracy', 'Sensitivity','Specificity')

byClass.f = byClass %>%
    filter(X %in% col2keep) %>%
    mutate(percIncorp = as.character(percIncorp),
           percTaxa = as.character(percTaxa))

byClass.f$percIncorp = factor(byClass.f$percIncorp, levels=sort(unique(as.numeric(byClass.f$percIncorp))))
byClass.f$percTaxa = factor(byClass.f$percTaxa, levels=sort(unique(as.numeric(byClass.f$percTaxa))))


ggplot(byClass.f, aes(percIncorp, byClass, fill=percTaxa)) +
    geom_boxplot(position='dodge') +
    labs(y='Value', x='Atom percent 13C') +
    facet_grid(X ~ .) +
    theme(
        text = element_text(size=16)
        )<jupyter_output><empty_output><jupyter_text>Rare taxa<jupyter_code>%%R -i workDir
setwd(workDir)

byClass = read.csv('./atomIncorp_taxaIncorp_abund/DESeq2-cMtx_rare_byClass.csv')
byClass %>% head
%%R -w 700 

col2keep = c('Balanced Accuracy', 'Sensitivity','Specificity')

byClass.f = byClass %>%
    filter(X %in% col2keep) %>%
    mutate(percIncorp = as.character(percIncorp),
           percTaxa = as.character(percTaxa))

byClass.f$percIncorp = factor(byClass.f$percIncorp, levels=sort(unique(as.numeric(byClass.f$percIncorp))))
byClass.f$percTaxa = factor(byClass.f$percTaxa, levels=sort(unique(as.numeric(byClass.f$percTaxa))))


ggplot(byClass.f, aes(percIncorp, byClass, fill=percTaxa)) +
    geom_boxplot(position='dodge') +
    labs(y='Value', x='Atom percent 13C') +
    facet_grid(X ~ .) +
    theme(
        text = element_text(size=16)
        )<jupyter_output><empty_output><jupyter_text>__Notes:__* Dominants can produce higher sensitivity, but (likely) smaller sample size, so crap shoot on whether dominant taxon is in the optimal GC window for the DESeq2 analysis. * Rares essentially match whole-community accuracy.  * Likely due to sample size of rares (most of community for log-normal abundance distribution)__Questions:__* What is the sample size of dominants and rares?  Sample size of dominants and rares* Method:    * Get number of taxa from dominant/rare-parsed DESeq2 object files<jupyter_code>%%R -i workDir
setwd(workDir)

dom.files = grep('DESeq2_dom$', list.files('./atomIncorp_taxaIncorp_abund/', recursive=TRUE), value=TRUE)
rare.files = grep('DESeq2_rare$', list.files('./atomIncorp_taxaIncorp_abund/', recursive=TRUE), value=TRUE)

read_taxa = function(file.name, dir='atomIncorp_taxaIncorp_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('percIncorp','percTaxa','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$percTaxa = factor(tbl.dom.rare$percTaxa, levels=sort(unique(as.numeric(tbl.dom.rare$percTaxa))))

ggplot(tbl.dom.rare, aes(percIncorp, N.taxa, fill=percTaxa)) +
    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)
        )<jupyter_output><empty_output><jupyter_text>__Notes:__* Figure is showing number of taxa in each group* Dominant n-taxa very low (mean values) with a high variance  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<jupyter_code>%%bash -s "$R_dir"

find atomIncorp_taxaIncorp_abund/ -name "*DESeq2" | \
    perl -pe 's/(.+\/).+/\1/' | \
    xargs -P 8 -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_taxaIncorp_abund/', recursive=TRUE), value=TRUE)


tbl.l = list()
for(f in deseq.files){
    f = paste(c('./atomIncorp_taxaIncorp_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('percIncorp','percTaxa','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(percIncorp ~ percTaxa) +
    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, percTaxa) %>%
    summarize(cor.pearson = cor(mean_rel_abund_perc, TP_FN))

tbl.s.c %>% head
%%R -w 600 -h 350
ggplot(tbl.s.c, aes(percIncorp, cor.pearson, fill=percTaxa)) +
    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)
        )<jupyter_output><empty_output><jupyter_text>__Notes:__* Suggests that the correlation between sensitivity and taxon abundance is dependant on % atom incorp & percent of taxa that incorporated  * If only a few incorporator taxa, sensitivity will be poor regardless of abundance.  * Higher % atom produces weak positive correlation between taxon abundance and correct classification.* Tecnically, values are point-biserial correlation coefficients (not pearson) ***  SANDBOX<jupyter_code>%%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()<jupyter_output><empty_output><jupyter_text>***<jupyter_code>%%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<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.15493674567398574, 'fraction_numerical': 0.013668750908826523, 'mean_word_length': 3.792682926829268, 'pattern_counts': {'":': 0, '<': 40, '<?xml version=': 0, '>': 84, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 1, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '9', 'dir': 'jupyter-structured-clean-dedup', 'id': '244785', 'license': 'MIT', 'n_tokens_mistral': 5677, 'n_tokens_neox': 5353, 'n_words': 971, 'path': 'ipynb/bac_genome/n1210/.ipynb_checkpoints/atomIncorp_taxaIncorp_abund-checkpoint.ipynb', 'repo_name': 'arischwartz/test', 'size': '13754'}