03-Scaling-coverage-data.Rmd

# Scaling coverage data {#scaling-coverage-data}

In this example we will look at the two scaling approaches that can be utilized in the SIPmg package. The user can decide either 

1.	A robust linear regression model
2.	An ordinary least squares linear regression model
We will also briefly discuss data filtering methods when there are heavily influencing outliers, as a result of upstream methods.

## Basic approach for scaling

Sequencing effort on each fraction and subsequently obtained coverage is typically used to estimate relative abundance for a feature of interest. However, as we discussed in the introductory vignette, this data is compositional and a quantitative estimate of concentrations can be highly informative. For this purpose, synthetic spike-in standards (sequins) with known concentrations are added before recovery of fractionated DNA, from the CsCl medium into TE buffer or water. For details on sequins please refer to [Hardwick et al. (2018)](https://www.nature.com/articles/s41467-018-05555-0). These reference standards in each fraction, can now be used to estimate the concentration of the feature of interest in each fraction.

The approach used for scaling is detailed in the following steps:

1.	For each fraction, sequin coverages and known concentrations are obtained to make a standard curve.

2.	Limit of detection of sequins, i.e., the lowest concentration where at least one sequin has detectable coverage, is determined for each group of sequin concentrations.

3.	For each group of sequin concentrations, mean, standard deviation, and the coefficient of variation of coverage are determined. Groups of sequins with coefficient of variation greater than the set threshold will be flagged.

4.	Sequin groups with coverage values above the limit of detection and below the coefficient of variation threshold are filtered in preparation for regression.

5.	Linear regression or robust linear regression, based on the user choice, is performed on sequin coverage values as the independent variable and concentration of the sequin group as the dependent variable. The user can also decide on performing log scaling of coverage and concentration values for regression.

6.	The regression parameters are extracted and plots with regression parameters and best fit line are saved for inspection.

7.	Coverage values above the limit of detection and below the coefficient of variation threshold are filtered for further analysis. The rest of the values are flagged and reported.

8.	The filtered values are scaled based on regression parameters to estimate absolute concentrations in each library.

9.	Absolute concentrations of biological features are saved as a fresh dataset for the subsequent isotope incorporator identification and AFE estimation pipeline.

## Choice of regression model

Ordinary least squares regression can be very sensitive to outliers, resulting a poor model performance when the data is contaminated with outliers. Although there are methods, like Cook's distance based filtering to address outliers, it may not always be the best idea to remove these data points.

See [these](https://statisticsbyjim.com/basics/remove-outliers/) [discussions](https://stats.stackexchange.com/questions/175/how-should-outliers-be-dealt-with-in-linear-regression-analysis) and more for better insights into handling outliers. 

Robust linear regression is an alternative to this situation. Robust linear regression assigns appropriate weights to the data points, minimizing the influence of outliers on the model performance, without deleting data. In this pipeline, a Huber loss function from the [`MASS` R package](https://cran.r-project.org/web/packages/MASS/MASS.pdf) was used for robust linear regression.

The user is free to choose between the two methods for regression. 

## Scaling the data

For illustrating the differences in the regression models, a different dataset, that was impacted by one or many upstream methods, will be used than the one in the introductory vignette. 

This dataset can be found here

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

```

For both scaling functions the following data and parameters are required

1.	Coverage data. The output of [checkm coverage](https://github.com/Ecogenomics/CheckM/wiki/Utility-Commands#coverage) command can be directly used.

2.	Sequin dilutions and metadata.

3.	Whether or not log10 scaling to be performed (a BOOLEAN value of TRUE or FALSE).

4.	Coefficient of variation. In this dataset, a higher coefficient of variation compared to the value used in the example workflow of the [introduction vignette](#quick-example) is used to allow data scaling to occur.

If the robust linear regression is chosen, the function is `scale_features_rlm` which is performed the function file `sequin_scaling_rlm.R`.

```{r robust linear regression, warning = FALSE, message = FALSE}
mag_tab_scaled_rlm <- SIPmg::scale_features_rlm(f_tibble, sequins, seq_dil, log_scale, coe_of_variation = coe_of_variation, save_plots = FALSE)
```


If a linear regression is chosen, the function is `scale_features_lm` which is performed by the function file `sequin_scaling_lm.R`. 

If the user chooses linear regression, they can choose to perform outlier filtering based on the traditional Cook’s distance threshold of 4/n (n being the sample size). Only the outliers in the sequin coverage data are flagged as outliers with the Cook’s distance threshold and are filtered out. Later, the remaining data is subjected to linear regression to obtain model parameters and to scale coverage values of feature of interest.


```{r linear regression, warning = FALSE, message = FALSE}
mag_tab_scaled_lm <- SIPmg::scale_features_lm(f_tibble, sequins, seq_dil, log_scale, coe_of_variation = coe_of_variation, cook_filtering = TRUE, save_plots = FALSE)
```

## Comparison of the fits

```{r load images, echo = F, warning = FALSE, message = FALSE}
rlm_example = knitr::include_graphics("rlm-example.png")

lm_example = knitr::include_graphics("lm-example.png")
```

Robust linear regression provides a better model fit compared to the linear regression. Additionally, robust linear regression has less variability compared to the linear regression model.

As discussed previously, removing the outliers may not always be the best idea. However, the outliers can be visualized to assess how "far out" are the outliers and how influential the data points are, negatively influencing the linear regression method.

The `mag_tab_scaled_lm` function provides the plots for the visualization of Cook's distance for the sequin data in each fraction. In the plot below, an outlier is visualized. For this data, the Cook's distance threshold was 0.09, and the outlier had a Cook's distance of 10.8. The highly influential data point could well be the reason the linear regression model had a poor fit.

```{r, outlier sequin coverages, warning = FALSE, message = FALSE}

#Cook's distance threshold of the data set
4/(length(mag_tab_scaled_lm$scale_fac$cooksd[[3]]))

#Cook's distance of the outlier before filtering the data
max(mag_tab_scaled_lm$scale_fac$cooksd[[3]])

#Before filtration
cooksd_example = knitr::include_graphics("cooksd-example.png")
```

Upon filtering the data to remove the outlier, the fit becomes much better for the linear regression model

```{r better fit, warning = FALSE, message = FALSE}
#Linear regression plot with filtered outliers in sequin data
filtered_lm_example = knitr::include_graphics("filtered_lm-example.png")
```

## What if there were no sequins added in the study
In our study, we realized that the sequin-scaled data provided the best correlation (Spearman correlation coefficient = 0.85), compared to methods without the use of sequins. However, not all studies would have access to sequins or could choose not to add sequins. We tested qSIP methods with absolute concentration of MAGs obtained as a product of fraction DNA concentration and either with MAG relative abundance (as performed by [Greenlon et al.](https://journals.asm.org/doi/10.1128/mSystems.00151-20)) or relative coverage (as performed by [Starr et al.](https://journals.asm.org/doi/10.1128/mSphere.00085-21)), or simply with relative coverage, i.e., without converting the MAG coverages to absolute concentration. We found that without sequins, simply using relative coverage provided better specificity (0.99) and good correlation (Spearman correlation coefficient = 0.76), compared to the other two methods. The user could choose any of the three approaches to obtain normalized coverages using the equation `coverage_normalization()` where the parameter approach would be chosen accordingly (either "greenlon", "starr", or "relative_coverage"). In this example, simple relative coverage will be estimated without converting to absolute concentration based on fraction DNA concentrations. The relative coverage table can then be converted to a phyloseq object to run incorporator identification and/or AFE estimation. For more details, please refer to the [incorporator identification section][SIP metagenomics and identification of isotope incorporators]

```{r, estimating relative coverage, warning = FALSE, message = FALSE}
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 %>% tibble::column_to_rownames(var = "Feature")), taxa_are_rows = TRUE) #Phyloseq OTU table

```



## Decision on the method to scale coverages

It can be seen from the above plots that the regression parameters differ based on the method of regression and whether or not data filtering is used. These regression parameters directly influence the estimation of concentrations of features of interest in the microbiome. Thus, the method for regression must be a careful choice. 

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