EPITOME_CS_integration.Rmd

---
title: "Epigenomic memory of infection: Data integration"
subtitle : "H3K4me2, transcriptome time courses and chromatin state"
author: "Claudia Chica"
date: "02.2021"
output: 
  html_document: 
    keep_md: yes
    number_sections: yes
    smart: no
    toc: yes
    toc_float: yes
editor_options: 
  chunk_output_type: console
---

```{r setup, message = FALSE, warning = FALSE, echo = FALSE}
knitr::opts_chunk$set(message=FALSE,warning=FALSE,echo=FALSE)#, cache=TRUE, cache.lazy = FALSE)
knitr::opts_knit$set(progress = FALSE, verbose = FALSE)
#knitr::opts_chunk$set(dev = "svg")

WDIR="/Volumes/bioit/12891_Chevalier_EpiMemStrep"

library(ggplot2)
library(plyranges)
library(RColorBrewer)
library(knitr)
library(pheatmap)
library(patchwork)
library(FactoMineR)
library(org.Hs.eg.db)
library(msigdbr)


```

# Methylation profile

```{r loadEPIGresults}
# Load epigenomic counts
mark="H3K4me2"
fileOut=paste0(WDIR,"/results_project12891/scripts/Robjects/countNormLength_",mark,".RData")
load(fileOut)

# Discard ONE peak with negative normalised values in a SINGLE replicate (artefact from cycliclowess normalisation)
countNormLength_H3K4me2=countNormLength_H3K4me2[countNormLength_H3K4me2$norm.H3K4me2_UI_REP1>0,]

# Get list of DM peaks
files=list.files(path = paste0(WDIR,"/results_project12891/tables/tables_v3_UCSC/",mark,"/"),pattern = "resAnDif_with_annotation_sorted.final.xls", full.names = T)
DM_results=list() ; i=1
for (fileDM in files) {
  no_col <- max(count.fields(fileDM, sep = "\t"))
  tmp=read.table(fileDM,header = F,sep="\t", fill = TRUE,skip = 1,
                 col.names=1:no_col,colClasses = c("character",rep("numeric",16),rep("character",no_col-17)))
  row.names(tmp) <- tmp[,1]
  tmp=tmp[,14:17]
  colnames(tmp) <- c("baseMean", "log2FoldChange", "pvalue", "padj")
  tmp=tmp[row.names(countNormLength_H3K4me2),]
  DM_results[[i]]=tmp ; i=i+1
}
rm(tmp)
names(DM_results) <- c("D7vsH3","D7vsUI","H3vsUI")

# DM profiles
fileOut=paste0(WDIR,"/results_project12891/scripts/Robjects/DMpeaks_",mark,".RData")
load(fileOut)

# Add not DEGs
DMpeaks_H3K4me2=intersect(DMpeaks_H3K4me2,row.names(DM_results$D7vsUI))
tmp=setdiff(row.names(countNormLength_H3K4me2),DMpeaks_H3K4me2)
DMs=data.frame(peakID=DMpeaks_H3K4me2,
               Kind="DM",
Sign=ifelse(DM_results$D7vsUI[DMpeaks_H3K4me2,]$log2FoldChange>=0,"Gain","Loss"))
notDMs=data.frame(peakID=tmp,
                  Kind="notDM",
Sign=ifelse(DM_results$D7vsUI[tmp,]$log2FoldChange>=0,"Gain","Loss"))
DMs=rbind(DMs,notDMs)
DMs$peakID=as.character(DMs$peakID)
rm(tmp,notDMs)

kable(table(DMs$Kind,DMs$Sign))
```

# Transcriptional profile
```{r loadTOMEresults}
# Expression data normalised and batch corrected 
fileOut=paste0(WDIR,"/results_project12891/scripts/Robjects/Data_RMAandNoBatch_Annot.RData")
load(fileOut)

# DA results, same object used for the functional analysis
fileOut=paste0(WDIR,"/results_project12891/scripts/Robjects/DE_results_Annot.RData")
load(fileOut)

SYMBOL2ENTREZ=unique(data.frame(SYMBOL=as.character(DE_results$H3vsUI$SYMBOL),
                         ENTREZID=as.character(DE_results$H3vsUI$ENTREZID),stringsAsFactors = FALSE))

# DEGs padjTh=.05 ; lfcTh=0
fileOut=paste0(WDIR,"/results_project12891/tables/transcriptome/DEGs_H3K4me2profiles_CompleteAnnotationPeaks_nearestGene.txt")
DEGs=read.table(fileOut,header = T)[,7:10]
DEGs=unique(na.omit(DEGs))
DEGs$SYMBOL=as.character(DEGs$SYMBOL)
DEGs$ENTREZID=as.character(DEGs$ENTREZID)
# Add not DEGs
tmp=setdiff(SYMBOL2ENTREZ$ENTREZID,DEGs$ENTREZID)
H3vsUI=DE_results$H3vsUI[!duplicated(DE_results$H3vsUI$ENTREZID),]
DIvsD7=DE_results$DIvsD7[!duplicated(DE_results$DIvsD7$ENTREZID),]
notDEGs=data.frame(ENTREZID=SYMBOL2ENTREZ[SYMBOL2ENTREZ$ENTREZID %in% tmp,]$ENTREZID,
                   SYMBOL=SYMBOL2ENTREZ[SYMBOL2ENTREZ$ENTREZID %in% tmp,]$SYMBOL,
                   Kind="notDE",
# Sign=ifelse(H3vsUI[H3vsUI$ENTREZID %in% tmp,]$logFC>=0,
#      ifelse(DIvsD7[DIvsD7$ENTREZID %in% tmp,]$logFC>=0,"Up","UpDown"),
#      ifelse(DIvsD7[DIvsD7$ENTREZID %in% tmp,]$logFC<0,"Down","DownUp")))
Sign=ifelse(H3vsUI[H3vsUI$ENTREZID %in% tmp,]$logFC>=0,"Up","Down"))
DEGs=rbind(DEGs,notDEGs)
DEGs$SYMBOL=as.character(DEGs$SYMBOL)
DEGs$ENTREZID=as.character(DEGs$ENTREZID)
rm(tmp,H3vsUI,DIvsD7,notDEGs)

kable(table(DEGs$Kind,DEGs$Sign))

```

