---
title: "MOFA2 analysis of Corces-Buenrostro AML dataset, Part 2"
author: "Delphine Potier & Carl Herrmann"
date: "`r date()`"
output: 
  html_document:
    toc: true
    toc_float: true
    toc_depth: 4
    number_sections: false
    code_folding: show
    highlight: zenburn
params:
  work_dir: "~/mydatalocal/"
editor_options: 
  chunk_output_type: console
---

**Data pre-processing is crucial for MOFA analysis.**

**We will go back to the different pre-processing steps (and to the various MOFA options) to investigate how they can impact the results.**

The following script allow to obtain the multiview matrix we used to perform MOFA analysis. At different points, the possibility to use another option will be proposed, you can choose to modify if and investigate how it impact the results.

```{r setup, include=FALSE}
knitr::opts_knit$set(root.dir=params$work_dir)
knitr::opts_chunk$set(echo = TRUE)
options(knitr.table.format="html")
setwd(params$work_dir)
```

# Library installation/loading

```{r libraries_loading, warning=FALSE, message=FALSE, results='hide'}
BiocManager::install("gage", update = FALSE)
# Load libraries
library(data.table)
library(ggplot2)
library(DESeq2) 
library(MOFA2)
library(FactoMineR)
library(factoextra)
library(ComplexHeatmap)
library(viridis)
library(DT)
library(msigdbr)
library(GGally)
library(dplyr)
library(gage)
```



# Data loading and preprocessing

## RNA-seq

**RNA-seq raw data loading**

```{r RNAseq-data_loading}
##----------------------------------------------------------------------------##
##                             Download counts                                ##
##----------------------------------------------------------------------------##
# set ftp url to RNA-seq data
ftp_url <- file.path("ftp://ftp.ncbi.nlm.nih.gov/geo/series/GSE74nnn/GSE74246",
                     "suppl/GSE74246_RNAseq_All_Counts.txt.gz")

read_delim_gz <- function(file_url) {
  con <- gzcon(url(file_url))
  txt <- readLines(con)
  return(read.delim(textConnection(txt), row.names = 1))
}

# read in data matrix
corces_rna_counts <- read_delim_gz(ftp_url)
```



**Data cleaning**

```{r RNAseq-data_cleaning, fig.width = 9, fig.height = 8}
##----------------------------------------------------------------------------##
##                      Data selection and sample QC                          ##
##----------------------------------------------------------------------------##

#Remove leukemic and erythroblast samples:
corces_rna_counts <- corces_rna_counts[,-grep("Ery|rHSC|LSC|Blast", colnames(corces_rna_counts))]

#Inspect correlation matrix:
cor_dm <- cor(corces_rna_counts)
Heatmap(cor_dm, col = magma(100), name = "Correlation") # need library(ComplexHeatmap)

#X5852.GMP is an outlier and will be removed, has much smaller library size as other GMPS:
corces_rna_counts <- corces_rna_counts[,-grep("X5852.GMP", colnames(corces_rna_counts))]

#Remove rows with rowSum==0:
corces_rna_counts <- corces_rna_counts[!rowSums(corces_rna_counts) == 0,]

rm(cor_dm)
dim(corces_rna_counts)
```

We obtain a matrix with **`r dim(corces_rna_counts)[2]`** samples and **`r dim(corces_rna_counts)[1]`** genes

$~$
$~$
$~$

**Annotation formatting**

```{r RNA-seq_annotation}
####################################################### 
##----------------------------------------------------------------------------##
##                              Annotation                                    ##
##----------------------------------------------------------------------------##
# extract celltypes from colnames
col.anno <- gsub(".*\\.", "", colnames(corces_rna_counts))
#Use short names
col.anno[grep("NK", col.anno)] <- "NK"
col.anno[grep("CD4", col.anno)] <- "CD4"
col.anno[grep("CD8", col.anno)] <- "CD8"


# Define color vector
type.color <- setNames(c("#771155", "#AA4488", "#CC99BB", "#114477", "#4477AA", "#77AADD",
                         "#117777", "#44AAAA", "#77CCCC", "#777711", "#AAAA44", "#DDDD77"),
                       c("HSC", "MPP", "LMPP", "CMP", "GMP", "MEP",
                         "CLP", "CD4", "CD8", "NK", "Bcell", "Mono"))

# type.color <- setNames(c("#356345", "#75A962", "#5CB4A1", "#F3CD99", "#F5BE44", "#D75B60",
#                          "#B6E1E9", "#4083B3", "#2E368B", "#B18DBA", "#55286D", "#EE7C32"),
#                        c("HSC", "MPP", "LMPP", "CMP", "GMP", "MEP",
#                          "CLP", "CD4", "CD8", "NK", "Bcell", "Mono"))

# Annotation data frame
corces_rna_annot <- data.frame(sample = colnames(corces_rna_counts),
                               Celltype = as.factor(col.anno),
                               color    = type.color[match(col.anno, names(type.color))],
                               row.names = colnames(corces_rna_counts),
                               stringsAsFactors = FALSE)
```


