<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 Goal
Simulating fullCyc Day1 control gradients
Not simulating incorporation (all 0% isotope incorp.)
Don't know how much true incorporatation for emperical data
Richness = genome reference pool size (n=1147)
Simulating 10 replicate gradients to assess simulation stochasticity.
Using a relative large bandwidth to create fairly smooth distributions
Init
Step1:  Nestly
assuming fragments already simulated
Step2:  BD min/max
what is the min/max BD that we care about?
Step3:  Loading data
Emperical
SIP data
Step4:  bulk soil samples
Step5:  Simulated
Step6:  'bulk soil' community files
Step7:  BD span of just overlapping taxa (redundant; but consistent wiht other notebooks)
Taxa overlapping between emperical data and genomes in dataset
These taxa should have the same relative abundances in both datasets.
The comm file was created from the emperical dataset phyloseq file.
Step8:  Check
Are all target (overlapping) taxa the same relative abundances in both datasets?
Step9:  Correlation between relative abundance and BD_range diff
Are low abundant taxa more variable in their BD span
Step10:  Notes
between Day1_rep10, Day1_richFromTarget_rep10, and Day1_add_Rich_rep10
Step11:  Check
Step12:  Notes
The abundance distributions of the overlapping OTUs look pretty similar between 'absolute' and 'relative' (post-PCR & post-sequencing simulation). 
The difference between absolute and relative are probably due to the PCR simulation
Calculating center of mass for overlapping taxa
weighted mean BD, where weights are relative abundances
Step13:  Notes
Error bars are stdev of simulation reps
The center of mass for most OTUs is shifted heavier in the simulations vs the emperical data
Also, most have a mean simulated BD of ~1.73
This is approximately the middle of the BD range (1.725)
It suggests that many of the taxa are equally dispersed across the gradient
Step14:  R^2 for each SIM_rep
Step15:  Notes
Step16:  Notes
The center of mass seems to be correct
Step17:  Notes
Step18:  Notes
Odd... these outlier OTUs are simulated from genomes with high G+C.
I would have expected the OTUs to be shifted right in the emperical data.
Are the simulations over-shifting the distributions of heavy G+C taxa?
Plotting rep genome G+C of all overlapping taxa
<ASSISTANT_TASK:>
Python Code:
import os
import glob
import re
import nestly
%load_ext rpy2.ipython
%load_ext pushnote
%%R
library(ggplot2)
library(dplyr)
library(tidyr)
library(gridExtra)
library(phyloseq)
## BD for G+C of 0 or 100
BD.GCp0 = 0 * 0.098 + 1.66
BD.GCp100 = 1 * 0.098 + 1.66
Explanation: Goal
Simulating fullCyc Day1 control gradients
Not simulating incorporation (all 0% isotope incorp.)
Don't know how much true incorporatation for emperical data
Richness = genome reference pool size (n=1147)
Simulating 10 replicate gradients to assess simulation stochasticity.
Using a relative large bandwidth to create fairly smooth distributions
Init
End of explanation
workDir = '/home/nick/notebook/SIPSim/dev/fullCyc/n1147_frag_norm_9_2.5_n5/'
buildDir = os.path.join(workDir, 'Day1_rep10_justOverlap')
R_dir = '/home/nick/notebook/SIPSim/lib/R/'
fragFile= '/home/nick/notebook/SIPSim/dev/bac_genome1147/validation/ampFrags.pkl'
targetFile = '/home/nick/notebook/SIPSim/dev/fullCyc/CD-HIT/target_taxa.txt'
physeqDir = '/var/seq_data/fullCyc/MiSeq_16SrRNA/515f-806r/lib1-7/phyloseq/'
physeq_bulkCore = 'bulk-core'
physeq_SIP_core = 'SIP-core_unk'
nreps = 10
prefrac_comm_abundance = '1e9'
seq_per_fraction = ['lognormal', 9.432, 0.5, 10000, 30000] # dist, mean, scale, min, max
bulk_days = [1]
nprocs = 14
# building tree structure
nest = nestly.Nest()
    ## varying params
