In [1]:
import os
import glob
import re
import nestly
In [2]:
%load_ext rpy2.ipython
%load_ext pushnote
In [3]:
%%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
In [4]:
workDir = '/home/nick/notebook/SIPSim/dev/fullCyc/n1147_frag_norm_9_2.5_n5/'
buildDir = os.path.join(workDir, 'Day1_rep10')
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 = 12
In [5]:
# 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)
### 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')
In [6]:
%%writefile $bashFile
#!/bin/bash
export PATH={R_dir}:$PATH
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
In [7]:
!chmod 777 $bashFile
!cd $workDir; \
nestrun --template-file $bashFile -d Day1_rep10 --log-file log.txt -j 2
In [9]:
%%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')
In [102]:
%%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
In [103]:
%%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)
In [104]:
%%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
In [105]:
%%R
physeq.bulk.n = transform_sample_counts(physeq.bulk, function(x) x/sum(x))
physeq.bulk.n
In [106]:
%%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)
In [107]:
%%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)
In [108]:
OTU_files = !find $buildDir -name "OTU_abs1e9_sub.txt"
#OTU_files = !find $buildDir -name "OTU_abs1e9.txt"
OTU_files
Out[108]:
In [114]:
%%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/', '', x)
#SIM_rep = gsub('/OTU_abs1e9_sub.txt', '', SIM_rep)
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
In [115]:
%%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
In [116]:
# loading comm files
comm_files = !find $buildDir -name "bulk-core_comm_target.txt"
comm_files
Out[116]:
In [117]:
%%R -i comm_files
df.comm = list()
for (f in comm_files){
rep = gsub('.+/Day1_rep10/([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)
In [118]:
%%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)
In [119]:
%%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)
In [120]:
%%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'
)
In [121]:
%%R -i targetFile
df.target = read.delim(targetFile, sep='\t')
df.target %>% nrow %>% print
df.target %>% head(n=3)
In [122]:
%%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'
)
In [123]:
%%R
# 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)
In [124]:
%%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'
)
In [125]:
%%R
join_abund_dists = function(df.EMP.j, df.SIM.j, df.target){
## emperical
df.EMP.j.f = df.EMP.j %>%
filter(abundance > 0) %>%
#filter(!OTU %in% c('OTU.32', 'OTU.2', 'OTU.4')) %>% # TEST
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) %>%
#filter(!taxon %in% c('OTU.32', 'OTU.2', 'OTU.4')) %>% # TEST
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
In [126]:
%%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
In [127]:
%%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)
# status
cat('Number of overlapping taxa:', df.j$OTU %>% unique %>% length, '\n')
# make figure
grid.arrange(p.EMP, p.SIM, ncol=1, heights=c(1,5))
In [128]:
%%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),
median_rel_abund_c = median(rel_abund_c)) %>%
ungroup()
return(df)
}
df.j.cm = center_mass(df.j)
In [129]:
%%R -w 650
# 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),
median_rel_abund_c = first(median_rel_abund_c)) %>%
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),
median_rel_abund_c = first(median_rel_abund_c)) %>%
ungroup()
# check
cat('Number of OTUs:', df.j.cm.s$OTU %>% unique %>% length, '\n')
# plotting
ggplot(df.j.cm.s, aes(emperical, simulated, color=median_rel_abund_c,
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.705, 1.74)) +
scale_color_gradient(trans='log') +
labs(title='Center of mass') +
theme_bw() +
theme(
text = element_text(size=16)
)
In [130]:
%%R
df.j.cm.s.f = df.j.cm.s %>%
mutate(CM_diff = emperical - simulated)
ggplot(df.j.cm.s.f, aes(median_rel_abund_c, CM_diff)) +
geom_point() +
scale_x_log10() +
labs(x='Relative abundance', y='Center of mass (Emperical - Simulated)', title='Center of mass') +
theme_bw() +
theme(
text = element_text(size=16)
)
possibilities: