RESTORBIOME2: Exploration and normalisation

Experiment with BA

Authors
Affiliation

Olivier Rué

Migale bioinformatics facility

Christelle Hennequet-Antier

Migale bioinformatics facility

Published

May 2, 2023

Modified

February 19, 2025

Show the code 
#Use packages and functions
library(tidyverse)
library(data.table)
library(SummarizedExperiment) # manipulation RNASeq assays
library(rtracklayer) ## for annotation
#library("DEFormats")
library(edgeR)
library(kableExtra)
library(dplyr)
library(purrr)
library(stringr)
library(ggplot2)
library(reshape2)
library(mixOmics) #PCA

Read counts

The matrix of raw counts looks like that.

Show the code 
#cat("EXP = ", params$my.interest)

### Get Raw Count data in edgeR object
# At gene level
#inputFile <-  "count_allSamples_BA.csv"
inputFile <- list.files("../data/counts") |> str_subset(params$my.interest)

#cat("File raw = ",inputFile,"\n")
rawdata <- fread(file.path("../data/counts",inputFile))
# rename first column gene ID to be homogenous between assays
colnames(rawdata)[1] <- "Gene_ID"
nameEXP <- params$my.interest

# genecount: count matrix
geneCount <- rawdata[,-c(1)] |> as.matrix()
row.names(geneCount) <- rawdata %>% dplyr::select(1) %>% unlist() %>% as.character()
#row.names(geneCount) %>% head()
#print(dim(geneCount))
head(geneCount) %>% kbl() %>% kable_styling() %>% 
  scroll_box(width = "100%", height = "200px")
BA_G_R1 BA_G_R2 BA_G_R3 BA_R_R1 BA_R_R2 BA_R_R3
gene-BAD_RS00005 6201 4936 4267 3441 2897 3652
gene-BAD_RS00010 6815 5705 4232 8147 5974 7703
gene-BAD_RS00015 1956 1488 1078 1824 1403 1736
gene-BAD_RS00020 533 441 273 635 401 525
gene-BAD_RS00025 8473 5992 4906 3805 2754 3120
gene-BAD_RS00030 17965 12906 11256 5707 3939 4261

This matrix contains raw counts from BA experiment with 1729 genes and 6 samples.

This is the experimental design.

Show the code 
# create Design from sample's names in file
sampleInfo.all <- colnames(rawdata)[-c(1)] %>%
  stringr::str_split(., "[_]", simplify = FALSE) %>%
  transpose() %>%
  purrr::simplify_all() 
names(sampleInfo.all) <- c("Bacteria","Treatment","Replicate")
#sampleInfo.all$sample <- paste(sampleInfo.all$Treatment, sampleInfo.all$Rep, sep="_")
sampleInfo.all$sample <- colnames(rawdata)[-c(1)]

Design <- as.data.frame(sampleInfo.all)
Design <- Design %>% mutate( group = paste(Bacteria, Treatment, sep="_"))
Design %>% kbl() %>% kable_styling() %>% 
  scroll_box(width = "100%", height = "200px")
Bacteria Treatment Replicate sample group
BA G R1 BA_G_R1 BA_G
BA G R2 BA_G_R2 BA_G
BA G R3 BA_G_R3 BA_G
BA R R1 BA_R_R1 BA_R
BA R R2 BA_R_R2 BA_R
BA R R3 BA_R_R3 BA_R
Show the code 
# Design_file <- paste("../output/Design", nameEXP, ".csv", sep="")
# write.table(Design, Design_file, row.names=FALSE, sep=",")

We used gene annotation produced using DBCAN .

Show the code 
# GFF annotation file
filegff <- case_when(
  nameEXP == 'BA' ~ "../data/GFF/NC_008618_bifido_adolescentis_spikes_cazymes_dbcan.gff3",
  nameEXP == 'BU' ~ "../data/GFF/CP102263_1_Bacteroides_uniformis_spikes_cazymes_dbcan_pulpred_cc.gff3",
  nameEXP == 'BC' ~ "../data/GFF/NZ_AP012325_bifido_catenulatum_spikes_cazymes_dbcan.gff3",
  nameEXP == 'ER' ~ "../data/GFF/NC_012781_eubacterium_rectale_spikes_cazymes_dbcan.gff3",
)

# use gff produced by DBCAN
gff <- readGFF(filegff)


df_annot <- gff %>% as_tibble() %>% unnest_longer(Parent)  


# remove comulmn with all NA
df_annot <- df_annot %>% select_if(~sum(!is.na(.)) > 0) 
# Remark: duplicated gene annotation
# which(duplicated(df_annot[,"Parent"]))
# keeep the first line when multiple annotation
df_annot <- df_annot%>%
  group_by(seqid, source, type, ID, Parent) %>%
  filter(row_number() <= 1)