nest.add('rep', [x + 1 for x in xrange(nreps)])
## set params
nest.add('bulk_day', bulk_days, create_dir=False)
nest.add('abs', [prefrac_comm_abundance], create_dir=False)
nest.add('percIncorp', [0], create_dir=False)
nest.add('percTaxa', [0], create_dir=False)
nest.add('np', [nprocs], create_dir=False)
nest.add('subsample_dist', [seq_per_fraction[0]], create_dir=False)
nest.add('subsample_mean', [seq_per_fraction[1]], create_dir=False)
nest.add('subsample_scale', [seq_per_fraction[2]], create_dir=False)
nest.add('subsample_min', [seq_per_fraction[3]], create_dir=False)
nest.add('subsample_max', [seq_per_fraction[4]], create_dir=False)
nest.add('bandwidth', [0.6], create_dir=False)
nest.add('cmd', ['"print unless /NA$/"'], create_dir=False)
### input/output files
nest.add('buildDir', [buildDir], create_dir=False)
nest.add('R_dir', [R_dir], create_dir=False)
nest.add('fragFile', [fragFile], create_dir=False)
nest.add('targetFile', [targetFile], create_dir=False)
nest.add('physeqDir', [physeqDir], create_dir=False)
nest.add('physeq_bulkCore', [physeq_bulkCore], create_dir=False)
# building directory tree
nest.build(buildDir)
# bash file to run
bashFile = os.path.join(buildDir, 'SIPSimRun.sh')
%%writefile $bashFile
#!/bin/bash
export PATH={R_dir}:$PATH
#-- making DNA pool similar to gradient of interest
echo '# Creating comm file from phyloseq'
phyloseq2comm.r {physeqDir}{physeq_bulkCore} -s 12C-Con -d {bulk_day} > {physeq_bulkCore}_comm.txt
printf 'Number of lines: '; wc -l {physeq_bulkCore}_comm.txt
echo '## Adding target taxa to comm file'
comm_add_target.r {physeq_bulkCore}_comm.txt {targetFile} > {physeq_bulkCore}_comm_target.txt
printf 'Number of lines: '; wc -l {physeq_bulkCore}_comm_target.txt
echo '## Selecting just target taxa'
perl -ne {cmd} {physeq_bulkCore}_comm_target.txt | comm_set_abund.r - > tmp.txt
rm -f {physeq_bulkCore}_comm_target.txt
mv tmp.txt {physeq_bulkCore}_comm_target.txt
echo '## parsing out genome fragments to make simulated DNA pool resembling the gradient of interest'
## all OTUs without an associated reference genome will be assigned a random reference (of the reference genome pool)
### this is done through --NA-random
SIPSim fragment_KDE_parse {fragFile} {physeq_bulkCore}_comm_target.txt \
    --rename taxon_name --NA-random > fragsParsed.pkl
echo '#-- SIPSim pipeline --#'
echo '# converting fragments to KDE'
SIPSim fragment_KDE \
    fragsParsed.pkl \
    > fragsParsed_KDE.pkl
    
echo '# adding diffusion'
SIPSim diffusion \
    fragsParsed_KDE.pkl \
    --bw {bandwidth} \
    --np {np} \
    > fragsParsed_KDE_dif.pkl    
echo '# adding DBL contamination'
SIPSim DBL \
    fragsParsed_KDE_dif.pkl \
    --bw {bandwidth} \
    --np {np} \
    > fragsParsed_KDE_dif_DBL.pkl      
    
echo '# making incorp file'
SIPSim incorpConfigExample \
  --percTaxa {percTaxa} \
  --percIncorpUnif {percIncorp} \
  > {percTaxa}_{percIncorp}.config
echo '# adding isotope incorporation to BD distribution'
SIPSim isotope_incorp \
    fragsParsed_KDE_dif_DBL.pkl \
    {percTaxa}_{percIncorp}.config \
    --comm {physeq_bulkCore}_comm_target.txt \
    --bw {bandwidth} \
    --np {np} \
    > fragsParsed_KDE_dif_DBL_inc.pkl
    
echo '# simulating gradient fractions'
SIPSim gradient_fractions \
    {physeq_bulkCore}_comm_target.txt \
    > fracs.txt 
echo '# simulating an OTU table'
SIPSim OTU_table \
    fragsParsed_KDE_dif_DBL_inc.pkl \
    {physeq_bulkCore}_comm_target.txt \
    fracs.txt \
    --abs {abs} \
    --np {np} \
    > OTU_abs{abs}.txt
    
#echo '# simulating PCR'
#SIPSim OTU_PCR \
#    OTU_abs{abs}.txt \
#    > OTU_abs{abs}_PCR.txt    
    
