In [55]:
workDir = '/ebio/abt3_projects/methanogen_host_evo/SIPSim_pt2/data/bac_genome1147/microBetaDiv/'
In [56]:
library(dplyr)
library(tidyr)
library(ggplot2)
as.Num = function(x) x %>% as.character %>% as.numeric
In [57]:
# classifying true positives, neg, ...
clsfy = function(guess,known){
if(is.na(guess) | is.na(known)){
return(NA)
}
if(guess == TRUE){
if(guess == known){
return('True positive')
} else {
return('False positive')
}
} else
if(guess == FALSE){
if(guess == known){
return('True negative')
} else {
return('False negative')
}
} else {
stop('Error: true or false needed')
}
}
In [58]:
# files on simulation accuracy
files = list.files(path=workDir, pattern='*-cMtx_byClass.txt', full.names=TRUE)
files
In [59]:
# combining files
df_byClass = list()
for (f in files){
ff = strsplit(f, '/') %>% unlist
fff = ff[length(ff)]
df_byClass[[fff]] = read.delim(f, sep='\t')
}
df_byClass = do.call(rbind, df_byClass)
df_byClass$file = gsub('\\.[0-9]+$', '', rownames(df_byClass))
df_byClass$method = gsub('-.+', '', df_byClass$file)
rownames(df_byClass) = 1:nrow(df_byClass)
df_byClass %>% head(n=3)
In [60]:
# renaming methods
rename = data.frame(method = c('DESeq2', 'DESeq2_multi', 'heavy', 'qSIP'),
method_new = c('HR-SIP', 'MW-HR-SIP', 'Heavy-SIP', 'q-SIP'))
df_byClass = inner_join(df_byClass, rename, c('method'='method')) %>%
select(-method) %>%
rename('method' = method_new)
df_byClass %>% head(n=3)
In [61]:
# summarize by SIPSim rep & library rep
df_byClass.s = df_byClass %>%
group_by(method, shared_perc, perm_perc, variables) %>%
summarize(mean_value = mean(values, na.rm=TRUE),
sd_value = sd(values, na.rm=TRUE))
# plotting
options(repr.plot.width=8, repr.plot.height=5)
p = ggplot(df_byClass.s, aes(variables, mean_value, color=method,
ymin=mean_value-sd_value,
ymax=mean_value+sd_value)) +
geom_pointrange(alpha=0.8, size=0.2) +
labs(y='Value') +
facet_grid(perm_perc ~ shared_perc) +
theme_bw() +
theme(
axis.title.x = element_blank(),
axis.text.x = element_text(angle=65, hjust=1)
)
plot(p)
In [62]:
# summarize by SIPSim rep & library rep
vars = c('Balanced Accuracy', 'Sensitivity', 'Specificity')
df_byClass.s.f = df_byClass.s %>%
filter(variables %in% vars) %>%
ungroup() %>%
mutate(perm_perc = perm_perc %>% as.character,
perm_perc = perm_perc %>% reorder(perm_perc %>% as.numeric))
# plotting
options(repr.plot.width=9, repr.plot.height=5)
p.pnt = ggplot(df_byClass.s.f, aes(shared_perc, mean_value,
color=perm_perc,
group=perm_perc,
ymin=mean_value-sd_value,
ymax=mean_value+sd_value)) +
geom_point(alpha=0.8) +
geom_linerange(alpha=0.8, size=0.5) +
geom_line() +
scale_color_discrete('% of rank\nabundances\npermuted') +
labs(x='% taxa shared among pre-fractionation communities') +
facet_grid(variables ~ method) +
theme_bw() +
theme(
axis.title.y = element_blank()
)
plot(p.pnt)
In [63]:
outF = file.path(workDir, 'microBetaDiv_acc.pdf')
ggsave(outF, p.pnt, width=10, height=6)
cat('File written:', outF, '\n')
In [64]:
outF = file.path(workDir, 'microBetaDiv_acc.jpeg')
ggsave(outF, p.pnt, width=10, height=6)
cat('File written:', outF, '\n')
In [65]:
BDshift_files = list.files(path=workDir, pattern='BD-shift_stats.txt', full.names=TRUE, recursive=TRUE)
BDshift_files %>% length %>% print
In [66]:
df_shift = list()
for(F in BDshift_files){
tmp = read.delim(F, sep='\t')
FF = strsplit(F, '/') %>% unlist
FFl = length(FF)
tmp$shared_perc = FF[FFl-3]
tmp$perm_perc = FF[FFl-2]
tmp$rep = FF[FFl-1]
tmp$file = FF[FFl]
df_shift[[F]] = tmp
}
df_shift = do.call(rbind, df_shift)
rownames(df_shift) = 1:nrow(df_shift)
df_shift = df_shift %>%
filter(library %in% c(2,4,6)) %>%
group_by(taxon, shared_perc, perm_perc, rep) %>%
summarize(median = median(median)) %>%
ungroup() %>%
rename('median_true_BD_shift' = median) %>%
mutate(true_incorporator = ifelse(median_true_BD_shift > 0.002, TRUE, FALSE))
# status
df_shift %>% nrow %>% print
df_shift %>% head(n=3)
In [67]:
comm_files = list.files(path=workDir, pattern='comm.txt', full.names=TRUE, recursive=TRUE)
comm_files %>% length %>% print
In [68]:
df_comm = list()
for(F in comm_files){
tmp = read.delim(F, sep='\t')
FF = strsplit(F, '/') %>% unlist
FFl = length(FF)
tmp$shared_perc = FF[FFl-3]
tmp$perm_perc = FF[FFl-2]
tmp$rep = FF[FFl-1]
tmp$file = FF[FFl]
df_comm[[F]] = tmp
}
df_comm = do.call(rbind, df_comm)
rownames(df_comm) = 1:nrow(df_comm)
# mean abund for libraries
df_comm = df_comm %>%
rename('taxon' = taxon_name) %>%
group_by(taxon, shared_perc, perm_perc, rep) %>%
summarize(mean_rel_abund_perc = mean(rel_abund_perc),
mean_rank_abund = mean(rank)) %>%
ungroup()
# status
df_comm %>% nrow %>% print
df_comm %>% head(n=3)
In [69]:
MW_files = list.files(path=workDir, pattern='_MW_DESeq2_incorp.txt', full.names=TRUE, recursive=TRUE)
MW_files %>% length %>% print
In [70]:
df_MW = list()
#for(F in MW_files){
for(F in MW_files){
tmp = read.delim(F, sep='\t')
FF = strsplit(F, '/') %>% unlist
FFl = length(FF)
tmp$shared_perc = FF[FFl-3]
tmp$perm_perc = FF[FFl-2]
tmp$rep = FF[FFl-1]
tmp$file = FF[FFl]
df_MW[[F]] = tmp
}
df_MW = do.call(rbind, df_MW)
rownames(df_MW) = 1:nrow(df_MW)
# status
df_MW %>% nrow %>% print
df_MW %>% head(n=3)
In [71]:
join_vars = c('taxon' = 'taxon',
'shared_perc'='shared_perc',
'perm_perc'='perm_perc',
'rep'='rep')
# joining tables
df_MW %>% nrow %>% print
df.j = df_MW %>%
left_join(df_shift, join_vars)
df.j %>% nrow %>% print
df.j = df.j %>%
left_join(df_comm, join_vars)
# status
df.j %>% nrow %>% print
df.j %>% head(n=3)
In [72]:
# calling true_pos + false_neg
df.j = df.j %>%
mutate(incorp_cls = mapply(clsfy, incorp, true_incorporator))
# status
df.j %>% nrow %>% print
df.j$incorp_cls %>% table %>% print
In [73]:
# function for calculating sensitivity
calc_sensitivity = function(incorp_cls){
tp = sum(incorp_cls == 'True positive')
fn = sum(incorp_cls == 'False negative')
x = tp / (tp + fn)
ifelse(is.na(x), 0, x)
}
# grouping by abundance and calculating sensitivity
df.j.s = df.j %>%
mutate(n_group = ntile(log10(mean_rel_abund_perc), 10)) %>%
group_by(n_group, shared_perc, perm_perc, rep) %>%
summarize(min_abund = min(mean_rel_abund_perc, na.rm=TRUE),
mean_abund = mean(mean_rel_abund_perc, na.rm=TRUE),
max_abund = max(mean_rel_abund_perc, na.rm=TRUE),
sensitivity = calc_sensitivity(incorp_cls)) %>%
group_by(n_group, shared_perc, perm_perc) %>%
summarize(mean_abund = mean(mean_abund),
mean_sensitivity = mean(sensitivity),
sd_sensitivity = sd(sensitivity)) %>%
ungroup() %>%
mutate(shared_perc = shared_perc %>% as.character,
shared_perc = shared_perc %>% reorder(shared_perc %>% as.numeric),
perm_perc = perm_perc %>% as.character,
perm_perc = perm_perc %>% reorder(perm_perc %>% as.numeric))
# status
df.j.s %>% head(n=3)
In [74]:
# plotting
options(repr.plot.width=8, repr.plot.height=4)
p_sens_abund = ggplot(df.j.s, aes(mean_abund, mean_sensitivity,
color=perm_perc,
ymin=mean_sensitivity-sd_sensitivity,
ymax=mean_sensitivity+sd_sensitivity)) +
geom_line(alpha=0.7) +
geom_linerange(alpha=0.7) +
geom_point(alpha=0.7) +
scale_x_log10(breaks=c(1e-3, 1e-2, 1e-1, 1e0)) +
scale_color_discrete('% of rank\nabundances\npermuted') +
labs(x='Mean % abundance',
y='Sensitivity') +
facet_wrap(~ shared_perc) +
theme_bw()
p_sens_abund
In [75]:
outF = file.path(workDir, 'microBetaDiv_sens-abund.pdf')
ggsave(outF, p_sens_abund, width=10, height=5)
cat('File written:', outF, '\n')
In [76]:
outF = file.path(workDir, 'microBetaDiv_sens-abund.jpeg')
ggsave(outF, p_sens_abund, width=10, height=5)
cat('File written:', outF, '\n')
In [77]:
betaDiv_files = list.files(path=workDir, pattern='comm_betaDiv.txt', full.names=TRUE, recursive=TRUE)
betaDiv_files %>% length %>% print
In [78]:
df_beta = list()
for(F in betaDiv_files){
tmp = read.delim(F, sep='\t')
FF = strsplit(F, '/') %>% unlist
FFl = length(FF)
tmp$shared_perc = FF[FFl-3]
tmp$perm_perc = FF[FFl-2]
tmp$rep = FF[FFl-1]
df_beta[[F]] = tmp
}
df_beta = do.call(rbind, df_beta)
rownames(df_beta) = 1:nrow(df_beta)
# status
df_beta %>% head(n=3)
In [79]:
# mean bray per parameter set
df_beta_s = df_beta %>%
group_by(shared_perc, perm_perc, rep) %>%
summarize(mean_BC = mean(bray)) %>%
ungroup() %>%
mutate(shared_perc = shared_perc %>% as.Num,
perm_perc = perm_perc %>% as.Num,
rep = rep %>% as.Num)
# status
df_beta_s %>% head(n=3)
In [80]:
# confusion matrix data
vars = c('Balanced Accuracy', 'Sensitivity', 'Specificity')
df_byClass.f = df_byClass %>%
filter(variables %in% vars) %>%
mutate(shared_perc = shared_perc %>% as.Num,
perm_perc = perm_perc %>% as.Num,
rep = rep %>% as.Num)
df_byClass.f.j = inner_join(df_byClass.f, df_beta_s,
c('shared_perc'='shared_perc',
'perm_perc'='perm_perc',
'rep'='rep')) %>%
#filter(variables == 'Balanced Accuracy') %>%
mutate(perm_perc = perm_perc %>% as.character,
perm_perc = perm_perc %>% reorder(perm_perc %>% as.numeric),
shared_perc = shared_perc %>% as.character,
shared_perc = shared_perc %>% reorder(shared_perc %>% as.numeric))
# status
df_byClass.f %>% head(n=3)
In [81]:
# plotting
options(repr.plot.width=9, repr.plot.height=5)
p.betadiv = ggplot(df_byClass.f.j %>% filter(library==6), aes(mean_BC, values)) +
stat_smooth(level=0.99, n=10) +
labs(x='Mean Bray-Curtis distance among pre-fractionation communities') +
facet_grid(variables ~ method) +
theme_bw() +
theme(
axis.title.y = element_blank()
)
p.betadiv
In [82]:
outF = file.path(workDir, 'microBetaDiv_bc-smooth.pdf')
ggsave(outF, p.betadiv, width=10, height=6)
cat('File written:', outF, '\n')
In [83]:
outF = file.path(workDir, 'microBetaDiv_bc-smooth.jpeg')
ggsave(outF, p.betadiv, width=10, height=6)
cat('File written:', outF, '\n')
In [84]:
BDshift_files = list.files(path=workDir, pattern='BD-shift_stats.txt', full.names=TRUE, recursive=TRUE)
BDshift_files %>% length %>% print
In [85]:
df_shift = list()
for(F in BDshift_files){
tmp = read.delim(F, sep='\t')
FF = strsplit(F, '/') %>% unlist
FFl = length(FF)
tmp$shared_perc = FF[FFl-3]
tmp$perm_perc = FF[FFl-2]
tmp$rep = FF[FFl-1]
tmp$file = FF[FFl]
df_shift[[F]] = tmp
}
df_shift = do.call(rbind, df_shift)
rownames(df_shift) = 1:nrow(df_shift)
df_shift = df_shift %>%
filter(library %in% c(2,4,6)) %>%
group_by(taxon, shared_perc, perm_perc, rep) %>%
summarize(median = median(median)) %>%
ungroup() %>%
rename('median_true_BD_shift' = median)
# status
df_shift %>% nrow %>% print
df_shift %>% head(n=3)
In [86]:
incorp_files = list.files(path=workDir, pattern='OTU_abs1e9_PCR_sub_filt_MW_DESeq2_incorp.txt', full.names=TRUE, recursive=TRUE)
incorp_files %>% length %>% print
In [87]:
df_incorp = list()
for(F in incorp_files){
tmp = read.delim(F, sep='\t')
FF = strsplit(F, '/') %>% unlist
FFl = length(FF)
tmp$shared_perc = FF[FFl-3]
tmp$perm_perc = FF[FFl-2]
tmp$rep = FF[FFl-1]
tmp$file = FF[FFl]
df_incorp[[F]] = tmp
}
df_incorp = do.call(rbind, df_incorp)
rownames(df_incorp) = 1:nrow(df_incorp)
df_incorp %>% head(n=3) %>% print
In [88]:
# just incorporators
df_incorp = df_incorp %>%
filter(incorp == TRUE) %>%
dplyr::distinct(taxon, incorp, shared_perc, perm_perc, rep) %>%
rename('HWHRSIP_incorp' = incorp)
df_incorp %>% nrow %>% print
df_incorp %>% head(n=3) %>% print
In [89]:
atomX_files = list.files(path=workDir, pattern='*_qSIP_atom.txt', full.names=TRUE, recursive=TRUE)
atomX_files %>% length %>% print
In [90]:
df_atomX = list()
for(F in atomX_files){
tmp = read.delim(F, sep='\t')
FF = strsplit(F, '/') %>% unlist
FFl = length(FF)
tmp$shared_perc = FF[FFl-3]
tmp$perm_perc = FF[FFl-2]
tmp$rep = FF[FFl-1]
tmp$file = FF[FFl]
df_atomX[[F]] = tmp
}
df_atomX = do.call(rbind, df_atomX)
rownames(df_atomX) = 1:nrow(df_atomX)
df_atomX %>% head(n=3) %>% print
In [91]:
# table join
df_atomX %>% nrow %>% print
df.j = left_join(df_atomX, df_shift, c('taxon' = 'taxon',
'shared_perc'='shared_perc',
'perm_perc'='perm_perc',
'rep'='rep')) %>%
filter(!is.na(BD_diff)) %>%
mutate(true_incorporator = ifelse(median_true_BD_shift > 0.002, TRUE, FALSE),
true_atom_fraction_excess = median_true_BD_shift / 0.036,
atom_fraction_excess = ifelse(is.na(atom_CI_low), 0, atom_fraction_excess))
df.j %>% nrow %>% print
df.j %>% head(n=3) %>% print
In [92]:
df.j$true_incorporator %>% summary
In [93]:
df.j = left_join(df.j, df_incorp, c('taxon' = 'taxon',
'shared_perc'='shared_perc',
'perm_perc'='perm_perc',
'rep'='rep')) %>%
mutate(HWHRSIP_incorp = ifelse(is.na(HWHRSIP_incorp), FALSE, TRUE))
df.j %>% nrow %>% print
df.j %>% head(n=3) %>% print
In [94]:
# difference between true and estimated
## q-SIP incorporators
df.j.dis.qSIP = df.j %>%
filter(atom_CI_low > 0) %>% # just incorporators identified by q-SIP
mutate(delta_excess = atom_fraction_excess * 100 - true_atom_fraction_excess * 100) %>%
group_by(shared_perc, perm_perc) %>%
summarize(mean_delta_excess = mean(delta_excess),
sd_delta_excess = sd(delta_excess)) %>%
ungroup() %>%
mutate(shared_perc = shared_perc %>% reorder(shared_perc %>% as.numeric),
perm_perc = perm_perc %>% reorder(perm_perc %>% as.numeric))
# plotting
options(repr.plot.width=6, repr.plot.height=3)
p_qSIP = ggplot(df.j.dis.qSIP, aes(shared_perc, mean_delta_excess,
color=perm_perc, group=perm_perc,
ymin=mean_delta_excess-sd_delta_excess,
ymax=mean_delta_excess+sd_delta_excess)) +
geom_linerange(alpha=0.4, size=1) +
geom_point() +
geom_line() +
#facet_grid(true_incorporator ~ .) +
scale_color_discrete('% incorp-\norators') +
labs(x='% taxa shared among pre-fractionation communities',
y='13C atom % excess\n(truth - estimate)') +
theme_bw()
p_qSIP
In [95]:
# difference between true and estimated
## q-SIP incorporators
df.j.dis.qSIP = df.j %>%
filter(atom_CI_low > 0) %>% # just incorporators identified by q-SIP
filter(HWHRSIP_incorp == TRUE) %>% # just MW-HR-SIP incorporators
mutate(delta_excess = atom_fraction_excess * 100 - true_atom_fraction_excess * 100) %>%
group_by(shared_perc, perm_perc) %>%
summarize(mean_delta_excess = mean(delta_excess),
sd_delta_excess = sd(delta_excess)) %>%
ungroup() %>%
mutate(shared_perc = shared_perc %>% reorder(shared_perc %>% as.numeric),
perm_perc = perm_perc %>% reorder(perm_perc %>% as.numeric))
# plotting
options(repr.plot.width=6, repr.plot.height=2.5)
p_qSIP = ggplot(df.j.dis.qSIP, aes(shared_perc, mean_delta_excess,
color=perm_perc, group=perm_perc,
ymin=mean_delta_excess-sd_delta_excess,
ymax=mean_delta_excess+sd_delta_excess)) +
geom_linerange(alpha=0.4, size=1) +
geom_point() +
geom_line() +
#facet_grid(true_incorporator ~ .) +
scale_color_discrete('% of rank\nabundances\npermuted') +
labs(x='% taxa shared among pre-fractionation communities',
y='13C atom % excess\n(truth - estimate)') +
theme_bw()
p_qSIP
In [96]:
# difference between true and estimated
## q-SIP incorporators
tmp1 = df.j %>%
filter(atom_CI_low > 0) %>% # just incorporators identified by q-SIP
mutate(delta_excess = atom_fraction_excess * 100 - true_atom_fraction_excess * 100) %>%
group_by(shared_perc, perm_perc) %>%
summarize(mean_delta_excess = mean(delta_excess),
sd_delta_excess = sd(delta_excess)) %>%
ungroup() %>%
mutate(shared_perc = shared_perc %>% reorder(shared_perc %>% as.numeric),
perm_perc = perm_perc %>% reorder(perm_perc %>% as.numeric),
incorp_called = 'No filter')
## MW-HR-SIP incorporators
tmp2 = df.j %>%
filter(atom_CI_low > 0) %>% # just incorporators identified by q-SIP
filter(HWHRSIP_incorp == TRUE) %>% # just MW-HR-SIP incorporators
mutate(delta_excess = atom_fraction_excess * 100 - true_atom_fraction_excess * 100) %>%
group_by(shared_perc, perm_perc) %>%
summarize(mean_delta_excess = mean(delta_excess),
sd_delta_excess = sd(delta_excess)) %>%
ungroup() %>%
mutate(shared_perc = shared_perc %>% reorder(shared_perc %>% as.numeric),
perm_perc = perm_perc %>% reorder(perm_perc %>% as.numeric),
incorp_called = 'MW-HR-SIP filter')
# combining tables
df.j.dis.qSIP = rbind(tmp1, tmp2) %>%
mutate(incorp_called = factor(incorp_called, levels=c('No filter', 'MW-HR-SIP filter')))
# plotting
options(repr.plot.width=6, repr.plot.height=4)
p_qSIP = ggplot(df.j.dis.qSIP, aes(shared_perc, mean_delta_excess,
color=perm_perc, group=perm_perc,
ymin=mean_delta_excess-sd_delta_excess,
ymax=mean_delta_excess+sd_delta_excess)) +
geom_line() +
geom_linerange(alpha=0.4, size=1) +
geom_point() +
facet_grid(incorp_called ~ .) +
scale_color_discrete('% of rank\nabundances\npermuted') +
labs(x='% taxa shared among pre-fractionation communities',
y='13C atom % excess\n(truth - estimate)') +
theme_bw()
p_qSIP
In [97]:
dBD_files = list.files(path=workDir, pattern='*_dBD.txt', full.names=TRUE, recursive=TRUE)
dBD_files %>% length %>% print
In [98]:
df_dBD = list()
for(F in dBD_files){
tmp = read.delim(F, sep='\t')
FF = strsplit(F, '/') %>% unlist
FFl = length(FF)
tmp$shared_perc = FF[FFl-3]
tmp$perm_perc = FF[FFl-2]
tmp$rep = FF[FFl-1]
tmp$file = FF[FFl]
df_dBD[[F]] = tmp
}
df_dBD = do.call(rbind, df_dBD)
rownames(df_dBD) = 1:nrow(df_dBD)
df_dBD %>% head(n=3) %>% print
In [99]:
df.j = inner_join(df_dBD, df_shift, c('taxon' = 'taxon',
'shared_perc'='shared_perc',
'perm_perc'='perm_perc',
'rep'='rep')) %>%
filter(!is.na(delta_BD)) %>%
mutate(true_incorporator = ifelse(median_true_BD_shift > 0.002, TRUE, FALSE),
true_atom_fraction_excess = median_true_BD_shift / 0.036,
atom_fraction_excess = delta_BD / 0.036)
df.j %>% head(n=3)
In [100]:
df.j = left_join(df.j, df_incorp, c('taxon' = 'taxon',
'shared_perc'='shared_perc',
'perm_perc'='perm_perc',
'rep'='rep')) %>%
mutate(HWHRSIP_incorp = ifelse(is.na(HWHRSIP_incorp), FALSE, TRUE))
df.j %>% nrow %>% print
df.j %>% head(n=3)
In [101]:
# difference between true and estimated
tmp1 = df.j %>%
mutate(delta_excess = atom_fraction_excess * 100 - true_atom_fraction_excess * 100) %>%
group_by(shared_perc, perm_perc) %>%
summarize(mean_delta_excess = mean(delta_excess),
sd_delta_excess = sd(delta_excess)) %>%
ungroup() %>%
mutate(shared_perc = shared_perc %>% reorder(shared_perc %>% as.numeric),
perm_perc = perm_perc %>% reorder(perm_perc %>% as.numeric),
incorp_called = 'No filter')
tmp2 = df.j %>%
filter(HWHRSIP_incorp == TRUE) %>% # just MW-HR-SIP incorporators
mutate(delta_excess = atom_fraction_excess * 100 - true_atom_fraction_excess * 100) %>%
group_by(shared_perc, perm_perc) %>%
summarize(mean_delta_excess = mean(delta_excess),
sd_delta_excess = sd(delta_excess)) %>%
ungroup() %>%
mutate(shared_perc = shared_perc %>% reorder(shared_perc %>% as.numeric),
perm_perc = perm_perc %>% reorder(perm_perc %>% as.numeric),
incorp_called = 'MW-HR-SIP filter')
# combining tables
df.j.dis.dBD = rbind(tmp1, tmp2) %>%
mutate(incorp_called = factor(incorp_called, levels=c('No filter', 'MW-HR-SIP filter')))
# plotting
options(repr.plot.width=8, repr.plot.height=4)
p_dBD = ggplot(df.j.dis.dBD, aes(shared_perc, mean_delta_excess,
color=perm_perc, group=perm_perc,
ymin=mean_delta_excess-sd_delta_excess,
ymax=mean_delta_excess+sd_delta_excess)) +
geom_line() +
geom_linerange(alpha=0.4, size=1) +
geom_point() +
facet_grid(incorp_called ~ .) +
scale_color_discrete('% of rank\nabundances\npermuted') +
labs(x='% taxa shared among pre-fractionation communities',
y='13C atom % excess\n(truth - estimate)') +
theme_bw()
p_dBD
In [102]:
df.jj = rbind(df.j.dis.qSIP %>% mutate(method='qSIP'),
df.j.dis.dBD %>% mutate(method='Delta BD')) %>%
mutate(method = gsub('qSIP', 'q-SIP', method))
p.comb = ggplot(df.jj, aes(shared_perc, mean_delta_excess,
color=perm_perc, group=perm_perc,
ymin=mean_delta_excess-sd_delta_excess,
ymax=mean_delta_excess+sd_delta_excess)) +
geom_line() +
geom_linerange(alpha=0.5, size=1) +
geom_point() +
scale_color_discrete('% of rank\nabundances\npermuted') +
labs(x='% taxa shared among pre-fractionation communities',
y='atom % excess 13C\n(truth - estimate)') +
facet_grid(incorp_called ~ method) +
theme_bw()
p.comb
In [103]:
outF = file.path(workDir, 'microBetaDiv_BDshift.pdf')
ggsave(outF, p.comb, width=8, height=4)
cat('File written:', outF, '\n')
In [104]:
outF = file.path(workDir, 'microBetaDiv_BDshift.jpeg')
ggsave(outF, p.comb, width=8, height=4)
cat('File written:', outF, '\n')
In [ ]: