---
title: "Quality Control Report"
output: 
    html_document:
        theme: cosmo
        number_sections: true
        toc: true
        toc_float:
            collapsed: true
toc-title: "rRMSAnalyzer"
    
params:
    library_col: "library"
    condition_col: "condition"
    ribo_name: "ribo"
    project_name: "Unnamed project"
---
```{r include = FALSE}
knitr::opts_chunk$set(echo=FALSE)

all_sites <- length(ribo_to_check[[1]][[1]][[1]])

# Get the number of valid sites
ribom <- extract_data(ribo_to_check)
valid_sites <- sum(complete.cases(ribom))

```

```{=html}
<style type="text/css">
  body{
 font-family: Helvetica;    
}

details{
background-color: #E4F3F3;
}

</style>
```
<details>

<summary>

<b> ℹ️ Summary of input data </b>

</summary>
Project name : **`r params$project_name`**<br><br>
User name : **`r Sys.getenv("USERNAME")`** <br>
Date : **`r Sys.time()`**

RiboClass Name : **`r params$ribo_name`**. <br> Number of uploaded samples: **`r length(ribo_to_check[[1]])`** samples. <br> RNA(s) used by the coverage tool: **`r as.character(ribo_to_check[["rna_names"]][["original_name"]])`** <br> There are **`r as.character(all_sites)`** genomic positions.


__**C-score**__ <br>
Method of computation : **`r as.character(ribo_to_check["cscore_method"])`** <br>
Window length : **`r as.character(ribo_to_check["cscore_window"])`**<br>
**`r as.character(valid_sites)`** genomic positions with a valid
C-score in all samples.
<br>
<br>
__**ComBat-seq**__<br>
ComBat-seq correction:
**`r ifelse(is.null(ribo_to_check[["combatSeq_count"]]),"No","Yes")`** <br>
Name of column in metadata used for ComBat-seq adjustment:
**`r ifelse(is.null(ribo_to_check[["col_used_combatSeq"]]),"None",as.character(ribo_to_check["col_used_combatSeq"]))`**

</details>

# Introduction

## Content

This quality control (QC) reports different metrics to ensure that subsequent analyses are free from bias related to the riboMethSeq technology. 

In this report, the end read count data is used to ensure the quality of the sequencing coverage on the samples. Then, this measure is used to ensure the technical homogeneity of the cohort by comparing samples with each other. QC verifies that the coverage is uniform and reproducible between samples, eliminating the possibility of biases due to sequencing and outlier samples. Finally, a cohort check using the C-score values is performed to verify that subsequent analyses will be free of biases that may affect biological conclusions.

## Project

This QC report is generated from the analysis of the project **`r params$project_name`**, which is loaded inside the RiboClass object **`r params$ribo_name`**. This **`r params$ribo_name`** object was generated from **`r length(ribo_to_check[[1]])`** samples and **`r as.character(all_sites)`** genomic positions.

The following (r)RNA were used by the coverage tool: **`r as.character(ribo_to_check[["rna_names"]][["original_name"]])`**

## Method

**C-score computation** : In the following analyses, the C-Score was computed using the **`r as.character(ribo_to_check["cscore_method"])`** of the **`r as.character(ribo_to_check["cscore_window"])`** neighboring positions (Birkedal et al, 2015; Marcel et al, 2020). In the end, there are **`r as.character(valid_sites)`** genomic positions with a valid C-score in all samples. Of note, C-score for genomic positions on RNA’s extremities cannot be computed because of the lack of neighbors in the local window.

**Sequencing biases correction** : The possible technical biases induced by the sequencing step are corrected using the **ComBat-seq** method (Zhang et al, 2020; Paraqindes et al, 2023). In this analysis, the correction is **`r ifelse(is.null(ribo_to_check[["combatSeq_count"]]),"unused","used")`**. Here we used the variable "**`r ifelse(is.null(ribo_to_check[["col_used_combatSeq"]]),"None",as.character(ribo_to_check["col_used_combatSeq"]))`**" to perform this correction.

## Metadata

The following metadata are used in the project **`r params$project_name`** :

```{r, echo=FALSE}
knitr::kable(ribo_to_check[[2]])
```
 
# Global coverage

In this section, we present QC measures that are derived from the end read count data on all sequenced positions (i.e., global coverage data) to ensure consistent and reproducible sample coverage. A good level of coverage is essential, because the level of methylation is estimated from its variations.

##  Coverage distribution

The coverage distribution of absolute end read counts at each genomic positions is visuallized here using a boxplot for each individual samples, controlling for coverage variation. 

Samples with a median coverage below 100 (the blue line), are highlighted in red and should be interpreted with caution or excluded from further analysis. This threshold ensures a level of precision within the percentage range for calculating the C-score, which is targeted for good quality analyses.

```{r boxplot2,warning=FALSE}
rRMSAnalyzer::boxplot_count(ribo_to_check) + ggplot2::ggtitle("End read counts distribution by sample")
```

## Relative Log Coverage

The **Relative Log Coverage** (RLC)  is a measure to ensure  the reproducibility of the coverage across all samples.

$$\log_2(\frac{pos_{sample} + 1}{pos_{all}})$$


Where :