echo '# subsampling from the OTU table (simulating sequencing of the DNA pool)'
SIPSim OTU_subsample \
    --dist {subsample_dist} \
    --dist_params mean:{subsample_mean},sigma:{subsample_scale} \
    --min_size {subsample_min} \
    --max_size {subsample_max} \
    OTU_abs{abs}.txt \
    > OTU_abs{abs}_sub.txt
        
echo '# making a wide-formatted table'
SIPSim OTU_wideLong -w \
    OTU_abs{abs}_sub.txt \
    > OTU_abs{abs}_sub_w.txt
    
echo '# making metadata (phyloseq: sample_data)'
SIPSim OTU_sampleData \
    OTU_abs{abs}_sub.txt \
    > OTU_abs{abs}_sub_meta.txt
!chmod 777 $bashFile
!cd $workDir; \
    nestrun  --template-file $bashFile -d Day1_rep10_justOverlap --log-file log.txt -j 2
%pushnote Day1_rep10_justOverlap complete
Explanation: Nestly
assuming fragments already simulated
End of explanation
%%R
## min G+C cutoff
min_GC = 13.5
## max G+C cutoff
max_GC = 80
## max G+C shift
max_13C_shift_in_BD = 0.036
min_BD = min_GC/100.0 * 0.098 + 1.66    
max_BD = max_GC/100.0 * 0.098 + 1.66    
max_BD = max_BD + max_13C_shift_in_BD
cat('Min BD:', min_BD, '\n')
cat('Max BD:', max_BD, '\n')
Explanation: BD min/max
what is the min/max BD that we care about?
End of explanation
%%R 
# simulated OTU table file
OTU.table.dir = '/home/nick/notebook/SIPSim/dev/fullCyc/frag_norm_9_2.5_n5/Day1_default_run/1e9/'
OTU.table.file = 'OTU_abs1e9_PCR_sub.txt'
#OTU.table.file = 'OTU_abs1e9_sub.txt'
#OTU.table.file = 'OTU_abs1e9.txt'
%%R -i physeqDir -i physeq_SIP_core -i bulk_days
# bulk core samples
F = file.path(physeqDir, physeq_SIP_core)
physeq.SIP.core = readRDS(F) 
physeq.SIP.core.m = physeq.SIP.core %>% sample_data
physeq.SIP.core = prune_samples(physeq.SIP.core.m$Substrate == '12C-Con' & 
                                physeq.SIP.core.m$Day %in% bulk_days, 
                                physeq.SIP.core) %>%
    filter_taxa(function(x) sum(x) > 0, TRUE)
physeq.SIP.core.m = physeq.SIP.core %>% sample_data        
physeq.SIP.core
%%R 
## dataframe
df.EMP = physeq.SIP.core %>% otu_table %>%
    as.matrix %>% as.data.frame
df.EMP$OTU = rownames(df.EMP)
df.EMP = df.EMP %>%    
    gather(sample, abundance, 1:(ncol(df.EMP)-1)) 
df.EMP = inner_join(df.EMP, physeq.SIP.core.m, c('sample' = 'X.Sample')) 
df.EMP.nt = df.EMP %>%
    group_by(sample) %>%
    mutate(n_taxa = sum(abundance > 0)) %>%
    ungroup() %>%
    distinct(sample) %>%
    filter(Buoyant_density >= min_BD, 
           Buoyant_density <= max_BD)
    
df.EMP.nt %>% head(n=3)
Explanation: Loading data
Emperical
SIP data
End of explanation
%%R
physeq.dir = '/var/seq_data/fullCyc/MiSeq_16SrRNA/515f-806r/lib1-7/phyloseq/'
physeq.bulk = 'bulk-core'
physeq.file = file.path(physeq.dir, physeq.bulk)
physeq.bulk = readRDS(physeq.file)
physeq.bulk.m = physeq.bulk %>% sample_data
physeq.bulk = prune_samples(physeq.bulk.m$Exp_type == 'microcosm_bulk' &
                            physeq.bulk.m$Day %in% bulk_days, physeq.bulk)
