“Publication Figure 3”
百度云盘链接: https://pan.baidu.com/s/15g7caZp354zIWktpnWzWhQ
提取码: 4sh7
Libraries
Standard Import
library(tidyverse)
library(cowplot)
library(scales)
library(ggpubr)
Special
library(ggridges)
library(grid)
# library(Hmisc)
# library(ggExtra)
Paths
bin_dir = '../bin/'
data_dir = '../data/'
results_dir = '../results/'
Custom Scripts / Theme
source(paste0(bin_dir, 'theme_settings.R'))
source (paste0(bin_dir, 'utilities.R'))
source(paste0(bin_dir, 'species_competition_functions.R'))
source(paste0(bin_dir, 'distmat_utils.R'))
source(paste0(bin_dir, 'analysis_metadata.R'))
date <- format(Sys.time(), "%d%m%y")
Import Tables
Metadata
samples.metadata <- read_tsv(paste0(data_dir, 'samples.metadata.tsv'))
microbe.taxonomy <- read_tsv(paste0(data_dir, 'microbe.taxonomy.tsv'))
microbe.metadata <- microbe.taxonomy %>%
right_join(., read_tsv(paste0(data_dir, 'microbe.metadata.tsv')))
Taxonomy
MetaPhlAn2 tables
combine with taxonomy
mp_species.long <- microbe.taxonomy %>%
right_join(., read_tsv(paste0(data_dir, 'samples.mp_profiles.species.tsv')),
by = 'species') %>%
left_join(., samples.metadata, by = 'Name')
rCDI subset
annotate with metadata
mp_species.rcdi.long <-
mp_species.long %>%
filter(Study %in% c(rcdi.studies))
SameStr Data
SameStr Case-Summary table
samples.sstr_data <- read_tsv(paste0(data_dir, 'samples.sstr_data_wmetadata.tsv')) %>%
left_join(., microbe.taxonomy, by = 'species')
Strain Finder Data
Strain Finder Tables
Sys.setenv("VROOM_CONNECTION_SIZE" = 10485760)
sf_strain_table <- read_tsv(paste0(data_dir, 'ALM.Strain_Table.tsv')) %>%
rename(Sample = '...1') # X1
sf_species_table <- read_tsv(paste0(data_dir, 'ALM.MG_OTU_Bacteria.tsv')) %>%
rename(bact = `BACT ID`) %>%
rename_all(.funs = funs(str_to_lower(.)))
Mock Data
Multi-Strain Species Population
mock.sstr_data <- read_tsv(paste0(data_dir, 'mock.sstr_data_wmetadata.tsv')) %>%
left_join(., microbe.taxonomy, by = 'species')
Fast ANI
FastAni Comparison Data
refs <- read_tsv(paste0(data_dir, 'references.mp_fastani.tsv'))
Set Figure
fig = paste0('Fig_3.')
Specificity of Shared Taxa
Comparing distributions of shared (genera, species, and) strains in related and unrelated sample pairs in order to assess the specificity with which methods call shared strains.
Strain Finder
Format
Taxa per sample matrix
# mg-OTU strain-level matrix
sf_strain_table.matrix <-
sf_strain_table %>%
mutate_at(.vars = vars(everything(), -Sample),
.funs = funs(as.numeric(. > 0))) %>%
column_to_rownames('Sample') %>%
as.matrix()
colnames(sf_strain_table.matrix) <-
tibble(tax = colnames(sf_strain_table.matrix)) %>%
separate(col = tax, into = c('bact', 'gdna'), sep = ';') %>%
group_by(bact) %>%
mutate(strain_id = paste0(bact, '.', dense_rank(gdna))) %>%
pull(strain_id)
# mg-OTU species-level matrix
sf_species_table.matrix <-
sf_strain_table %>%
mutate_at(.vars = vars(everything(), -Sample),
.funs = funs(as.numeric(. > 0))) %>%
pivot_longer(names_to = 'tax', values_to = 'rel_abund', cols = starts_with('BACT')) %>%
separate(col = tax, into = c('bact', 'gdna'), sep = ';') %>%
dplyr::select(-gdna) %>%
group_by(Sample, bact) %>%
summarize(rel_abund = as.numeric(sum(rel_abund, na.rm = T) > 0), .groups = 'drop') %>%
pivot_wider(id_cols = 'Sample', names_from = 'bact', values_from = 'rel_abund') %>%
column_to_rownames('Sample') %>%
as.matrix()
# mg-OTU genus-level matrix
sf_genus_table.matrix <-
sf_strain_table %>%
mutate_at(.vars = vars(everything(), -Sample),
.funs = funs(as.numeric(. > 0))) %>%
pivot_longer(names_to = 'tax', values_to = 'rel_abund', cols = starts_with('BACT')) %>%
separate(col = tax, into = c('bact', 'gdna'), sep = ';') %>%
dplyr::select(-gdna) %>%
group_by(Sample, bact) %>%
summarize(rel_abund = as.numeric(sum(rel_abund, na.rm = T) > 0), .groups = 'drop') %>%
left_join(., sf_species_table %>% dplyr::select(bact, genus), by = 'bact') %>%
group_by(Sample, genus) %>%
summarize(rel_abund = as.numeric(sum(rel_abund, na.rm = T) > 0), .groups = 'drop') %>%
pivot_wider(id_cols = 'Sample', names_from = 'genus', values_from = 'rel_abund') %>%
column_to_rownames('Sample') %>%
as.matrix()
Sample/Sample co-occurrence matrix at mg-OTU species- and strain-level
# genus
sf_genus_table.dist <- sf_genus_table.matrix %*% t(sf_genus_table.matrix)
sf_genus_table.long <-
distmat_to_long(distmat = as.data.frame(sf_genus_table.dist),
value_name = 'n.shared_genus',
rm_diag = T)
sf_genus_counts <- diag(sf_genus_table.dist) %>%
as.data.frame() %>%
rownames_to_column('Sample') %>%
rename(n.genus = '.')
# species
sf_species_table.dist <- sf_species_table.matrix %*% t(sf_species_table.matrix)
sf_species_table.long <-
distmat_to_long(distmat = as.data.frame(sf_species_table.dist),
value_name = 'n.shared_species',
rm_diag = T)
sf_species_counts <- diag(sf_species_table.dist) %>%
as.data.frame() %>%
rownames_to_column('Sample') %>%
rename(n.species = '.')
# strain
sf_strain_table.dist <- sf_strain_table.matrix %*% t(sf_strain_table.matrix)
sf_strain_table.long <-
distmat_to_long(distmat = as.data.frame(sf_strain_table.dist),
value_name = 'n.shared_strains',
rm_diag = T)
sf_strain_counts <- diag(sf_strain_table.dist) %>%
as.data.frame() %>%
rownames_to_column('Sample') %>%
rename(n.strains = '.')
join data
sf <-
full_join(sf_species_table.long, sf_strain_table.long, by = c('row','col')) %>%
full_join(., sf_genus_table.long, by = c('row','col')) %>%
left_join(., sf_genus_counts, by = c('row' = 'Sample')) %>%
left_join(., sf_genus_counts, by = c('col' = 'Sample'), suffix = c('.row', '.col')) %>%
left_join(., sf_species_counts, by = c('row' = 'Sample')) %>%
left_join(., sf_species_counts, by = c('col' = 'Sample'), suffix = c('.row', '.col')) %>%
left_join(., sf_strain_counts, by = c('row' = 'Sample')) %>%
left_join(., sf_strain_counts, by = c('col' = 'Sample'), suffix = c('.row', '.col')) %>%
rename(Sample.row = row,
Sample.col = col,
n.shared_species.sf = n.shared_species,
n.shared_strains.sf = n.shared_strains,
n.shared_genus.sf = n.shared_genus) %>%
rename_at(.vars = vars(matches('n.*col'), matches('n.*row')),
.funs = funs(gsub(gsub(., pattern = '.row', replacement = '.sf.row'),
pattern = '.col', replacement = '.sf.col'))) %>%
left_join(., samples.metadata, by = c('Sample.row' = 'Sample')) %>%
left_join(., samples.metadata, by = c('Sample.col' = 'Sample'), suffix = c('.row', '.col')) %>%
filter(!is.na(Name.row), !is.na(Name.col)) %>%
rowwise() %>%
mutate(Name.paste = paste0(sort(c(Name.row, Name.col)), collapse = ',')) %>%
ungroup() %>%
mutate(n.not_shared_strains.sf = (n.strains.sf.row - n.shared_strains.sf) + (n.strains.sf.col - n.shared_strains.sf),
n.not_shared_species.sf = (n.species.sf.row - n.shared_species.sf) + (n.species.sf.col - n.shared_species.sf),
n.not_shared_genus.sf = (n.genus.sf.row - n.shared_genus.sf) + (n.genus.sf.col - n.shared_genus.sf),
n.total_strains.sf = n.shared_strains.sf + n.not_shared_strains.sf,
n.total_species.sf = n.shared_species.sf + n.not_shared_species.sf,
n.total_genus.sf = n.shared_genus.sf + n.not_shared_genus.sf,
f.shared_strains.sf = n.shared_strains.sf / n.total_strains.sf,
f.shared_species.sf = n.shared_species.sf / n.total_species.sf,
f.shared_genus.sf = n.shared_genus.sf / n.total_genus.sf)
Shared Taxa
Density Ridge Plot based on Strain Finder Data published with the Alm dataset, showing co-occurrence of genera, species, and strains in samples which are common in unrelated but more common in related samples.
sf %>%
mutate(related = ifelse(replace_na(Donor.Subject.row == Donor.Subject.col, F), 'Related', 'Unrelated')) %>%
pivot_longer(names_to = 'metric', values_to = 'value',
cols = c('n.shared_genus.sf', 'n.shared_species.sf', 'n.shared_strains.sf')) %>%
mutate(metric = case_when(
metric == 'n.shared_strains.sf' ~ 'StrainFinder/\nmg-OTUs (strain)',
metric == 'n.shared_species.sf' ~ 'StrainFinder/\nmg-OTUs (species)',
metric == 'n.shared_genus.sf' ~ 'StrainFinder/\nmg-OTUs (genus)',
T ~ metric
)) %>%
ggplot(aes(value, fct_relevel(metric,
'StrainFinder/\nmg-OTUs (strain)',
'StrainFinder/\nmg-OTUs (species)',
'StrainFinder/\nmg-OTUs (genus)'), fill = related)) +
geom_density_ridges(panel_scaling = F, alpha = 0.75, scale = 5) +
scale_fill_manual(values = c('black','white')) +
scale_x_continuous(trans = pseudo_log_trans(), breaks = c(0, 3, 10, 30, 100, 300)) +
theme_cowplot() +
labs(x = 'Shared Bacteria') +
theme(legend.key.size = unit(c(.75, .5, .75, .5), "lines"),
legend.title = element_blank(),
axis.title.y = element_blank(),
)
SameStr
Format
Generate MetaPhlAn co-occurrence table, taxa counts
cooccurrence <- mp_cooccurrence(mp_species.long)
taxa_counts <- cooccurrence$mp_taxa_counts
taxa_cooccurrence <- cooccurrence$mp_taxa_table.long
Merge with SameStr strain-level data
sstr <-
samples.sstr_data %>%
group_by(Study.pair, Name.row, Name.col,
Case_Name.pair, Case_Name.row, Case_Name.col, Same_Case.label,
Related.label, Related.bool, Pair_Type) %>%
summarize(.groups = 'drop',
n.analyzed_strains = sum(overlap > min_overlap, na.rm = T),
n.shared_strains.mvs = sum(distance.mvs > min_similarity & overlap > min_overlap, na.rm = T),
n.shared_strains.cvs = sum(distance.cvs > min_similarity & overlap > min_overlap, na.rm = T)) %>%
# tag successful fmt cases
mutate(fmt_success = !Case_Name.pair %in% cases_failed) %>%
mutate(Study.row = str_split_fixed(Name.row, pattern = '_', 2)[,1],
Study.col = str_split_fixed(Name.col, pattern = '_', 2)[,1]) %>%
mutate(Related.text = ifelse(Related.bool, 'Related', 'Unrelated')) %>%
rowwise() %>%
mutate(Sample_Pair = paste0(sort(c(Name.row, Name.col)), collapse = ',')) %>%
ungroup() %>%
left_join(., taxa_cooccurrence, by = 'Sample_Pair') %>%
left_join(., taxa_counts, by = c('Name.row' = 'Sample')) %>%
left_join(., taxa_counts, by = c('Name.row' = 'Sample'), suffix = c('.row','.col')) %>%
mutate(n.not_shared_species = (n.species.row - n.shared_species) + (n.species.col - n.shared_species),
n.not_shared_genus = (n.genus.row - n.shared_genus) + (n.genus.col - n.shared_genus),
n.not_shared_family = (n.family.row - n.shared_family) + (n.family.col - n.shared_family),
n.total_species = n.shared_species + n.not_shared_species,
n.total_genus = n.shared_genus + n.not_shared_genus,
n.total_family = n.shared_family + n.not_shared_family,
f.shared_strains.mvs = replace_na(n.shared_strains.mvs / n.analyzed_strains, 0),
f.shared_strains.cvs = replace_na(n.shared_strains.cvs / n.analyzed_strains, 0),
f.shared_species = n.shared_species / n.total_species,
f.shared_genus = n.shared_genus / n.total_genus,
f.shared_family = n.shared_family / n.total_family)
Subset rCDI Alm data
sstr.alm <-
sstr %>%
# ALM only, same as above
filter(Study.pair == 'ALM')
sstr.control <-
sstr %>%
# Control only
filter(Study.row %in% control.studies &
Study.col %in% control.studies) %>%
mutate(Related.text = ifelse(Related.bool, 'Related', 'Unrelated'))
Shared Taxa
Density Ridge Plot based on SameStr Data for the same Alm dataset shown above, showing co-occurrence of genera, species, and strains in samples. Co-occurrence is common in unrelated genera and species, but more separable compared to SF. For shared strains, distributions divide related from unrelated sample pairs which mostly don’t share any strains.
Alm
sstr.alm %>%
pivot_longer(names_to = 'metric', values_to = 'value',
cols = c('n.shared_genus', 'n.shared_species', 'n.shared_strains.mvs')) %>%
mutate(metric = case_when(
metric == 'n.shared_strains.mvs' ~ 'SameStr/\nMetaPhlAn2 (strain)',
metric == 'n.shared_species' ~ 'SameStr/\nMetaPhlAn2 (species)',
metric == 'n.shared_genus' ~ 'SameStr/\nMetaPhlAn2 (genus)',
T ~ metric
)) %>%
ggplot(aes(value, fct_relevel(metric,
'SameStr/\nMetaPhlAn2 (strain)',
'SameStr/\nMetaPhlAn2 (species)',
'SameStr/\nMetaPhlAn2 (genus)',), fill = Related.bool)) +
geom_density_ridges(alpha = 0.75, scale = 4, stat = 'binline', bins = 18) + # scale = 4
scale_x_continuous(trans = pseudo_log_trans(), breaks = c(0, 3, 10, 30, 100, 300)) +
scale_fill_manual(values = c('black','white')) +
theme_cowplot() +
labs(x = 'Shared Bacteria') +
theme(legend.key.size = unit(c(.75, .5, .75, .5), "lines"),
legend.title = element_blank(),
axis.title.y = element_blank(),
)
Most unrelated strains
samples.sstr_data %>%
filter(Study.pair == 'ALM') %>%
filter(overlap > min_overlap) %>%
group_by(genus, species, Related.bool) %>%
summarize(count = sum(Shared_Strain.bool == T, na.rm = T),
total = sum(!is.na(Shared_Strain.bool))) %>%
mutate(fraction = count / total) %>%
filter(total > 20) %>%
arrange(Related.bool, desc(fraction))
Control
plot.control <-
sstr %>%
filter(Study.row %in% control.studies &
Study.col %in% control.studies) %>%
mutate(Related.text = ifelse(Related.bool, 'Related', 'Unrelated')) %>%
pivot_longer(names_to = 'metric', values_to = 'value',
cols = c('f.shared_strains.mvs', 'f.shared_strains.cvs')) %>%
mutate(metric = case_when(
metric == 'f.shared_strains.mvs' ~ 'SameStr',
metric == 'f.shared_strains.cvs' ~ 'Consensus',
)) %>%
ggplot(aes(value, Related.text, fill = metric)) +
geom_density_ridges(alpha = 0.75, scale = 4, stat = 'binline', bins = 18) + # scale = 4
scale_x_continuous(labels = percent_format(), expand = c(0,0), breaks = c(0, 0.5, 1)) +
scale_fill_manual(values = c('black','white')) +
theme_cowplot() +
labs(x = 'as [%] of analyzed Spp.') +
theme(legend.title = element_blank(),
axis.title.y = element_blank(),
)
plot.control + theme(legend.position = 'none')
Exporting Plot
output_name = 'SharedStrains.Control'
ggsave(plot.control + theme(legend.position = 'none'),
device = 'pdf', dpi = 300, width = 2.85, height = 3.25,
filename = paste0(results_dir, fig, output_name, '.RidgePlot.pdf'))
sstr.control %>%
dplyr::select(Related.text, n.shared_strains.mvs) %>%
group_by(Related.text) %>%
summarize_all(.funs = funs(n(), min, max, median, mean, sd))
Tool Comparison
Format
Merge Tool Data
alm_table <-
sf %>%
left_join(., sstr.alm %>%
rename_at(.vars = vars(n.shared_strains.mvs,
n.shared_species,
n.shared_genus,
n.shared_family,
f.shared_strains.mvs,
f.shared_species,
f.shared_genus,
f.shared_family),
.funs = funs(paste0(gsub(., pattern = '.mvs', replacement = ''), '.sstr'))) %>%
rowwise() %>%
mutate(Name.paste = paste0(sort(c(Name.row, Name.col)), collapse = ',')) %>%
ungroup() %>%
dplyr::select(-Name.row, -Name.col, -Case_Name.row, -Case_Name.col)) %>%
# label study, case
mutate(Case_Split.row = str_split(Case_Name.row, ';'),
Case_Split.col = str_split(Case_Name.col, ';')) %>%
rowwise() %>%
mutate(Study.same = Study.row == Study.col,
Study.pair = ifelse(Study.same, Study.row, 'Unrelated'),
Case_Name.same = (Case_Name.col %in% Case_Split.row) | (Case_Name.row %in% Case_Split.col),
Case_Name.pair = ifelse(Case_Name.same, intersect(Case_Split.row, Case_Split.col), 'Unrelated')) %>%
dplyr::select(-starts_with('Case_Split.')) %>%
ungroup() %>%
# label comparison type + study, case
mutate(
Pair_Type = 'No-Pair',
Pair_Type = case_when(
(Sample_Type.row == 'recipient' & Sample_Type.col == 'donor') |
(Sample_Type.col == 'recipient' & Sample_Type.row == 'donor') ~ 'Pre-FMT/Donor',
(Sample_Type.row == 'post' & Sample_Type.col == 'recipient') |
(Sample_Type.col == 'post' & Sample_Type.row == 'recipient') ~ 'Pre-FMT/Post-FMT',
(Sample_Type.row == 'post' & Sample_Type.col == 'donor') |
(Sample_Type.col == 'post' & Sample_Type.row == 'donor') ~ 'Donor/Post-FMT',
Sample_Type.row == 'recipient' & Sample_Type.col == 'recipient' ~ 'Pre-FMT/Pre-FMT',
Sample_Type.row == 'post' & Sample_Type.col == 'post' ~ 'Post-FMT/Post-FMT',
Sample_Type.row == 'donor' & Sample_Type.col == 'donor' ~ 'Donor/Donor')) %>%
# label family relationship
mutate(Same_Case.label = ifelse(Case_Name.same, 'Same Case', 'Distinct Case'),
Same_Study.label = ifelse(Study.same, 'Same Study', 'Distinct Study'),
) %>%
dplyr::select(Name.paste, contains('n.shared'), starts_with('f.'), ends_with('.pair'),
Related.label, Related.bool, Pair_Type) %>%
mutate(Related.text = ifelse(Related.bool, 'Related', 'Unrelated'))
alm_table.n <-
alm_table %>%
pivot_longer(names_to = 'metric', values_to = 'value',
cols = starts_with('n.shared')) %>%
separate(metric, into = c('metric','taxonomic_level','tool'), sep = '[.]') %>%
mutate(taxonomic_level = gsub(taxonomic_level, pattern = 'shared_', replacement = ''),
taxonomic_level = gsub(taxonomic_level, pattern = 'strains', replacement = 'strain')) %>%
dplyr::select(-metric)
alm_table.f <-
alm_table %>%
pivot_longer(names_to = 'metric', values_to = 'value',
cols = starts_with('f.shared')) %>%
separate(metric, into = c('metric','taxonomic_level','tool'), sep = '[.]') %>%
mutate(taxonomic_level = gsub(taxonomic_level, pattern = 'shared_', replacement = ''),
taxonomic_level = gsub(taxonomic_level, pattern = 'strains', replacement = 'strain')) %>%
dplyr::select(-metric)
Shared Strains n
Number of shared strains between related and unrelated samples of the ALM rCDI cohort based on CVS, MVS of MetaPhlAn2 alignments and Strain Finder on AMPHORA markers
plot.strain_count <-
alm_table.n %>%
filter(taxonomic_level == 'strain') %>%
mutate(tool = case_when(
tool == 'sstr' ~ 'MVS',
tool == 'sf' ~ 'SF',
tool == 'cvs' ~ 'CVS',
T ~ tool)) %>%
mutate(tool = fct_relevel(tool,'SF','CVS','MVS')) %>%
ggplot(aes(value, tool, fill = Related.text)) +
geom_density_ridges(alpha = 0.75, scale = 4, stat = 'binline', bins = 18) +
scale_x_continuous(expand = c(0,0),
trans = pseudo_log_trans(1, base = 10),
breaks = c(0, 1, 10, 30, 100, 300)
) +
theme_cowplot() +
labs(x = 'Count [n]') +
scale_fill_manual(values = c('black', 'white')) +
theme(strip.background = element_blank(),
axis.title.y = element_blank(),
legend.title = element_blank(),
legend.position = c(0.6, 0.9))
legend <- cowplot::get_legend(plot.strain_count)
plot.strain_count + theme(legend.position = 'none')
grid.newpage()
grid.draw(legend)
Whereas Strain Finder is finding a lot of co-occurring strains in unrelated samples, the consensus approach is very conservative. Compared to SameStr MVS, consensus finds slightly fewer potential false-positives (means ~.45 vs. .22 strains in unrelated samples) but also misses on a large fraction of potential true-positives (means ~15 vs. 8 in related samples)
alm_table.n %>%
filter(taxonomic_level == 'strain') %>%
ungroup() %>%
mutate(tool = case_when(
tool == 'sstr' ~ 'MVS',
tool == 'sf' ~ 'SF',
tool == 'cvs' ~ 'CVS',
T ~ tool)) %>%
dplyr::select(Name.paste, Related.text, tool, value) %>%
pivot_wider(names_from = 'tool', values_from = 'value') %>%
mutate(`SF-MVS` = SF - MVS) %>%
dplyr::select(-Name.paste) %>%
group_by(Related.text) %>%
summarize_all(.funs = funs(n(), min, max, median, mean, sd)) %>%
group_by(Related.text) %>%
mutate_all(.funs = funs(round(., 2))) %>%
pivot_longer(names_to = 'names', values_to = 'value', cols = matches('.*_.*')) %>%
separate(names, into = c('tool','metric'), sep = '_') %>%
pivot_wider(names_from = 'metric', values_from = 'value')
Related
CVS 7.99±6.83
MVS 14.77±11.34
SF 78.35±73.50
DIFF 93.13±82.59
Unrelated
CVS 0.22±0.55
MVS 0.45±0.75
SF 34.71±37.52
DIFF 35.16±37.68
Shared Strains %
Fraction of shared strains betweent related and unrelated samples of the ALM rCDI cohort
plot.strain_fraction <-
alm_table.f %>%
filter(taxonomic_level == 'strain') %>%
mutate(tool = case_when(
tool == 'sstr' ~ 'SameStr MVS/\nMetaPhlAn2',
tool == 'sf' ~ 'StrainFinder/\nmg-OTUs',
tool == 'cvs' ~ 'Consensus/\nMetaPhlAn2',
T ~ tool)) %>%
mutate(tool = fct_relevel(tool,
'StrainFinder/\nmg-OTUs',
'Consensus/\nMetaPhlAn2',
'SameStr MVS/\nMetaPhlAn2',
)) %>%
ggplot(aes(value, tool, fill = Related.text)) +
geom_density_ridges(alpha = 0.75, scale = 4, stat = 'binline', bins = 18) +
scale_x_continuous(expand = c(0,0),
trans = pseudo_log_trans(base = 10),
labels = percent_format()
) +
theme_cowplot() +
labs(x = 'Taxa Shared between Samples [%]') +
scale_fill_manual(values = c('black', 'white')) +
theme(strip.background = element_blank(),
axis.title.y = element_blank(),
legend.title = element_blank(),
legend.position = c(0.6, 0.9))
plot.strain_fraction + theme(legend.position = 'none')
Exporting Plot
output_name = 'SharedStrains.Alm'
ggsave(plot.strain_count + theme(legend.position = 'none'),
device = 'pdf', dpi = 300, width = 3.5, height = 3.25,
filename = paste0(results_dir, fig, output_name, '.Counts', '.RidgePlot.pdf'))
ggsave(plot.strain_fraction + theme(legend.position = 'none', axis.line.y = element_blank()),
device = 'pdf', dpi = 300, width = 5, height = 4,
filename = paste0(results_dir, fig, output_name, '.Fraction', '.RidgePlot.pdf'))
ggsave(legend,
device = 'pdf', dpi = 300, width = 5, height = 5,
filename = paste0(results_dir, fig, output_name, '.Legend.pdf'))
Shared F>G>Spp. %
alm_table %>%
pivot_longer(names_to = 'metric', values_to = 'value',
cols = starts_with('f.shared')) %>%
separate(metric, into = c('metric','taxonomic_level','tool'), sep = '[.]') %>%
dplyr::select(-metric) %>%
filter(tool != 'cvs') %>%
mutate(taxonomic_level = str_to_title(gsub(taxonomic_level, pattern = 'shared_', replacement = ''))) %>%
mutate(tool = case_when(tool == 'sstr' ~ 'SameStr/\nMetaPhlAn2', tool == 'sf' ~ 'StrainFinder/\nmg-OTUs', T ~ tool)) %>%
ggplot(aes(value, tool, fill = Related.text)) +
geom_density_ridges(alpha = 0.75, scale = 3) +
scale_x_continuous(limits = c(0, 1), labels = percent_format(), expand = c(0,0)) +
theme_cowplot() +
labs(x = 'Fraction of Taxa Shared between Samples') +
scale_fill_manual(values = c('black', 'white')) +
facet_wrap(~taxonomic_level, nrow = 1, scales = 'free_x') +
theme(strip.background = element_blank(),
axis.title.y = element_blank(),
legend.title = element_blank(),
legend.position = c(0.8, 0.9))
Sensitivity of Taxa Detection
Format
Comparison of genera found with MP2/Amphora
counts <-
full_join(sf_genus_counts, sf_species_counts, by = 'Sample') %>%
mutate(Name = paste0('ALM_', Sample, '.pair')) %>%
full_join(., taxa_counts %>%
rename(Name = 'Sample') %>% dplyr::select(-n.family),
by = 'Name', suffix = c('.sf','.mp')) %>%
filter(!is.na(n.species.mp), !is.na(n.species.sf))
lm_eqn <- function(df){
m <- lm(mp ~ sf, df);
eq <- substitute(italic(y) == a + b %.% italic(x)*","~~italic(r)~"="~rv,
list(a = format(unname(coef(m)[1]), digits = 2),
b = format(unname(coef(m)[2]), digits = 2),
rv = format(sqrt(summary(m)$r.squared), digits = 2)))
as.character(as.expression(eq));
}
counts.data <- counts %>%
pivot_longer(names_to = 'metric', values_to = 'n', cols = starts_with('n.')) %>%
separate(metric, into = c('metric','taxonomic_level','tool'), sep = '[.]') %>%
dplyr::select(-metric) %>%
pivot_wider(names_from = 'tool', values_from = 'n') %>%
mutate(difference = mp - sf)
counts.data.genus <-
counts.data %>%
filter(taxonomic_level == 'genus')
counts.data.species <-
counts.data %>%
filter(taxonomic_level == 'species')
get_stats <- function(df, digits = 2) {
df.mean <- mean(df, na.rm = T)
df.sd <- sd(df, na.rm = T)
df.format <- paste('µ=', round(df.mean, digits = digits), '±', round(df.sd, digits = digits), sep = '')
return(
list(
'mean' = df.mean,
'sd' = df.sd,
'format' = df.format))
}
digits = 1
stats.g.diff <- get_stats(counts.data.genus$difference, digits = digits)
stats.g.mp <- get_stats(counts.data.genus$mp, digits = digits)
stats.g.sf <- get_stats(counts.data.genus$sf, digits = digits)
stats.s.diff <- get_stats(counts.data.species$difference, digits = digits)
stats.s.mp <- get_stats(counts.data.species$mp, digits = digits)
stats.s.sf <- get_stats(counts.data.species$sf, digits = digits)
Statistics
counts.data %>%
dplyr::select(-Sample, -Name) %>%
group_by(taxonomic_level) %>%
summarize_all(.funs = funs(mean, sd)) %>%
dplyr::select(taxonomic_level, starts_with('sf_'), starts_with('mp_'), starts_with('dif')) %>%
mutate_at(.vars = vars(everything(), -taxonomic_level),
.funs = funs(round(., digits = 2)))
genus:
sf = 23.78 ± 16.67
mp = 50.54 ± 15.0
species:
sf = 48.48 ± 33.88
mp = 97.62 ± 39.6
Detected Taxa per Sample
yy.s = 80
yy.g = 30
nudge = 70
plot <-
counts.data %>%
ggplot(aes(y = mp, x = sf, fill = taxonomic_level)) +
geom_point(shape = 21, col = 'black', size = 5) +
annotate(geom = 'text', y = yy.s, x = yy.s + nudge + 5, label = stats.s.diff$format, size = 5, parse = F) +
annotate(geom = 'text', y = yy.g, x = yy.g + nudge, label = stats.g.diff$format, size = 5, parse = F) +
annotate("segment", show.legend = F,
arrow = arrow(length = unit(0.3, "cm"), type = "closed"), size = 1, col = 'black',
y=yy.s, yend=yy.s + stats.s.diff$mean, x=yy.s, xend=yy.s) +
annotate("segment", show.legend = F,
arrow = arrow(length = unit(0.3, "cm"), type = "closed"), size = 1, col = 'black',
y=yy.g, yend=yy.g + stats.g.diff$mean, x=yy.g, xend=yy.g) +
geom_abline(slope = 1, intercept = 0, linetype = 'dashed') +
theme_cowplot() +
scale_x_continuous(limits = c(0, 200)) +
scale_y_continuous(limits = c(0, 200)) +
scale_fill_manual(values = c('grey25', 'white')) +
scale_color_manual(values = c('black', 'black')) +
theme(aspect.ratio = 1,
plot.subtitle = element_text(hjust = 0.5),
legend.title = element_blank(),
legend.position = c(0.02, 0.9)) +
labs(y = 'SameStr / MetaPhlAn2', x = 'Strain Finder / mg-OTUs', subtitle = 'Detected Taxa per Sample')
plot
Exporting Plot
output_name = 'DetectedTaxa.Species.Alm'
ggsave(plot,
device = 'pdf', dpi = 300, width = 3.75, height = 3.75,
filename = paste0(results_dir, fig, output_name, '.CorrPlot.pdf'))
Taxa Differences
Number of detected mg-OTUs:
species == 306
genus == 116
sf_strain_table %>%
mutate_at(.vars = vars(everything(), -Sample),
.funs = funs(as.numeric(. > 0))) %>%
pivot_longer(names_to = 'tax', values_to = 'rel_abund', cols = starts_with('BACT')) %>%
separate(col = tax, into = c('bact', 'gdna'), sep = ';') %>%
dplyr::select(-gdna) %>%
group_by(bact) %>%
summarize(present = as.numeric(sum(rel_abund, na.rm = T) > 0)) %>%
summarize(count = sum(present, na.rm = T))
sf_strain_table %>%
mutate_at(.vars = vars(everything(), -Sample),
.funs = funs(as.numeric(. > 0))) %>%
pivot_longer(names_to = 'tax', values_to = 'rel_abund', cols = starts_with('BACT')) %>%
separate(col = tax, into = c('bact', 'gdna'), sep = ';') %>%
dplyr::select(-gdna) %>%
group_by(Sample, bact) %>%
summarize(rel_abund = as.numeric(sum(rel_abund, na.rm = T) > 0)) %>%
left_join(., sf_species_table %>% dplyr::select(bact, genus), by = 'bact') %>%
group_by(genus) %>%
summarize(present = as.numeric(sum(rel_abund, na.rm = T) > 0)) %>%
summarize(count = sum(present, na.rm = T))
Number of detected MetaPhlAn2 species:
species == 399
genus == 154
mp_species.long %>%
dplyr::select(species, Name, rel_abund) %>%
filter(grepl(Name, pattern = '^ALM'), rel_abund > 0) %>%
filter(!grepl(species, pattern = 'unclassified')) %>%
group_by(species) %>%
summarize(present = sum(rel_abund > 0, na.rm = T)) %>%
summarize(count = sum(present > 0, na.rm = T))
mp_species.long %>%
dplyr::select(genus, Name, rel_abund) %>%
filter(grepl(Name, pattern = '^ALM'), rel_abund > 0) %>%
filter(!grepl(genus, pattern = 'unclassified')) %>%
group_by(genus) %>%
summarize(present = sum(rel_abund > 0, na.rm = T)) %>%
summarize(count = sum(present > 0, na.rm = T))
Format
found_sf <-
sf_strain_table %>%
mutate_at(.vars = vars(everything(), -Sample),
.funs = funs(as.numeric(. > 0))) %>%
pivot_longer(names_to = 'tax', values_to = 'rel_abund', cols = starts_with('BACT')) %>%
separate(col = tax, into = c('bact', 'gdna'), sep = ';') %>%
dplyr::select(-gdna) %>%
group_by(Sample, bact) %>%
summarize(rel_abund = as.numeric(sum(rel_abund, na.rm = T) > 0), .groups = 'drop') %>%
left_join(., sf_species_table %>% dplyr::select(bact, genome, genus), by = 'bact') %>%
group_by(Sample, genus) %>%
summarize(species_count.sf = sum(rel_abund, na.rm = T), .groups = 'drop')
found_mp <-
mp_species.long %>%
dplyr::select(genus, species, Name, rel_abund) %>%
filter(grepl(Name, pattern = '^ALM'), rel_abund > 0) %>%
mutate(genus = gsub(genus, pattern = '_noname', replacement = '')) %>%
group_by(genus, Name) %>%
summarize(species_count.mp = sum(rel_abund > 0, na.rm =T), .groups = 'drop') %>%
mutate(Sample = gsub(str_split_fixed(Name, pattern = '_', 2)[,2], pattern = '.pair', replacement = '')) %>%
dplyr::select(-Name)
taxa.genus.important <-
full_join(found_sf, found_mp) %>%
mutate_at(.vars = vars(starts_with('species_count')),
.funs = funs(ifelse(is.na(.), 0, .))) %>%
group_by(genus) %>%
summarize(species_count.sf = sum(species_count.sf, na.rm = T),
species_count.mp = sum(species_count.mp, na.rm = T), .groups = 'drop') %>%
mutate(total = species_count.mp + species_count.sf) %>%
arrange(desc(total)) %>%
top_n(30, total) %>%
pull(genus)
found <-
full_join(found_sf, found_mp) %>%
mutate_at(.vars = vars(starts_with('species_count')),
.funs = funs(ifelse(is.na(.), 0, .))) %>%
filter(genus %in% taxa.genus.important)
Ridge Plot
plot <-
found %>%
pivot_longer(names_to = 'metric', values_to = 'count', cols = c('species_count.sf', 'species_count.mp')) %>%
mutate(metric = case_when(
metric == 'species_count.sf' ~ 'mg-OTUs',
metric == 'species_count.mp' ~ 'MetaPhlAn2',
)) %>%
ggplot(aes(count, genus, fill = metric)) +
geom_density_ridges(panel_scaling = F, alpha = 0.75, scale = 2) +
scale_x_continuous(trans = pseudo_log_trans()) +
labs(x = 'Detected Species') +
theme_cowplot() +
theme(axis.title.y = element_blank()) +
scale_fill_manual(values = c('black','white'))
plot + theme(legend.position = 'none')
legend <- cowplot::get_legend(plot)
grid.newpage()
grid.draw(legend)
Exporting Plot
output_name = 'DetectedTaxa.Species.Alm'
ggsave(plot + theme(legend.position = 'none'),
device = 'pdf', dpi = 300, width = 4, height = 6,
filename = paste0(results_dir, fig, output_name, '.RidgePlot.pdf'))
ggsave(legend,
device = 'pdf', dpi = 300, width = 1.25, height = 1.5,
filename = paste0(results_dir, fig, output_name, '.Legend.pdf'))
found %>%
ungroup() %>%
dplyr::select(-Sample) %>%
group_by(genus) %>%
summarize_all(.funs = funs(mean, sd)) %>%
mutate(diff = species_count.sf_mean - species_count.mp_mean) %>%
group_by(genus) %>%
mutate_all(.funs = funs(round(., 2))) %>%
ungroup() %>%
rename_at(.vars = vars(starts_with('species_')),
.funs = funs(gsub(., pattern = 'species_count.', replacement = ''))) %>%
dplyr::select(genus, starts_with('sf'), starts_with('mp'), diff) %>%
arrange(diff) %>%
head(10)
Genus SF vs. MP
Streptococcus 2.68±4.02 vs. 8.15±2.90 (-5.47)
Lachnospiraceae 0.00±0.00 vs. 4.83±2.54 (-4.83)
Bacteroides 3.87±4.70 vs. 8.44±4.70 (-4.57)
Lactobacillus 1.41±2.37 vs. 5.89±2.69 (-4.48)
Clostridium 2.43±3.25 vs. 6.18±3.57 (-3.75)
SF Majority Strain
sf_min_abund <- function(min_abund) {
sf_strain_table.matrix <-
sf_strain_table %>%
mutate_at(.vars = vars(everything(), -Sample),
.funs = funs(as.numeric(. > min_abund))) %>% # set > .5 filter here to compare SF majority strains
column_to_rownames('Sample') %>%
as.matrix()
colnames(sf_strain_table.matrix) <-
tibble(tax = colnames(sf_strain_table.matrix)) %>%
separate(col = tax, into = c('bact', 'gdna'), sep = ';') %>%
group_by(bact) %>%
mutate(strain_id = paste0(bact, '.', dense_rank(gdna))) %>%
pull(strain_id)
# strain
sf_strain_table.dist <- sf_strain_table.matrix %*% t(sf_strain_table.matrix)
sf_strain_table.long <-
distmat_to_long(distmat = as.data.frame(sf_strain_table.dist),
value_name = paste0('n.shared_strains.', 'min_abund_', min_abund),
rm_diag = T)
return(sf_strain_table.long)
}
df <- sf_min_abund(0)
for (m in c(5e-3, 5e-2, 5e-1)) {
df <- full_join(df, sf_min_abund(m), by = c('row','col'))
}
related.v <-
alm_table.n %>%
filter(taxonomic_level == 'strain', tool == 'sf') %>%
dplyr::select(Name.paste, Related.text)
sf <-
df %>%
rename(Sample.row = row,
Sample.col = col,
) %>%
rename_at(.vars = vars(matches('n.*col'), matches('n.*row')),
.funs = funs(gsub(gsub(., pattern = '.row', replacement = '.sf.row'),
pattern = '.col', replacement = '.sf.col'))) %>%
left_join(., samples.metadata, by = c('Sample.row' = 'Sample')) %>%
left_join(., samples.metadata, by = c('Sample.col' = 'Sample'), suffix = c('.row', '.col')) %>%
filter(!is.na(Name.row), !is.na(Name.col)) %>%
rowwise() %>%
mutate(Name.paste = paste0(sort(c(Name.row, Name.col)), collapse = ',')) %>%
ungroup() %>%
left_join(., related.v) %>%
# mutate(related = ifelse(replace_na(Donor.Subject.row == Donor.Subject.col, F), 'Related', 'Unrelated')) %>%
pivot_longer(names_to = 'metric', values_to = 'value',
cols = starts_with('n.shared_strains')) %>%
separate(metric, into = c('metric','min_abund'), sep = '.abund_') %>%
mutate(min_abund = as.numeric(min_abund))
Strain Finder shows a high number of cooccurring strains in unrelated sample pairs. Cooccurrence distributions of related / unrelated sample pairs overlap widely, even when considering only strains that were detected at >=50% strain abundance in the species population - a problem that might be due to the amphora markers.
plot <-
sf %>%
mutate(metric = paste0('>=', round(min_abund * 100, 5), '%')) %>%
ggplot(aes(value, fct_reorder(metric, min_abund, .desc = T), fill = Related.text)) +
geom_density_ridges(alpha = 0.75, scale = 2) +
scale_fill_manual(values = c('black','white')) +
scale_x_continuous(expand = c(0,0),
trans = pseudo_log_trans(1, base = 10),
breaks = c(0, 1, 10, 30, 100, 300)
) +
theme_cowplot() +
theme(aspect.ratio = 1) +
labs(x = 'Shared Strains', y = 'Strain Finder\nStrain Abundance') +
theme(legend.title = element_blank(),
legend.position = 'top')
plot
Exporting plot
output_name = 'StrainFinder.StrainCooccurrence'
ggsave(plot,
device = 'pdf', dpi = 300, width = 4, height = 4,
filename = paste0(results_dir, fig, output_name, '.Plot.pdf'))
Mock Multi-Strain Species Population
Detecting strains hidden in Mock Multi-Strain Species Population with varying numbers of strains, target coverage, noise coverage and total coverage.
Format
percent_format_signif <-
function(x, label.precision=3, label.percent=T) {
if (label.percent == T) {
paste0(signif(100 * x, label.precision), "%")
} else {
signif(digits = label.precision)
}
}
Calculate metrics
df <-
mock.sstr_data %>%
pivot_longer(names_to = 'metric',
values_to = 'distance',
cols = c('distance.mvs', 'distance.cvs')) %>%
mutate(shared = distance > min_similarity & overlap > min_overlap,
should = tsc > 0,
TP = shared & should,
FP = shared & !should,
TN = !shared & !should,
FN = !shared & should,
event = case_when(
TP ~ 'TP',
FP ~ 'FP',
TN ~ 'TN',
FN ~ 'FN',
T ~ 'Unknown'
)) %>%
mutate(metric = case_when(
metric == 'distance.cvs' ~ 'Consensus (CVS)',
metric == 'distance.mvs' ~ 'SameStr (MVS)',
T ~ 'Unknown'
)) %>%
group_by(metric, tsc, ngc, total_coverage, rel_abund.target) %>%
summarize(TP = sum(TP, na.rm = T),
FP = sum(FP, na.rm = T),
TN = sum(TN, na.rm = T),
FN = sum(FN, na.rm = T),
.groups = 'drop') %>%
mutate(Precision = TP / (TP + FP),
Recall = TP / (TP + FN),
`Precision/Recall` = Precision / Recall,
Accuracy = (TP+TN)/(TP+FP+FN+TN),
`Error Rate` = 1 - Accuracy)
N Mock Observations
nrow(mock.sstr_data)
Accuracy Heatmap
plot <-
df %>%
mutate(`Error Rate` =
Hmisc::cut2(`Error Rate`,
cuts = seq(0, 1, 0.2),
minmax = T,
formatfun = scales::percent,
oneval = F) %>%
as.character() %>%
gsub(pattern = '\\[', replacement = '') %>%
gsub(pattern = '\\]', replacement = '') %>%
gsub(pattern = '\\)', replacement = '') %>%
gsub(pattern = ',', replacement = '-')) %>%
filter(ngc > 0) %>%
ggplot(aes(
tsc,
as.factor(ngc))
) +
geom_tile(aes(fill = `Error Rate`), col = 'black') +
scale_x_continuous(expand = c(0,0),
breaks = seq(0, 20, 4)) +
scale_y_discrete(expand = c(0,0)) +
facet_wrap(~ metric, scales = 'free_y') +
theme_cowplot() +
theme(aspect.ratio = 1,
strip.background = element_blank(),
panel.spacing.x = unit(1, "lines"),
axis.text = element_text(size = 9),
axis.title = element_text(size = 12),
axis.title.x = element_text(hjust = 0.5),
legend.position = 'right',
legend.title = element_text(size = 12),
legend.key = element_rect(colour = "black")) +
labs(x = 'Mean Fold-Coverage of the Target Strain',
y = 'Noise Coverage') +
scale_fill_manual(values = c('white', 'gray75', 'gray50', 'gray25', 'black'))
plot + theme(legend.position = 'none')
legend <- cowplot::get_legend(plot)
grid.newpage()
grid.draw(legend)
Exporting Plot
output_name = 'AccuracyHeatmap.Mock'
ggsave(plot + theme(legend.position = 'none'),
device = 'pdf', dpi = 300, width = 5.5, height = 3,
filename = paste0(results_dir, fig, output_name, '.Plot.pdf'))
ggsave(legend,
device = 'pdf', dpi = 300, width = 1.25, height = 1.5,
filename = paste0(results_dir, fig, output_name, '.Legend.pdf'))
Comparison of average detection accuracy using Consensus vs. SameStr MVS with at least 5x coverage when target strain is dominant (>50% species population) or sub-dominant (15-50% species population)
df %>%
mutate(miss_rate = TP / (FN + TP)) %>%
filter(rel_abund.target > .5, total_coverage > 5) %>%
group_by(metric) %>%
summarize(Accuracy = 1 - mean(`Error Rate`, na.rm = T))
df %>%
filter(rel_abund.target > .15) %>%
filter(rel_abund.target <= .5, total_coverage > 5) %>%
group_by(metric) %>%
summarize(Accuracy = 1 - mean(`Error Rate`, na.rm = T))
Phylogenetic Resolution of MP Markers
Comparing full-genome based Fast-ANI scores to MetaPhlAn2 based similarity
Reference Genomes
Similarity Distribution of Reference Genomes
percent_format_signif <- function(x, label.precision=3, label.percent=T) {
if (label.percent == T) {
paste0(signif(100 * x, label.precision), "%")
} else {
signif(digits = label.precision)
}
}
species_ridge <-
refs %>%
filter(distance > .97 | ani > .97) %>%
pivot_longer(values_to = 'value', names_to = 'metric', cols = c('distance', 'ani')) %>%
mutate(metric = case_when(
metric == 'distance' ~ 'Nucleotide Similarity of\nMetaPhlAn2 Markers',
metric == 'ani' ~ 'fastANI Identity of\nFull Sequences',
T ~ metric
)) %>%
ggplot(aes(x = value,
y = fct_reorder(genus, genus_avg_marker_len),
group = species,
fill = species)) +
geom_density_ridges(scale = 3,
alpha = 0.5, stat = 'binline', show.legend = F,
) +
scale_x_reverse(labels = percent_format_signif,
breaks=c(min_similarity, 0.99, 0.98, 0.97),
limits = c(NA, 0.97),
expand = c(0.1, 0)) +
theme_cowplot() +
geom_vline(xintercept = min_similarity, show.legend = T, linetype = 'dashed') +
guides(fill = guide_legend(ncol = 1)) +
labs(x = 'Sequence Identity') +
facet_grid(cols = vars(metric)) +
theme(strip.background = element_blank(),
axis.title.y = element_blank())
species_ridge
refs %>%
pivot_longer(values_to = 'value', names_to = 'metric', cols = c('distance', 'ani')) %>%
group_by(metric) %>%
summarize(not_shared = sum(value < min_similarity, na.rm = T),
shared = sum(value > min_similarity, na.rm = T)) %>%
mutate(total = not_shared + shared,
rate = shared / total * 100)
Correlation of MP marker vs. fastANI based similarity
min = .97
refs.cor <-
cor.test(~ ani + distance,
data = refs %>%
filter(!is.na(ani), !is.na(distance)) %>%
filter(ani > min & distance > min),
method = 'spearman') %>%
broom::tidy(.) %>%
mutate(estimate.round = round(estimate, 2),
sig = cut(p.value,
breaks = c(-Inf, 0.001, 0.01, 0.05, Inf),
labels = c("***", "**", "*", "")))
refs.cor
plot.scatter <-
ggscatter(data = refs %>%
filter(!is.na(ani), !is.na(distance)) %>%
filter(ani > min & distance > min),
'ani', 'distance',
color = 'black',
fill = 'black',
alpha = 0.025,
xlab = 'Full-Sequence fastANI', ylab = 'MetaPhlAn2 Markers',
ggtheme = theme_cowplot()) +
theme(aspect.ratio = 1) +
scale_x_continuous(labels = percent_format_signif, breaks=c(min_similarity, 0.99, 0.98, 0.97)) +
scale_y_continuous(labels = percent_format_signif, breaks=c(min_similarity, 0.99, 0.98, 0.97)) +
annotate('text', label = paste0('r=',
refs.cor$estimate.round,
' ',
refs.cor$sig),
hjust = 0,
x = .97, y = .999, size = 6) +
geom_density_2d(col = 'black')
plot <-
ggExtra::ggMarginal(
plot.scatter,
type = 'density',
color = 'black',
fill = 'grey75',
size = 10)
plot
Exporting Plots
output_name = 'MetaPhlAn2_FastANI.References'
ggsave(species_ridge,
device = 'pdf', dpi = 300, width = 6, height = 8,
filename = paste0(results_dir, fig, output_name, '.RidgePlot.pdf'))
ggsave(plot,
device = 'pdf', dpi = 300, width = 4, height = 4,
filename = paste0(results_dir, fig, output_name, '.CorrPlot.pdf'))