```{r load_peaks2Gene}
mark="H3K4me2"
### Nearest gene
fileOut=paste0(WDIR,"/results_project12891/scripts/Robjects/mergedPeaks_annotNearest_",mark,".RData")
load(fileOut)
peaks2GeneNearest_H3K4me2=data.frame(
  peakId=mergedPeaks_annotNearest_H3K4me2@ranges@NAMES,
  SYMBOL=as.character(mergedPeaks_annotNearest_H3K4me2@elementMetadata@listData[["mcols.GENESYMBOL"]]),
  distance=mergedPeaks_annotNearest_H3K4me2@elementMetadata@listData[["distance"]],
  stringsAsFactors=FALSE)
peaks2GeneNearest_H3K4me2=merge(peaks2GeneNearest_H3K4me2,SYMBOL2ENTREZ,by.x=2,by.y=1,all.x=TRUE)[,c(2,1,3,4)]

### Target gene from MEME
# fileOut=paste0(WDIR,"/results_project12891/scripts/Robjects/mergedPeaks_annotTgene_",mark,".RData")
# load(fileOut)
# peaks2GeneTgene_H3K4me2=data.frame(peakId=mergedPeaks_annotTgene_H3K4me2$peakId,
#                   SYMBOL=as.character(mergedPeaks_annotTgene_H3K4me2$Gene_Name),
#                   ENTREZID=mergedPeaks_annotTgene_H3K4me2$ENTREZID,
#                   distance=0,stringsAsFactors=FALSE)

# Load ENTREZ ENSEMBL ids mapping to cross Tgene and EPI-TOME information
ENSEMBLEids <- org.Hs.egENSEMBL 
mapped_genes <- mappedkeys(ENSEMBLEids)
ENTREZtoENSEMBL <- as.data.frame(ENSEMBLEids[mapped_genes])
colnames(ENTREZtoENSEMBL)=c("ENTREZID","ENSEMBLEID")

# Load Tgene MEME results
fileIn=paste0(WDIR,"/results_project12891/MEME_Tgene/H2K4me2/links.tsv")
mergedPeaks_Tgene=read.table(fileIn,header = T,sep = "\t")
## Get only significant associations peak-gene
pvalTh=.05 #; FDRth=.5
mergedPeaks_TgeneSig=subset(mergedPeaks_Tgene,mergedPeaks_Tgene$CnD_P_Value < pvalTh)
mergedPeaks_TgeneSig$chr=gsub("chr","",gsub(":[0-9-]*","",mergedPeaks_TgeneSig$RE_Locus))
mergedPeaks_TgeneSig$start=as.character(as.integer(gsub("chr[0-9XY]*:","",gsub("-[0-9]*","",mergedPeaks_TgeneSig$RE_Locus)))-1)
mergedPeaks_TgeneSig$end=gsub(".*-","",mergedPeaks_TgeneSig$RE_Locus)
mergedPeaks_TgeneSig$peakID=apply(mergedPeaks_TgeneSig[,18:20],1,function(x){paste0(x[1],"_",x[2],"_",x[3])})
mergedPeaks_TgeneSig$ENSEMBL=gsub("\\.[0-9-]*","",mergedPeaks_TgeneSig$Gene_ID)
peaks2GeneTgene_H3K4me2=na.omit(merge(mergedPeaks_TgeneSig,
                           ENTREZtoENSEMBL,by.x=22,by.y=2,all.x=TRUE)[,c(22,23,10,14,17)])

# Select unique gene-region associations 
# select association minimun Tgene pvalue 
peaks2GeneTgene_H3K4me2$peak2gene=paste0(peaks2GeneTgene_H3K4me2$peakID,"_",
  peaks2GeneTgene_H3K4me2$ENTREZID)
peaks2GeneTgene_H3K4me2$tmp=paste0(peaks2GeneTgene_H3K4me2$peakID,"_",
  peaks2GeneTgene_H3K4me2$ENTREZID,"_",as.character(peaks2GeneTgene_H3K4me2$CnD_P_Value))
peaks2GeneTgene_H3K4me2=peaks2GeneTgene_H3K4me2[!duplicated(peaks2GeneTgene_H3K4me2$tmp),]

tmp=aggregate(peaks2GeneTgene_H3K4me2$CnD_P_Value,list(peaks2GeneTgene_H3K4me2$peak2gene),min)
tmp$tmp=paste0(tmp$Group.1,"_",as.character(tmp$x))
peaks2GeneTgene_H3K4me2=peaks2GeneTgene_H3K4me2[peaks2GeneTgene_H3K4me2$tmp %in% tmp$tmp,-7]

rm(mergedPeaks_TgeneSig,ENSEMBLEids,mapped_genes,ENTREZtoENSEMBL)
```

# Chromatin states

```{r selectPeaks, eval=FALSE, include=FALSE}
# Select peaks using counts
logFCthUp_H3K4me2=.005 #quantile .75
logFCthDown_H3K4me2=-.005 #quantile .25
tmp=countNormLength_H3K4me2
selPeaksAllRep=row.names(
  subset(tmp,(
  (log(tmp$norm.H3K4me2_D7_REP1/tmp$norm.H3K4me2_UI_REP1,2) > logFCthUp_H3K4me2 & 
   log(tmp$norm.H3K4me2_D7_REP2/tmp$norm.H3K4me2_UI_REP2,2) > logFCthUp_H3K4me2)  |
  (log(tmp$norm.H3K4me2_D7_REP1/tmp$norm.H3K4me2_UI_REP1,2) < logFCthDown_H3K4me2 &  
   log(tmp$norm.H3K4me2_D7_REP2/tmp$norm.H3K4me2_UI_REP2,2) < logFCthDown_H3K4me2))))

selPeaksAll=row.names(tmp)

# Select peaks related to DEGs 
selPeaks_DEGs_nearest=unique(merge(DEGs[DEGs$Kind!="notDE",],peaks2GeneNearest_H3K4me2[,c(4,1)],by=1)[,5])
selPeaks_DEGs_target=unique(merge(DEGs[DEGs$Kind!="notDE",],peaks2GeneTgene_H3K4me2[,c(3,1)],by=1)[,5])

print(c("Number of:"))
print(c("Total peaks",dim(tmp)[1]))
print(c("Selected peaks - methylation logFold",length(selPeaksAll)))
print(c("Selected peaks - reproducible methylation logFold",length(selPeaksAllRep)))
print(c("Selected peaks - nearest to DEG",length(selPeaks_DEGs_nearest)))
print(c("Selected peaks - with DEG target",length(selPeaks_DEGs_target)))

print(c("Methylation logFold nearest to DEG",length(intersect(selPeaksAll,selPeaks_DEGs_nearest))))
print(c("Methylation logFold with DEG target",length(intersect(selPeaksAll,selPeaks_DEGs_target))))
print(c("Reproducible methylation logFold nearest to DEG",length(intersect(selPeaksAllRep,selPeaks_DEGs_nearest))))
print(c("Reproducible methylation logFold with DEG target",length(intersect(selPeaksAllRep,selPeaks_DEGs_target))))

ggplotdf=data.frame(log2FC=c(
  DM_results$D7vsUI[selPeaksAll,]$log2FoldChange,
  DM_results$D7vsUI[selPeaksAllRep,]$log2FoldChange,
  DM_results$D7vsUI[selPeaks_DEGs_nearest,]$log2FoldChange,
  DM_results$D7vsUI[selPeaks_DEGs_target,]$log2FoldChange),
  peakPopulation=factor(c(rep("All peaks",length(selPeaksAll)),
                   rep("Reproducible Log2FC",length(selPeaksAllRep)),
                   rep("DEGs nearest",length(selPeaks_DEGs_nearest)),
                   rep("DEGs target",length(selPeaks_DEGs_target))),
                   levels = c("All peaks","Reproducible Log2FC","DEGs nearest","DEGs target")))

ggplot(ggplotdf, aes(x=log2FC,fill=peakPopulation)) +
    geom_density(alpha=0.5) + facet_grid(~ peakPopulation)
```

```{r get_mergedPeaks_annotTXS}
# Load all merged peaks coords
fileIn=paste0(WDIR,"/results_project12891/tables/tables_v3_UCSC/",mark,"/",mark,"_mergedPeaks.bed")
mergedPeaks=read_bed(fileIn)

# Annotate merged peaks
fileAnnot=paste0(WDIR,"/results_project12891/igv/igv_session_v3/annotation/hg19_transcripts.geneSymbol.reformat.bed")
txs_gr <- read_bed(fileAnnot)
tss_gr <- GRanges(seqnames = txs_gr@seqnames,
                   ranges = IRanges(start = txs_gr@ranges@start-1e3, end = txs_gr@ranges@start+1e3),
                   strand = txs_gr@strand,
                   mcols = data.frame(GENESYMBOL=mcols(txs_gr)$name))
txs_gr <- GRanges(seqnames = txs_gr@seqnames,
                   ranges = txs_gr@ranges,
                   strand = txs_gr@strand,
                   mcols = data.frame(GENESYMBOL=mcols(txs_gr)$name))

mergedPeaks_annotTXS=array(NA,dim=length(mergedPeaks)) ; names(mergedPeaks_annotTXS)=mergedPeaks$name
ovtss=unique(findOverlaps(mergedPeaks,tss_gr)@from)
mergedPeaks_annotTXS[ovtss]="TSS"
ovintraG=setdiff(unique(findOverlaps(mergedPeaks,txs_gr)@from),ovtss)
mergedPeaks_annotTXS[ovintraG]="intraG"
mergedPeaks_annotTXS[is.na(mergedPeaks_annotTXS)]="interG"
```