physeq.bulk.m = physeq.bulk %>% sample_data
physeq.bulk
%%R
physeq.bulk.n = transform_sample_counts(physeq.bulk, function(x) x/sum(x))
physeq.bulk.n
%%R
# making long format of each bulk table
bulk.otu = physeq.bulk.n %>% otu_table %>% as.data.frame
ncol = ncol(bulk.otu)
bulk.otu$OTU = rownames(bulk.otu)
bulk.otu = bulk.otu %>%
    gather(sample, abundance, 1:ncol) 
bulk.otu = inner_join(physeq.bulk.m, bulk.otu, c('X.Sample' = 'sample')) %>%
    dplyr::select(OTU, abundance) %>%
    rename('bulk_abund' = abundance)
bulk.otu %>% head(n=3)
%%R
# joining tables
df.EMP.j = inner_join(df.EMP, bulk.otu, c('OTU' = 'OTU')) %>%
    filter(Buoyant_density >= min_BD, 
           Buoyant_density <= max_BD) 
    
df.EMP.j %>% head(n=3)
Explanation: bulk soil samples
End of explanation
OTU_files = !find $buildDir -name "OTU_abs1e9_sub.txt"
OTU_files
%%R -i OTU_files
# loading files
df.SIM = list()
for (x in OTU_files){
    SIM_rep = gsub('/home/nick/notebook/SIPSim/dev/fullCyc/n1147_frag_norm_9_2.5_n5/Day1_rep10_justOverlap/', '', x)
    SIM_rep = gsub('/OTU_abs1e9_sub.txt', '', SIM_rep)
    df.SIM[[SIM_rep]] = read.delim(x, sep='\t') 
    }
df.SIM = do.call('rbind', df.SIM)
df.SIM$SIM_rep = gsub('\\.[0-9]+$', '', rownames(df.SIM))
rownames(df.SIM) = 1:nrow(df.SIM)
df.SIM %>% head
%%R
## edit table
df.SIM.nt = df.SIM %>%
    filter(count > 0) %>%
    group_by(SIM_rep, library, BD_mid) %>%
    summarize(n_taxa = n()) %>%
    filter(BD_mid >= min_BD, 
           BD_mid <= max_BD)
df.SIM.nt %>% head    
Explanation: Simulated
End of explanation
# loading comm files
comm_files = !find $buildDir -name "bulk-core_comm_target.txt"
comm_files
%%R -i comm_files
df.comm = list()
for (f in comm_files){
    rep = gsub('.+/Day1_rep10_justOverlap/([0-9]+)/.+', '\\1', f)
    df.comm[[rep]] = read.delim(f, sep='\t') %>%
        dplyr::select(library, taxon_name, rel_abund_perc) %>%
        rename('bulk_abund' = rel_abund_perc) %>%
        mutate(bulk_abund = bulk_abund / 100)
}
df.comm = do.call('rbind', df.comm)
df.comm$SIM_rep = gsub('\\.[0-9]+$', '', rownames(df.comm))
rownames(df.comm) = 1:nrow(df.comm)
df.comm %>% head(n=3)
%%R
## joining tables
df.SIM.j = inner_join(df.SIM, df.comm, c('SIM_rep' = 'SIM_rep',
                                         'library' = 'library',
                                         'taxon' = 'taxon_name')) %>%
    filter(BD_mid >= min_BD, 
           BD_mid <= max_BD)
    
df.SIM.j %>% head(n=3)
%%R 
# filtering & combining emperical w/ simulated data
## emperical 
max_BD_range = max(df.EMP.j$Buoyant_density) - min(df.EMP.j$Buoyant_density)
df.EMP.j.f = df.EMP.j %>%
    filter(abundance > 0) %>%
    group_by(OTU) %>%
    summarize(mean_rel_abund = mean(bulk_abund),
              min_BD = min(Buoyant_density),
              max_BD = max(Buoyant_density),
              BD_range = max_BD - min_BD,
              BD_range_perc = BD_range / max_BD_range * 100) %>%
    ungroup() %>%
    mutate(dataset = 'emperical',
           SIM_rep = NA)
## simulated
max_BD_range = max(df.SIM.j$BD_mid) - min(df.SIM.j$BD_mid)
df.SIM.j.f = df.SIM.j %>%
    filter(count > 0) %>%
    group_by(SIM_rep, taxon) %>%
    summarize(mean_rel_abund = mean(bulk_abund),
              min_BD = min(BD_mid),
              max_BD = max(BD_mid),
              BD_range = max_BD - min_BD,
              BD_range_perc = BD_range / max_BD_range * 100) %>%
    ungroup() %>%
    rename('OTU' = taxon) %>%
    mutate(dataset = 'simulated')