$~$
$~$
$~$


## ATAC-seq

**Raw data loading**

```{r ATACseq-data_loading}
# Download and read RNAseq data.
##––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––##
##                                 Set paths                                  ##
##––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––##
atac.data.path  <- file.path("data/atacseq/")
GSE74912.path <- file.path(atac.data.path, "GSE74912_ATACseq_All_Counts.txt.gz")
dir.create(atac.data.path, recursive = TRUE, showWarnings = FALSE)

##––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––##
##                               Download files                              ##
##––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––##
# download data
if(!file.exists(GSE74912.path)){
  download.file(url = "ftp://ftp.ncbi.nlm.nih.gov/geo/series/GSE74nnn/GSE74912/suppl/GSE74912_ATACseq_All_Counts.txt.gz", GSE74912.path)
}

##––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––##
##                                Data loading                                ##
##––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––##
# read in data matrix
atac.counts <- read.delim(gzfile(GSE74912.path), stringsAsFactors = FALSE)
#atac.counts[1:5,1:5]
```



**Data cleaning**

```{r ATACseq-data_cleaning, fig.width = 9, fig.height = 8}
##----------------------------------------------------------------------------##
##                   Data selection and preprocessing                         ##
##----------------------------------------------------------------------------##

# Separate annotation column from data.
atac.row.anno <- atac.counts[,1:3]
atac.counts <- atac.counts[,-c(1:3)]

# Remove rows with rowSums < 2000 and leukemic and erythroblast samples.
rownames(atac.counts) <- do.call(paste, c(as.list(atac.row.anno), sep = "_"))

# remove rows with rowSums < 2000
atac.row.anno <- atac.row.anno[rowSums(atac.counts) > 2000,]
atac.counts <- atac.counts[rowSums(atac.counts) > 2000,]

# remove leukemic and erythroblast samples
atac.counts <- atac.counts[,-grep("Ery|LSC|pHSC|Leuk|CD34", colnames(atac.counts))]


#Remove X6792.7A, due to low coverage.
atac.counts <- atac.counts[,-grep("X6792.7A", colnames(atac.counts))]

saveRDS(atac.counts, file = "data/atac_counts.RDS")
```

We obtain a matrix with **`r dim(atac.counts)[2]`** samples and **`r dim(atac.counts)[1]`** genes

$~$
$~$
$~$

**Annotation formatting**

```{r ATACseq_annotation}

##––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––##
##                              Annotation                                    ##
##––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––##
# extract celltypes from colnames
col.anno <- colnames(atac.counts)
col.anno[grep("CD4", col.anno)] <- "CD4"
col.anno[grep("CD8", col.anno)] <- "CD8"
col.anno[grep("NK", col.anno)] <- "NK"
col.anno[grep("Nkcell", col.anno)] <- "NK"
col.anno[grep("Bcell", col.anno)] <- "Bcell"
col.anno[grep("CLP", col.anno)] <- "CLP"
col.anno[grep("1(A|B)", col.anno)] <- "HSC"
col.anno[grep("2(A|B)", col.anno)] <- "MPP"
col.anno[grep("3(A|B)", col.anno)] <- "LMPP"
col.anno[grep("4(A|B)", col.anno)] <- "CMP"
col.anno[grep("5(A|B)", col.anno)] <- "GMP"
col.anno[grep("6(A|B)", col.anno)] <- "MEP"
col.anno[grep("7(A|B)", col.anno)] <- "Mono"


# Annotation data frame
atac.annot <- data.frame(sample = colnames(atac.counts),
                         Celltype = as.factor(col.anno),
                         color    = type.color[match(col.anno, names(type.color))],
                         row.names = colnames(atac.counts),
                         stringsAsFactors = FALSE)
```


$~$
$~$
$~$


## Match RNAseq and ATACseq  

Clean to keep only RNA-seq and ATAC-seq matching samples

*Of note, it is possible in MOFA to work with partially unmatched data : you can choose to modify the script an run MOFA with some unmatched datasets*