```{r annotationQuantification}
DIR=paste0(WDIR,"/results_project12891/Encode_data/bigWig/")
fileMetadata=paste0(DIR,"ENCODE_metadata_bigWig_foldChangeOverControl.released.txt")
Encode_metadata <- read.table(file=fileMetadata,header=TRUE,sep="\t")

marks=c("H3K4me3", "H3K4me2", "H3K27ac", "H3K4me1", "H3K79me2", "H3K27me3")
AnnotationCols <- data.frame(
  Mark=factor(sub("-human","",Encode_metadata$Experiment.target),levels = marks),
  row.names = Encode_metadata$File.accession)

colsG=brewer.pal(n = 3, name = "Greys")
colsM=brewer.pal(n = 4, name = "Greens")
colsS=brewer.pal(n = 3, name = "Blues")[2:3]
nClust=8 ; colsC=brewer.pal(n = nClust, name = "Set2")
colsAnnot <- list(Mark=c("H3K4me3"=colsS[1], "H3K4me2"=colsS[2], "H3K27ac"=colsM[1], "H3K4me1"=colsM[2], "H3K79me2"=colsM[3], "H3K27me3"=colsM[4]),
  DM_sign=c("Gain"="#FF0000", "Loss"="#FF000030"),
  DM_kind=c("DM"="#000000", "notDM"="transparent"),
  Localization=c("TSS"=colsG[1], "intraG"=colsG[2], "interG"=colsG[3]),
  CS=c("1"=colsC[1],"2"=colsC[2],"3"=colsC[3],"4"=colsC[4],"5"=colsC[5],"6"=colsC[6],"7"=colsC[7],"8"=colsC[8]))
```

```{bash getMtx_ENCODElogchange, eval=FALSE, include=FALSE}
WDIR=/pasteur/projets/policy01/BioIT/12891_Chevalier_EpiMemStrep/results_project12891/
module load UCSC-tools/v373
module load bedtools/2.25.0

# 1. Download bigWig files
# 2. Get bedGraph files
cd Encode_data/bigWig/
for file in *.bigWig ; do 
if [ ! -f ${file/bigWig/bedGraph} ] ; then echo $file ; 
srun bigWigToBedGraph $file ${file/bigWig/bedGraph} &  fi ; done

# 3. Overlap with mergedPeaks
mark=H3K4me2
filePeaks=$WDIR/tables/tables_v3_UCSC/*/${mark}_mergedPeaks.bed
cd $WDIR/Encode_data/bigWig

# 3.1. Small files
emacs overlapBedGraph_mergedPeaks.sh
#!/bin/sh
i=$(head -n ${SLURM_ARRAY_TASK_ID} $1 | tail -n 1) 
WDIR=$2
filePeaks=$3
mark=$4

if [ ! -f $WDIR/${i/.bedGraph/_mergedPeaks_${mark}.ov} ] ; then 
srun -o $i.ovb -e $i.e intersectBed -a $filePeaks -b $i -wb
awk '{s[$4]=s[$4]+$NF ; n[$4]=n[$4]+1} END {for (p in n) {print p"\t"s[p]/n[p]}}' $i.ovb | 
sort -k1,1 > $WDIR/${i/.bedGraph/_mergedPeaks_${mark}.ov} ; fi

ls *bedGraph > fileList.txt
n=`wc -l fileList.txt | awk '{print $1}'`
sbatch -p hubbioit --qos hubbioit -o ov.o -e ov.e --mem=200G --array=1-$n%$n overlapBedGraph_mergedPeaks.sh fileList.txt $WDIR/Encode_data/bigWig/ $filePeaks $mark

# 4 Assemble matrix
echo "Peak_name" > header0 ; i=0
for file in *$mark*ov ; do let j=$i+1
echo ${file/_mergedPeaks_${mark}.ov} | paste header$i - > header$j ; ((i++)) ; done
paste *$mark*ov | awk '{line=$1 ; for (i=2;i<=NF;i=i+2) {line=line"\t"$i} ; print line}' |
cat header$j - > ENCODE_mergedPeaks_${mark}.ov

# Clean
rm $WDIR/Encode_data/bigWig/*.ovb $WDIR/Encode_data/bigWig/*.bedGraph
rm $WDIR/Encode_data/bigWig/*.e $WDIR/Encode_data/bigWig/ov.o
rm header*
```

```{r ENCODE_bigWig}
DIR=paste0(WDIR,"/results_project12891/Encode_data/bigWig/")
fileSignal=paste0(DIR,"ENCODE_mergedPeaks_",mark,".ov")
matSignal <- read.table(file=fileSignal,header=TRUE,row.names = 1,sep="\t")

markOrder=c(#H3K4me3
  "ENCFF628ANV", "ENCFF510LTC", "ENCFF556OVF",
  #"ENCFF313GRK", "ENCFF709KWY", "ENCFF635VEW",
  "ENCFF250OYQ", "ENCFF189JCB", "ENCFF421ZNH",
  #H3K4me2
  "ENCFF876DTJ", "ENCFF975YWP", "ENCFF419LFZ",
  #H3K4me1
  "ENCFF165ZPD", "ENCFF569ZJQ", "ENCFF613NHX",
  #H3K27ac
  #"ENCFF216GUZ", "ENCFF597SRK", "ENCFF562LEK",
  "ENCFF137KNW", "ENCFF103BLQ", "ENCFF663ILW",
  #H3K79me2
  "ENCFF892LYD", "ENCFF154HXA", "ENCFF931LYX",
  #H3K27me3
  "ENCFF819QSK", "ENCFF134YLO", "ENCFF336AWS")
```

```{r ENCODE_stats}
ENCODE_stats <- function(matSignalSel,selPeaks,mergedPeaks_annotTXS,DMs) {
allZero=which(apply(matSignalSel,1,sum)==0)
if (length(allZero)>0) matSignalSel=matSignalSel[-allZero,]

mtx=t(matSignalSel)
d=dist(t(scale(mtx)))
hc=hclust(d,method = "ward.D2")
nClust=8
clust=cutree(hc,k=nClust)

selPeaks=intersect(selPeaks,row.names(matSignalSel))

# Complete AnnotationRows
AnnotationRows <- data.frame(
  Localization=factor(mergedPeaks_annotTXS[selPeaks],levels = c("TSS","intraG","interG")),
  CS=factor(clust[selPeaks]),
  row.names = selPeaks
)
AnnotationRows=merge(AnnotationRows,DMs,by.x=0,by.y=1)
row.names(AnnotationRows)=AnnotationRows$Row.names ; AnnotationRows=AnnotationRows[,-1]
colnames(AnnotationRows)[3:4]=c("DM_kind","DM_sign")

return(AnnotationRows)
} 
```

```{r ENCODE_heatmap, fig.height=10, fig.width=8}
ENCODE_heatmap <- function(matSignalSel,AnnotationRows,colsAnnot) {
breaks=seq(-4,4,.1) ; cols=colorRampPalette(c("purple","black","yellow"))(length(breaks)+1)

allZero=which(apply(matSignalSel,1,sum)==0)
if (length(allZero)>0) matSignalSel=matSignalSel[-allZero,]

pheatmap(matSignalSel, 
         scale = "row", cluster_rows = T, cluster_cols = F,
         clustering_method = "ward.D2",
         show_rownames = F, show_colnames = F,
         annotation_col = subset(AnnotationCols,row.names(AnnotationCols) %in% markOrder), 
         annotation_row = AnnotationRows, 
         annotation_colors = colsAnnot,
         color = cols, border = NA,
         #main = paste(mark,"DM peaks",length(DMpeaks)),
         fontsize = 12)
}
```

