```{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/tcga-bic_',
  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}
#### Load required R libraries ####

## CRAN libraries
requiredLib <- c("knitr", 
                 "FactoMineR", 
                 "factoextra",
                 "openxlsx" ## to write excel tables
                 )
for (lib in requiredLib) {
  if (!require(lib, character.only = TRUE)) {
    install.packages(lib, )
  }
  require(lib, character.only = TRUE)
}

## Bioconductor libraries
if (!requireNamespace("BiocManager", quietly = TRUE))
    install.packages("BiocManager")

## biomaRt is used to get  information (name, description) from Ensembl gene IDs
if (!require("biomaRt", character.only = TRUE)) {
  BiocManager::install("biomaRt")
}
require("biomaRt", character.only = TRUE)


```


## Downloading and loading the data and metadata files

The preprocessed datasets are available on the course github repository. 

### A convenient function to download files only once

We will use the function `download_only_once()` that we defined in a previous course on data exploration. 

This function takes as input a base URL, a file name and a local folder. 

- It checks if the local folder already exists, and if not  creates it. 
- It then checks if the file is already present in this local folder, and if not downloads it. 

This facilitates the downloading of the different files required for the practical. Open the code box below and run the code in your R console.


```{r function_download_only_once}
#' @title Download a file only if it is not yet here
#' @author Jacques van Helden email{Jacques.van-Helden@@france-bioinformatique.fr}
#' @param url_base base of the URL, that will be prepended to the file name
#' @param file_name name of the file (should not contain any path)
#' @param local_folder path of a local folder where the file should be stored
#' @return the function returns the path of the local file, built from local_folder and file_name
#' @export©
download_only_once <- function(
  url_base, 
  file_name,
  local_folder) {

  ## Define the source URL  
  url <- file.path(url_base, file_name)
  message("Source URL\n\t",  url)

  ## Define the local file
  local_file <- file.path(local_folder, file_name)
  
  ## Create the local data folder if it does not exist
  dir.create(local_folder, showWarnings = FALSE, recursive = TRUE)
  
  ## Download the file ONLY if it is not already there
  if (!file.exists(local_file)) {
    message("Downloading file from source URL to local file\n\t", 
            local_file)
    download.file(url = url, destfile = local_file)
  } else {
    message("Local file already exists, no need to download\n\t", 
            local_file)
  }
  
  return(local_file)
}
```


### Exercise: data download

1. Use the function `download_only_once()` to download the BIC  transciptome data and metadata from the github web site and store it in a local folder (for example `~/m3-stat-R/TCGA-BIC_analysis`). 

Base URL: <https://github.com/DU-Bii/module-3-Stat-R/raw/master/stat-R_2021/data/TCGA_BIC_subset/>

Files : 

- Expression table (1000 top-ranking genes from differential analysis): 

    - File name: `BIC_log2-norm-counts_edgeR_DEG_top_1000.tsv.gz`

    - This file contains log2-transformed and standardised counts, with 1000 genes (rows) x 819 samples (columns)