```{r}
##––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––##
##                            ATACseq annotations                             ##
##––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––##

# replace numbered celltypes by character names
atac.anno.cellID <- col.anno
rm(col.anno)


# Paste donor ID and cell type
atac.anno <- paste0(sapply(strsplit(colnames(atac.counts), "\\."), "[[", 1), ".", atac.anno.cellID)

##––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––##
##                             RNAseq annotations                             ##
##––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––##
# Keep only Donor ID
rnaseqIDs <- setNames(sapply(strsplit(colnames(corces_rna_counts), "\\."), "[[", 1), colnames(corces_rna_counts))
rnaseqIDs <- sub("^X", "", rnaseqIDs)

##––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––##
##                         Match RNAseq and ATACseq                           ##
##––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––##
rna.atac.matched.samples <- lapply(setNames(1:length(rnaseqIDs), names(rnaseqIDs)), function(i) {
  # Find same donor in RNAseq and ATACseq
  atac.matched <- grep(rnaseqIDs[i], atac.anno, value = TRUE)
  if (length(atac.matched) > 0) {
    # if same donor in both omics search if same cell type
    # extract cell tupe fron RNAseq colname
    anno <- sapply(strsplit(names(rnaseqIDs)[i], "\\."), "[[", 2)
    anno <- sub("Tcell", "", anno)
    anno <- sub("cell", "", anno)
    # find in ATACseq
    anno.matched <- grep(anno, atac.matched, value = TRUE)
    if (anno == "MPP") {
      anno.matched <- grep("LMPP", anno.matched, value = TRUE, invert = TRUE)
    }
    if (length(anno.matched) > 0) {
      data.frame(rnaID  = names(rnaseqIDs)[i],
                 atacID = colnames(atac.counts)[atac.anno %in% anno.matched],
                 cellID = atac.anno.cellID[atac.anno %in% anno.matched],
                 atac.anno = anno.matched, 
                 row.names = colnames(atac.counts)[atac.anno %in% anno.matched])
    }
  }
} )

# Keep only matched samples
rna.atac.matched.samples <- rna.atac.matched.samples[!sapply(rna.atac.matched.samples, is.null)]

##––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––##
##                             Format annotations                             ##
##––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––##
# Keep only one ATACseq replicate
rna.atac.annot <- do.call(rbind, rna.atac.matched.samples)
rna.atac.annot <- rna.atac.annot[!duplicated(rna.atac.annot$rnaID),]

# Format annotation
rna.atac.annot$original.atacID <- rna.atac.annot$atacID
rna.atac.annot$atacID <- rna.atac.annot$rnaID
rownames(rna.atac.annot) <- rna.atac.annot$rnaID 


rna.atac.annot <- rna.atac.annot[,c(2,3,1,4,5)]
colnames(rna.atac.annot) <- c("sample", "Celltype", "rna.sample", "atac.sample", "original.atacID")

rna.atac.annot$color <- type.color[match(rna.atac.annot$Celltype, names(type.color))]

#Save annotations
saveRDS(rna.atac.annot, file = "data/atac_annotations.RDS")

##––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––##
##                          Removing unmatched samples                        ##
##––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––##

## RNA-seq ##
#Keep only samples matching with ATAC-seq
corces_rna_counts <- corces_rna_counts[,rna.atac.annot$rna.sample]

#Remove rows with rowSum==0:
corces_rna_counts <- corces_rna_counts[!rowSums(corces_rna_counts) == 0,]


### ATAC-seq ##
#Keep only samples matching with ATAC-seq
atac.counts <- atac.counts[,rna.atac.annot$original.atacID]

#Remove rows with rowSum==0:
atac.counts <- atac.counts[!rowSums(atac.counts) == 0,] #useless the size stay the same
```

We obtain a matrix with **`r dim(corces_rna_counts)[2]`** samples and **`r dim(corces_rna_counts)[1]`** genes for RNA-seq and a matrix with **`r dim(atac.counts)[2]`** samples and **`r dim(atac.counts)[1]`** genes for ATAC-seq.

## Data normalization {.tabset}

MOFA requires a proper normalization for the model to work.
In our case (RNA-seq / ATAC-seq, which are count based data), it is recommended to perform size factor normalisation (to remove library size effects) + variance stabilisation (data need to be continuous and normal-ish distributed; the closer it looks to a gaussian  the better, but this is not necessary).


If this step is not done correctly, the model might learn a very strong Factor 1 that will capture differences in the total expression per sample, and more subtle sources of variation will be downweighted.

*Note 1 : We followed MOFA recommendation but you can choose to use an other normalization method (e.g. use log2 normalization that was used with ButchR in the 1 part of the workshop or rlog). The script to perform a log2 normalization of the data is already here, you can check the histograms to see count distribution before and after normalization with vst or log2. You are free to test any other method you like or check what it gives without normalization.*


*Note 2 : It could be interesting to filter out genes display a very low expression in to few samples, if you want to try the command lines are in comments.*