## All H3K4me2 peaks

```{r ENCODE_counts, fig.height=10, fig.width=8, cache=TRUE}
selPeaksAll=row.names(countNormLength_H3K4me2)
selPeaks=selPeaksAll
matSignalSel=na.omit(matSignal[selPeaks,markOrder])
fileOut=paste0(WDIR,"/results_project12891/scripts/Robjects/AnnotationRows_allPeaks.RData")
if (file.exists(fileOut)) {load(fileOut)} else {
AnnotationRows=ENCODE_stats(matSignalSel,selPeaks,mergedPeaks_annotTXS,DMs)
AnnotationRows_allPeaks=AnnotationRows
save(AnnotationRows_allPeaks,file=fileOut)}

# Tables
AnnotationRows=AnnotationRows_allPeaks
t=table(AnnotationRows$CS)
kable(ceiling(t/sum(t)*100),col.names = c("CS","Freq"))

t1=table(AnnotationRows$Localization)
tt=table(AnnotationRows$CS,AnnotationRows$Localization)
kable(ceiling(t1/sum(t1)*100),col.names = c("Localization","Freq"))
kable(ceiling(tt/as.numeric(t)*100),row.names = T)
# 
# t1=table(AnnotationRows$DM_sign)
# kable(ceiling(t1/sum(t1)*100),col.names = c("H3K4me2 Profile","Freq"))
# tt=table(AnnotationRows$CS,AnnotationRows$DM_sign)
# kable(ceiling(tt/as.numeric(t)*100),row.names = T)
# 
# t1=table(AnnotationRows$DM_kind)
# kable(ceiling(t1/sum(t1)*100),col.names = c("H3K4me2 Dynamic","Freq"))
# tt=table(AnnotationRows$CS,AnnotationRows$DM_kind)
# kable(ceiling(tt/as.numeric(t)*100),row.names = T)

ENCODE_heatmap(matSignalSel,AnnotationRows,colsAnnot)
```

# Multi factorial analysis

Multi factorial analysis (MFA) (Escofier & Pages, Computational Statistics & Data Analysis, 1994) studies several groups of variables (numerical and/or categorical) defined on the same set of individuals. MFA is a factor analysis applied to the array including all groups of variables. Roughly, the behaviour of the method is equivalent to PCA (concerning quantitative variables) or to MCA (concerning qualitative variables).


```{r MFA_functions}
#http://juliejosse.com/wp-content/uploads/2019/03/MFA_staf2015.pdf

individualContributionPlot <- function(mfa.res,nPeaks) {
contribs=data.frame()
for (i in 1:5) {
rs=mfa.res[["ind"]][["contrib"]][order(mfa.res[["ind"]][["contrib"]][,i],decreasing = TRUE),i] %>% cumsum() 
peakId=names(rs)
contribs=rbind(contribs,data.frame(peakId,seq(1,length(rs),1),as.vector(rs),rep(i,length(rs))))}
colnames(contribs)=c("peakId","Rank","RankedSum","Dimension")
contribs$Dimension=factor(contribs$Dimension)

p <- ggplot(contribs,aes(x=Rank,y=RankedSum)) + geom_point(aes(color=Dimension))

print(p)

return(contribs)
}

MFAstats <- function(mfa.res,xlim,ylim) {
# Eigen values of separate analysis
# Describe dimensionality of each variable group
#print(kable(sapply(mfa.res[["separate.analyses"]],"[[", 1),
#            caption = "Eigen values of separate analysis"))
# Eigen values of global analysis
# How many components are necessary to describe the total inertia of the dataset
print(kable(mfa.res[["eig"]][1:10,],
            caption = "Eigen values of global analysis"))
# Correlation coeffs
# How much of the stt defined by component C depends on groups(s) G : G -> C
# Fundamental to see if there are structures COMMON to groups -> useful to integrate 
print(kable(mfa.res[["group"]][["correlation"]],caption = "Group correlation with factor"))
# Inertia or group contribution 
# How much of the inertia of group G is explained by component C : C -> G
print(kable(mfa.res[["group"]][["cos2"]],caption = "Group inertia per factor"))
# Group contribution 
plotGrouprepresentation(mfa.res)
# Invidual contribution to each component
# Determine wheter an axis is due to only some individuals 
#mfa.res[["ind"]][["contrib"]]
#individualContributions=individualContributionPlot(mfa.res,nPeaks)
# Categorical variable representation (centroid of inds )
# Determine how are 
#dimdesc(mfa.res)$Dim.2
#splotIndrepresentation(mfa.res,xlim,ylim)
# Variable contribution 
plotVarrepresentation(mfa.res)
}

plotGrouprepresentation <- function(mfa.res) {
p1 <- plot.MFA(mfa.res,choix = "group",axes = c(1,2), title = "")
p2 <- plot.MFA(mfa.res,choix = "group",axes = c(1,3), title = "")
p3 <- plot.MFA(mfa.res,choix = "group",axes = c(1,4), title = "")
p4 <- plot.MFA(mfa.res,choix = "group",axes = c(1,5), title = "")
print((p1 | p2) / (p3 | p4))
}

plotIndrepresentation <- function(mfa.res,xlim,ylim) {
  if (ylim!="" & xlim!="") {
p1 <- plot.MFA(mfa.res,choix = "ind", invisible = "ind", axes = c(1,2), 
               title = "",ylim=ylim,xlim=xlim)
p2 <- plot.MFA(mfa.res,choix = "ind", invisible = "ind", axes = c(1,3),
               title = "",ylim=ylim,xlim=xlim)
p3 <- plot.MFA(mfa.res,choix = "ind", invisible = "ind", axes = c(1,4),
               title = "",ylim=ylim,xlim=xlim)
p4 <- plot.MFA(mfa.res,choix = "ind", invisible = "ind", axes = c(1,5),           
               title = "",ylim=ylim,xlim=xlim) } else  
{
p1 <- plot.MFA(mfa.res,choix = "ind", invisible = "ind", axes = c(1,2), title = "")
p2 <- plot.MFA(mfa.res,choix = "ind", invisible = "ind", axes = c(1,3), title = "")
p3 <- plot.MFA(mfa.res,choix = "ind", invisible = "ind", axes = c(1,4), title = "")
p4 <- plot.MFA(mfa.res,choix = "ind", invisible = "ind", axes = c(1,5), title = "")
  
}
print((p1 | p2) / (p3 | p4))
}

plotVarrepresentation <- function(mfa.res) {
  title=""
  p1 <- plot.MFA(mfa.res,choix="var", axes=c(1,2),title=title,cex=.7)
  p2 <- plot.MFA(mfa.res,choix="var", axes=c(1,3),title=title,cex=.7)
  p3 <- plot.MFA(mfa.res,choix="var", axes=c(1,4),title=title,cex=.7)
  p4 <- plot.MFA(mfa.res,choix="var", axes=c(1,5),title=title,cex=.7)
  
print((p1 | p2) / (p3 | p4) + plot_layout(guides = "collect"))

}

clusterStats <- function(hcpc.res) {
### Cluster description
## By variables: 
# Test if variable (OMIC samples) values distribution per cluster is random or not 
# For quantitative variables
print(kable(hcpc.res[["desc.var"]]$quanti,caption = "Description by quantitative variables"))
# For qualitative variables
print(kable(hcpc.res[["desc.var"]][["category"]],caption = "Description by qualitative variables"))
## By components 
# Test if individuals coordinates distribution per cluster is random or not 
# For quantitative variables
print(kable(hcpc.res[["desc.axes"]][["quanti"]],
            caption = "Description by factors (quantitative variables)"))
# For qualitative variables
print(kable(hcpc.res[["desc.var"]][["test.chi2"]],
            caption = "Description by factors (qualitative variables)"))
## By individuals (genomic regions)
# Nearest individuals to the cluster's centroid
print("Description by individuals (nearest):")
print(hcpc.res[["desc.ind"]]$para)
# Furthest individuals to the cluster's centroid
print("Description by individuals (furthest):")
print(hcpc.res[["desc.ind"]]$dist)
}

```