## join
df.j = rbind(df.EMP.j.f, df.SIM.j.f) %>%
    filter(BD_range_perc > 0,
            mean_rel_abund > 0)
df.j$SIM_rep = reorder(df.j$SIM_rep, df.j$SIM_rep %>% as.numeric)
df.j %>% head(n=3)
%%R -h 400
## plotting
ggplot(df.j, aes(mean_rel_abund, BD_range_perc, color=SIM_rep)) +
    geom_point(alpha=0.3) +
    scale_x_log10() +
    scale_y_continuous() +
    labs(x='Pre-fractionation abundance', y='% of total BD range') +
    facet_grid(dataset ~ .) +
    theme_bw() +
    theme(
        text = element_text(size=16),
        panel.grid = element_blank()#,
        #legend.position = 'none'
        )
Explanation: 'bulk soil' community files
End of explanation
%%R -i targetFile
df.target = read.delim(targetFile, sep='\t')
df.target %>% nrow %>% print
df.target %>% head(n=3)
%%R
# filtering to just target taxa
df.j.t = df.j %>% 
    filter(OTU %in% df.target$OTU) 
df.j %>% nrow %>% print
df.j.t %>% nrow %>% print
## plotting
ggplot(df.j.t, aes(mean_rel_abund, BD_range_perc, color=SIM_rep)) +
    geom_point(alpha=0.5, shape='O') +
    scale_x_log10() +
    scale_y_continuous() +
    #scale_color_manual(values=c('blue', 'red')) +
    labs(x='Pre-fractionation abundance', y='% of total BD range') +
    facet_grid(dataset ~ .) +
    theme_bw() +
    theme(
        text = element_text(size=16),
        panel.grid = element_blank()#,
        #legend.position = 'none'
        )
Explanation: BD span of just overlapping taxa (redundant; but consistent wiht other notebooks)
Taxa overlapping between emperical data and genomes in dataset
These taxa should have the same relative abundances in both datasets.
The comm file was created from the emperical dataset phyloseq file.
End of explanation
%%R -w 600 -h 500
# formatting data
df.1 = df.j.t %>% 
    filter(dataset == 'simulated') %>%
    select(SIM_rep, OTU, mean_rel_abund, BD_range, BD_range_perc)
df.2 = df.j.t %>%
    filter(dataset == 'emperical') %>%
    select(SIM_rep, OTU, mean_rel_abund, BD_range, BD_range_perc)
df.12 = inner_join(df.1, df.2, c('OTU' = 'OTU')) %>%
    mutate(BD_diff_perc = BD_range_perc.y - BD_range_perc.x)
df.12$SIM_rep.x = reorder(df.12$SIM_rep.x, df.12$SIM_rep.x  %>% as.numeric)
## plotting
p1 = ggplot(df.12, aes(mean_rel_abund.x, mean_rel_abund.y)) +
    geom_point(alpha=0.5) +
    scale_x_log10() +
    scale_y_log10() +
    labs(x='Relative abundance (simulated)', y='Relative abundance (emperical)') +
    facet_wrap(~ SIM_rep.x)
    theme_bw() +
    theme(
        text = element_text(size=16),
        panel.grid = element_blank(),
        legend.position = 'none'
        )
p1
Explanation: Check
Are all target (overlapping) taxa the same relative abundances in both datasets?
End of explanation
%%R -w 800 -h 500
ggplot(df.12, aes(mean_rel_abund.x, BD_diff_perc)) +
    geom_point(alpha=0.5) +
    scale_x_log10() +
    labs(x='Pre-fractionation relative abundance', 
         y='Difference in % of gradient spanned\n(emperical - simulated)',
         title='Overlapping taxa') +
    facet_wrap(~ SIM_rep.x) +
    theme_bw() +
    theme(
        text = element_text(size=16),
        panel.grid = element_blank(),
        legend.position = 'none'
        )