### RNA-seq {.tabset}
```{r RNA-seq_data_preprocessing}
##----------------------------------------------------------------------------##
##                              Normalize counts                              ##
##----------------------------------------------------------------------------#
dds <- DESeqDataSetFromMatrix(countData = corces_rna_counts, colData = corces_rna_annot[rna.atac.annot$rna.sample,], design = ~ Celltype)

# do size factor normalization
dds <- estimateSizeFactors(dds)
#sizeFactors(dds)

# # It could be interesting to filter for low expression
# # Code to filter out genes where there are less than 1 sample (for B cell we have only one sample) with normalized counts greater than or equal to 5.
# idx <- rowSums( counts(dds, normalized=TRUE) >= 5 ) >= 1
# dds <- dds[idx,]


## Using VST
###############################

# do vst  (recommanded by MOFA)
vsd <- vst(dds, blind = FALSE) # blind = FALSE to take in account the experimental design

corces_rna_norm <- assay(vsd)

## Or logs
###############################
# This is how where normalized counts in the first part of this course
# do +1 log2 transformation
corces_rna_norm2 <- apply(counts(dds, normalized=TRUE) + 1, 2, log2)

#corces_rna_norm <- corces_rna_norm2
rm(ftp_url, dds)
```

$~$


#### Histogram of counts before normalisation

```{r RNA-seq_hist_raw}
hist(sapply(corces_rna_counts, as.numeric), breaks = 100, main = "Histogram of counts before normalisation")
```

#### Histogram of vst counts

Normalization recommanded by MOFA
```{r RNA-seq_hist_vst-norm}
hist(corces_rna_norm, breaks = 100, main = "Histogram of VST normalized counts")

hist(corces_rna_norm - rowMeans(corces_rna_norm), breaks = 100, main = "Histogram of VST normalized and centered counts")
```


#### Histogram of log2 counts

Normalization used in multivariate analyis part I

```{r RNA-seq_hist_log-norm}
hist(corces_rna_norm2, breaks = 100, main = "Histogram of log2 normalized counts")

hist(corces_rna_norm2 - rowMeans(corces_rna_norm2), breaks = 100, main = "Histogram of VST normalized and centered counts")
```

$~$
$~$
$~$

### ATAC-seq {.tabset}

```{r ATACseq_data_preprocessing}
##----------------------------------------------------------------------------##
##                              Normalize counts                              ##
##----------------------------------------------------------------------------##
## Using VST
###############################
dds <- DESeqDataSetFromMatrix(countData = atac.counts, colData = atac.annot[rna.atac.annot$original.atacID,], design = ~ Celltype)

# do size factor normalization
dds <- estimateSizeFactors(dds)
#sizeFactors(dds)

# # It could be interesting to filter for too low accessibility
# # Code to filter out genes where there are less than 1 sample (for B cell we have only one sample) with normalized counts greater than or equal to 5.
# idx <- rowSums( counts(dds, normalized=TRUE) >= 5 ) >= 1
# dds <- dds[idx,]

# do vst
vsd <- vst(dds, blind = FALSE) # blind = FALSE to take in account the experimental design

atac.norm.mat <- assay(vsd)

## Or logs
###############################
# This is how where normalized counts in the first part of this course
# do +1 log2 transformation
atac.norm.mat2 <- apply(counts(dds, normalized=TRUE) + 1, 2, log2)

#atac.norm.mat <- atac.norm.mat2
rm(vsd, GSE74912.path, dds)
```

$~$

#### Histogram of counts before normalisation

```{r ATACseq_hist_raw}
hist(sapply(atac.counts, as.numeric), breaks = 100, main = "Histogram of counts before normalisation")
```

#### Histogram of vst counts

Normalization recommanded by MOFA
```{r ATACseq_hist_vst-norm}
hist(atac.norm.mat, breaks = 100, main = "Histogram of VST normalized counts")

hist(atac.norm.mat - rowMeans(atac.norm.mat), breaks = 100, main = "Histogram of VST normalized and centered counts")
```


#### Histogram of log2 counts

Normalization used in multivariate analyis part I

```{r ATACseq_hist_log-norm}
hist(atac.norm.mat2, breaks = 100, main = "Histogram of log2 normalized counts")

hist(atac.norm.mat2 - rowMeans(atac.norm.mat2), breaks = 100, main = "Histogram of VST normalized and centered counts")
```

$~$
$~$
$~$


# Features selection {.tabset}

Features (genes/regions) displaying no or low variation should be removed as they can cause numerical issues in the model.

MOFA recommends to select the top most variable features for each assay (of course other methods can be used).
Moreover our datasets have different sizes; this has to be taken in account by selecting similar amount of features for both as bigger data modalities will tend to be over-represented in the MOFA model and the model might miss sources of variation unique to the small data set).

*This is another crucial step of the data pre-processing you can play with.*

Here we will select the top 5000 most variable features for both RNA-seq and ATAC-seq.

*This number is arbitrary and can be changed, one can also compare with an imbalanced number of selected features depending on the downstream analysis to be run (e.g. higher number of ATAC-seq regions to perform motif enrichment analysis of differently accessible DNA regions).*

*Moreover it is possible to use other selection methods (e.g. PCA, using the top X features participating to the first Y PCs or taking in account or with a more advanced method) *


## RNA-seq {.tabset}

### HVG selection