```{r MFA_prepareMtx}
# Differential analysis info
AnnotationRows=AnnotationRows_allPeaks[,2:4]
#DEGs2peaks=merge(DEGs,peaks2GeneTgene_H3K4me2[,c(3,1)],by=1)
#AnnotationRows=merge(AnnotationRows,DEGs2peaks[,c(5,4,3,1)],by.x=0,by.y=1)
DEGs2peaks=merge(DEGs,peaks2GeneTgene_H3K4me2,by.x=1,by.y=2)
AnnotationRows=merge(AnnotationRows,DEGs2peaks[,c(5,3,4,1,2,6:8)],by.x=0,by.y=1)
AnnotationRows$Distance=log(abs(AnnotationRows$Distance)+1,10)
#AnnotationRows=AnnotationRows[!duplicated(AnnotationRows[,1]),]
#AnnotationRows <- AnnotationRows[,-c(1,8)] ; 
colnames(AnnotationRows)[5:6] <- c("DEG_kind","DEG_sign")
AnnotationRows$DEG_kind=factor(AnnotationRows$DEG_kind,levels = c("common","only_H3vsUI","only_DIvsD7","notDE"))
#AnnotationRows$DEG_sign=factor(AnnotationRows$DEG_sign,levels = c("Up","Down","UpDown","DownUp"))
AnnotationRows$DEG_sign=factor(AnnotationRows$DEG_sign,levels = c("Up","Down"))
#Re-order columns
AnnotationRows=AnnotationRows[,c(1,7:11,2:6)]

# Merge methylation, transcription and epigenome signal
METHOME=merge(t(scale(t(countNormLength_H3K4me2))),
              DM_results$D7vsUI[,1:2],by=0)
colnames(METHOME)[8:9]=c("H3K4me2_log2baseMean","H3K4me2_D7vsUI_log2FC")
METHOME$H3K4me2_log2baseMean=log(METHOME$H3K4me2_log2baseMean,2)

# Merge TOME data
Data=Data_RMA_NoBatchComBat_Annot
t=na.omit(Data[!duplicated(Data$ENTREZID),])[,c(1,5,9,2,6,10,3,7,11,4,8,12)]
TOME=log(t,2)
row.names(TOME)=na.omit(Data[!duplicated(Data$ENTREZID),14])

# Merge ENCODE data
samples=AnnotationCols[markOrder,]
names(samples)=markOrder
ENCODEmean=t(scale(t(
  data.frame(
  H3K4me3=apply(matSignal[,names(samples[samples=="H3K4me3"])],1,function(x) mean((x))),
  H3K4me2=apply(matSignal[,names(samples[samples=="H3K4me2"])],1,function(x) mean((x))),
  H3K4me1=apply(matSignal[,names(samples[samples=="H3K4me1"])],1,function(x) mean((x))),
  H3K27ac=apply(matSignal[,names(samples[samples=="H3K27ac"])],1,function(x) mean((x))),
  H3K79me2=apply(matSignal[,names(samples[samples=="H3K79me2"])],1,function(x) mean((x))),
  H3K27me3=apply(matSignal[,names(samples[samples=="H3K27me3"])],1,function(x) mean((x))),
  row.names = row.names(matSignal))
)))

# EPITOME all 
EPITOME=merge(AnnotationRows,METHOME,by=1)
EPITOME=merge(EPITOME,ENCODEmean,by.x=1,by.y=0)
EPITOME=merge(EPITOME,TOME,by.x=2,by.y=0)
row.names(EPITOME)=paste0(EPITOME$Row.names,"_",EPITOME$ENTREZID)
colnames(EPITOME)[2]="peakID"
```

```{r MFA_getMtx1}

#### Get INPUT matrix

# Select all peak2gene associations where there is a DMR or a DEG
mfa.EPITOME=subset(EPITOME,EPITOME$DM_kind=="DM" | EPITOME$DEG_kind!="notDE")
corrTh=.4
mfa.EPITOME_corrTh=subset(EPITOME,(EPITOME$DM_kind=="DM" | EPITOME$DEG_kind!="notDE") &
                    abs(EPITOME$Correlation)>corrTh)
distTh=20e3
mfa.EPITOME_distTh=subset(EPITOME,(EPITOME$DM_kind=="DM" | EPITOME$DEG_kind!="notDE") &
                   EPITOME$Distance<log(distTh+1,10))
# Save INPUT matix objects
fileOut=paste0(WDIR,"/results_project12891/scripts/Robjects/mfa.EPITOME.RData")
if (!file.exists(fileOut)) {save(mfa.EPITOME,file = fileOut)}
fileOut=paste0(WDIR,"/results_project12891/scripts/Robjects/mfa.EPITOME_distTh.RData")
if (!file.exists(fileOut)) {save(mfa.EPITOME_distTh,file = fileOut)}
fileOut=paste0(WDIR,"/results_project12891/scripts/Robjects/mfa.EPITOME_corrTh.RData")
if (!file.exists(fileOut)) {save(mfa.EPITOME_corrTh,file = fileOut)}

```

## Input matrix

MFA is performed on the `r dim(mfa.EPITOME)[1]` peak-gene associations that include a DMR or a DEG. This is the distribution of quantitative and qualitative group variables.

```{r MFA_getMtx2}
# Check matrix behaviour
boxplot(mfa.EPITOME[,c(26:37,18:19,4,5,20:25)],las=2,pch=19,cex=.5,
        ylab="Scaled values by group",cex.axis=.8,
        col=c(rep(palette()[8],12),rep(palette()[6],2),
              rep(palette()[3],1),rep(palette()[4],1),rep(palette()[5],6)),
        main="Distribution of quantitative variables for MFA")

kable(table(mfa.EPITOME$DEG_kind,mfa.EPITOME$DEG_sign),caption = "DEGs for MFA")
kable(table(mfa.EPITOME$DM_kind,mfa.EPITOME$DM_sign),caption = "DMRs for MFA")
```

## MFA components description