Explanation: Correlation between relative abundance and BD_range diff
Are low abundant taxa more variable in their BD span
End of explanation
%%R
join_abund_dists = function(df.EMP.j, df.SIM.j, df.target){
    ## emperical 
    df.EMP.j.f = df.EMP.j %>%
        filter(abundance > 0) %>%
        dplyr::select(OTU, sample, abundance, Buoyant_density, bulk_abund) %>%
        mutate(dataset = 'emperical', SIM_rep = NA) %>%
        filter(OTU %in% df.target$OTU) 
    
    ## simulated
    df.SIM.j.f = df.SIM.j %>%
        filter(count > 0) %>%
        dplyr::select(taxon, fraction, count, BD_mid, bulk_abund, SIM_rep) %>%
        rename('OTU' = taxon,
               'sample' = fraction,
               'Buoyant_density' = BD_mid,
               'abundance' = count) %>%
        mutate(dataset = 'simulated') %>%
        filter(OTU %in% df.target$OTU) 
    
    ## getting just intersecting OTUs
    OTUs.int = intersect(df.EMP.j.f$OTU, df.SIM.j.f$OTU)
    
    df.j = rbind(df.EMP.j.f, df.SIM.j.f) %>%
        filter(OTU %in% OTUs.int) %>%
        group_by(sample) %>%
        mutate(rel_abund = abundance / sum(abundance))
    
    cat('Number of overlapping OTUs between emperical & simulated:', 
            df.j$OTU %>% unique %>% length, '\n\n')
    return(df.j)
    }
df.j = join_abund_dists(df.EMP.j, df.SIM.j, df.target)
df.j %>% head(n=3) %>% as.data.frame                  
%%R
# closure operation
df.j = df.j %>%
    ungroup() %>%
    mutate(SIM_rep = SIM_rep %>% as.numeric) %>%
    group_by(dataset, SIM_rep, sample) %>%
    mutate(rel_abund_c = rel_abund / sum(rel_abund)) %>%
    ungroup()
df.j %>% head(n=3) %>% as.data.frame
%%R -h 1500 -w 800
# plotting 
plot_abunds = function(df){
    p = ggplot(df, aes(Buoyant_density, rel_abund_c, fill=OTU)) +
        geom_area(stat='identity', position='dodge', alpha=0.5) +
        labs(x='Buoyant density', 
             y='Subsampled community\n(relative abundance for subset taxa)') +
        theme_bw() +
        theme( 
            text = element_text(size=16),
            legend.position = 'none',
              axis.title.y = element_text(vjust=1),        
              axis.title.x = element_blank(),
              plot.margin=unit(c(0.1,1,0.1,1), "cm")
             )
    return(p)
    }
# simulations
df.j.f = df.j %>%
    filter(dataset == 'simulated')
p.SIM = plot_abunds(df.j.f)
p.SIM = p.SIM + facet_grid(SIM_rep ~ .)
# emperical
df.j.f = df.j %>%
    filter(dataset == 'emperical')
p.EMP = plot_abunds(df.j.f)
# make figure
grid.arrange(p.EMP, p.SIM, ncol=1, heights=c(1,5))
Explanation: Notes
between Day1_rep10, Day1_richFromTarget_rep10, and Day1_add_Rich_rep10:
Day1_rep10 has the most accurate representation of BD span (% of gradient spanned by taxa).
Accuracy drops at ~1e-3 to ~5e-4, but this is caused by detection limits (veil-line effect).
Comparing abundance distributions of overlapping taxa
End of explanation
OTU_files = !find $buildDir -name "OTU_abs1e9.txt"
OTU_files
%%R -i OTU_files
# loading files
df.SIM.abs = list()
for (x in OTU_files){
    SIM_rep = gsub('/home/nick/notebook/SIPSim/dev/fullCyc/n1147_frag_norm_9_2.5_n5/Day1_rep10/', '', x)
    SIM_rep = gsub('/OTU_abs1e9.txt', '', SIM_rep)
    df.SIM.abs[[SIM_rep]] = read.delim(x, sep='\t') 
    }
df.SIM.abs = do.call('rbind', df.SIM.abs)
df.SIM.abs$SIM_rep = gsub('\\.[0-9]+$', '', rownames(df.SIM.abs))
rownames(df.SIM.abs) = 1:nrow(df.SIM.abs)
df.SIM.abs %>% head
%%R
# subset just overlapping taxa
# & closure operation
df.SIM.abs.t = df.SIM.abs %>%
    filter(taxon %in% df.target$OTU) %>%
    group_by(SIM_rep, fraction) %>%
    mutate(rel_abund_c = count / sum(count)) %>%
    rename('Buoyant_density' = BD_mid,
           'OTU' = taxon)