```{r RNAseq_HVG_selection}
# Select highly variable genes
##############
topVarGenes <- head(order(rowVars(corces_rna_norm), decreasing = TRUE), 5000)
matTopVarGenes <-corces_rna_norm[topVarGenes,]
matTopVarGenes <- matTopVarGenes - rowMeans(matTopVarGenes)
```


#### Heatmap

```{r RNAseq_HVG_heatmap, fig.width = 9}
pheatmap::pheatmap(matTopVarGenes, annotation_col = corces_rna_annot[,2, drop = FALSE], annotation_colors = list(Celltype = c(HSC = "#771155", MPP = "#AA4488", LMPP = "#CC99BB", CMP = "#114477", GMP = "#4477AA", MEP = "#77AADD", Mono = "#DDDD77", CD4 = "#44AAAA", CD8 = "#77CCCC", NK = "#777711", Bcell = "#AAAA44", CLP = "#117777")), show_rownames = FALSE )
```


#### Checking expected markers 

Historically, cluster differentiation genes (usually named CD[0-9]*) are genes used to separate and characterize the different hematopoietic populations by FACS. 

Among them classical markers are :

* Bcells : CD19, CD79A, CD79B
* Tcells : CD3D, CD3E, CD3G, CD4, CD8A, CD8B 
* Monocytes : CD14, CD300E, CD163
* NK cells : CD160 (activation)
* Stage markers: CD24, CD5, CD52, CD53, CD34, CD84, CD33, (+ some markers have different gene/protein names : CD62L (=SELL), CD117 (=KIT),CD45 (=PTPRC), CD90 (=THY1), CD10 (=MME))

$~$
$~$

**Heatmap of some CD genes of interest**

```{r RNAseq_CD_genes_visualization, fig.width = 9, fig.height = 7}
CDgenes <- c("CD19","CD79A","CD3D","CD3E","CD3G","CD4","CD8A","CD8B","CD24","CD5","CD52","CD53","CD34","CD14","CD300E","CD163","CD160","CD74","SELL","CD83","CD58","CD84","CD22","CD300A","KIT","PTPRC","THY1","MME","CD33")

pheatmap::pheatmap(corces_rna_norm[CDgenes,], annotation_col = corces_rna_annot[,2, drop = FALSE], annotation_colors = list(Celltype = c(HSC = "#771155", MPP = "#AA4488", LMPP = "#CC99BB", CMP = "#114477", GMP = "#4477AA", MEP = "#77AADD", Mono = "#DDDD77", CD4 = "#44AAAA", CD8 = "#77CCCC", NK = "#777711", Bcell = "#AAAA44", CLP = "#117777")), show_rownames = TRUE )
```

$~$
$~$

It can be interesting to check which ones we get among HVG if classical genes are found.

```{r RNAseq_selected_CD_genes_visualization, fig.width = 9, fig.height = 10}
CDgenes <- rownames(matTopVarGenes)[grepl("^CD[0-9]+", rownames(matTopVarGenes))]
pheatmap::pheatmap(corces_rna_norm[CDgenes,], annotation_col = corces_rna_annot[,2, drop = FALSE], annotation_colors = list(Celltype = c(HSC = "#771155", MPP = "#AA4488", LMPP = "#CC99BB", CMP = "#114477", GMP = "#4477AA", MEP = "#77AADD", Mono = "#DDDD77", CD4 = "#44AAAA", CD8 = "#77CCCC", NK = "#777711", Bcell = "#AAAA44", CLP = "#117777")), show_rownames = TRUE )
```

**Many CD genes are selected among selected HVG**



### PCA informative genes selection {.tabset}

Alternatively we could select informative genes using PCA

```{r RNAseq_runPCA, warning=FALSE}
# Find informative genes using PCA
##############

# Run PCA
pca_results <- PCA(t(corces_rna_norm), scale.unit=TRUE, ncp=25, graph=F) 

# Check PC1 and PC2
fviz_pca_ind(pca_results, axes = c(1,2))
```

PC1 and PC2 separate immature cell types from the more mature ones

```{r RNAseq_barplotPCA, warning=FALSE}
# Check PC barplot
eig.val <- pca_results$eig
barplot(eig.val[, 2], 
        names.arg = 1:nrow(eig.val), 
        main = "Variances Explained by PCs (%)",
        xlab = "Principal Components",
        ylab = "Percentage of variances",
        col ="steelblue")
# Add connected line segments to the plot
lines(x = 1:nrow(eig.val), eig.val[, 2], 
      type = "b", pch = 19, col = "red")
```

We can see an elbow between 6 and 13 components. 
We can further check the information given by each PC plotting them before choosing the number of PC/genes to select.

#### PC1/2

```{r RNAseq_PC1_PC2, fig.width = 7, fig.height = 7}
fviz_pca_ind(pca_results, axes = c(1,2))+ 
  geom_point(aes(colour = factor(corces_rna_annot[rownames(pca_results$ind$coord),]$color))) +
  guides(colour = guide_legend(title = "color")) + theme(legend.position = "none")
```