```{r MFA, fig.height=10, fig.width=10, cache=TRUE}

#### MFA
colInds=c(26:37,18:19,4,5,20:25,7,10:11,8:9)

# fileOut=paste0(WDIR,"/results_project12891/scripts/Robjects/mfa.res_TOMEMETHDISTCORREPIGDEGDMR_CS_distTh.RData")
# if (file.exists(fileOut)) {load(fileOut)} else {
#   mfa.res_TOMEMETHDISTCORREPIGDEGDMR_CS_distTh  <- MFA(mfa.EPITOME_distTh[,colInds], ncp=5, 
#     group=c(12,2,1,1,6,1,2,2), type=c(rep("c",5),rep("n",3)),
#     name.group=c("TOME","METH","DIST","CORR","EPIG","CS","DEG","DMR"),num.group.sup=c(6))
# save(mfa.res_TOMEMETHDISTCORREPIGDEGDMR_CS_distTh,file=fileOut)}
# 
# fileOut=paste0(WDIR,"/results_project12891/scripts/Robjects/mfa.res_TOMEMETHDISTCORREPIGDEGDMR_CS_corrTh.RData")
# if (file.exists(fileOut)) {load(fileOut)} else {
#   mfa.res_TOMEMETHDISTCORREPIGDEGDMR_CS_corrTh  <- MFA(mfa.EPITOME_corrTh[,colInds], ncp=5, 
#     group=c(12,2,1,1,6,1,2,2), type=c(rep("c",5),rep("n",3)),
#     name.group=c("TOME","METH","DIST","CORR","EPIG","CS","DEG","DMR"),num.group.sup=c(6))
# save(mfa.res_TOMEMETHDISTCORREPIGDEGDMR_CS_corrTh,file=fileOut)}

fileOut=paste0(WDIR,"/results_project12891/scripts/Robjects/mfa.res_TOMEMETHDISTCORREPIGDEGDMR_CS.RData")
if (file.exists(fileOut)) {load(fileOut)} else {
  mfa.res_TOMEMETHDISTCORREPIGDEGDMR_CS  <- MFA(mfa.EPITOME[,colInds], ncp=5, 
    group=c(12,2,1,1,6,1,2,2), type=c(rep("c",5),rep("n",3)),
    name.group=c("TOME","METH","DIST","CORR","EPIG","CS","DEG","DMR"),num.group.sup=c(6))
save(mfa.res_TOMEMETHDISTCORREPIGDEGDMR_CS,file=fileOut)}

#### MFA stats
mfa.res=mfa.res_TOMEMETHDISTCORREPIGDEGDMR_CS
title="" ; xlim=c(-1,1) ; ylim=c(-1,1)
MFAstats(mfa.res,"","")
#print(plot.MFA(mfa.res,choix="axes", axes=c(1,2),title=title))
plot.MFA(mfa.res,choix="ind", invisible="ind", axes=c(1,2), ylim=c(-1,2))
p1 <- plot.MFA(mfa.res,choix="ind", invisible="ind", axes=c(1,2), xlim=xlim, ylim=ylim)
p2 <- plot.MFA(mfa.res,choix="ind", invisible="ind", axes=c(3,2), xlim=xlim, ylim=ylim)
p3 <- plot.MFA(mfa.res,choix="ind", invisible="ind", axes=c(3,4), xlim=xlim, ylim=ylim)
p4 <- plot.MFA(mfa.res,choix="ind", invisible="ind", axes=c(5,4), xlim=xlim, ylim=ylim)
print((p1 | p2) / (p3 | p4))
```

# Hierarchical clustering on MFA components

Clustering is performed on the first 5 components obtained with the MFA.

```{r HCPCclustering, fig.height=10, fig.width=10, cache=TRUE}

#### Clustering
nClust=12 ; maxClustSize=1400
hcpc.res <- HCPC(mfa.res, nb.clust=nClust, conso=0, min=3, max=10,
                 metric = "euclidean", method="ward",graph = FALSE,nb.par = maxClustSize)

#clusterStats(hcpc.res)

# Cluster description in term of qualitative variables
lapply(hcpc.res[["desc.var"]][["category"]],function(x) {kable(subset(x,x[,"v.test"]>5.3))})

# Cluster description
tmp=hcpc.res[["data.clust"]][,23:28]
table=merge(mfa.EPITOME[,1:6],tmp,by=0)
colnames(table)[1]="peakIDgene_association"
colnames(table)[13]="MFAclust"

layout(matrix(1:2,nrow = 2))
boxplot(table$Distance ~ table$MFAclust,pch=19,cex=.5,
        xlab="",ylab="Distance peak - gene")
boxplot(table$Correlation ~ table$MFAclust,pch=19,cex=.5,
        xlab="MFA cluster",ylab="Correlation peak_hist-gene_exp")

# Table of associations with clusters
fileOut=paste0(WDIR,"/results_project12891/tables/integration/mfa_clust_TOMEMETHDISTCORREPIGDEGDMR_CS.txt")
if (!file.exists(fileOut)) {
# Add individual data
for (i in 1:nClust) {
  table=merge(table,as.data.frame(hcpc.res[["desc.ind"]][["para"]][[i]]),by.x=1,by.y=0,all.x=TRUE)
  colnames(table)[dim(table)[2]]=paste0("distFromCentroid_MFAclust",i)
}

write.table(table,file = fileOut,sep="\t",quote=FALSE,row.names = FALSE, col.names = TRUE)
}

```

## Statistics on clusters

Number of total genes/peaks per cluster (Ngenes,Npeaks). Number of genes/peaks belonging only to a given cluster (NgeneOnlyClust, NpeakOnlyClust). Nearest cluster in terms of shared genes/peaks (nearestClustGene,nearestClustPeak) and corresponding number of shared genes/peaks (nearestClustNgene,nearestClustNpeak).

```{r}
geneFreq_allAss=table(table$ENTREZID)
singleClustGenes=names(subset(geneFreq_allAss,geneFreq_allAss==1))
peakFreq_allAss=table(table$peakID)
singleClustPeaks=names(subset(peakFreq_allAss,peakFreq_allAss==1))
clusts=split(table,table$MFAclust)
geneClust=unique(table[,c(2,13)])
peakClust=unique(table[,c(3,13)])

clustStats=data.frame(Counts=c("Ngenes","NgeneOnlyClust","nearestClustGene","nearestClustNgene",
                               "Npeaks","NpeaksOnlyClust","nearestClustPeak","nearestClustNpeak")) ; i=1
for (clust in clusts) {
  Ngenes=length(unique(clust$ENTREZID))
  NgeneOnlyClust=length(setdiff(unique(clust$ENTREZID),
                        subset(geneClust,geneClust$MFAclust!=i)$ENTREZID))
  nearestClustGene=names(sort(table(geneClust[geneClust$ENTREZID %in% unique(clust$ENTREZID),]$MFAclust),decreasing = T)[2])
  nearestClustNgene=as.numeric(sort(table(geneClust[geneClust$ENTREZID %in% unique(clust$ENTREZID),]$MFAclust),decreasing = T)[2])
  Npeaks=length(unique(clust$peakID))
  NpeaksOnlyClust=length(setdiff(unique(clust$peakID),
                         subset(peakClust,peakClust$MFAclust!=i)$peakID))
  nearestClustPeak=names(sort(table(peakClust[peakClust$peakID %in% unique(clust$peakID),]$MFAclust),decreasing = T)[2])
  nearestClustNpeak=as.numeric(sort(table(peakClust[peakClust$peakID %in% unique(clust$peakID),]$MFAclust),decreasing = T)[2])
  clustStats=cbind(clustStats,c(Ngenes,NgeneOnlyClust,nearestClustGene,nearestClustNgene,Npeaks,NpeaksOnlyClust,nearestClustPeak,nearestClustNpeak)) ; i=i+1
}
colnames(clustStats)=c("Counts",seq(1,nClust,1))
kable(clustStats)
```


# Explore gene sets from functional analysis 

Check representation of selected gene sets (GS) among genes of every MFA cluster. Values correspond to the ratio between the observed number of genes belonging to a given GS and the expected number per cluster. 

```{r functions_genesetClusterRepresentation, cache=TRUE}
getClustGSenrich <- function(table,selectedGSs,m_t2g_lists) {
  tmp=unique(table[,c(2,13)])
  t=table(tmp$MFAclust)
  expected=t/sum(t)
  overRep=data.frame(matrix(nrow = length(selectedGSs),ncol = nClust))
  colnames(overRep)=seq(1,nClust,1) ; row.names(overRep)=selectedGSs
  for (selectedGS in selectedGSs) {
    ENTREZ_IDs=unique(c(m_t2g_lists[[selectedGS]]))
    subTable=tmp[tmp$ENTREZID %in% ENTREZ_IDs,]
  
    # Clustering distribution
    obs=table(subTable$MFAclust)/sum(table(subTable$MFAclust))
  
    overRep[selectedGS,] <- obs/expected
  }
    return(overRep)

}

m_t2g_df <- msigdbr(species = "Homo sapiens", category = "C2") %>% 
  dplyr::select(gs_subcat, gs_name, entrez_gene)
tmp=m_t2g_df[m_t2g_df$gs_subcat=="CP:REACTOME",]
m_t2g_list_reactome=split(as.character(tmp$entrez_gene),gsub("REACTOME_","",tmp$gs_name))
# tmp=m_t2g_df[m_t2g_df$gs_subcat=="CP:KEGG",]
# m_t2g_list_kegg=split(as.character(tmp$entrez_gene),gsub("KEGG_","",tmp$gs_name))

# m_t2g_lists=c(m_t2g_list_reactome,m_t2g_list_kegg)
m_t2g_lists=m_t2g_list_reactome
```

