RESTORBIOME2: Exploration and normalisation

Experiment with BC

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")
BC_G_R1 BC_G_R2 BC_G_R3 BC_R_R1 BC_R_R2 BC_R_R3
gene-BBCT_RS00005 2244 1589 2602 1720 1491 1500
gene-BBCT_RS00010 3481 2690 3314 944 883 786
gene-BBCT_RS00015 293 224 348 86 91 60
gene-BBCT_RS00020 105 82 119 16 34 19
gene-BBCT_RS00025 5255 3925 5994 10416 9780 8155
gene-BBCT_RS00030 5678 4422 6988 3117 3020 2712

This matrix contains raw counts from BC experiment with 1719 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
BC G R1 BC_G_R1 BC_G
BC G R2 BC_G_R2 BC_G
BC G R3 BC_G_R3 BC_G
BC R R1 BC_R_R1 BC_R
BC R R2 BC_R_R2 BC_R
BC R R3 BC_R_R3 BC_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 BC 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
BC_G_R1 14335635
BC_G_R2 10399918
BC_G_R3 13878411
BC_R_R1 13452356
BC_R_R2 12861188
BC_R_R3 12069280

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 1693 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
BC_G_R2  BC_G 10399748    1.3293540       BC         G        R2 BC_G_R2
BC_R_R3  BC_R 12069105    0.7091508       BC         R        R3 BC_R_R3
BC_R_R2  BC_R 12860993    0.7463302       BC         R        R2 BC_R_R2
BC_R_R1  BC_R 13452165    0.7623289       BC         R        R1 BC_R_R1
BC_G_R3  BC_G 13878114    1.4185006       BC         G        R3 BC_G_R3
BC_G_R1  BC_G 14335421    1.3143705       BC         G        R1 BC_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