---
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

(**<font color="red">Denis, please give the parameters</font>**)


### 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: <https://google.github.io/styleguide/Rguide.xml>
    - 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