-   pos~sample~ : the end read count at a specific position for a specific sample.
-   pos~all~ : the end read count **median of all samples** for the same specific position.

The following image shows the **Relative Log Coverage boxplot**. The blue lines give the median +/- 2mad variation. Median outside the blue lines identifies putative outlier sample (highlighted in red). Sample that presents a lower or higher coverage compared to other samples might be considered carefully or removed in the following biological analyses as an outlier sample.

```{r}
rRMSAnalyzer::plot_rle(ribo_to_check,show_outlier = T) + ggplot2::ggtitle("Relative Log Coverage")
```

## RNA fraction per sample

The repartition of RNA end read counts is an indicator of the RNA(s) that present a lower or higher coverage compared to the other RNAs within the same sample. A sample that shows an unusual repartition might be considered with caution, as it could impact the C-score robustness. The following image shows a stacked barplots representing end read count fractions of total end read counts for each RNA per sample. In the case of rRNAs, most end read counts are assigned to the longest rRNA, namely the 18S (~ 0.25) and the 28S (~ 0.50).

```{r}
rRMSAnalyzer::plot_counts_fraction(ribo_to_check) + ggplot2::ggtitle("RNA fraction per sample")
```

# Profiles comparison

In this section, we explore the difference between the samples using global coverage and C-scores. By using such a global information, technical variabilities are more visible than biological variabilities, allowing us to highlight possible biases due to the experiment. 

A sample with an atypical coverage pattern should be treated with caution or excluded from further analysis.

It is possible for samples to cluster in sequencing batches. In this case, ComBat-seq correcting algorithm should be used to avoid biaise in subsequent analyses. 

Of note, rRNA sequences are enriched in GC-stretches and in nucleotide repeats which introduces a bias during read alignment. The coverage profile is therefore not uniform, however, it will be typical of the rRNA sequence and reproducible across samples.

## Correlation matrix

In order to identify distinct coverage profiles, coverage profiles are compared with one other within the cohort. Pairwise sample correlations of coverage profile are calculated using a Pearson correlation based distance. A heatmap is then generated that summarizes the correlation values for the entire series. These values goes from 0 (no correlation, indicated in red) to 1 (perfect correlation, indicated in blue). 

```{r heatmap_corr}
rRMSAnalyzer::plot_heatmap_corr(ribo_to_check,"Count")
```

## Correspondance analysis

To identify distinct coverage profile, we compare the coverage profiles of all samples in the series. The Correspondence Analysis (COA) plot shows the distance between each sample based on its coverage profile. Samples that are far appart are considered putative outliers. At this stage, it is possible to visualise the bias linked to the sequencing batch or to evaluate the performance of the ComBat-seq correction algorithm if it has been applied (Zhang et al, 2020; Paraqindes et al, in revision).

```{r}
rRMSAnalyzer::plot_coa(ribo_to_check,color_col = params$library_col)
```

# C-scores analysis

## Principal component analysis of C-scores

In this section, Qcs are made from C-score measurements, which will be used in the subsequent analysis. This final QC will therefore check that the dataset is free of biases that could affect the biological conclusions.

Since more than 1.5% of the C-score corresponds to non-biological noise, normalization during C-score calculation should limit the dispersion of the samples based on their C-score at all genomic positions in the whole serie. Thus, Principal Component Analysis (PCA) on C-scores helps to identify putative batch effect that can be corrected using the inter-normalization option (ComBat-seq tool) (Zhang et al, 2020; Paraqindes et al, 2023). 

The image shows **a Principal Component Analysis (PCA) plot illustrating the distance between each sample based on its C-score at all genomic positions**. The samples are colored by library. 

The PCA cannot use genomic positions with an invalid C-score.
Thus, it uses **`r as.character(valid_sites)`** genomic positions instead of all **`r as.character(all_sites)`**.

```{r not corrected cscore PCA,echo=FALSE}
rRMSAnalyzer::plot_pca(ribo_to_check,params$library_col)
```

# About

## Bibliography

-   Yuqing Zhang, Giovanni Parmigiani, W Evan Johnson, ComBat-seq: batch effect adjustment for RNA-seq count data, NAR Genomics and Bioinformatics, Volume 2, Issue 3, 1 September 2020, lqaa078, [https://doi.org/10.1093/nargab/lqaa078](https://doi.org/10.1093/nargab/lqaa078)

-   <mark> Marcel et al 2020, </mark>

-   Birkedal, U., Christensen-Dalsgaard, M., Krogh, N., Sabarinathan, R., Gorodkin, J. and Nielsen, H. (2015), Profiling of Ribose Methylations in RNA by High-Throughput Sequencing. Angew. Chem. Int. Ed., 54: 451-455. https://doi.org/10.1002/anie.201408362

## Citation
**When using this report, please cite :**


## License

rRMSAnalyzer<br>
Copyright (C) 2023 Centre de Recherche en Cancérologie de Lyon

This program is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.

This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
GNU General Public License for more details.

You should have received a copy of the GNU General Public License
along with this program.  If not, see <http://www.gnu.org/licenses/>.

## SessionInfo

```{r,echo=TRUE}
sessionInfo()
```


-----------------------
This report has been automatically generated with [rRMSAnalyzer](https://github.com/RibosomeCRCL/rRMSAnalyzer) 2.0.0 (April 2023)