#glimpse(df_annot)


rawdata_annot <- left_join(rawdata[, 1], df_annot, by=c("Gene_ID" = "Parent"), keep=FALSE) %>%
  as.data.frame()

#glimpse(rawdata_annot)

In this BA experiment, the effect of each treatment (G, R) on gene expression will be explore. The number of biological repetitions is 3, 3 for G, R treatment, respectively.

Show the code 
# create a SummarizedExperiment object
se <- SummarizedExperiment(assays=list(counts=geneCount),
                     rowData=rawdata_annot, colData=Design)
#print(se)

#counts matrix = assay(se)
#size of library
#libsize <- assay(se) %>% colSums()

#exemple to remove sample BT_G_Rep1
#se <- se[, se$sample !="BT_G_Rep1"]
#colData(se)
#dim(se)
#print(se)

#counts matrix = assay(se)
libsize <- assay(se) %>% colSums()
libsize %>% kbl(col.names = c("libsize")) %>% kable_styling(full_width = FALSE)
libsize
BA_G_R1 18392341
BA_G_R2 13758350
BA_G_R3 12109038
BA_R_R1 13997179
BA_R_R2 11289032
BA_R_R3 13912504

Genes with very low counts across all libraries may be filtered.

Show the code 
# Filter gene with low count, at least number of biological replicates with cpm > valfilt 
#summary(libsize)
#valfilt <- round(10*1e06 / min(libsize), 1)
valfilt <- 1
nbrepbio <- min(table(colData(se)$group))
keep <- rowSums(edgeR::cpm(se) > valfilt) >= nbrepbio
#summary(keep)
se <- se[keep,]
#dim(se)

The filtering on genes is based on count per million (cpm) greater than 1 in at least 3 samples corresponding to the minimum number of biological replicates. We kept 1713 expressed genes for further analyses from 6 samples.

Normalisation

A normalization factor is calculated to take into account the different sizes of the sequencing banks (i.e. the total read count) and the distribution of reads per sample on sequencing run, as discussed [1]. Normalization by trimmed mean of M values (TMM) [2] is performed by using the calcNormFactors function from edgeR R package. It calculates a set of normalization factors, one for each sample, to eliminate composition biases between libraries.

Here, the table contains library size and normalisation factors using TMM method ordered by library size. These graphs are another way of verifying the quality control carried out during the sequencing and bioinformatics analysis steps.

Show the code 
dge <- calcNormFactors(se)
#cat("Ordered library size and normalisation factors TMM's method \n")
# dge$samples
dge$samples[order(dge$samples$lib.size), ]
        group lib.size norm.factors Bacteria Treatment Replicate  sample
BA_R_R2  BA_R 11288923    0.9468454       BA         R        R2 BA_R_R2
BA_G_R3  BA_G 12108903    1.0737004       BA         G        R3 BA_G_R3
BA_G_R2  BA_G 13758228    1.0889577       BA         G        R2 BA_G_R2
BA_R_R3  BA_R 13912355    0.9092578       BA         R        R3 BA_R_R3
BA_R_R1  BA_R 13997036    0.9886040       BA         R        R1 BA_R_R1
BA_G_R1  BA_G 18392123    1.0048875       BA         G        R1 BA_G_R1
Show the code 
# barplot
p<-ggplot(dge$samples, aes(x=dge$samples$sample, y=dge$samples$norm.factors, fill=dge$samples$group)) +
    geom_col() + xlab("Samples") + ylab("TMM factors") + labs(fill = "group") + theme_bw() +
    theme(axis.text.x = element_text(angle=90)) + ggtitle("TMM size factors \n")
p

Show the code 
# Library size with group
p<-ggplot(dge$samples, aes(x=dge$samples$sample, y=dge$samples$lib.size, fill=dge$samples$group)) +
    geom_col() + xlab("Samples") + ylab("Library size") + labs(fill = "group") + theme_bw() +
    theme(axis.text.x = element_text(angle=90)) + ggtitle("Library size \n")
p

Show the code 
# bp <- ggplot(dge$samples, aes(x=dge$samples$group, y=dge$samples$lib.size, fill=dge$samples$group)) + geom_boxplot() + 
#       geom_jitter(position=position_jitter(0.1)) + theme_bw() + 
#     theme(axis.text.x = element_text(angle=90)) +
#       xlab("group") + ylab("Library size") + labs(fill = "group")
# bp

Multidimensional scaling plot

Multidimensional scaling (MDS) plot shows the relationships between the samples. The top (500 genes) are used to calculate the distance between expression profiles of samples. The distance approximate the log2 fold change between the samples.

As expected, samples are grouped by treatment.

Show the code 
limma::plotMDS(dge, col = as.numeric(dge$samples$group), labels = dge$samples$sample, cex=0.6, top = 500, main = "MDS plot top 500 genes")

