.Rhistory

)
#Plot abundance vs BD and save the plots in the current path
incorporator_abs_abundance = incorporator_abs_abundance %>%
dplyr::mutate(
plots = purrr::map(.data$summary_coverage, ~ ggplot2::ggplot(data = .x, aes(x = mean_BD, y = mean_abs_abundance)) +
geom_point(aes(color = Isotope)) +
geom_line(aes(color = Isotope)) +
geom_errorbar(aes(ymin = mean_abs_abundance - sd_abs_abundance,
ymax = mean_abs_abundance + sd_abs_abundance)) +
ylab("Mean absolute \n abundance (attamole/uL)") +
xlab("Mean buoyant \n density (g/mL)") +
theme_bw()))
View(incorporator_abs_abundance$plots[[1]])
incorporator_abs_abundance$plots[[1]]
incorporator_abs_abundance$plots[[2]]
library(xfun)
detach("package:xfun", unload = TRUE)
remove.packages("xfun", lib="~/R/win-library/4.1")
install.packages("xfun")
install.packages("xfun")
library(qSIPmg)
library(qSIPmg)
library(qSIPmg)
library(qSIPmg)
usethis::use_vignette("qSIPmg-vignette")
library(qSIPmg)
library(qSIPmg)
EBImage::readImage("rlm-example.png")
EBImage::readImage("rlm-example.png")
install.packages("~/GitHub/qSIPmg_1.0.0.tar.gz", repos = NULL, type = "source")
browseVignettes("mypackage")
browseVignettes("qSIPmg-vignette")
library(qSIPmg)
library(qSIPmg)
library(qSIPmg)
remove.packages("qSIPmg")
library(qSIPmg)
?qSIPmg::scale_features_lm()
dir
dir()
setwd("~/GitHub/qSIPmg")
library(qSIPmg)
?qSIPmg::scale_features_lm
install.packages("~/GitHub/qSIPmg/qSIPmg_1.0.0.tar.gz", repos = NULL, type = "source")
detach("package:qSIPmg", unload = TRUE)
remove.packages("qSIPmg", lib="~/R/win-library/4.1")
install.packages("~/GitHub/qSIPmg/qSIPmg_1.0.0.tar.gz", repos = NULL, type = "source")
libarary(qSIPmg)
library(qSIPmg)
library(qSIPmg)
library(qSIPmg)
install.packages("roxygen2")
library(SIPmg)
install.packages("~/GitHub/SIPmg_1.0.0.tar.gz", repos = NULL, type = "source")
install.packages("bookdown")
library(SIPmg)
install.packages("EBImage")
#Load required libraries
library(tidyverse)
#Load required libraries
library(tidyverse)
library(phyloseq)
library(HTSSIP)
library(ggpubr)
library(SIPmg)
set.seed(seed = 1000)
## Load data
#Coverage metadata
#Uncomment if your coverage data is in the format mentioned above for this file. Remains commented if you are using the output from `checkm coverage`
f_tibble <- readr::read_csv("mock_input_data/coverage_metadata.csv")
#Sequins metadata
sequins <- readr::read_csv(file="mock_input_data/sequins_metadata.csv")
#Dilutions data
seq_dil = readr::read_csv(file = "mock_input_data/dilutions_data.csv")
#Log scale BOOLEAN. True or False depending on how you would want the MAG coverages to be scaled. Select TRUE if you need MAG concentrations scaled on the log scale
log_scale = TRUE
#coe_of_variation. Acceptable coefficient of variation for coverage and detection (eg. 20 - for 20 % threshold of coefficient of variation) (Coverages above the threshold value will be flagged in the plots)
coe_of_variation = 20
#Taxonomy
gtdbtk_bac_summary = readr::read_delim("mock_input_data/gtdbtk.bac120.summary.tsv",
"\t", escape_double = FALSE, trim_ws = TRUE)
gtdbtk_archaea = readr::read_delim("mock_input_data/gtdbtk.ar122.summary.tsv",
"\t", escape_double = FALSE, trim_ws = TRUE)
#GC content
GC_content <- readr::read_csv(file = "mock_input_data/GC_content.csv")
#Fractions
fractions = readr::read_csv("mock_input_data/fractions.csv")
taxonomy_tibble = dplyr::bind_rows(gtdbtk_bac_summary, gtdbtk_archaea) #Combine bacteria and archaea taxonomy files if it has not been done yet
#mag_tab is a tibble with absolute concentrations of MAGs obtained by scaling MAG coverages using linear regression models on sequin coverages and concentration
##Scale MAG coverages to obtain MAG absolute concentrations and save scaling plots in the working directory
#For rlm scaling using scale_features_rlm
#For rlm scaling using scale_features_lm
mag_tab_scaled <- SIPmg::scale_features_rlm(f_tibble, sequins, seq_dil, log_scale, coe_of_variation = coe_of_variation)
mag_tab = as.matrix(mag_tab_scaled$mag_tab) #Extract absolute abundances as a matrix
##Load data
#Coverage data
f_tibble <- readr::read_csv("mock_input_data/coverages_outliers.csv")
#Sequins metadata
sequins <- readr::read_csv(file="mock_input_data/sequin_meta_outliers.csv")
#Dilutions data
seq_dil = readr::read_csv(file = "mock_input_data/seq_dilution_outliers.csv")
#Log scale BOOLEAN. True or False depending on how you would want the MAG coverages to be scaled. Select TRUE if you need MAG concentrations scaled on the log scale
log_scale = TRUE
#coe_of_variation. Acceptable coefficient of variation for coverage and detection (eg. 20 - for 20 % threshold of coefficient of variation) (Coverages above the threshold value will be flagged in the plots)
coe_of_variation = 50
mag_tab_scaled_rlm <- SIPmg::scale_features_rlm(f_tibble, sequins, seq_dil, log_scale, coe_of_variation = coe_of_variation)
mag_tab_scaled_lm <- SIPmg::scale_features_lm(f_tibble, sequins, seq_dil, log_scale, coe_of_variation = coe_of_variation, cook_filtering = TRUE)
rlm_example = EBImage::readImage("rlm-example.png")
rlm_example = EBImage::resize(rlm_example,dim(rlm_example)[1]/2)
lm_example = EBImage::readImage("lm-example.png")
lm_example = EBImage::resize(lm_example,dim(lm_example)[1]/2)
filtered_lm_example = EBImage::readImage("filtered_lm-example.png")
fitlered_lm_example = EBImage::resize(filtered_lm_example,dim(filtered_lm_example)[1]/2)
cooksd_example = EBImage::readImage("cooksd-example.png")
cooksd_example = EBImage::resize(cooksd_example,dim(cooksd_example)[1]/2)
# Robust linear regression plot
EBImage::display(rlm_example)
#Linear regression plot without filtering sequin data
EBImage::display(lm_example)
taxonomy_tibble = dplyr::bind_rows(gtdbtk_bac_summary, gtdbtk_archaea) #Combine bacteria and archaea taxonomy files if it has not been done yet
#mag_tab is a tibble with absolute concentrations of MAGs obtained by scaling MAG coverages using linear regression models on sequin coverages and concentration
##Scale MAG coverages to obtain MAG absolute concentrations and save scaling plots in the working directory
#For rlm scaling using scale_features_rlm
#For rlm scaling using scale_features_lm
mag_tab_scaled_lm <- SIPmg::scale_features_lm(f_tibble, sequins, seq_dil, log_scale, coe_of_variation = coe_of_variation, cook_filtering = TRUE)
mag_tab = as.matrix(mag_tab_scaled_lm$mag_tab) #Extract absolute abundances as a matrix
mag.table = phyloseq::otu_table(mag_tab, taxa_are_rows = TRUE) #Phyloseq OTU table
taxonomy.object = SIPmg::tax.table(taxonomy_tibble) # Create a taxonomy phyloseq object
samples.object = SIPmg::sample.table(fractions) # Create a samples phyloseq object
phylo.qSIP = SIPmg::phylo.table(mag.table, taxonomy.object, samples.object) # Make a phyloseq table for downstream qSIP analysis
taxonomy_tibble = dplyr::bind_rows(gtdbtk_bac_summary, gtdbtk_archaea) #Combine bacteria and archaea taxonomy files if it has not been done yet
taxonomy.object = SIPmg::tax.table(taxonomy_tibble) # Create a taxonomy phyloseq object
mag.table = phyloseq::otu_table(mag_tab, taxa_are_rows = TRUE) #Phyloseq OTU table
samples.object = SIPmg::sample.table(fractions) # Create a samples phyloseq object
phylo.qSIP = SIPmg::phylo.table(mag.table, taxonomy.object, samples.object) # Make a phyloseq table for downstream qSIP analysis
taxa_names(phylo.qSIP)
samples.object
taxonomy.object
mag.table = phyloseq::otu_table(mag_tab, taxa_are_rows = TRUE) #Phyloseq OTU table
taxonomy.object = SIPmg::tax.table(taxonomy_tibble)
samples.object = SIPmg::sample.table(fractions) # Create a samples phyloseq object
taxa_names(mag.table)
taxa_names(taxonomy.object)
## Load data
#Coverage metadata
#Uncomment if your coverage data is in the format mentioned above for this file. Remains commented if you are using the output from `checkm coverage`
f_tibble <- readr::read_csv("mock_input_data/coverage_metadata.csv")
#Sequins metadata
sequins <- readr::read_csv(file="mock_input_data/sequins_metadata.csv")
#Dilutions data
seq_dil = readr::read_csv(file = "mock_input_data/dilutions_data.csv")
#Log scale BOOLEAN. True or False depending on how you would want the MAG coverages to be scaled. Select TRUE if you need MAG concentrations scaled on the log scale
log_scale = TRUE
#coe_of_variation. Acceptable coefficient of variation for coverage and detection (eg. 20 - for 20 % threshold of coefficient of variation) (Coverages above the threshold value will be flagged in the plots)
coe_of_variation = 20
#Taxonomy
gtdbtk_bac_summary = readr::read_delim("mock_input_data/gtdbtk.bac120.summary.tsv",
"\t", escape_double = FALSE, trim_ws = TRUE)
gtdbtk_archaea = readr::read_delim("mock_input_data/gtdbtk.ar122.summary.tsv",
"\t", escape_double = FALSE, trim_ws = TRUE)
#GC content
GC_content <- readr::read_csv(file = "mock_input_data/GC_content.csv")
#Fractions
fractions = readr::read_csv("mock_input_data/fractions.csv")
taxonomy_tibble = dplyr::bind_rows(gtdbtk_bac_summary, gtdbtk_archaea) #Combine bacteria and archaea taxonomy files if it has not been done yet
#mag_tab is a tibble with absolute concentrations of MAGs obtained by scaling MAG coverages using linear regression models on sequin coverages and concentration
##Scale MAG coverages to obtain MAG absolute concentrations and save scaling plots in the working directory
#For rlm scaling using scale_features_rlm
#For rlm scaling using scale_features_lm
mag_tab_scaled_lm <- SIPmg::scale_features_lm(f_tibble, sequins, seq_dil, log_scale, coe_of_variation = coe_of_variation, cook_filtering = TRUE)
mag_tab = as.matrix(mag_tab_scaled_lm$mag_tab) #Extract absolute abundances as a matrix
mag.table = phyloseq::otu_table(mag_tab, taxa_are_rows = TRUE) #Phyloseq OTU table
taxonomy.object = SIPmg::tax.table(taxonomy_tibble) # Create a taxonomy phyloseq object
samples.object = SIPmg::sample.table(fractions) # Create a samples phyloseq object
phylo.qSIP = SIPmg::phylo.table(mag.table, taxonomy.object, samples.object) # Make a phyloseq table for downstream qSIP analysis
atomX = SIPmg::qSIP_atom_excess_MAGs(phylo.qSIP,
control_expr='Isotope=="12C"',
treatment_rep='Replicate',
Gi = GC_content)
#Bootstrap confidence intervals
ncores = 6
doParallel::registerDoParallel(ncores)
df_atomX_boot = SIPmg::qSIP_bootstrap_fcr(atomX, n_boot=10, parallel = T) #Change "parallel = FALSE" to compute using a single-core
#GC content
GC_content <- readr::read_csv(file = "mock_input_data/GC_content.csv")
#Fractions
fractions = readr::read_csv("mock_input_data/fractions.csv")
df_atomX_boot = SIPmg::qSIP_bootstrap_fcr(atomX, n_boot=10, parallel = T) #Change "parallel = FALSE" to compute using a single-core
atomX = SIPmg::qSIP_atom_excess_MAGs(phylo.qSIP,
control_expr='Isotope=="12C"',
treatment_rep='Replicate',
Gi = GC_content)
df_atomX_boot = SIPmg::qSIP_bootstrap_fcr(atomX, n_boot=10) #Change "parallel = FALSE" to compute using a single-core
, parallel = T
df_atomX_boot = SIPmg::qSIP_bootstrap_fcr(atomX, n_boot=10, parallel = T) #Change "parallel = FALSE" to compute using a single-core
df_atomX_boot = SIPmg::qSIP_bootstrap_fcr(atomX, n_boot=10) #Change "parallel = FALSE" to compute using a single-core
CI_threshold = 0
df_atomX_boot = df_atomX_boot %>%
dplyr::mutate(Incorporator_qSIP = A_CI_fcr_low > CI_threshold,
Incorporator_delbd = A_delbd - A_delbd_sd > 0,
OTU = reorder(OTU, -A))
View(df_atomX_boot)
windows = data.frame(density_min=c(1.71,1.72, 1.73),
density_max=c(1.74,1.75,1.76))
padj_cutoff = 0.05
mw.hr.sip = qSIPmg::HRSIP(physeq = phylo.qSIP, design = ~Isotope,
density_windows = windows,
sparsity_threshold = seq(0, 0.3, 0.05),
padj_cutoff = padj_cutoff,
parallel = T)
mw.hr.sip = qSIPmg::HRSIP(physeq = phylo.qSIP, design = ~Isotope,
density_windows = windows,
sparsity_threshold = seq(0, 0.3, 0.05),
padj_cutoff = padj_cutoff)
mw.hr.sip = mw.hr.sip %>%
mutate(incorp = padj < padj_cutoff)
View(mw.hr.sip)
#Get incorporator info
n_qSIP_incorp = df_atomX_boot %>%
dplyr::filter(Incorporator_qSIP == TRUE) %>%
nrow
n_qSIP_incorp
df_atomX_boot = df_atomX_boot %>%
inner_join(taxonomy_tibble)
View(taxonomy_tibble)
df_atomX_boot = df_atomX_boot %>%
inner_join(taxonomy_tibble %>%
select(user_genome, classification) %>%
rename(OTU = user_genome))
??qSIPmg::incorporators_taxonomy()
?qSIPmg::incorporators_taxonomy()
SIPmg::incorporators_taxonomy(df_atomX_boot)
SIPmg::incorporators_taxonomy(bootstrapped_AFE_table = df_atomX_boot, taxonomy = taxonomy_tibble)
qSIP_incorp = df_atomX_boot %>%
select(OTU, classification, A, A_sd, Incorporator_qSIP) %>%
filter(Incorporator_qSIP == TRUE)
qSIP_incorp
View(qSIP_incorp)
n_qSIP_incorp = nrow(qSIP_incorp)
n_qSIP_incorp
delbd_incorp = df_atomX_boot %>%
select(OTU, classification, A_delbd, A_delbd_sd, Incorporator_delbd) %>%
filter(Incorporator_delbd == TRUE)
mw.hr.sip_incorp = mw.hr.sip %>%
filter(incorp == TRUE)
mw.hr.sip_incorp = mw.hr.sip %>%
select(OTU, taxa, incorp) %>%
filter(incorp == TRUE)
n_mw.hr.sip_incorp = nrow(n_mw.hr.sip_incorp)
n_mw.hr.sip_incorp = nrow(mw.hr.sip_incorp)
all_incorp_tibble = full_join(qSIP_incorp, full_join(delbd_incorp, mw.hr.sip_incorp, by = "OTU"), by = "OTU")
View(all_incorp_tibble)
#Get incorporator info
qSIP_incorp = df_atomX_boot %>%
select(OTU, classification, A, A_sd, Incorporator_qSIP) %>%
filter(Incorporator_qSIP == TRUE) %>%
select(-classification)
n_qSIP_incorp = nrow(qSIP_incorp)
delbd_incorp = df_atomX_boot %>%
select(OTU, classification, A_delbd, A_delbd_sd, Incorporator_delbd) %>%
filter(Incorporator_delbd == TRUE) %>%
select(-classification)
n_delbd_incorp = nrow(delbd_incorp)
mw.hr.sip_incorp = mw.hr.sip %>%
select(OTU, taxa, incorp) %>%
filter(incorp == TRUE) %>%
select(-taxa)
n_mw.hr.sip_incorp = nrow(mw.hr.sip_incorp)
all_incorp_tibble = full_join(qSIP_incorp, full_join(delbd_incorp, mw.hr.sip_incorp, by = "OTU"), by = "OTU")
View(all_incorp_tibble)
mw.hr.sip_incorp = mw.hr.sip %>%
select(OTU, taxa, incorp) %>%
filter(incorp == TRUE) %>%
rename("Incorporator_mw_hr.sip" = "incorp") %>%
select(-taxa)
View(mw.hr.sip_incorp)
n_mw.hr.sip_incorp = nrow(mw.hr.sip_incorp)
all_incorp_tibble = full_join(qSIP_incorp, full_join(delbd_incorp, mw.hr.sip_incorp, by = "OTU"), by = "OTU")
install.packages("devtools")
devtools::build_vignettes()
devtools::build_vignettes()
install.packages("knitr")
devtools::build_vignettes()
tools::pkgVignettes()
tools::pkgVignettes(package = ".")
tools::pkgVignettes(package = "SIPmg")
library(SIPmg)
devtools::build_vignettes()
devtools::build_vignettes(pkg = ".")
devtools::build_vignettes(pkg = ".")
devtools::build_vignettes(pkg = ".")
devtools::build_vignettes(pkg = ".")
?SIPmg::coverage_normalization()
rel.cov = SIPmg::coverage_normalization(f_tibble = f_tibble, approach = "relative_coverage")
f_tibble <- readr::read_csv("mock_input_data/coverage_metadata.csv")
f_tibble <- readr::read_csv("mock_input_data/coverage_metadata.csv")
rel.cov = SIPmg::coverage_normalization(f_tibble = f_tibble, approach = "relative_coverage")
mag.table = phyloseq::otu_table(as.matrix(rel.cov), taxa_are_rows = TRUE) #Phyloseq OTU table
View(rel.cov)
library(tidyverse)
mag.table = phyloseq::otu_table(as.matrix(rel.cov %>% column_to_rownames(var = "Feature")), taxa_are_rows = TRUE) #Phyloseq OTU table
library(SIPmg)
devtools::build_vignettes()
devtools::build_vignettes()
library(SIPmg)
2+2
?rename()
?reorder()
library(SIPmg)
?str_remove
?unique
devtools::load_all(".")
library(SIPmg)
?as_tibble()
?as_tibble
library(SIPmg)
?column_to_rownames()
library(SIPmg)
library(SIPmg)
library(SIPmg)
library(SIPmg)
library(SIPmg)
library(readr)
contains
?contains()
?if_else()
if (!requireNamespace("ggplot2", quietly = TRUE))
if (!requireNamespace("ggplot2", quietly = TRUE))
if (!requireNamespace("ggplot2", quietly = TRUE))
list.of.packages <- c("tidyverse", "phyloseq", "HTSSIP", "ggpubr")
new.packages <- list.of.packages[!(list.of.packages %in% installed.packages()[,"Package"])]
if(length(new.packages)) install.packages(new.packages)
R CMD check ./SIPmg_1.2.tar.gz
checking PDF version of manual without index ... ERROR
tinytex::install_tinytex()
library(SIPmg)
R CMD build ./SIPmg_1.2
devtools::check_win_devel()
devtools::build_manual()
rlang::last_error()
rlang::last_trace()
devtools::build_manual(manual = T)
tools::texi2pdf()
tools::texi2pdf(file = "vignettes_temp/SIPmg-vignette.html")
devtools::check_rhub()
tools::texi2pdf()
tinytex::install_tinytex()
devtools::build_manual()
devtools::build_manual(".")
devtools::document(roclets = c('rd', 'collate', 'namespace', 'vignette'))
devtools::manual(roclets = c('rd', 'collate', 'namespace', 'vignette'))
devtools::build_manual(roclets = c('rd', 'collate', 'namespace', 'vignette'))
devtools::build_manual(path = ".")
devtools::build_manual(path = ".")
rlang::last_error()`
rlang::last_error()
tinytex::tlmgr_install("makeindex")
## Load data
#Coverage metadata
#Uncomment if your coverage data is in the format mentioned above for this file. Remains commented if you are using the output from `checkm coverage`
f_tibble <- readr::read_csv("mock_input_data/coverage_metadata.csv")
#Sequins metadata
sequins <- readr::read_csv(file="mock_input_data/sequins_metadata.csv")
#Dilutions data
seq_dil = readr::read_csv(file = "mock_input_data/dilutions_data.csv")
#Log scale BOOLEAN. True or False depending on how you would want the MAG coverages to be scaled. Select TRUE if you need MAG concentrations scaled on the log scale
log_scale = TRUE
#coe_of_variation. Acceptable coefficient of variation for coverage and detection (eg. 20 - for 20 % threshold of coefficient of variation) (Coverages above the threshold value will be flagged in the plots)
coe_of_variation = 20
#Taxonomy
gtdbtk_bac_summary = readr::read_delim("mock_input_data/gtdbtk.bac120.summary.tsv",
"\t", escape_double = FALSE, trim_ws = TRUE)
gtdbtk_archaea = readr::read_delim("mock_input_data/gtdbtk.ar122.summary.tsv",
"\t", escape_double = FALSE, trim_ws = TRUE)
#GC content
GC_content <- readr::read_csv(file = "mock_input_data/GC_content.csv")
#Fractions
fractions = readr::read_csv("mock_input_data/fractions.csv")
taxonomy_tibble = dplyr::bind_rows(gtdbtk_bac_summary, gtdbtk_archaea) #Combine bacteria and archaea taxonomy files if it has not been done yet
#mag_tab is a tibble with absolute concentrations of MAGs obtained by scaling MAG coverages using linear regression models on sequin coverages and concentration
##Scale MAG coverages to obtain MAG absolute concentrations and save scaling plots in the working directory
#For rlm scaling using scale_features_rlm
#For rlm scaling using scale_features_lm
mag_tab_scaled <- SIPmg::scale_features_rlm(f_tibble, sequins, seq_dil, log_scale, coe_of_variation = coe_of_variation, save_plots = FALSE)
#Load required libraries
list.of.packages <- c("tidyverse", "phyloseq", "HTSSIP", "ggpubr","data.table")
new.packages <- list.of.packages[!(list.of.packages %in% installed.packages()[,"Package"])]
if(length(new.packages)) install.packages(new.packages, quiet = T, dependencies = T, repos = "http://cran.us.r-project.org")
if (!require("BiocManager", quietly = TRUE))
install.packages("BiocManager")
BiocManager::install("EBImage")
library(tidyverse)
library(phyloseq)
library(HTSSIP)
library(ggpubr)
library(SIPmg)
set.seed(seed = 1000)
library(SIPmg)
library(tidyverse)
library(phyloseq)
library(HTSSIP)
library(ggpubr)
library(SIPmg)
library(tidyverse)
library(phyloseq)
library(HTSSIP)
library(ggpubr)
library(SIPmg)
##Scale MAG coverages to obtain MAG absolute concentrations and save scaling plots in the working directory
#For rlm scaling using scale_features_rlm
#For rlm scaling using scale_features_lm
mag_tab_scaled <- SIPmg::scale_features_rlm(f_tibble, sequins, seq_dil, log_scale, coe_of_variation = coe_of_variation, save_plots = FALSE)
##Scale MAG coverages to obtain MAG absolute concentrations and save scaling plots in the working directory
#For rlm scaling using scale_features_rlm
#For rlm scaling using scale_features_lm
mag_tab_scaled <- SIPmg::scale_features_rlm(f_tibble, sequins, seq_dil, log_scale, coe_of_variation = coe_of_variation, save_plots = FALSE)
mag_tab = as.matrix(mag_tab_scaled$mag_tab) #Extract absolute abundances as a matrix
mag.table = phyloseq::otu_table(mag_tab, taxa_are_rows = TRUE) #Phyloseq OTU table
taxonomy.object = SIPmg::tax.table(taxonomy_tibble) # Create a taxonomy phyloseq object
samples.object = SIPmg::sample.table(fractions) # Create a samples phyloseq object
phylo.qSIP = SIPmg::phylo.table(mag.table, taxonomy.object, samples.object) # Make a phyloseq table for downstream qSIP analysis
atomX = SIPmg::qSIP_atom_excess_MAGs(phylo.qSIP,
control_expr='Isotope=="12C"',
treatment_rep='Replicate',
Gi = GC_content)
#Bootstrap confidence intervals
df_atomX_boot = SIPmg::qSIP_bootstrap_fcr(atomX, n_boot=10) #Change "parallel = FALSE" to compute using a single-core
CI_threshold = 0
df_atomX_boot = df_atomX_boot %>%
dplyr::mutate(Incorporator = A_CI_fcr_low > CI_threshold,
OTU = reorder(OTU, -A))
(atom_f_excess_plot = ggplot2::ggplot(df_atomX_boot, aes(OTU, A, ymin=A_CI_low, ymax=A_CI_high, color=Incorporator)) +
geom_pointrange(size=0.25) +
geom_linerange() +
geom_hline(yintercept=0, linetype='dashed', alpha=0.5) +
labs(x='MAGs', y='Atom fraction excess') +
theme_bw() +
coord_flip() +
ggtitle("Isotope incorporating MAGs"))
#Get incorporator info
n_incorp = df_atomX_boot %>%
dplyr::filter(Incorporator == TRUE) %>%
nrow
#Get incorporator list
incorporator_list = SIPmg::incorporators_taxonomy(taxonomy = taxonomy_tibble, bootstrapped_AFE_table = df_atomX_boot)
#Print incorporator information
cat('Number of incorporators:', n_incorp, '\n')
print(incorporator_list, n = nrow(incorporator_list))
#Load function for abundance plots
data_vis_table = SIPmg::sip_vis(fractions = fractions, mag_tab = mag_tab, afe_tab = df_atomX_boot, taxonomy = taxonomy_tibble)
library(SIPmg)
#Load required libraries
list.of.packages <- c("tidyverse", "HTSSIP", "ggpubr","data.table")
new.packages <- list.of.packages[!(list.of.packages %in% installed.packages()[,"Package"])]
if(length(new.packages)) install.packages(new.packages, quiet = T, dependencies = T, repos = "http://cran.us.r-project.org")
if (!require("BiocManager", quietly = TRUE))
install.packages("BiocManager")
BiocManager::install("phyloseq")
BiocManager::install("EBImage")
library(tidyverse)
library(phyloseq)
library(HTSSIP)
library(ggpubr)
library(SIPmg)
set.seed(seed = 1000)
## Load data
#Coverage metadata
#Uncomment if your coverage data is in the format mentioned above for this file. Remains commented if you are using the output from `checkm coverage`
f_tibble <- readr::read_csv("mock_input_data/coverage_metadata.csv")
#Sequins metadata
sequins <- readr::read_csv(file="mock_input_data/sequins_metadata.csv")
#Dilutions data
seq_dil = readr::read_csv(file = "mock_input_data/dilutions_data.csv")
#Log scale BOOLEAN. True or False depending on how you would want the MAG coverages to be scaled. Select TRUE if you need MAG concentrations scaled on the log scale
log_scale = TRUE
#coe_of_variation. Acceptable coefficient of variation for coverage and detection (eg. 20 - for 20 % threshold of coefficient of variation) (Coverages above the threshold value will be flagged in the plots)
coe_of_variation = 20
#Taxonomy
gtdbtk_bac_summary = readr::read_delim("mock_input_data/gtdbtk.bac120.summary.tsv",
"\t", escape_double = FALSE, trim_ws = TRUE)
gtdbtk_archaea = readr::read_delim("mock_input_data/gtdbtk.ar122.summary.tsv",
"\t", escape_double = FALSE, trim_ws = TRUE)
#GC content
GC_content <- readr::read_csv(file = "mock_input_data/GC_content.csv")
#Fractions
fractions = readr::read_csv("mock_input_data/fractions.csv")
taxonomy_tibble = dplyr::bind_rows(gtdbtk_bac_summary, gtdbtk_archaea) #Combine bacteria and archaea taxonomy files if it has not been done yet
#mag_tab is a tibble with absolute concentrations of MAGs obtained by scaling MAG coverages using linear regression models on sequin coverages and concentration
##Scale MAG coverages to obtain MAG absolute concentrations and save scaling plots in the working directory
#For rlm scaling using scale_features_rlm
#For rlm scaling using scale_features_lm
mag_tab_scaled <- SIPmg::scale_features_rlm(f_tibble, sequins, seq_dil, log_scale, coe_of_variation = coe_of_variation, save_plots = FALSE)
mag_tab = as.matrix(mag_tab_scaled$mag_tab) #Extract absolute abundances as a matrix
mag.table = phyloseq::otu_table(mag_tab, taxa_are_rows = TRUE) #Phyloseq OTU table
taxonomy.object = SIPmg::tax.table(taxonomy_tibble) # Create a taxonomy phyloseq object
samples.object = SIPmg::sample.table(fractions) # Create a samples phyloseq object
phylo.qSIP = SIPmg::phylo.table(mag.table, taxonomy.object, samples.object) # Make a phyloseq table for downstream qSIP analysis
atomX = SIPmg::qSIP_atom_excess_MAGs(phylo.qSIP,
control_expr='Isotope=="12C"',
treatment_rep='Replicate',
Gi = GC_content)
#Bootstrap confidence intervals
df_atomX_boot = SIPmg::qSIP_bootstrap_fcr(atomX, n_boot=10) #Change "parallel = FALSE" to compute using a single-core
library(SIPmg)
library(SIPmg)
#Bootstrap confidence intervals
df_atomX_boot = SIPmg::qSIP_bootstrap_fcr(atomX, n_boot=10) #Change "parallel = FALSE" to compute using a single-core
library(SIPmg)
library(SIPmg)
library(SIPmg)
library(SIPmg)
library(SIPmg)