--- title: 'RNA-seq: analysis of Psoriasis data (Li et al, 2014)' author: "Denis Puthier and Jacques van Helden" date: '`r Sys.Date()`' output: html_document: fig_caption: yes highlight: zenburn self_contained: yes theme: cerulean toc: yes toc_depth: 3 toc_float: yes slidy_presentation: fig_caption: yes fig_height: 4 fig_width: 6 highlight: tango incremental: no keep_md: no self_contained: yes slide_level: 2 smaller: yes theme: cerulean toc: yes widescreen: yes ioslides_presentation: colortheme: dolphin fig_caption: yes fig_height: 4 fig_width: 6 fonttheme: structurebold highlight: tango incremental: no keep_md: no slide_level: 2 smaller: yes theme: cerulean toc: yes widescreen: yes pdf_document: fig_caption: yes highlight: zenburn toc: yes toc_depth: 3 beamer_presentation: colortheme: dolphin fig_caption: yes fig_height: 4 fig_width: 6 fonttheme: structurebold highlight: tango incremental: no keep_tex: no slide_level: 2 theme: Montpellier toc: yes font-import: http://fonts.googleapis.com/css?family=Risque font-family: Garamond transition: linear bibliography: "RNA-seq_psoriasis_analysis_biblio.bib" #css: ../../html/course.css --- ```{r setup, include=FALSE} knitr::opts_chunk$set(echo = FALSE, eval = TRUE, cache = TRUE, message = FALSE, warning = FALSE, comment = "") ## Options to display pretty numbers library(knitr) knit_hooks$set(inline = function(x) { prettyNum(x, big.mark=" ") }) options(scipen = 6, digits = 3) ``` # Exercise Li and coworkers [-@Li:2014dp] used RNA-seq to characterize the transcriptome of 92 skin samples of people suffering from psoriasis, and 82 control samples. We pre-processed the raw sequences (NGS reads) to obtain a data table containing the counts per gene (row) for each sample (columns). ## Goals of this practical 1. **Detection and interpretation of differentially expressed genes:** run a statistical analysis of the read counts, in order to detect differentially expressed genes (**DEG**) between psoriasis and control samples [@Li:2014dp], and to study the functional associations of these genes, and provide a biological interpretation of these results. 2. **Methodological evaluation:** a. **Robusntess analysis:** analyse the impact of individual particularities of the patients/controls, and measure the impact of the sample size by applying a sub-sampling approach. b. **Negative control:** check the reliability of the DEG detection method (DESeq2) by running the same analysis with two subset of samples belonging to the same group. For example, select 10 WT samples and 10 other WT samples, and run differential analysis between them. In principle, the program should return a *negative* result, i.e. declare that no gene is differentially expressed, since the samples come from the same group. Consequently, if the negative control declares genes as differentially expressed, these genes should be considered as false positives. The good answer would thus be to declare not a single DEG, or, if some genes are declared differentially expressed, it should be with a barely significant p-value. If the negative control returns many DEG and/or associates genes with very low p-values, it means that we have a problem either with the method used, or with the dataset (for example it does not comply with the underlying assumptions for the test). More precisely, the negative control is an empirical way to measure if the actual rate of false positive coresponds to the expected rate(i.e. if the p-value, or derived statistics, can be considered as reliable indications of the significance). 3. **Usage of the good practices**. An mportant goal of this course -- and the report -- is to learn how to use good practices for the analysis of NGS data. This includes a. **Tractability:** you and other people should be able to track the origin of all your results. For this, you need to keep a trace of each step of each analysis. b. **Reproducibility:** other people should be able not only to trace the origin of your results, but also to reproduce them by themselves. NGS data anlaysis lends itself particularly well to reproducibility, since everything is done on computers and managed via software (tools and scripts). c. **Portability:** the analysis done on your computer should be reproducible on other computers as well. For this you need to ensure for isntance that all paths are defined relative rather than absolute, and that the path definitions rely on platform-independent methods. **************************************************************** ## Dataset The raw datasets were downloaded from Gene Expression Omnibus series **[GSE54456](https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE54456)**. ### Processing (**Denis, please give the parameters**) ### Access to the count tbale ```{r} ## If required, download the count table from the Web site (this is done only once). url.counts <- "https://dputhier.github.io/ASG/practicals/rnaseq_psoriasis_Li2014/data/counts_all_sample_header.txt.gz" ## Local paths dir.psoriasis <- ("~/ASG/practicals/rnaseq_psoriasis_Li2014/") dir.counts <- file.path(dir.psoriasis, "data") file.counts <- file.path("data", "counts_all_sample_header.txt.gz") ## Create a directory to download the dataset if it does not exist dir.create(dir.counts, showWarnings = FALSE, recursive = TRUE) setwd(dir.psoriasis) ## Download the psoriasis file if required if (!file.exists(file.counts)) { download.file(url=url.counts, destfile = file.counts) } counts <- read.delim(file.counts, row.names = 1) names(counts) <- sub(pattern="_count.txt", replacement="", x=names(counts)) names(counts) <- sub(pattern="__", replacement="_", x=names(counts)) # View(counts) # head(names(counts)) kable(counts[101:108, 1:3], caption = "**Table 1.** Small piece of the count table. ") ``` The [count table](`r file.counts`) contains `r nrow(counts)` rows (one per gene) and `r ncol(counts)` columns (one per sample). ```{r pheno_table} ## Generate a "pheno" table (description of each sample) from the headers of the count table. This is a bit tricky, but since we are often led to use such tricks it is good to see how to proceed. ## Load the limma library which implement strplit2, muech more convenient than strsplit if (!require(limma)) { source("http://bioconductor.org/biocLite.R") biocLite("limma") } library(limma) psoriasis.pheno <- data.frame(strsplit2(names(counts), split = "_")) names(psoriasis.pheno) <- c("GSM_ID", "M_ID", "Group","tissue", "SRX_ID") row.names(psoriasis.pheno)<- names(counts) # head(psoriasis.pheno) # View(psoriasis.pheno) ## Write a file with the pheno setwd(dir.psoriasis) file.pheno <- file.path("data", "counts_all_sample_pheno.tab") write.table(psoriasis.pheno, file=file.pheno, quote=FALSE, row.names = TRUE, col.names = TRUE, sep="\t") ``` The [pheno table](`r file.pheno`) contains `r nrow(psoriasis.pheno)` rows (one per sample) and `r ncol(psoriasis.pheno)` columns (one per sample attribute). ```{r} kable(head(psoriasis.pheno), row.names = FALSE, caption = "**Table 2. First rows of the pheno table. ** Each row describes one sample. ") kable(as.data.frame(table(psoriasis.pheno$Group)), col.names = c("Group", "Number of samples"), caption = "**Table 3. Number of samples per group.** ") ``` ## Analyses 1. **Differential analysis**. Analyse the full dataset to detect differentially expressed genes. This analysis should include the following steps. a. Log2-transformation of the counts (with an epsilon). b. Graphical description of the data (historgrams, barplots, boxplots, ...). c. Computation of summary statistics per sample (min, max, mean, median, quartiles, number of zeros, ...) for raw counts and log2-transformed counts. d. Detection of differentially expressed genes (including graphical representations: volcano plot, p-value histogram, ...). e. Functional enrichment of the differentially expressed genes. f. ... any other type of analysis, figure, table that you might find useful to intepret the data. 2. **Sub-sampling**. Run the same analysis on randomly selected subsets of the samples, with various subset sizes (n=2,3,4,5,10,20,40). 3. **Negative control. ** Run the same analysis with two subsets of samples belonging to the same group (psoriasis versus psoriasis, control versus control). **************************************************************** # Report Les rapports peuvent être rédigés en français ou en anglais / reports can be written in either English or French. ## Format of the report 1. **Source document in Rmd**. The primary report is an R markdown document (.Rmd extension) which **must** contain all the code used to run the anlayses, as well as the main tables and figures produced by the analysis, and a text structured according to the common practice for scientific articles. - The R code should be compliant with the following guidelines: - The R code should be properly documented. - This document should enable us to reproduce your analysis on our own computer (avoid absolute paths). 2. **Compiled report (html or pdf)**. This report should look like a small scientific article (see structure hereafter) with Figures, Tables, and interpretation of the results. The **R** code should not be displayed in the compiled report (set the knitr option `echo=FALSE` when generating the last version). Think of your report as a document written for a biologists who want to understand the approach and the results, but is not interested by the technical details of the R programming. ## Structure of the report In total, the report should not exceed 5-6 pages, including figures, but without counting the bibliographic references and appendices (for which there is no limit). 1. **Introduction:** a brief summary (5-10 lines) of the biological context (the disease, the transcriptome), the biological question addressed in the report, and the general approach envisaged to answer these questions. 2. **Material and Methods:** a summary of the bioinformatics / statistical methods and libraries used for the analysis (1/2 to 1 page), with a brief explanation about which tool was used to do what, and links to the official web page or publication about the tool. If required, additional details about the methods and parameters can be provided in appendix. 3. **Data description:** a brief description of the data source (with link to the GEO record), its content (how many samples, how many groups, ...) and the data type (paired-ends or single-end, ...). 4. **Results and discussion:** the results should be presented and discussed together in this section. This section contains the main figures and tables that are used to interpret the results. Additional figures and tables can be provided as appendices, or as separate files (for example full result tables with all the genes, ...). 5. **Conclusion and perspective** (~1/2 page): summarize the results, show in what they did -- or did not -- enable you to answer the initial questions of the introduction, and add some perspectives about possible future extensions of the work presented here. 6. **Appendices:* any additional information or result that might be helpful to consult in order to get a deeper understanding of your results. **Format:** the report can be submitted in either html or pdf format. The original Rmd file that was used to generate the report must be submitted together with the report. This file should allow the teachers to reprodue the analysis on their computers. ## Evaluation criteria The evaluation of your report will be based on multiple criteria. 1. Can we reproduce your results using your Rmd file? 2. Did you clearly formulate the biological question, and the relationship between these questions and the bioinformatics approaches used to answer them? 3. Relevance of the results and your interpretation. 4. Clarity of the text. 5. Clarity of the code. **************************************************************** # References