PC1 and 2 separate lymphoid cells from myeloid and immature cells

#### PC3/4
```{r RNAseq_PC3_PC4, fig.width = 7, fig.height = 7}
fviz_pca_ind(pca_results, axes = c(3,4))+ 
  geom_point(aes(colour = factor(corces_rna_annot[rownames(pca_results$ind$coord),]$color))) +
  guides(colour = guide_legend(title = "color")) + theme(legend.position = "none")
```

PC3 separate T cells from monocytes and CLP

PC4 separate MEP, LMPP, GMP, and NK cells from the rest

#### PC5/6
```{r RNAseq_PC5_PC6, fig.width = 7, fig.height = 7}
fviz_pca_ind(pca_results, axes = c(5,6))+ 
  geom_point(aes(colour = factor(corces_rna_annot[rownames(pca_results$ind$coord),]$color))) +
  guides(colour = guide_legend(title = "color")) + theme(legend.position = "none")
```

PC5/6 distinguish NK cells from other matuer lymphoid cells (T and B cells) 


#### PC7/8
```{r RNAseq_PC7_PC8, fig.width = 7, fig.height = 7}
fviz_pca_ind(pca_results, axes = c(7,8))+ 
  geom_point(aes(colour = factor(corces_rna_annot[rownames(pca_results$ind$coord),]$color))) +
  guides(colour = guide_legend(title = "color")) + theme(legend.position = "none")
```


PC7 separate B cells fom the rest 

PC8 separate seem noisy (CLP at both ends)

#### PC9/10
```{r RNAseq_PC9_PC10, fig.width = 7, fig.height = 7}
fviz_pca_ind(pca_results, axes = c(9,10))+ 
  geom_point(aes(colour = factor(corces_rna_annot[rownames(pca_results$ind$coord),]$color))) +
  guides(colour = guide_legend(title = "color")) + theme(legend.position = "none")
```

PC9 separate MEP from GMP/LMPP

PC10 seems noisy

#### PC11/12
```{r RNAseq_PC11_PC12, fig.width = 7, fig.height = 7}
fviz_pca_ind(pca_results, axes = c(11,12)) + 
  geom_point(aes(colour = factor(corces_rna_annot[rownames(pca_results$ind$coord),]$color))) +
  guides(colour = guide_legend(title = "color")) + theme(legend.position = "none")
```
Noisy PCs

#### Genes contribution to each PCs

```{r RNAseq_PCAtable, warning=FALSE}
# Get genes participation to each PC
DT::datatable(pca_results$var$contrib,filter = 'top',extensions = 'Buttons', options = list(dom = 'Blfrtip', buttons = c('excel', "csv"), fixedHeader = TRUE))  %>% formatSignif(c(1,c(1:23)))
```

```{r RNAseq_PCA_selection, fig.width = 9, fig.height = 5}
fviz_contrib(pca_results, choice = "var", axes = 1, top = 50)
fviz_contrib(pca_results, choice = "var", axes = 2, top = 50)
fviz_contrib(pca_results, choice = "var", axes = 3, top = 50)
fviz_contrib(pca_results, choice = "var", axes = 4, top = 50)
fviz_contrib(pca_results, choice = "var", axes = 5, top = 50)
fviz_contrib(pca_results, choice = "var", axes = 6, top = 50)
fviz_contrib(pca_results, choice = "var", axes = 7, top = 50)
fviz_contrib(pca_results, choice = "var", axes = 8, top = 50)
fviz_contrib(pca_results, choice = "var", axes = 9, top = 50)
fviz_contrib(pca_results, choice = "var", axes = 10, top = 50)
fviz_contrib(pca_results, choice = "var", axes = 11, top = 50)
fviz_contrib(pca_results, choice = "var", axes = 12, top = 50)
fviz_contrib(pca_results, choice = "var", axes = 13, top = 50)
```

#### Gene selection

We choose to select the 9 first PC only (from PC10, we see mostly noise) and exclude a noisy PC8. Then the top 500 genes for each PC is selected. This is arbitrary and if you like you could do something morre elaborated like taking all genes that explains 50% of variance of each PCs 