## GS specific to 1st infection
```{r genesetClusterRepresentation_only1st}
  
selectedGSs_only1st=c("NUCLEAR_RECEPTOR_TRANSCRIPTION_PATHWAY",
              "CHEMOKINE_RECEPTORS_BIND_CHEMOKINES",
              "CHROMATIN_MODIFYING_ENZYMES",
              "PKMTS_METHYLATE_HISTONE_LYSINES",
              "HATS_ACETYLATE_HISTONES",
              "VESICLE_MEDIATED_TRANSPORT",
              "MEMBRANE_TRAFFICKING",
              "DNA_DOUBLE_STRAND_BREAK_RESPONSE",
              "ANTIGEN_PROCESSING_UBIQUITINATION_PROTEASOME_DEGRADATION",
              "PROCESSING_OF_DNA_DOUBLE_STRAND_BREAK_ENDS",
              "DNA_DAMAGE_TELOMERE_STRESS_INDUCED_SENESCENCE")
tab=getClustGSenrich(table,selectedGSs_only1st,m_t2g_lists)
kable(tab,digits = rep(1,nClust))
```              

## GS specific to 2nd infection
```{r genesetClusterRepresentation_only2nd}
selectedGSs_only2nd=c("INFLUENZA_INFECTION",
              "CELLULAR_RESPONSES_TO_EXTERNAL_STIMULI",
              "INFECTIOUS_DISEASE","DISEASE",
              "SIGNALING_BY_ROBO_RECEPTORS","SIGNALING_BY_NUCLEAR_RECEPTORS")
 tab=getClustGSenrich(table,selectedGSs_only2nd,m_t2g_lists)
kable(tab,digits = rep(1,nClust))
```

## GS common to both infections
```{r genesetClusterRepresentation_common}
selectedGSs_common=c("INTERLEUKIN_4_AND_INTERLEUKIN_13_SIGNALING",
              "SIGNALING_BY_INTERLEUKINS","DNA_DOUBLE_STRAND_BREAK_REPAIR","DNA_REPAIR",
              "FOXO_MEDIATED_TRANSCRIPTION",
              "NEGATIVE_REGULATION_OF_MAPK_PATHWAY",
              "CYTOKINE_SIGNALING_IN_IMMUNE_SYSTEM",
              "INTERLEUKIN_10_SIGNALING",
              "HDR_THROUGH_SINGLE_STRAND_ANNEALING_SSA_",
              "FORMATION_OF_SENESCENCE_ASSOCIATED_HETEROCHROMATIN_FOCI_SAHF_")
tab=getClustGSenrich(table,selectedGSs_common,m_t2g_lists)
kable(tab,digits = rep(1,nClust))

```

```{bash searchSpecificGenes, include=FALSE, eval=FALSE}
cd /Volumes/bioit/12891_Chevalier_EpiMemStrep/results_project12891/tables/
# Count genes only 2nd in all clusters
awk '$(NF-2)>.4 {print $2}' transcriptome/Sig_H3vsUI_D7vsUI_DIvsD7_DIvsH3_Annot.txt> id ; 
grep -f id integration/mfa_clust_TOMEMETHDISTCORREPIGDEGDMR_CS.txt|awk '{print $4,$(NF-12)}' | 
sort -u | awk '{print $2}'|sort |uniq -c  

# Find genes only 2nd in cluster 7
awk '$(NF-2)>.4 {print $2}' transcriptome/Sig_H3vsUI_D7vsUI_DIvsD7_DIvsH3_Annot.txt> id
grep -f id integration/mfa_clust_TOMEMETHDISTCORREPIGDEGDMR_CS.txt|
awk '$(NF-12)==7 {print $4}' | sort -u 

rm id

```


# References

```{r}
sessionInfo()
```


