In [88]:
import os
import glob
import re
import nestly
In [89]:
%load_ext rpy2.ipython
%load_ext pushnote
In [90]:
%%R
library(ggplot2)
library(dplyr)
library(tidyr)
library(gridExtra)
library(phyloseq)
In [91]:
## 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
print 'Min BD: {}'.format(min_BD)
print 'Max BD: {}'.format(max_BD)
In [92]:
workDir = '/home/nick/notebook/SIPSim/dev/fullCyc/n1147_frag_norm_9_2.5_n5/'
buildDir = os.path.join(workDir, 'rep4_DBL-comm_bw_ALL4')
R_dir = '/home/nick/notebook/SIPSim/lib/R/'
fragFile = '/home/nick/notebook/SIPSim/dev/bac_genome1147/validation/ampFrags_kde.pkl'
commFile = '/home/nick/notebook/SIPSim/dev/fullCyc/fullCyc_12C-Con_trm_comm.txt'
# emperical data for validation
emp_shan_file = '/home/nick/notebook/SIPSim/dev/fullCyc_trim/SIP-core_unk_shan.txt'
emp_BDspan_file = '/home/nick/notebook/SIPSim/dev/fullCyc_trim/SIP-core_unk_trm_BD-span.txt'
emp_corr_file = '/home/nick/notebook/SIPSim/dev/fullCyc_trim/SIP-core_unk_trm_corr.txt'
nreps = 4
In [19]:
# building tree structure
nest = nestly.Nest()
# varying params
nest.add('DBL_scaling', [0, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8])
nest.add('bandwidth', [0.002, 0.02, 0.2, 0.4, 0.6, 0.8])
nest.add('rep', [x + 1 for x in xrange(nreps)])
## set params
nest.add('abs', ['1e9'], create_dir=False)
nest.add('percIncorp', [0], create_dir=False)
nest.add('percTaxa', [0], create_dir=False)
nest.add('np', [8], create_dir=False)
nest.add('subsample_dist', ['lognormal'], create_dir=False)
nest.add('subsample_mean', [9.432], create_dir=False)
nest.add('subsample_scale', [0.5], create_dir=False)
nest.add('subsample_min', [10000], create_dir=False)
nest.add('subsample_max', [30000], create_dir=False)
nest.add('min_BD', [min_BD], create_dir=False)
nest.add('max_BD', [max_BD], 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('commFile', [commFile], create_dir=False)
# building directory tree
nest.build(buildDir)
# bash file to run
bashFile = os.path.join(buildDir, 'SIPSimRun.sh')
In [20]:
%%writefile $bashFile
#!/bin/bash
export PATH={R_dir}:$PATH
echo '#-- SIPSim pipeline --#'
echo '# shuffling taxa in comm file'
comm_shuffle_taxa.r {commFile} > comm.txt
echo '# simulating gradient fractions'
SIPSim gradient_fractions \
--BD_min {min_BD} \
--BD_max {max_BD} \
comm.txt \
> fracs.txt
echo '# adding diffusion'
SIPSim diffusion \
{fragFile} \
-n 100000 \
--bw {bandwidth} \
--np {np} \
> ampFrags_KDE_dif.pkl
echo '# adding DBL contamination; abundance-weighted smearing'
SIPSim DBL \
ampFrags_KDE_dif.pkl \
-n 100000 \
--comm comm.txt \
--commx {DBL_scaling} \
--np {np} \
> ampFrags_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 \
ampFrags_KDE_dif_DBL.pkl \
{percTaxa}_{percIncorp}.config \
-n 100000 \
--comm comm.txt \
--np {np} \
> ampFrags_KDE_dif_DBL_inc.pkl
echo '# simulating an OTU table'
SIPSim OTU_table \
ampFrags_KDE_dif_DBL_inc.pkl \
comm.txt \
fracs.txt \
--abs {abs} \
--np {np} \
> OTU_abs{abs}.txt
#-- w/ PCR simulation --#
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}_PCR.txt \
> OTU_abs{abs}_PCR_sub.txt
echo '# making a wide-formatted table'
SIPSim OTU_wideLong -w \
OTU_abs{abs}_PCR_sub.txt \
> OTU_abs{abs}_PCR_sub_w.txt
echo '# making metadata (phyloseq: sample_data)'
SIPSim OTU_sampleData \
OTU_abs{abs}_PCR_sub.txt \
> OTU_abs{abs}_PCR_sub_meta.txt
#-- w/out PCR simulation --#
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
#-- making summary tables --#
# PCR
shannon_calc.r OTU_abs{abs}_PCR_sub.txt > OTU_abs{abs}_PCR_sub_shan.txt
BD_span_calc.r OTU_abs{abs}_PCR_sub.txt comm.txt > OTU_abs{abs}_PCR_sub_BD-span.txt
correlogram_make.r OTU_abs{abs}_PCR_sub.txt > OTU_abs{abs}_PCR_sub_corr.txt
# no PCR
shannon_calc.r OTU_abs{abs}_sub.txt > OTU_abs{abs}_sub_shan.txt
BD_span_calc.r OTU_abs{abs}_sub.txt comm.txt > OTU_abs{abs}_sub_BD-span.txt
correlogram_make.r OTU_abs{abs}_sub.txt > OTU_abs{abs}_sub_corr.txt
#-- removing large intermediate files --#
rm -f ampFrags_KDE_dif.pkl
rm -f ampFrags_KDE_dif_DBL.pkl
rm -f ampFrags_KDE_dif_DBL_inc.pkl
In [ ]:
!chmod 777 $bashFile
!cd $workDir; \
nestrun --template-file $bashFile -d rep4_DBL-comm_bw_ALL4 --log-file log.txt -j 2
In [ ]:
%pushnote rep4_DBL-comm_bw_ALL4 complete
In [23]:
%%R
# function for loading dataset files
load.data.files = function(sim.files, emp.file){
# loading
## simulations
df = list()
for(x in sim.files){
# simulation
tmp = read.delim(x, sep='\t')
xx = strsplit(x, '/')[[1]]
tmp$DBL_scale = xx[10] %>% as.numeric
tmp$bw = xx[11] %>% as.numeric
tmp$SIM_rep = xx[12] %>% as.numeric
tmp$dataset = 'Simulation'
df[[x]] = tmp
# emperical (matched for each simulation)
if(xx[12] %>% as.numeric == 1){
tmp = read.delim(emp.file, sep='\t')
tmp$DBL_scale = xx[10] %>% as.numeric
tmp$bw = xx[11] %>% as.numeric
tmp$SIM_rep = 1
tmp$dataset = 'Emperical'
xy = paste0(x, '_EMP')
df[[xy]] = tmp
}
}
df = do.call(rbind, df) %>% as.data.frame
rownames(df) = 1:nrow(df)
# return
return(df)
}
In [93]:
sim_shan_files = !find $buildDir -name "OTU_abs1e9_PCR_sub_shan.txt"
print len(sim_shan_files)
In [94]:
# checking for empty files
for x in sim_shan_files:
ret = !ls -thlc $x
if ret[0].split(' ')[4] == '0':
print ret
In [95]:
%%R -i sim_shan_files -i emp_shan_file
df.shan = load.data.files(sim_shan_files, emp_shan_file)
df.shan %>% tail(n=3)
In [96]:
%%R -w 800 -h 600
# summarizing
df.shan.s = df.shan %>%
group_by(dataset, bw, DBL_scale, BD_bin = ntile(Buoyant_density, 24)) %>%
summarize(mean_shannon = mean(shannon),
sd_shannon = sd(shannon),
mean_BD = mean(Buoyant_density))
ggplot(df.shan.s, aes(mean_BD, mean_shannon, color=dataset,
ymin=mean_shannon-sd_shannon, ymax=mean_shannon+sd_shannon)) +
geom_pointrange() +
facet_grid(DBL_scale ~ bw) +
labs(x='Buoyant density (binned; 24 bins)', y='Shannon index') +
theme_bw() +
theme(
text = element_text(size=16),
axis.text.x = element_text(angle=45, hjust=1)
)
In [213]:
%%R -w 650 -h 600
# pairwise correlations for each dataset
df.shan.bin = df.shan %>%
group_by(BD_bin = ntile(Buoyant_density, 24))
calc.pearson = function(x){
cor(x[,'shannon.x'], x['shannon.y'], method='pearson')[1,1]
}
df.shan.corr = inner_join(df.shan.bin, df.shan.bin, c('BD_bin' = 'BD_bin',
'bw' = 'bw',
'DBL_scale' = 'DBL_scale')) %>%
group_by(bw, DBL_scale, dataset.x, dataset.y) %>%
nest() %>%
mutate(model = purrr::map(data, calc.pearson)) %>%
unnest(pearson = model %>% purrr::map(function(x) x)) %>%
ungroup() %>%
select(-data, -model) %>%
mutate(pearson_txt = round(pearson, 2))
# plotting
ggplot(df.shan.corr, aes(dataset.x, dataset.y, fill=pearson)) +
geom_tile() +
geom_text(aes(label=pearson_txt), color='white', size=5) +
scale_fill_gradient(low='black', high='red') +
labs(title='Shannon index') +
facet_grid(DBL_scale ~ bw) +
theme(
text = element_text(size=16),
axis.text.x = element_text(angle=45, hjust=1)
)
In [214]:
%%R -w 500 -h 250
# getting emperical-emperical corr
emp.val = df.shan.corr %>%
filter((dataset.x == 'Emperical' &
dataset.y == 'Emperical')) %>%
group_by() %>%
summarize(max_value = max(pearson)) %>%
ungroup() %>%
select(max_value) %>% as.matrix %>% as.vector
emp.val = emp.val[1]
# filtering
df.shan.corr.f = df.shan.corr %>%
filter((dataset.x == 'Simulation' &
dataset.y == 'Emperical')) %>%
mutate(DBL_scale = DBL_scale %>% as.character,
bw = bw %>% as.character,
pearson_txt = ifelse(round(pearson,2) >= round(emp.val, 2),
gsub('$', '*', pearson_txt),
pearson_txt)) %>%
complete(DBL_scale, bw)
df.shan.corr.f %>% head(n=3)
# plotting
p.shan = ggplot(df.shan.corr.f, aes(DBL_scale,bw, fill=pearson)) +
geom_tile() +
geom_text(aes(label=pearson_txt), color='white', size=5.5) +
scale_x_discrete(expand=c(0,0)) +
scale_y_discrete(expand=c(0,0)) +
scale_color_manual(values=c('white', 'black')) +
scale_fill_gradient('Pearson', low='black', high='red') +
labs(x='gamma', y='KDE Bandwidth', title='Shannon index') +
theme(
text = element_text(size=16),
plot.title = element_text(size=16)
)
p.shan
In [215]:
%%R -w 650 -h 600
# pairwise correlations for each dataset
# df.shan.bin = df.shan %>%
# group_by(BD_bin = ntile(Buoyant_density, 24))
# calc.pearson.p = function(x){
# cor.test(x$shannon.x, x$shannon.y, method='pearson', alternative='greater')$p.value
# }
# df.shan.corr = inner_join(df.shan.bin, df.shan.bin, c('BD_bin' = 'BD_bin',
# 'bw' = 'bw',
# 'DBL_scale' = 'DBL_scale')) %>%
# group_by(bw, DBL_scale, dataset.x, dataset.y) %>%
# nest() %>%
# mutate(model = purrr::map(data, calc.pearson.p)) %>%
# unnest(pearson = model %>% purrr::map(function(x) x)) %>%
# ungroup() %>%
# select(-data, -model) %>%
# mutate(pearson = p.adjust(pearson, method='BH'),
# pearson_txt = round(pearson, 3))
# # plotting
# ggplot(df.shan.corr, aes(dataset.x, dataset.y, fill=pearson)) +
# geom_tile() +
# geom_text(aes(label=pearson_txt), color='white', size=5) +
# scale_fill_gradient(low='black', high='red') +
# labs(title='Shannon index') +
# facet_grid(DBL_scale ~ bw) +
# theme(
# text = element_text(size=16),
# axis.text.x = element_text(angle=45, hjust=1)
# )
In [220]:
%%R -w 650 -h 600
# pairwise correlations for each dataset
df.shan.bin = df.shan %>%
group_by(BD_bin = ntile(Buoyant_density, 24))
calc.pearson = function(x){
cor(x[,'shannon.x'], x['shannon.y'], method='pearson')[1,1]
}
# dataset
df.shan.j = inner_join(df.shan.bin, df.shan.bin, c('BD_bin' = 'BD_bin',
'bw' = 'bw',
'DBL_scale' = 'DBL_scale')) %>%
group_by(bw, DBL_scale, dataset.x, dataset.y)
# correlations
df.shan.j.corr = df.shan.j %>%
nest() %>%
mutate(model = purrr::map(data, calc.pearson)) %>%
unnest(pearson = model %>% purrr::map(function(x) x)) %>%
ungroup() %>%
select(-data, -model) %>%
filter((dataset.x == 'Simulation' &
dataset.y == 'Emperical')) %>%
mutate(boot_rep_id = NA)
df.shan.j.corr %>% head
In [ ]:
%%R
# bootstrapping
subsamp = function(x,y){
base::sample(c(x,y), size=length(x))
}
cor.boot = function(rep_id, df){
# shuffling among the 2 groups
df = df %>%
mutate(shannon.x = subsamp(shannon.x, shannon.y),
shannon.y = subsamp(shannon.x, shannon.y))
# correlations
df %>%
nest() %>%
mutate(model = purrr::map(data, calc.pearson.p)) %>%
unnest(pearson = model %>% purrr::map(function(x) x)) %>%
ungroup() %>%
select(-data, -model) %>%
filter((dataset.x == 'Simulation' &
dataset.y == 'Emperical')) %>%
mutate(boot_rep_id = rep_id)
}
# making bootstrapped tables
params = data.frame(rep_id = 1:100)
df.shan.j.b.corr = plyr::mdply(params, cor.boot, df=df.shan.j)
df.shan.j.b.corr %>% head(n=3)
In [ ]:
%%R
inner_join(df.shan.j.corr, df.shan.j.b.corr,
c('bw'='bw', 'DBL_scale'='DBL_scale')) %>%
head %>% as.data.frame
In [ ]:
%%R
#df.shan.j.corr.j =
inner_join(df.shan.j.corr, df.shan.j.b.corr,
c('bw'='bw', 'DBL_scale'='DBL_scale')) %>%
mutate(gt = pearson.x > pearson.y) %>% # actual is > shuffled
group_by(bw, DBL_scale) %>%
summarize(p = (length(gt) - sum(gt)) / length(gt))
In [30]:
sim_BDspan_files = !find $buildDir -name "OTU_abs1e9_PCR_sub_BD-span.txt"
print len(sim_BDspan_files)
print emp_BDspan_file
In [31]:
%%R -i sim_BDspan_files -i emp_BDspan_file
df.BDspan = load.data.files(sim_BDspan_files, emp_BDspan_file)
df.BDspan %>% head
In [32]:
%%R -w 700 -h 600
# plotting
ggplot(df.BDspan, aes(mean_preFrac_abund, BD_range_perc, fill=dataset)) +
geom_hex(alpha=0.5) +
scale_x_log10() +
facet_grid(DBL_scale ~ bw) +
labs(x='Pre-fractionation abundance', y='BD span') +
theme_bw() +
theme(
text = element_text(size=16),
axis.text.x = element_text(angle=45, hjust=1)
)
In [33]:
%%R -i sim_BDspan_files -i emp_BDspan_file
# binning by pre-fractionation abundances
n.tile = 20
df.BDspan = df.BDspan %>%
group_by(dataset, library, DBL_scale, bw, preFrac_abund_bin = ntile(mean_preFrac_abund, n.tile)) %>%
summarize(mean_preFrac_abund = mean(mean_preFrac_abund),
var_BD_range = var(BD_range),
sd_BD_range = sd(BD_range))
df.BDspan %>% tail(n=3)
In [34]:
%%R -w 675 -h 600
calc.spearman = function(x){
cor(x[,'var_BD_range.x'], x['var_BD_range.y'], method='spearman')[1,1]
}
df.BDspan.corr = inner_join(df.BDspan, df.BDspan, c('preFrac_abund_bin' = 'preFrac_abund_bin',
'DBL_scale' = 'DBL_scale',
'bw' = 'bw')) %>%
group_by(DBL_scale, bw, dataset.x, dataset.y) %>%
nest() %>%
mutate(model = purrr::map(data, calc.spearman)) %>%
unnest(spearman = model %>% purrr::map(function(x) x)) %>%
ungroup() %>%
select(-data, -model) %>%
mutate(spearman_txt = round(spearman, 2))
# plotting
ggplot(df.BDspan.corr, aes(dataset.x, dataset.y, fill=spearman)) +
geom_tile() +
geom_text(aes(label=spearman_txt), color='white', size=5) +
scale_fill_gradient(low='black', high='red') +
labs(title='BD span') +
facet_grid(DBL_scale ~ bw) +
theme(
text = element_text(size=16),
axis.text.x = element_text(angle=45, hjust=1)
)
In [76]:
%%R -w 500 -h 250
# getting emperical-emperical corr
emp.val = df.BDspan.corr %>%
filter((dataset.x == 'Emperical' &
dataset.y == 'Emperical')) %>%
group_by() %>%
summarize(max_value = max(spearman, na.rm=TRUE)) %>%
ungroup() %>%
select(max_value) %>% as.matrix %>% as.vector
emp.val = emp.val[1]
# filtering
df.BDspan.corr.f = df.BDspan.corr %>%
filter((dataset.x == 'Simulation' &
dataset.y == 'Emperical')) %>%
mutate(DBL_scale = DBL_scale %>% as.character,
bw = bw %>% as.character,
gt_emp = ifelse(round(spearman,2) >= round(emp.val,2), 'bold.italic', 'plain')) %>%
complete(DBL_scale, bw)
# plotting
p.span = ggplot(df.BDspan.corr.f, aes(DBL_scale, bw, fill=spearman)) +
geom_tile() +
scale_x_discrete(expand=c(0,0)) +
scale_y_discrete(expand=c(0,0)) +
geom_text(aes(label=spearman_txt, fontface=gt_emp), color='white', size=5.5) +
scale_color_manual(values=c('white', 'black')) +
scale_fill_gradient('Spearman', low='black', high='red') +
labs(x='gamma', y='KDE Bandwidth', title='OTU BD range') +
theme(
text = element_text(size=16),
plot.title = element_text(size=16)
)
p.span
In [36]:
sim_corr_files = !find $buildDir -name "OTU_abs1e9_PCR_sub_corr.txt"
print len(sim_corr_files)
print emp_corr_file
In [42]:
# checking for empty files
for x in sim_corr_files:
ret = !ls -thlc $x
if ret[0].split(' ')[4] == '0':
print ret
In [43]:
%%R -i sim_corr_files -i emp_corr_file
df.corr = load.data.files(sim_corr_files, emp_corr_file)
# binning
df.corr = df.corr %>%
filter(!is.na(Mantel.corr)) %>%
group_by(DBL_scale, bw, dataset, library, class.index.bin = ntile(class.index, 12))
df.corr %>% tail(n=3) %>% as.data.frame
In [44]:
%%R -w 800 -h 600
# plotting
df.corr.s = df.corr %>%
group_by(DBL_scale, bw, dataset, class.index.bin) %>%
summarize(mean_Mantel.corr = mean(Mantel.corr),
sd_Mantel.corr = sd(Mantel.corr),
mean_class.index = mean(class.index))
ggplot(df.corr.s, aes(mean_class.index, mean_Mantel.corr, color=dataset,
ymin=mean_Mantel.corr-sd_Mantel.corr,
ymax=mean_Mantel.corr+sd_Mantel.corr)) +
geom_pointrange(size=0.2) +
labs(x='Class index (binned; 12 bins)', y='Mantel correlation coef.') +
facet_grid(DBL_scale ~ bw) +
theme_bw() +
theme(
text = element_text(size=16),
axis.text.x = element_text(angle=45, hjust=1)
)
In [45]:
%%R -w 700 -h 600
# pairwise correlations for each dataset
df.shan.bin = df.shan %>%
group_by(BD_bin = ntile(Buoyant_density, 24))
calc.pearson = function(x){
cor(x[,'Mantel.corr.x'], x['Mantel.corr.y'], method='pearson')[1,1]
}
df.corr.lm = inner_join(df.corr, df.corr, c('class.index.bin' = 'class.index.bin',
'bw' = 'bw',
'DBL_scale' = 'DBL_scale')) %>%
group_by(bw, DBL_scale, dataset.x, dataset.y) %>%
nest() %>%
mutate(model = purrr::map(data, calc.pearson)) %>%
unnest(pearson = model %>% purrr::map(function(x) x)) %>%
ungroup() %>%
select(-data, -model) %>%
mutate(pearson_txt = round(pearson, 2))
# plotting
ggplot(df.corr.lm, aes(dataset.x, dataset.y, fill=pearson)) +
geom_tile() +
geom_text(aes(label=pearson_txt), color='white', size=5) +
scale_fill_gradient(low='black', high='red') +
labs(title='Beta diversity correlogram') +
facet_grid(DBL_scale ~ bw) +
theme(
text = element_text(size=16),
axis.text.x = element_text(angle=45, hjust=1)
)
In [71]:
%%R -w 500 -h 250
# getting emperical-emperical corr
emp.val = df.corr.lm %>%
filter((dataset.x == 'Emperical' &
dataset.y == 'Emperical')) %>%
group_by() %>%
summarize(max_value = max(pearson)) %>%
ungroup() %>%
select(max_value) %>% as.matrix %>% as.vector
emp.val = emp.val[1]
print(emp.val)
# filtering
df.corr.lm.f = df.corr.lm %>%
filter((dataset.x == 'Simulation' &
dataset.y == 'Emperical')) %>%
mutate(DBL_scale = DBL_scale %>% as.character,
bw = bw %>% as.character,
gt_emp = ifelse(round(pearson,2) >= round(emp.val,2), 'bold.italic', 'plain')) %>%
complete(DBL_scale, bw)
df.corr.lm.f %>% head(n=3)
# plotting
p.corr = ggplot(df.corr.lm.f, aes(DBL_scale,bw, fill=pearson)) +
geom_tile() +
geom_text(aes(label=pearson_txt,fontface=gt_emp), color='white', size=5.5) +
scale_x_discrete(expand=c(0,0)) +
scale_y_discrete(expand=c(0,0)) +
scale_color_manual(values=c('white', 'black')) +
scale_fill_gradient('Pearson', low='black', high='red') +
labs(x='gamma', y='KDE Bandwidth', title='Correlograms of Jaccard dissimilarities') +
theme(
text = element_text(size=16),
plot.title = element_text(size=16)
)
p.corr
In [85]:
%%R -w 500 -h 650
p.comb = cowplot::ggdraw() +
geom_rect(aes(xmin=0, ymin=0, xmax=1, ymax=1), fill='white') +
cowplot::draw_plot(p.shan, 0, 0.66, 0.96, 0.33) +
cowplot::draw_plot(p.corr, 0, 0.33, 0.965, 0.33) +
cowplot::draw_plot(p.span, 0, 0.00, 1, 0.33) +
cowplot::draw_plot_label(c('A)', 'B)', 'C)'), c(0,0,0), c(0.99, 0.66, 0.33))
p.comb
In [87]:
%%R -i workDir
# saving figure
F = file.path(workDir, 'SIM-v-emp_param_optimize.pdf')
ggsave(F, p.comb, width=7.7, height=10)
cat('File written:', F, '\n')
In [ ]: