--- title: "DU-Bii Study cases - Mouse fibrotic kidney" author: "Olivier Sand and Jacques van Helden" date: '`r Sys.Date()`' output: html_document: self_contained: no fig_caption: yes highlight: zenburn theme: cerulean toc: yes toc_depth: 3 toc_float: yes code_folding: "hide" ioslides_presentation: slide_level: 2 self_contained: no colortheme: dolphin fig_caption: yes fig_height: 5 fig_width: 7 fonttheme: structurebold highlight: tango smaller: yes toc: yes widescreen: yes slidy_presentation: smart: no slide_level: 2 self_contained: yes fig_caption: yes fig_height: 5 fig_width: 7 highlight: tango incremental: no keep_md: yes smaller: yes theme: cerulean toc: yes widescreen: yes beamer_presentation: colortheme: dolphin fig_caption: yes fig_height: 5 fig_width: 7 fonttheme: structurebold highlight: tango incremental: no keep_tex: no slide_level: 2 theme: Montpellier toc: yes revealjs::revealjs_presentation: theme: night transition: none self_contained: true css: ../slides.css pdf_document: fig_caption: yes highlight: zenburn toc: yes toc_depth: 3 powerpoint_presentation: slide_level: 2 fig_caption: yes fig_height: 5 fig_width: 7 toc: yes font-import: http://fonts.googleapis.com/css?family=Risque subtitle: DUBii 2020 font-family: Garamond transition: linear editor_options: chunk_output_type: console --- ```{r settings, include=FALSE, echo=FALSE, eval=TRUE} options(width = 300) # options(encoding = 'UTF-8') knitr::opts_chunk$set( fig.width = 7, fig.height = 5, fig.path = 'figures/mouse-kidney_', fig.align = "center", size = "tiny", echo = TRUE, eval = TRUE, warning = FALSE, message = FALSE, results = TRUE, comment = "") options(scipen = 12) ## Max number of digits for non-scientific notation # knitr::asis_output("\\footnotesize") ``` ```{r libraries, echo=FALSE, eval=TRUE} requiredLib <- c("knitr", "stringr", "vioplot") for (lib in requiredLib) { if (!require(lib, character.only = TRUE)) { install.packages(lib, ) } require(lib, character.only = TRUE) } ``` ## Reference Pavkovic, M., Pantano, L., Gerlach, C.V. et al. Multi omics analysis of fibrotic kidneys in two mouse models. Sci Data 6, 92 (2019). - - - Mouse fibrotic kidney browser: - Data on Zenodo: ## Samples >Samples from two mouse models were collected. The first one is a reversible chemical-induced injury model (folic acid (FA) induced nephropathy). The second one is an irreversible surgically-induced fibrosis model (unilateral ureteral obstruction (UUO)). mRNA and small RNA sequencing, as well as 10-plex tandem mass tag (TMT) proteomics were performed with kidney samples from different time points over the course of fibrosis development. > In the FA model, mice were sacrificed before the treatment (day 0) and 1, 2, 7, and 14 days after a single injection of folic acid. > For the UUO model, mice were sacrificed before obstruction (day 0) and 3, 7, and 14 days after the ureter of the left kidney was obstructed via ligation. For both studies, kidneys were removed at each time point for total RNA isolation and protein sample preparation. We will first explore the UUO transcriptome data. ## Data sources | Doc | URL | |:-----------------------------|:--------------------------| | Total RNA for the experiment on Unilateral ureter obstruction (UUO) model | | ## Parameters ```{r parameters} #### Define parameters for the analysis #### ## Keep a trace of the original parameters par.ori <- par(no.readonly = TRUE) ## Analysis parameters parameters <- list( dataset = "uuo", ## Supported: uuo, fa datatype = "transcriptome", epsilon = 0.1, minCount = 10, forceDownload = FALSE) kable(as.data.frame(parameters)) ``` ## Output directories ```{r outdirs} #### Output directories #### outdirs <- list() # outdirs$main <- getwd() outdirs$main <- "." ## Data directory, where the data will be downloaded and uncompressed outdirs$data <- file.path(outdirs$main, "data") dir.create(outdirs$data, recursive = TRUE, showWarnings = FALSE) ## Main result directory outdirs$results <- file.path(outdirs$main, "results") # Transcriptome results outdirs$transcriptome <- file.path(outdirs$results, "transcriptome") dir.create(outdirs$transcriptome, recursive = TRUE, showWarnings = FALSE) ``` ## Download transcriptome data ```{r download_transcriptome} #### Download transcriptome data #### archiveFile <- "MouseKidneyFibrOmics-v1.0.zip" archiveURL <- file.path("https://zenodo.org/record/2592516/files/hbc", archiveFile) localDataArchive <- file.path(outdirs$data, archiveFile) if (file.exists(localDataArchive) & !parameters$forceDownload) { message("Data archive already downloaded:\n\t", localDataArchive) } else { message("Downloading data archive from zenodo: ", archiveURL) download.file(url = archiveURL, destfile = localDataArchive) ## Uncompess the archive message("Uncompressing data archive") unzip(zipfile = localDataArchive, exdir = outdirs$data) } ## Define destination directory # outdirs$csv <- file.path(outdirs$data, "CSV") # dir.create(outdirs$csv, showWarnings = FALSE, recursive = TRUE) ``` ## Load raw counts ```{r load_raw_counts} #### Load raw counts data table #### ## Note: the "raw" counts are decimal numbers, I suspect that they have been somewhat normalised. To check. rawCountFile <- file.path( outdirs$data, paste0("hbc-MouseKidneyFibrOmics-a39e55a/tables/", parameters$dataset, "/results/counts/raw_counts.csv.gz")) # "hbc-MouseKidneyFibrOmics-a39e55a/tables/fa/results/counts/raw_counts.csv.gz") rawValues <- read.csv(file = rawCountFile, header = 1, row.names = 1) ``` The RNA-seq transcriptome data was loaded as raw counts. This table contains `r nrow(rawValues)` rows (genes) and `r ncol(rawValues)` columns (samples). ## Build metadata table ```{r metadata} #### Build metadata table #### metadata <- data.frame( dataType = "transcriptome", sampleName = colnames(rawValues)) metadata[ , c("condition", "sampleNumber")] <- str_split_fixed(string = metadata$sampleName, pattern = "_", n = 2) ## Colors per condition colPerCondition <- c(normal = "#BBFFBB", day3 = "#FFFFDD", day7 = "#FFDD88", day14 = "#FF4400") metadata$color <- colPerCondition[metadata$condition] ``` ## Compute sample-wise statistics ```{r sample_wise_stat} sampleStat <- metadata sampleStat$mean <- apply(X = rawValues, 2, mean) sampleStat$median <- apply(X = rawValues, 2, median) sampleStat$sd <- apply(X = rawValues, 2, sd) sampleStat$min <- apply(X = rawValues, 2, min) sampleStat$perc05 <- apply(X = rawValues, 2, quantile, prob = 0.05) sampleStat$Q1 <- apply(X = rawValues, 2, quantile, prob = 0.25) sampleStat$median <- apply(X = rawValues, 2, quantile, prob = 0.5) sampleStat$Q3 <- apply(X = rawValues, 2, quantile, prob = 0.75) sampleStat$perc95 <- apply(X = rawValues, 2, quantile, prob = 0.95) sampleStat$max <- apply(X = rawValues, 2, max) sampleStat$iqr <- apply(X = rawValues, 2, IQR) ## Print statistics per sample kable(format(x = format(digits = 5, sampleStat))) ``` ## Distribution of raw counts ```{r RNA-seq_count_distrib, fig.width=12, fig.height=5, out.width="70%", fig.cap="Distribution of raw counts"} hist(unlist(rawValues), breaks = 1000, main = "Raw count distribution", xlab = "Raw counts", ylab = "Number of genes (all samples)") ``` The distribution of raw counts is not very informative, because the range is defined by some outlier, i.e. a gene having a huge number of reads. Even with strong zoom on the abcsissa range from 0 to 500, the histogram shows a steep drop in the first bins. ## Distribution of raw counts - truncated abscissa ```{r RNA-seq_count_distrib_truncated, fig.width=12, fig.height=5, out.width="70%", fig.cap="Distribution of raw counts"} #### Count distrib - truncated abscissa #### hist(unlist(rawValues), breaks = 500000, main = "Raw count distribution", xlab = "Raw counts (truncated abscissa", ylab = "Number of genes (all samples)", xlim = c(0, 500)) ``` ## Log2-transformed counts A typical approach is to normalise the counts by applying a log2 transformation . This however creates a problem when the counts of a given gene in a given sample is 0. To circumvent this, we can add an epsilon ($\epsilon = `r parameters$epsilon`$) before the log2 transformation. ```{r log2_counts, fig.width=12, fig.height=5, out.width="70%"} #### Log2 transformatiojn of the counts #### log2Values <- log2(rawValues + parameters$epsilon) hist(unlist(log2Values), breaks = seq(from = -5, to = 22, by = 0.1), main = "log2-counts distribution", xlab = "log2(Counts + epsilon)", ylab = "Number of genes", col = "#BBFFBB") ``` ## Two-columns test :::::::::::::: {.columns} ::: {.column} contents 1 ... ::: ::: {.column} pas l'air de marcher ... ::: :::::::::::::: ## Box plots We can now inspect the distribution of counts per sample with the `boxplot()` function. ```{r box plots, fig.width=5, fig.height=7, out.width="25%"} #### Box plots #### par(mar = c(4, 6, 5, 1)) boxplot(log2Values, col = metadata$color, horizontal = TRUE, las = 1, main = "log2-transformed", xlab = "log2(counts)") par(par.ori) ``` ## Box plot comment We notice an obvious problem: the vast majority of counts is very small. This can result from different causes, which will not be investigated in this context. ## Gene filtering ```{r gene_filtering, fig.width=5, fig.height=7, out.width="50%"} ## Filter out the genes with very low counts in all conditions undetectedGenes <- apply(rawValues, MARGIN = 1, FUN = sum) < parameters$minCount # table(undetectedGenes) log2Filtered <- log2Values[!undetectedGenes, ] ``` We filtered out all the genes whose maximal count value across all samples was lower than `r parameters$minCount`. Among the `r nrow(rawValues)` genes from the raw count table, `r sum(undetectedGenes)` were considered undetected according to this criterion. We use the remaining `r nrow(log2Filtered)` genes for the subsequent analyses. ## Boxplot after gene filtering ```{r boxplot_filtered, fig.width=5, fig.height=7, out.width="32%"} #### Box plots after filtering #### par(mar = c(4, 6, 5, 1)) boxplot(log2Filtered, col = metadata$color, horizontal = TRUE, las = 1, main = "Filtered genes", xlab = "log2(counts)") par(par.ori) ``` ## Normalisation (more precisely: scaling) Before going any further, it is important to ensure some normalisation of the counts, in order to correct for biases due to inter-sample differences in sequencing depth. For the sake of simplicity, we will use here a very simple criterion: median-based normalisation. The principle is to multiply the counts of each sample by a scaling factor in order to bring each sample to the same median count. ## Normalisation code ```{r median_normalisation, fig.width=5, fig.height=7, out.width="25%"} #### Median-based normalisation #### sampleMedians <- apply(log2Filtered, 2, median) seriesMedian <- median(sampleMedians) scalingFactors <- seriesMedian / sampleMedians log2Standardised <- data.frame(matrix( nrow = nrow(log2Filtered), ncol = ncol(log2Filtered))) colnames(log2Standardised) <- colnames(log2Filtered) rownames(log2Standardised) <- rownames(log2Filtered) for (j in 1:ncol(log2Filtered)) { log2Standardised[, j] <- log2Filtered[, j] * scalingFactors[j] } ## Check the remaining medians apply(log2Standardised, 2, median) ``` ## Normalized boxplot ```{r median_normalisation_boxplot, fig.width=5, fig.height=7, out.width="32%"} #### Box plots after scaling #### par(mar = c(4, 6, 5, 1)) boxplot(log2Standardised, col = metadata$color, horizontal = TRUE, las = 1, main = "Median-based scaled", xlab = "log2(counts") par(par.ori) ``` ## Violin plots We can also inspect the distribution of counts per sample with the `vioplot()` function. ```{r violin plots, fig.width=5, fig.height=7, out.width="25%"} #### Violin plots #### par(mar = c(4, 6, 5, 1)) vioplot::vioplot(log2Values, col = metadata$color, horizontal = TRUE, las = 1, main = "log2-transformed", xlab = "log2(counts)") par(par.ori) ``` ## Scatter plots We can also inspect the distribution of counts per sample with the `plot()` function. ```{r scatter plots, fig.width=5, fig.height=7, out.width="25%"} #### Scatter plots #### #par(mar = c(4, 6, 5, 1)) #plot(log2Values[,metadata$sampleName], # col = metadata$color) #par(par.ori) ``` ## Combinations ## Exporting the result We export the pre-processed data in separate tables for further reuse. ```{r save_tables} #### Save tables #### outfiles <- vector() ## Raw counts, all the variables outfiles["raw"] <- file.path(outdirs[parameters$datatype], paste0(parameters$dataset, "_", parameters$datatype, "_raw.tsv.gz")) write.table(x = format(digits = 3, big.mark = "", decimal.mark = ".", rawValues), dec = ".", file = gzfile(outfiles["raw"], "w"), quote = FALSE, sep = "\t") ## Log2-transformed counts, all the genes outfiles["log2"] <- file.path(outdirs[parameters$datatype], paste0(parameters$dataset, "_", parameters$datatype, "_log2.tsv.gz")) write.table(x = format(digits = 3, big.mark = "", decimal.mark = ".", log2Values), dec = ".", file = gzfile(outfiles["log2"], "w"), quote = FALSE, sep = "\t") ## Filtered genes only, log2-transformed counts outfiles["filtered"] <- file.path(outdirs[parameters$datatype], paste0(parameters$dataset, "_", parameters$datatype, "_log2_filtered.tsv.gz")) write.table(x = format(digits = 3, big.mark = "", decimal.mark = ".", log2Filtered), dec = ".", file = gzfile(outfiles["filtered"], "w"), quote = FALSE, sep = "\t") ## Filtered genes only, log2-transformed and standardized counts outfiles["standardised"] <- file.path(outdirs[parameters$datatype], paste0(parameters$dataset, "_", parameters$datatype, "_log2_norm.tsv.gz")) write.table(x = format(digits = 3, big.mark = "", decimal.mark = ".", log2Standardised), dec = ".", file = gzfile(outfiles["standardised"], "w"), quote = FALSE, sep = "\t") ## Metadata outfiles["metadata"] <- file.path(outdirs[parameters$datatype], paste0(parameters$dataset, "_", parameters$datatype, "_metadata.tsv")) write.table(x = metadata, file = outfiles["metadata"] , quote = FALSE, sep = "\t") ## Build a table to display the links in the report fileTable <- data.frame(outfiles) fileTable$basename <- basename(fileTable$outfiles) fileTable$dirname <- dirname(fileTable$outfiles) fileTable$link <- paste0("[", fileTable$basename, "]", "(", fileTable$outfiles, ")") ## Print the directories kable(t(as.data.frame.list(outdirs)), col.names = "Directories", caption = "Directories") ## Print the link table kable(cbind(rownames(fileTable), fileTable$link), row.names = FALSE, col.names = c("Contents", "Link"), caption = "Output files") ``` ## Session info ```{r session_info} sessionInfo() ```