02-Quick-example.Rmd

# Using SIPmg R package {#quick-example}

This R package was designed to facilitate statistical analyses using SIP metagenomic data. Particularly, coverage and taxonomic classification of metagenome-assembled genomes and metadata of fractionated DNA is used to identify incorporators and quantify isotopic enrichment. 


Overall, the package can perform the following:

1. Coverage normalization to absolute concentrations or relative coverage. Scaling of coverages to absolute abundances based on post-fractionation reference standards, called sequins.

2. Statistical analysis based on either of 
  + AFE mathematical model by [Hungate et al. (2015)](https://journals.asm.org/doi/10.1128/AEM.02280-15) or
  + Delta buoyant density $AFE_{ΔBD}$ = $\frac{W_{Lab}- W_{Light}}{I_{max}}$.        
  Where Imax is the maximum linear shift in DNA BD (upon 100% labeling), as        discussed by [Birnie and Rickwood (1978)](https://www.ncbi.nlm.nih.gov/pmc/articles/PMC2595786/)
  + Analyses as offered by [HTS-SIP](https://cran.r-project.org/web/packages/HTSSIP/index.html)
  
3. Data visualization for identification of isotope incorporators.

## Example workflow

qSIP analysis is typically performed using sequins, MAG coverage, and MAG taxonomy is performed using this R markdown. Briefly, coverages are normalized, MAG coverages are then scaled based on linear regression models from sequin coverage and concentrations, a phyloseq object is created from MAG absolute concentrations and taxonomy (GTDB taxonomy output format is required) data. \n
Note: Sequins that were spiked in the DNA-SIP fractions will be used in scaling and creating linear regression models for evaluating absolute MAG concentrations. Please see `sequin_scaling.R` \n
For normalizing coverage values please see `pooling_funs.R` \n

This R markdown primarily uses `Tidyverse`, `phyloseq`, and `HTSSIP` packages \n
This markdown uses functions from `functions_qSIP_MAGs.R` \n

For this data, MetaSIPSim (developed by [Barnett et al, 2020](https://bmcbioinformatics.biomedcentral.com/articles/10.1186/s12859-020-3372-6)) was used to simulate a data set from a hypothetical SIP workflow performed on the MBARC-26 mock community (as defined in [Singer et al, 2016](https://www.nature.com/articles/sdata201681)). The data set was further augmented with simulated post-fractionation sequin spike-ins and the data set was simulated to have all 26 organisms as incorporators. Some of the 26 organisms also had plasmids, which were assembled as distinct genomic units, adding the total genomic features in the simulated study to be 38

```{r Load libraries, warning = FALSE, message = FALSE}
#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")

library(tidyverse)
library(phyloseq)
library(HTSSIP)
library(ggpubr)
library(SIPmg)


set.seed(seed = 1000)
```

## Load required data 

The following files are required: \n

### MAG coverage data

Pooled coverages data as a comma separated file (.csv file) across samples for `Features` including sequins that were used as spike-ins. The followings columns are required: \n   

  1. Feature: A character string describing the `Feature` label
  
  2. Sample: The label for these n number of columns should be in the format of "'isotope'\_rep\_#\_fraction\_#". For instance, "12C\_rep\_1\_fraction\_1". The number of sample columns should be the same as the product of (replicates, fractions, isotopes)
    
### Sequin metadata

Load the sequins metadata as a comma separated file (.csv file) which has the following columns: \n

  1. Feature: As described above
    
  2. Length: Length of the sequin in bp
    
  3. GC_content: GC content of the sequin
  
  4. Sequence: Sequin nucleotide sequence
  
  5. Concentration: Concentration of the sequin in attamoles/uL
    
### Dilutions data

Load dilutions data as a comma separated file (.csv file) that contains the following columns \n

  1. Sample: Similar to the sample name as described above.
  
  2. Dilution: Dilution of sequins added to the fraction before sequencing.

### Fractions metadata

A fractions file as a comma separated file (.csv file) with the following columns: \n

  1. Replicate: Depends on how many replicates the study has \n
  
  2. Fractions: Typically in the range of 1-24 \n
  
  3. Buoyant_density: As calculated from the refractometer for each fraction and replicate \n
  
  4. Isotope - "12C", "13C", "14N", "15N" etc. \n
  
  5. DNA_concentration \n
  
  6. Sample - In the format "'isotope'\_rep\_#\_fraction\_#". For instance 12C\_rep\_1\_fraction\_1 \n

### GTDB style taxonomy data

A taxonomy file in the GTDB output format (.tsv format). Load the bacteria and archaea taxonomy outputs separately. The markdown requires loading the standard output files from GTDB-Tk separately for bacteria and archaea

### MAG absolute concentrations

MAG absolute concentrations obtained from scaling_sequins.R. \n
`mag_tab` object obtained from the above script is to be used as the input here

### GC content

GC content of the MAGs as a comma separated file (.csv file). The table should contain the following columns: \n

  1. OTU: MAG identifier such as the `Feature` label from the `sequin_scaling.R` script \n
  
  2. GC_content: GC content of the `Feature` in the range of 0-1 \n
  
### 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 \n

### Coefficient 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. Here a value of 20 is used. \n


```{r Load data, warning = FALSE, message = FALSE}
## 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")

```

## Coverage normalization and scaling

Coverage data of MAGs (or features of interest) and sequins is scaled using sequin concentration to obtain absolute concentration of MAGs (or features of interest). For this step, either linear regression or robust linear regression can be used. 

If linear regression is used, the function `scale_features_lm` is used. The functions in `sequin_scaling_lm.R` are used for this step.

If robust linear regression is used, the outliers impact on the regression models is minimized. For robust linear regression, `scale_features_rlm` is used. The functions in `sequin_scaling_rlm.R` are used for this step. In this example workflow, robust linear regression is used to scale coverage data in log10 scale. 

For more discussion on the choice of regression method and comparison of the two methods please see this [vignette](#scaling-coverage-data)

### This function uses the following datasets:

1. Coverage data (`f_tibble`)
2. Sequins metadata (`sequins`)
3. Dilution of sequins used to add into the fractions (`seq_dil`)

The regression scaling plots are saved in the project folder within a subfolder `sequin_scaling_plots`

```{r estimate absolute concentrations, warning = FALSE, message = FALSE}
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
```

#### Example of a scaling plot

An example of a scaling plot is as below

```{r example scaling plot, echo = FALSE, warning = FALSE, message = FALSE}
mag_tab_scaled$plots$plots[[1]]
```

## Preparation of scaled data for qSIP analysis

The scaled abundance data, taxonomy data, and metadata is then converted to phyloseq objects and a master phyloseq object is created for qSIP analysis.

```{r Make phyloseq objects, warning = FALSE, message = FALSE}

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
```

## Estimation of atom fraction excess

The following steps estimate the atom fraction excess (AFE) in the MAGs. This calculation is based on the mathematical model suggested by [Hungate et al. (2015)](https://journals.asm.org/doi/full/10.1128/AEM.02280-15). In the article, the authors suggest that GC contents of the individual biological features of interest should be accounted for to have a better estimate of AFE. Through the power of metagenomic analysis and the recovery of MAGs, one can estimate the GC content of features of interest and thereby obtain a better estimate of AFE. Thus, here the user would use the GC contents of MAGs in the AFE estimations and obtain isotopic enrichment of individual MAGs. 

An important consideration is that the AFE calculations are based on a mathematical model and are not an absolute estimate of isotopic enrichment. One of the goals of AFE estimation is to determine isotope incorporators. 

In this example workflow, the approach reported in [Hungate et al. (2015)](https://journals.asm.org/doi/full/10.1128/AEM.02280-15) is used to determine incorporators.

For this step, the following are used:

1. Master phyloseq object from the previous step
2. GC content

Thanks to [HTSSIP R package](https://journals.plos.org/plosone/article?id=10.1371/journal.pone.0189616), bootstrapping can be performed faster with multithreading.

```{r Calculate atom fraction excess, warning = FALSE, message = FALSE}
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

```


## Plot the atom fraction excess plot

Plot the atom fraction excess data with 5 % and 95 % quantiles as the confidence limits. In the same way as reported in [Hungate et al. (2015)](https://journals.asm.org/doi/full/10.1128/AEM.02280-15), a MAG is considered an incorporator if the lower confidence interval of its AFE is above zero.

```{r Plot atom fraction excess, warning = FALSE, message = FALSE}
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"))
```

## Check incorporator list

```{r incorporators, warning = FALSE, message = FALSE}

#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))
```

It appears that not all incorporators were identified with the bootstrapping method. This could be due to the model unable to detect the incorporators. Another way to detect incorporators, as mentioned above, is to test with the ΔBD method or with the HR-SIP or MW-HR-SIP methods. For a more detailed discussion, please refer this [vignette](#identifying-incorporators)

```{r, warning = FALSE, message = FALSE}
sessionInfo()
```