df.SIM.abs.t %>% head(n=3) %>% as.data.frame
%%R -w 800 -h 1200
# plotting
p.abs = plot_abunds(df.SIM.abs.t) 
p.abs + facet_grid(SIM_rep ~ .)
Explanation: Check: plotting closure abs-abunds of overlapping taxa
The overlapping taxa should have the same closure-transformed relative abundances for both:
absolute abundances (OTU table)
relative abundances (subsampled OTU table; as above)
Loading OTU table (abs abunds)
End of explanation
%%R
center_mass = function(df){
    df = df %>%
        group_by(dataset, SIM_rep, OTU) %>%
        summarize(center_mass = weighted.mean(Buoyant_density, rel_abund_c, na.rm=T)) %>%
        ungroup()
    return(df)
}
df.j.cm = center_mass(df.j) 
%%R
# getting mean cm for all SIM_reps
df.j.cm.s = df.j.cm %>%
    group_by(dataset, OTU) %>%
    summarize(mean_cm = mean(center_mass, na.rm=T),
              stdev_cm = sd(center_mass)) %>%
    ungroup() %>%
    spread(dataset, mean_cm) %>%
    group_by(OTU) %>%
    summarize(stdev_cm = mean(stdev_cm, na.rm=T),
              emperical = mean(emperical, na.rm=T),
              simulated = mean(simulated, na.rm=T)) %>%
    ungroup()
# check
cat('Number of OTUs:', df.j.cm.s$OTU %>% unique %>% length, '\n')
# plotting
ggplot(df.j.cm.s, aes(emperical, simulated,
                    ymin = simulated - stdev_cm,
                    ymax = simulated + stdev_cm)) +
    geom_pointrange() +
    stat_function(fun = function(x) x, linetype='dashed', alpha=0.5, color='red') +
    scale_x_continuous(limits=c(1.69, 1.74)) +
    scale_y_continuous(limits=c(1.7, 1.75)) +
    labs(title='Center of mass') +
    theme_bw() +
    theme(
        text = element_text(size=16)
    )
Explanation: Notes
The abundance distributions of the overlapping OTUs look pretty similar between 'absolute' and 'relative' (post-PCR & post-sequencing simulation). 
The difference between absolute and relative are probably due to the PCR simulation
Calculating center of mass for overlapping taxa
weighted mean BD, where weights are relative abundances
End of explanation
%%R
BD_MIN = df.j$Buoyant_density %>% min 
BD_MAX = df.j$Buoyant_density %>% max
BD_AVE = mean(c(BD_MIN, BD_MAX))
print(c(BD_MIN, BD_AVE, BD_MAX))
Explanation: Notes
Error bars are stdev of simulation reps
The center of mass for most OTUs is shifted heavier in the simulations vs the emperical data
Also, most have a mean simulated BD of ~1.73
This is approximately the middle of the BD range (1.725)
It suggests that many of the taxa are equally dispersed across the gradient
End of explanation
%%R
# formatting table
df.j.cm.j = inner_join(df.j.cm %>% 
                           filter(dataset == 'simulated') %>%
                           rename('cm_SIM' = center_mass),
                       df.j.cm %>% 
                           filter(dataset == 'emperical') %>%
                           rename('cm_EMP' = center_mass),
                       c('OTU' = 'OTU')) %>%
    select(-starts_with('dataset'))
df.j.cm.j %>% head
%%R -w 300 -h 400
# lm()
df.j.cm.j.lm = df.j.cm.j %>%
    group_by(SIM_rep.x) %>%
    do(fit = lm(cm_EMP ~ cm_SIM, data = .)) %>%
    mutate(R2 =  summary(fit)$coeff[2],
           data = 'simulated')
#df.j.cm.j.lm %>% head
# plotting
ggplot(df.j.cm.j.lm, aes(data, R2)) +
    geom_boxplot() +
    geom_jitter(height=0, width=0.2, color='red') +
    labs(y='R^2', title='simulated ~ emperical') +
    theme_bw() +
    theme(
        text = element_text(size=16)
    )
Explanation: R^2 for each SIM_rep
End of explanation
%%R -h 1100 -w 800
# cutoff on which OTU are major outliers (varying between simulated and emperical)
BD.diff.cut = 0.02
# which OTU to plot?
df.j.cm.s.f = df.j.cm.s %>% 
    mutate(cm_diff = abs(emperical - simulated)) %>%
    filter(cm_diff > BD.diff.cut, ! is.na(simulated)) 