```{r RNAseq_PCAgenes_selection, warning=FALSE}
genes_pc1 <- rownames(head(pca_results$var$contrib[order(pca_results$var$contrib[,1], decreasing = TRUE),], n = 500))
genes_pc2 <- rownames(head(pca_results$var$contrib[order(pca_results$var$contrib[,2], decreasing = TRUE),], n = 500))
genes_pc3 <- rownames(head(pca_results$var$contrib[order(pca_results$var$contrib[,3], decreasing = TRUE),], n = 500))
genes_pc4 <- rownames(head(pca_results$var$contrib[order(pca_results$var$contrib[,4], decreasing = TRUE),], n = 500))
genes_pc5 <- rownames(head(pca_results$var$contrib[order(pca_results$var$contrib[,5], decreasing = TRUE),], n = 500))
genes_pc6 <- rownames(head(pca_results$var$contrib[order(pca_results$var$contrib[,6], decreasing = TRUE),], n = 500))
genes_pc7 <- rownames(head(pca_results$var$contrib[order(pca_results$var$contrib[,7], decreasing = TRUE),], n = 500))
#genes_pc8 <- rownames(head(pca_results$var$contrib[order(pca_results$var$contrib[,8], decreasing = TRUE),], n = 500))
genes_pc9 <- rownames(head(pca_results$var$contrib[order(pca_results$var$contrib[,9], decreasing = TRUE),], n = 500))
#genes_pc10 <- rownames(head(pca_results$var$contrib[order(pca_results$var$contrib[,10], decreasing = TRUE),], n = 100))
#genes_pc11 <- rownames(head(pca_results$var$contrib[order(pca_results$var$contrib[,11], decreasing = TRUE),], n = 100))
#genes_pc12 <- rownames(head(pca_results$var$contrib[order(pca_results$var$contrib[,12], decreasing = TRUE),], n = 100))
#genes_pc13 <- rownames(head(pca_results$var$contrib[order(pca_results$var$contrib[,13], decreasing = TRUE),], n = 100))
#genes_pc14 <- rownames(head(pca_results$var$contrib[order(pca_results$var$contrib[,14], decreasing = TRUE),], n = 100))
#genes_pc15 <- rownames(head(pca_results$var$contrib[order(pca_results$var$contrib[,15], decreasing = TRUE),], n = 100))

topPCAGenes <- unique(c(genes_pc1, genes_pc2, genes_pc3, genes_pc4, genes_pc5, genes_pc6, genes_pc7, genes_pc9))#,genes_pc8, genes_pc10, genes_pc11, genes_pc12, genes_pc13, genes_pc14, genes_pc15))
matTopPCAGenes <-corces_rna_norm[topPCAGenes,]
matTopPCAGenes <- matTopPCAGenes - rowMeans(matTopPCAGenes)
```

Doing so we select **`r length(topPCAGenes)`** genes

```{r RNAseq_PCAgenes_selection_heatmap1, fig.width = 9, fig.height = 6}
pheatmap::pheatmap(matTopPCAGenes, annotation_col = corces_rna_annot[,2, drop = FALSE], annotation_colors = list(Celltype = c(HSC = "#771155", MPP = "#AA4488", LMPP = "#CC99BB", CMP = "#114477", GMP = "#4477AA", MEP = "#77AADD", Mono = "#DDDD77", CD4 = "#44AAAA", CD8 = "#77CCCC", NK = "#777711", Bcell = "#AAAA44", CLP = "#117777")), show_rownames = FALSE )

#Check CD genes we get
CDgenes2 <- rownames(matTopPCAGenes)[grepl("^CD[0-9]+", rownames(matTopPCAGenes))]
pheatmap::pheatmap(corces_rna_norm[CDgenes2,], annotation_col = corces_rna_annot[,2, drop = FALSE], annotation_colors = list(Celltype = c(HSC = "#771155", MPP = "#AA4488", LMPP = "#CC99BB", CMP = "#114477", GMP = "#4477AA", MEP = "#77AADD", Mono = "#DDDD77", CD4 = "#44AAAA", CD8 = "#77CCCC", NK = "#777711", Bcell = "#AAAA44", CLP = "#117777")), show_rownames = TRUE )
```


### Summary of genes selected through the different methods {.tabset}

**Intersection**

* **All genes**

**`r length(intersect(rownames(matTopVarGenes),rownames(matTopPCAGenes)))` genes are common to all selection methods**

**`r length(setdiff(rownames(matTopVarGenes),rownames(matTopPCAGenes)))`** genes are selected only in HVG.

**`r length(setdiff(rownames(matTopPCAGenes),rownames(matTopVarGenes)))`** genes are selected only in PCA (method 1).

$~$
$~$

* **Between CD genes**

**`r length(intersect(CDgenes,CDgenes2))` genes are common to all selection methods** (`r toString(intersect(CDgenes,CDgenes2))`)

**`r length(setdiff(CDgenes,CDgenes2))` genes are selected only in HVG :** `r toString(setdiff(CDgenes,CDgenes2))`**, while ** `r toString(setdiff(CDgenes2,CDgenes))` **are not found**

**`r length(setdiff(CDgenes2,CDgenes))` genes are selected only in PCA (method 1) **`r toString(setdiff(CDgenes2,CDgenes))`**, while ** `r toString(setdiff(CDgenes,CDgenes2))` **are not found**