Reproducibility token

Show the code 
sessioninfo::session_info(pkgs = "attached")
─ Session info ───────────────────────────────────────────────────────────────
 setting  value
 version  R version 4.4.2 (2024-10-31)
 os       Ubuntu 24.04.1 LTS
 system   x86_64, linux-gnu
 ui       X11
 language (EN)
 collate  fr_FR.UTF-8
 ctype    fr_FR.UTF-8
 tz       Europe/Paris
 date     2025-02-16
 pandoc   3.1.1 @ /usr/lib/rstudio/resources/app/bin/quarto/bin/tools/ (via rmarkdown)

─ Packages ───────────────────────────────────────────────────────────────────
 package              * version date (UTC) lib source
 Biobase              * 2.64.0  2024-04-30 [1] Bioconductor 3.19 (R 4.4.0)
 BiocGenerics         * 0.50.0  2024-04-30 [1] Bioconductor 3.19 (R 4.4.0)
 data.table           * 1.16.4  2024-12-06 [1] CRAN (R 4.4.2)
 dplyr                * 1.1.4   2023-11-17 [1] CRAN (R 4.4.0)
 edgeR                * 4.2.2   2024-10-13 [1] Bioconductor 3.19 (R 4.4.1)
 forcats              * 1.0.0   2023-01-29 [1] CRAN (R 4.4.0)
 GenomeInfoDb         * 1.40.1  2024-05-24 [1] Bioconductor 3.19 (R 4.4.0)
 GenomicRanges        * 1.56.2  2024-10-09 [1] Bioconductor 3.19 (R 4.4.1)
 ggplot2              * 3.5.1   2024-04-23 [1] CRAN (R 4.4.0)
 IRanges              * 2.38.1  2024-07-03 [1] Bioconductor 3.19 (R 4.4.1)
 kableExtra           * 1.4.0   2024-01-24 [1] CRAN (R 4.4.0)
 lattice              * 0.22-5  2023-10-24 [4] CRAN (R 4.3.1)
 limma                * 3.60.6  2024-10-02 [1] Bioconductor 3.19 (R 4.4.1)
 lubridate            * 1.9.4   2024-12-08 [1] CRAN (R 4.4.2)
 MASS                 * 7.3-61  2024-06-13 [4] CRAN (R 4.4.1)
 MatrixGenerics       * 1.16.0  2024-04-30 [1] Bioconductor 3.19 (R 4.4.0)
 matrixStats          * 1.5.0   2025-01-07 [1] CRAN (R 4.4.2)
 mixOmics             * 6.28.0  2024-04-30 [1] Bioconductor 3.19 (R 4.4.2)
 purrr                * 1.0.2   2023-08-10 [1] CRAN (R 4.4.0)
 readr                * 2.1.5   2024-01-10 [1] CRAN (R 4.4.0)
 reshape2             * 1.4.4   2020-04-09 [1] CRAN (R 4.4.0)
 rtracklayer          * 1.64.0  2024-04-30 [1] Bioconductor 3.19 (R 4.4.0)
 S4Vectors            * 0.42.1  2024-07-03 [1] Bioconductor 3.19 (R 4.4.1)
 stringr              * 1.5.1   2023-11-14 [1] CRAN (R 4.4.0)
 SummarizedExperiment * 1.34.0  2024-05-01 [1] Bioconductor 3.19 (R 4.4.0)
 tibble               * 3.2.1   2023-03-20 [1] CRAN (R 4.4.0)
 tidyr                * 1.3.1   2024-01-24 [1] CRAN (R 4.4.0)
 tidyverse            * 2.0.0   2023-02-22 [1] CRAN (R 4.4.0)

 [1] /home/orue/R/x86_64-pc-linux-gnu-library/4.4
 [2] /usr/local/lib/R/site-library
 [3] /usr/lib/R/site-library
 [4] /usr/lib/R/library

──────────────────────────────────────────────────────────────────────────────

References

1. Dillies M-A, Rau A, Aubert J, Hennequet-Antier C, Jeanmougin M, Servant N, et al. A comprehensive evaluation of normalization methods for illumina high-throughput RNA sequencing data analysis. Brief Bioinform. 2012;14:671–83.
2. Robinson MD, Oshlack A. A scaling normalization method for differential expression analysis of RNA-seq data. Genome Biology. 2010;11:R25. doi:10.1186/gb-2010-11-3-r25.

Reuse

This document will not be accessible without prior agreement of the partners

A work by Migale Bioinformatics Facility
Université Paris-Saclay, INRAE, MaIAGE, 78350, Jouy-en-Josas, France
Université Paris-Saclay, INRAE, BioinfOmics, MIGALE bioinformatics facility, 78350, Jouy-en-Josas, France