```{r MFA_inputMatrixBenchmark, include=FALSE, eval=FALSE}
ass=sample(row.names(mfa.EPITOME_distTh), 500, replace = FALSE, prob = NULL)
xlim=c(-.6,.6) ; ylim=c(-.6,.6)

title="TOME-METH_DIST-CORR-EPIG-CS-DEG-DMR"
colInds=c(26:37,18:19,4,5,20:25,7,10:11,8:9)
mfa.res_TOMEMETH_DISTCORREPIGCSDEGDMR_distTh <- MFA(mfa.EPITOME_distTh[ass,colInds], ncp=5, 
  group=c(12,2,1,1,6,1,2,2), type=c(rep("c",5),rep("n",3)),
  name.group=c("TOME","METH","DIST","CORR","EPIG","CS","DEG","DMR"),num.group.sup=c(3,4,5,6,7,8))

title="TOME-METH-DIST-CORR_EPIG-CS-DEG-DMR"
mfa.res_TOMEMETHDISTCORR_EPIGCSDEGDMR_distTh <- MFA(mfa.EPITOME_distTh[ass,colInds], ncp=5, 
  group=c(12,2,1,1,6,1,2,2), type=c(rep("c",5),rep("n",3)),
  name.group=c("TOME","METH","DIST","CORR","EPIG","CS","DEG","DMR"),
  num.group.sup=c(5,6,7,8))

title="TOME-METH-DIST-CORR-EPIG_CS-DEG-DMR"
mfa.res_TOMEMETHDISTCORREPIG_CSDEGDMR_distTh <- MFA(mfa.EPITOME_distTh[ass,colInds], ncp=5, 
  group=c(12,2,1,1,6,1,2,2), type=c(rep("c",5),rep("n",3)),
  name.group=c("TOME","METH","DIST","CORR","EPIG","CS","DEG","DMR"),num.group.sup=c(6,7,8))

title="TOME-METH-DIST-CORR-EPIG-CS_DEG-DMR"
mfa.res_TOMEMETHDISTCORREPIGCS_DEGDMR_distTh <- MFA(mfa.EPITOME_distTh[ass,colInds], ncp=5, 
  group=c(12,2,1,1,6,1,2,2), type=c(rep("c",5),rep("n",3)),
  name.group=c("TOME","METH","DIST","CORR","EPIG","CS","DEG","DMR"),num.group.sup=c(7,8))

title="TOME-METH-DIST-CORR-EPIG-CS-DEG_DMR"
mfa.res_TOMEMETHDISTCORREPIGCSDEG_DMR_distTh <- MFA(mfa.EPITOME_distTh[ass,colInds], ncp=5, 
  group=c(12,2,1,1,6,1,2,2), type=c(rep("c",5),rep("n",3)),
  name.group=c("TOME","METH","DIST","CORR","EPIG","CS","DEG","DMR"),num.group.sup=c(8))

title="TOME-METH-DIST-CORR-EPIG-CS-DEG-DMR"
mfa.res_TOMEMETHDISTCORREPIGCSDEGDMR_distTh <- MFA(mfa.EPITOME_distTh[ass,colInds], ncp=5, 
  group=c(12,2,1,1,6,1,2,2), type=c(rep("c",5),rep("n",3)),
  name.group=c("TOME","METH","DIST","CORR","EPIG","CS","DEG","DMR"))

title="TOME-METH-DMR-DEG-DIST-CORR-EPIG_CS"
mfa.res_TOMEMETHDISTCORREPIGDEGDMR_CS_distTh  <- MFA(mfa.EPITOME_distTh[ass,colInds], ncp=5, 
  group=c(12,2,1,1,6,1,2,2), type=c(rep("c",5),rep("n",3)),
  name.group=c("TOME","METH","DIST","CORR","EPIG","CS","DEG","DMR"),num.group.sup=c(6))

results=list(TOMEMETH=mfa.res_TOMEMETH_DISTCORREPIGCSDEGDMR_distTh,
     TOMEMETHDISTCORR=mfa.res_TOMEMETHDISTCORR_EPIGCSDEGDMR_distTh,
     TOMEMETHDISTCORREPIG=mfa.res_TOMEMETHDISTCORREPIG_CSDEGDMR_distTh,
     TOMEMETHDISTCORREPIGCS=mfa.res_TOMEMETHDISTCORREPIGCS_DEGDMR_distTh,
     TOMEMETHDISTCORREPIGCSDEG=mfa.res_TOMEMETHDISTCORREPIGCSDEG_DMR_distTh,
     TOMEMETHDISTCORREPIGCSDEGDMR=mfa.res_TOMEMETHDISTCORREPIGCSDEGDMR_distTh,
     TOMEMETHDISTCORREPIGDEGDMR=mfa.res_TOMEMETHDISTCORREPIGDEGDMR_CS_distTh
     )

xlim=c(-1,1) ; ylim=c(-1,1)
i=1
for (mfa.res in results) {
  title=names(results)[i] ; print(title)
  print(kable(mfa.res[["eig"]][1:10,],caption = "Eigen values of global analysis"))
  print(plot.MFA(mfa.res,choix="axes", axes=c(1,2),title=title))
  print(plot.MFA(mfa.res,choix="var", axes=c(1,2),title=title))
  print(plot.MFA(mfa.res,choix="var", axes=c(3,2),title=title))
  print(plot.MFA(mfa.res,choix="var", axes=c(3,4),title=title))
  plot.MFA(mfa.res,choix="ind", invisible="ind", axes=c(1,2))
  plot.MFA(mfa.res,choix="ind", invisible="ind", axes=c(1,2), xlim=xlim, ylim=ylim)
  plot.MFA(mfa.res,choix="ind", invisible="ind", axes=c(3,2), xlim=xlim, ylim=ylim)
  plot.MFA(mfa.res,choix="ind", invisible="ind", axes=c(3,4), xlim=xlim, ylim=ylim)
  i=i+1
}

MFAstats(mfa.res,xlim,ylim)

###################################
###################################
# OTHER INPUT matrices

mfa.EPITOME_notDM_distTh=subset(EPITOME,(EPITOME$DM_kind!="DM" & EPITOME$DEG_kind!="notDE") &
                   EPITOME$Distance<log(distTh+1,10))

fileOut=paste0(WDIR,"/results_project12891/scripts/Robjects/mfa.res_TOMEMETH_distTh.RData")
colInds=c(26:37,18:19,10:11,8:9)
if (file.exists(fileOut)) {load(fileOut)} else {
mfa.res_TOMEMETH_distTh2 <- MFA(mfa.EPITOME_distTh[,colInds], ncp=5, 
    group=c(12,2,2,2), type=c(rep("c",2),rep("n",2)),
    name.group=c("TOME","METH","DEG","DMR"),
    num.group.sup=c(3,4))
save(mfa.res_TOMEMETH_distTh2,file=fileOut)}

fileOut=paste0(WDIR,"/results_project12891/scripts/Robjects/mfa.res_TOMEMETH_EPIG.RData")
colInds=c(26:37,18:19,7,10:11,8:9)
if (file.exists(fileOut)) {load(fileOut)} else {
mfa.res_TOMEMETH_EPIG <- MFA(mfa.EPITOME[,colInds], ncp=5, 
    group=c(12,2,1,2,2), type=c(rep("c",2),rep("n",3)),
    name.group=c("TOME","METH","EPIG","DEG","DMR"),
    num.group.sup=c(3,4,5))
save(mfa.res_TOMEMETH_EPIG,file=fileOut)}

fileOut=paste0(WDIR,"/results_project12891/scripts/Robjects/mfa.res_TOMEMETH_EPIG_distTh.RData")
colInds=c(26:37,18:19,7,10:11,8:9)
if (file.exists(fileOut)) {load(fileOut)} else {
mfa.res_TOMEMETH_EPIG_distTh <- MFA(mfa.EPITOME_distTh[,colInds], ncp=5, 
    group=c(12,2,1,2,2), type=c(rep("c",2),rep("n",3)),
    name.group=c("TOME","METH","EPIG","DEG","DMR"),
    num.group.sup=c(3,4,5))
save(mfa.res_TOMEMETH_EPIG_distTh,file=fileOut)}

fileOut=paste0(WDIR,"/results_project12891/scripts/Robjects/mfa.res_TOMEMETH_EPIG_notDM_distTh.RData")
colInds=c(26:37,18:19,7,10:11,8:9)
if (file.exists(fileOut)) {load(fileOut)} else {
mfa.res_TOMEMETH_EPIG_notDM_distTh <- MFA(mfa.EPITOME_notDM_distTh[,colInds], ncp=5, 
    group=c(12,2,1,2,2), type=c(rep("c",2),rep("n",3)),
    name.group=c("TOME","METH","EPIG","DEG","DMR"),
    num.group.sup=c(3,4,5))
save(mfa.res_TOMEMETH_EPIG_notDM_distTh,file=fileOut)}


```

```{r MFAold, eval=FALSE, include=FALSE}

selPeaks=unique(c(as.character(DMs[DMs$Kind=="DM",]$peakID),selPeaks_DEGs_nearest))
selPeaks_nearest=unique(c(as.character(DMs[DMs$Kind=="DM",]$peakID),
                  intersect(selPeaks_DEGs_nearest,selPeaksAllRep)))
selPeaks_target=unique(c(as.character(DMs[DMs$Kind=="DM",]$peakID),
                  intersect(selPeaks_DEGs_target,selPeaksAllRep)))
nPeaks=length(selPeaks)

selPeaks_nearest=unique(c(as.character(DMs[DMs$Kind=="DM",]$peakID),
                  intersect(selPeaks_DEGs_nearest,selPeaksAllRep)))

#### MFA 
fileOut=paste0(WDIR,"/results_project12891/scripts/Robjects/mfa.res_TOMEMETH_EPIG_nearest.RData")
if (file.exists(fileOut)) {load(fileOut)} else {
mfa.res_TOMEMETH_EPIG_nearest <- MFA(na.omit(EPITOME[selPeaks_nearest,c(7:18,5,6,25:30)]), ncp=5, 
    group=c(12,2,2,2,2), type=c(rep("c",2),rep("n",3)),
    name.group=c("TOME","METH","EPIG","DM","DEG"),
    num.group.sup=c(4,5))
save(mfa.res_TOMEMETH_EPIG_nearest,file=fileOut)}

fileOut=paste0(WDIR,"/results_project12891/scripts/Robjects/mfa.res_TOMEMETH_EPIG_target.RData")
if (file.exists(fileOut)) {load(fileOut)} else {
mfa.res_TOMEMETH_EPIG_target <- MFA(na.omit(EPITOME[selPeaks_target,c(7:18,5,6,25:30)]), ncp=5, 
    group=c(12,2,2,2,2), type=c(rep("c",2),rep("n",3)),
    name.group=c("TOME","METH","EPIG","DM","DEG"),
    num.group.sup=c(4,5))
save(mfa.res_TOMEMETH_EPIG_target,file=fileOut)}

fileOut=paste0(WDIR,"/results_project12891/scripts/Robjects/mfa.res_EPIGMETH_nearest.RData")
if (file.exists(fileOut)) {load(fileOut)} else {
mfa.res_EPIGMETH_nearest <- MFA(na.omit(EPITOME[selPeaks_nearest,c(5,6,19:30)]), ncp=5, 
    group=c(2,6,1,1,2,2), type=c(rep("c",2),rep("n",4)),
    name.group=c("METH","EPIGENOME","LOC","CSkind","DM","DEG"),
    num.group.sup=c(4))
save(mfa.res_EPIGMETH_nearest,file=fileOut)}

fileOut=paste0(WDIR,"/results_project12891/scripts/Robjects/mfa.res_EPIGMETH_target.RData")