#### HVG All
```{r RNAseq_selected1, fig.width = 9, fig.height = 8}
pheatmap::pheatmap(matTopVarGenes, annotation_col = corces_rna_annot[,2, drop = FALSE], annotation_colors = list(Celltype = c(HSC = "#771155", MPP = "#AA4488", LMPP = "#CC99BB", CMP = "#114477", GMP = "#4477AA", MEP = "#77AADD", Mono = "#DDDD77", CD4 = "#44AAAA", CD8 = "#77CCCC", NK = "#777711", Bcell = "#AAAA44", CLP = "#117777")), show_rownames = FALSE )
```

#### PCA All
```{r RNAseq_selected2, fig.width = 9, fig.height = 8}
pheatmap::pheatmap(matTopPCAGenes, annotation_col = corces_rna_annot[,2, drop = FALSE], annotation_colors = list(Celltype = c(HSC = "#771155", MPP = "#AA4488", LMPP = "#CC99BB", CMP = "#114477", GMP = "#4477AA", MEP = "#77AADD", Mono = "#DDDD77", CD4 = "#44AAAA", CD8 = "#77CCCC", NK = "#777711", Bcell = "#AAAA44", CLP = "#117777")), show_rownames = FALSE )
```



#### HVG CD genes 
```{r RNAseq_selected1b, fig.width = 9, fig.height = 10}
CDgenes <- rownames(matTopVarGenes)[grepl("^CD[0-9]+", rownames(matTopVarGenes))]
pheatmap::pheatmap(corces_rna_norm[CDgenes,], annotation_col = corces_rna_annot[,2, drop = FALSE], annotation_colors = list(Celltype = c(HSC = "#771155", MPP = "#AA4488", LMPP = "#CC99BB", CMP = "#114477", GMP = "#4477AA", MEP = "#77AADD", Mono = "#DDDD77", CD4 = "#44AAAA", CD8 = "#77CCCC", NK = "#777711", Bcell = "#AAAA44", CLP = "#117777")), show_rownames = TRUE )
```

#### PCA CD genes 
```{r RNAseq_selected2b, fig.width = 9, fig.height = 6}
CDgenes2 <- rownames(matTopPCAGenes)[grepl("^CD[0-9]+", rownames(matTopPCAGenes))]
pheatmap::pheatmap(corces_rna_norm[CDgenes2,], annotation_col = corces_rna_annot[,2, drop = FALSE], annotation_colors = list(Celltype = c(HSC = "#771155", MPP = "#AA4488", LMPP = "#CC99BB", CMP = "#114477", GMP = "#4477AA", MEP = "#77AADD", Mono = "#DDDD77", CD4 = "#44AAAA", CD8 = "#77CCCC", NK = "#777711", Bcell = "#AAAA44", CLP = "#117777")), show_rownames = TRUE )
```






## ATAC-seq 

**Highly variable regions selection**

```{r ATACseq_HVG_selection}
# Select highly variable features
##############
topVarRegs <- head(order(rowVars(atac.norm.mat), decreasing = TRUE), 5000) #round(dim(atac.norm.mat)[1]/10)
matTopVarRegs <-atac.norm.mat[topVarRegs,]
matTopVarRegs <- matTopVarRegs - rowMeans(matTopVarRegs)
```

```{r ATACseq_HVR_heatmap, fig.width = 9}
# Heatmap of the most variable region
pheatmap::pheatmap(matTopVarRegs, annotation_col = atac.annot[,2, drop = FALSE], annotation_colors = list(Celltype = c(HSC = "#771155", MPP = "#AA4488", LMPP = "#CC99BB", CMP = "#114477", GMP = "#4477AA", MEP = "#77AADD", Mono = "#DDDD77", CD4 = "#44AAAA", CD8 = "#77CCCC", NK = "#777711", Bcell = "#AAAA44", CLP = "#117777")), show_rownames = FALSE)
```

# Build a multiview matrix (MOFA input)

This multiview matrix will contain RNA-seq and ATAC-seq cleaned and normalized matrix for the selected features.

```{r}
##––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––##
##                   Built multi view to easy access                          ##
##––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––##
multiview.norm.mat <- list(rna  = matTopVarGenes[, match(rna.atac.annot$rna.sample, colnames(matTopVarGenes))],
                           atac = matTopVarRegs[, match(rna.atac.annot$original.atacID, colnames(matTopVarRegs))])

##––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––##
##                   Format multi view to easy access                         ##
##––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––##
# Format Colnames
colnames(multiview.norm.mat$atac) <- colnames(multiview.norm.mat$rna)

##––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––##
##                    Save multi view to easy access                          ##
##––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––##

saveRDS(multiview.norm.mat, file = "data/multiview.RDS")
#multiview.norm.mat <- readRDS("Path/object.RDS")
```



Now you can go back to **Multivariate2_MOFA_part1.Rmd** to rerun MOFA and see how it affects the results compared to your initial run.