- Metadata: 

    - File name: `BIC_sample-classes.tsv.gz`

    - This file indicates the status of the 3 marker genes traditionnally used to diagnose the breast cancer type, as well as the cancer class derived from these 3 markers
    - For more information, see the pre-processing report:  [import_TCGA_from_Recount.html](https://du-bii.github.io/study-cases/Homo_sapiens/TCGA_study-case/import_TCGA_from_Recount.html)


2. After having downloaded these files, load them in variables having the following names (for the sake of consistency with the course material). 

- Expression table: bic_expr
- Metadata: bic_meta


### Solution: data download and load

```{r download_and_load}
## Define the remote URL and local folder
bic_url <- "https://github.com/DU-Bii/module-3-Stat-R/raw/master/stat-R_2021/data/TCGA_BIC_subset/"

bic_folder <- "~/m3-stat-R/TCGA-BIC_analysis"

## Download and load the expression data table
## Note: we use check.names=FALSE to avoid replacing hyphens by dots
## in sample names, because we want to keep them as in the 
## original data files. 
bic_expr_file <- download_only_once(
  url_base = bic_url, 
  file_name = "BIC_log2-norm-counts_edgeR_DEG_top_1000.tsv.gz",
  local_folder = bic_folder)
bic_expr <- read.delim(file = bic_expr_file, 
                       header = TRUE, 
                       row.names = 1,
                       check.names = FALSE)
# colnames(bic_expr)
# dim(bic_expr)
# View(head(bic_expr))

## Download the metadata file
bic_meta_file <- download_only_once(
  url_base = bic_url, 
  file_name = "BIC_sample-classes.tsv.gz",
  local_folder = bic_folder)
bic_meta <- read.delim(file = bic_meta_file, 
                       header = TRUE, 
                       row.names = 1,
                       check.names = FALSE)

```


## Exploring the metadata

### Exercise: metadata exploration

Check the content of the metadata file by looking at the first rows, and count the number of samples per class. 

### Solution: metadata exploration

```{r metadata_exploration}

## Show the head of the metadata table
kable(head(bic_meta, n = 10), caption = "First rows of the BIC metadata table")

## Number of samples per class
kable(sort(table(bic_meta$cancer.type), decreasing = TRUE), 
      caption = "Number of samples per class")

```


## Sorting samples by cancer class


Sort the expression and metadata tables so that the samples with the same class come together. This is a bit triccky so we provide immediately the solution, but you might attempt todo it in you way if you have time. 


### Some tips: sorting samples

- The expression table has one row per gene and one column per sample, whereas the metadata file has one row per sample. 

- You can use the function `order()` to obtain the indices of the metadata table by order of sample class. This will return a vector with all the indixes of the first class, then all the indices of the second class, etc. 

- You can then use the indices of the dataframes to re-order them. 

- We are not sure that the samples have the same order in the metadata and in the expression tables. So we will use the sample class to re-order the rows (samples) of the metadata file, but after this is done we will use the sample IDs of the metadata file (in the row names) in order to sort the columns of the expression table. This will guarantee that samples are in a consistent order between data and metadata tables.
    

### Solution: sorting samples

```{r sort_samples}
## Check that the row names of the metadata contain the same set 
## of IDs than the column names of the expression table
length(rownames(bic_meta)) ## the metadata contains 819 rows
length(colnames(bic_expr)) ## the data contains 819 columns
length(intersect(rownames(bic_meta), colnames(bic_expr))) ## Their intersection contains the same number of elements

## Get the sample order according to the cancer type column of the metadata
sample_order <- order(bic_meta$cancer.type)
print(bic_meta$cancer.type[sample_order])

## Sort the metadata rows (samples) according to this order
bic_meta <- bic_meta[sample_order, ]

## Sort the expression table to make sure that the samples (columns)
## come in the same order as the rows of the metadata table
bic_expr <- bic_expr[, row.names(bic_meta)]

```

## Generate readable sample labels

The TCGA samples have very long and structured IDs, whcih make them inconvenient to display on the type of graphs we would like to generate (box plots, PC plots, heatmaps, ...). 

We will thus generate readable labels, and add them as a separate column to the metadata table.

For this, we will use the class name and use the funciotn `make.names()` to avoid having duplicate labels (this will add a number besides the class name). 

```{r sample_labels}
#### Sample labels ####
bic_meta$label <- make.names(
  names = bic_meta$cancer.type, 
  unique = TRUE)

```


## Get gene IDs and names

The row names of the expression table indicate the Ensembl gene ID and the version of the genome annotation (a number after the dot). This ID is quite long, and we would prefer to have more readable labels for the genes. For this, we will get the gene name (also called "gene symbol"). 

```{r get_gene_names}
#### Get gene names ####

gene_info <- data.frame(
  ENSG = sub(x = row.names(bic_expr), 
             perl = TRUE, 
             replacement = "", 
             pattern = '\\..*'),
  row.names = row.names(bic_expr)
)
# head(gene_info)


## Get gene names from ENSEMBL using biomaRt package
library("biomaRt")
ensembl <- useMart("ensembl", dataset = "hsapiens_gene_ensembl")
ensembl_info <- getBM(attributes = c('ensembl_gene_id',
                                     'external_gene_name',
                                     'description'),
                      filters = 'ensembl_gene_id',
                      values = gene_info$ENSG,
                      mart = ensembl)
row.names(ensembl_info) <- ensembl_info$ensembl_gene_id

## Add gene name in a new column of the gene_info table
gene_info$name <- ensembl_info[gene_info$ENSG, "external_gene_name"]
gene_info$description <- ensembl_info[gene_info$ENSG, "description"]

## There is one gene with no name
# sum(is.na(gene_info$name))
## Replace undefined gene names by the gene ID
gene_info[is.na(gene_info$name), "name"] <-
  gene_info[is.na(gene_info$name), "ENSG"]

## Print the first rows of the gene info table
kable(head(gene_info, n = 10), caption = "First rows of the gene information table collected from BioMart")

## Export the gene info table in a TSV file
gene_info_file_tsv <- file.path(bic_folder, "BIC_top1000_gene_info.tsv")
message("Saving gene info table in TSV file\n\t", 
        gene_info_file_tsv)
write.table(x = gene_info, 
            file = gene_info_file_tsv, 
            quote = FALSE, 
            sep = "\t", 
            row.names = TRUE, 
            col.names = NA)

## Export the gene info table in Excel (xlsx) file
gene_info_file_xlsx <- file.path(bic_folder, "BIC_top1000_gene_info.xlsx")
message("Saving gene info table in xlsx file\n\t",
        gene_info_file_xlsx)
openxlsx::write.xlsx(x = gene_info, file = gene_info_file_xlsx)


```



## Assign sample colors

We add here a code to assign a color to each sample according to its class. 

```{r sample_colors}
#### Class and sample colors ####
class_names <- unique(bic_meta$cancer.type)

## Assign one color to each class
class_color <- c(
  "Basal.like" = "#DDEEFF",  
  "HER2pos" = "#88FF88",
  "Luminal.A" = "#FFBB55",
  "Luminal.B" = "#EE88FF",
  "Unclassified"  = "#8888FF"
)

## Assign a color to each sample according to its class
bic_meta$color <- class_color[bic_meta$cancer.type]
# table(bic_meta$color)

```

## Sample standardisation

We will standardise the samples using robust estimators, to avoid the impact of extreme values. 

```{r}
#### Sample standardisation ####

## Compute descriptive stats on the samples 
bic_meta$mean <- apply(bic_expr, 2, mean)
bic_meta$sd <- apply(bic_expr, 2, sd)
bic_meta$median <- apply(bic_expr, 2, median)
bic_meta$iqr <- apply(bic_expr, 2, IQR)

## Centering on the mean
expr_median <- median(unlist(bic_expr))
bic_expr_centered <- t(t(bic_expr) - bic_meta$median) + expr_median
# dim(bic_expr)
# dim(bic_expr_centered)
# apply(bic_expr_centered, 2, median)

## Compute the IQR of the whole expression table
expr_iqr <- IQR(unlist(bic_expr_centered)) 

## Scale each sample based  on its IQR
bic_expr_std <- t((t(bic_expr) - unlist(bic_meta$median)) 
                  /  unlist(bic_meta$iqr)) * expr_iqr + expr_median
# apply(bic_expr_std, 2, IQR)
# range( apply(bic_expr_std, 2, IQR))
# range( apply(bic_expr_std, 2, median))

```

## Dataset with readable row and column names

Create a dataset with 

- gene names as row.names 
- sample labels as column.names

```{r dataset_with_labels}
#### Data set with readable 
bic_expr_labels <- bic_expr_std
dim(bic_expr_labels)
colnames(bic_expr_labels) <- bic_meta$label
rownames(bic_expr_labels) <- make.names(gene_info$name) 

```


## Box plot

### Exercise: box plot

Select 30 samples at random, and draw a boxplot.

Sort the samples according to their index in the columns of the expression table (this will ensure that samples of the same class come together). 

### Solution: box plot

```{r box_plot_rand_samples, fig.height=8, fig.width=12, out.width="100%",  fig.cap="Box plot of a random selection of samples from the TCGA-BIC transcriptome dataset. "}
#### Box plot ####
sample_size <- 30

## select sample indices
selected_samples <- sort(sample(
  x = 1:ncol(bic_expr), 
  size = sample_size,
  replace = FALSE))

## Box plots
par.ori <- par(no.readonly = TRUE)
par(mfrow = c(1, 3))
par(mar = c(5, 6, 4, 1))
## Original data
boxplot(x = bic_expr[, selected_samples], 
        names = bic_meta$label[selected_samples],
        col = bic_meta$color[selected_samples],
        xlab = "log2(counts)",
        main = "Before normalisation",
        horizontal = TRUE, 
        las = 1, notch = TRUE,
        cex.axis = 1)
boxplot(x = bic_expr_centered[, selected_samples], 
        names = bic_meta$label[selected_samples],
        col = bic_meta$color[selected_samples],
        xlab = "log2(counts)",
        main = "Median-based centering",
        horizontal = TRUE, 
        las = 1, notch = TRUE,
        cex.axis = 0.7)
boxplot(x = bic_expr_std[, selected_samples], 
        names = bic_meta$label[selected_samples],
        col = bic_meta$color[selected_samples],
        xlab = "log2(counts)",
        main = "IQR-based scaling",
        horizontal = TRUE, 
        las = 1, notch = TRUE,
        cex.axis = 0.7)
par(par.ori)
```



## PCA

### Exercise

Generate a graph with the two first components and color samples according to their class. Do the sample classes segregate on this PC plot?


### Solution


```{r pca}
#### PC plot of the TAGC BIC samples ####
library(factoextra)


## Compute the principal components
## This is done on the transposed data table
bic_pca <- PCA(t(bic_expr_labels), 
               scale.unit = FALSE, 
               ncp = ncol(bic_expr_labels), 
               graph = FALSE)

## Table with the coordinates of each sample in the PC space
bic_pcs <- bic_pca$ind$coord

# Check the PC dimensions
dim(bic_pcs)


```

```{r pca_variance_lot, fig.width=7, fig.height=5, out.width="60%", fig.cap="Variance of the first components for the BIC dataset"}

#### Plot the variance per component ####
fviz_eig(bic_pca, addlabels = TRUE)

```

The plots of individuals show that the cancer types segregate partly onb the first and second PCs. 

```{r pca_ind_plot_PC1-2}
#### Plot individuals on the 2 first components ####
fviz_pca_ind(X = bic_pca, 
             axes = c(1, 2), 
             col.ind = bic_meta$cancer.type,
             label = "none") 
```

The 3rd and 4th PCs seem much less informative with respect to the cancer types: the individuals seem  intermingled irrespective of their cancer type. 

```{r pca_ind_plot_PC3-4}
#### Plot individuals on the 2 first components ####
fviz_pca_ind(X = bic_pca, 
             axes = c(3, 4), 
             col.ind = bic_meta$cancer.type,
             label = "none") 
```


## Variables

At the end of this tutorial, we dispose of the following data types. 

| Variable name | Data content |
|--------------|----------------------------------------------|
| `bic_meta` | metadata with a few added columns (sample color, estimators of central tendency and dispersion) |
| `gene_info`| ID, name and description of the 1000 genes used here |
| `bic_expr` | non-normalised expression table |
| `bic_expr_centered` | median-based expression table |
| `bic_expr_std` | sample-wise standardised expression table, all samples having the same median and IQR |
| `bic_expr_labels` | same content as bic_expr_std but with row names replaced by human-readable gene names, and column names by human-readable sample labels |

## Save a memory image

We can now save a memory image, which will enable us to reload the whole data with a single command for the next tutorials. 

```{r save_image}
#### Save memory image ####
mem_image <- file.path(bic_folder, "bic_data.Rdata")
message("Saving memory image in file\n\t", mem_image)
save.image(file = mem_image)
```

The memory image has been saved in a local folder, in the following file: `r mem_image`

## Command to reload this dataset

The memory image generated at the end of this tutorial is stored on the github repository, and can be reloaded as follows. 

```{r reload_mem_image, eval=FALSE}
#### Reload memory image from github repository ####
github_mem_img <- "https://github.com/DU-Bii/module-3-Stat-R/blob/master/stat-R_2021/data/TCGA_BIC_subset/bic_data.Rdata?raw=true"

## Define local destination folder
bic_folder <- "~/m3-stat-R/TCGA-BIC_analysis"
## Create it if required
dir.create(bic_folder, showWarnings = FALSE, recursive = TRUE)

## Define local destination for the memory image
mem_image <- file.path(bic_folder, "bic_data.Rdata")
if (file.exists(mem_image)) {
  message("Memory image already there, skipping download")
} else {
  message("Downloading memory image from\n", github_mem_img)
  download.file(url = github_mem_img, destfile = mem_image)
  message("Local memory image\t", mem_image)
}

## Load the memory image
message("Loading memory image", mem_image)
load(mem_image)

```



## Session info

As usual, we write the session info in the report for the sake of traceability. 

```{r session_info}
sessionInfo()

```