print('OTUs:')
print(df.j.cm.s.f$OTU)
# filtering to just target taxon
## Simulated
df.j.f = df.j %>%
    filter(dataset == 'simulated', 
           OTU %in% df.j.cm.s.f$OTU) 
p.SIM = plot_abunds(df.j.f)
p.SIM = p.SIM + facet_grid(SIM_rep ~ .)
## Emperical
df.j.f = df.j %>%
    filter(dataset == 'emperical', 
           OTU %in% df.j.cm.s.f$OTU)
p.EMP = plot_abunds(df.j.f)
# make figure
grid.arrange(p.EMP, p.SIM, ncol=1, heights=c(1,5))
Explanation: Notes:
Pretty low R^2, with 1 simulation rep at almost 0
Plotting abundances of some outlier taxa
why is simulated so different in center of mass vs emperical?
End of explanation
%%R
# subset outliers
df.SIM.abs.t = df.SIM.abs %>%
    filter(taxon %in% df.target$OTU) %>%
    group_by(SIM_rep, fraction) %>%
    mutate(rel_abund_c = count / sum(count)) %>%
    rename('Buoyant_density' = BD_mid,
           'OTU' = taxon) %>%
    filter(OTU %in% df.j.cm.s.f$OTU) 
df.SIM.abs.t %>% head(n=3) %>% as.data.frame
%%R -w 800 -h 1200
# plotting
p.abs = plot_abunds(df.SIM.abs.t) 
p.abs + facet_grid(SIM_rep ~ .)
Explanation: Notes
The center of mass seems to be correct: most taxa are shifted heavier relative to the emperical data
Check: which taxon is the highly abundant one?
why is it shifted so far to 'heavy'
the abundance distribution mode is ~1.73, which is a G+C of ~0.71
Abundant taxon: OTU.32
  * rep genome: Pseudonocardia_dioxanivorans_CB1190.fna
    * genome GC = 73.31
    * do the G+C contents of this amplicon region vary among taxa in the genus?
    * No, genome G+C for all 7 genomes in ncbi are ~70 
Check: plotting 'absolute' abudance distributions for major CM outliers
End of explanation
%%R -w 800
genomes = df.target %>% 
    filter(OTU %in% df.SIM.abs.t$OTU) 
df.genInfo = read.delim('/var/seq_data/ncbi_db/genome/Jan2016/bac_complete_spec-rep1_rn/genome_info.txt')
df.genInfo.f = df.genInfo %>% 
    filter(seq_name %in% genomes$genome_seqID) %>%
    mutate(genome_ID = gsub('\\.fna', '', file_name))
df.genInfo.f$genome_ID = reorder(df.genInfo.f$genome_ID, df.genInfo.f$total_GC)
# plotting
ggplot(df.genInfo.f, aes(genome_ID, total_GC)) +
    geom_point() +
    theme_bw() +
    theme(
        text = element_text(size=16),
        axis.text.x = element_text(angle=50, hjust=1)
    )
Explanation: Notes:
simulated absolute abundances seem similar to the simulated relative abundances: both shifted heavy vs emperical
What genomes are these outliers?
End of explanation
%%R -w 800
df.genInfo = read.delim('/var/seq_data/ncbi_db/genome/Jan2016/bac_complete_spec-rep1_rn/genome_info.txt')
df.genInfo.f = df.genInfo %>% 
    filter(seq_name %in% df.target$genome_seqID) %>%
    mutate(genome_ID = gsub('\\.fna', '', file_name))
df.genInfo.f$genome_ID = reorder(df.genInfo.f$genome_ID, df.genInfo.f$total_GC)
# plotting
ggplot(df.genInfo.f, aes(genome_ID, total_GC)) +
    geom_point() +
    theme_bw() +
    theme(
        text = element_text(size=16),
        axis.text.x = element_blank()
    )
Explanation: Notes
Odd... these outlier OTUs are simulated from genomes with high G+C.
I would have expected the OTUs to be shifted right in the emperical data.
Are the simulations over-shifting the distributions of heavy G+C taxa?
Plotting rep genome G+C of all overlapping taxa
End of explanation
<END_TASK>