---
title: "Visualization and gene set enrichment analysis of a gene cluster"
author: "Jim Zhang"
date: "`r Sys.Date()`"
output:
  html_document:
    fig_caption: yes
    toc: yes
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(dpi=300, fig.pos="H", fig.width=8, fig.height=6, dev=c('png', 'pdf'), echo=FALSE, warning=FALSE, message=FALSE);

library(DEGandMore);
library(awsomics);
library(gplots);
library(knitr);

if (!exists('name')) stop("cluster name not given\n"); # the name of the cluster
if (!exists('path')) stop("output file path not given\n"); # the path to output files
if (!exists('home')) stop("path to report home not given\n"); # the path to report home page
if (!exists('mtrx')) stop("data matrix not given\n"); # the data matrix of gene clusters
if (!exists('anno')) stop("gene annotation not given\n"); # the annotation of each row in the data matrix
if (!exists('grps')) stop("sample groups not given\n"); # the sample groups corresponding to the data matrix columns
if (!exists('univ')) stop("background genes not given\n"); # the background genes for enrichment analysis

tbl<-cbind(anno[rownames(mtrx), , drop=FALSE], round(mtrx, 4)); 
CreateDatatable(tbl, paste(path, 'data.html', sep='/')); 

#if (!file.exists(path)) dir.create(path, recursive = TRUE);
```

# `r name`

<div align="right">**_[Go back to parent page](`r home`)_**</div>

_This procedure is usually run as a subroutine of [ClReport.Rmd](https://raw.githubusercontent.com/zhezhangsh/DEGandMore/master/examples/MultiGroupCluster/ClReport.Rmd) procedure._

  - Total number of genes **`r length(univ)`**
  - Number of genes in cluster **`r nrow(mtrx)`**
  - Total number of samples **`r ncol(mtrx)`**
  - Number of sample groups **`r length(grps)`**

Click [here](data.html) to view gene list and data of the cluster.

## Visualization

<div align="right">**_[Go back to parent page](`r home`)_**</div>

```{r clustering_1, include=TRUE, out.width="600px",}
plot(hclust(as.dist(1-cor(mtrx))), main=name, xlab='Sample', sub=''); 
```

**Figure 1** Hierarchical clustering of samples using all genes of the cluster

```{r heatmap_1, include=TRUE, out.width="800px", fig.height=8, fig.width=max(6, 0.05*nrow(mtrx))}
sortColumn<-function(ms, g) {
  corr<-as.vector(cor(rowMeans(ms[, g[[2]], drop=FALSE]), ms[, g[[1]], drop=FALSE]));
  g[[1]]<-g[[1]][order(corr)];
  for (i in 2:length(g)) {
    corr<-as.vector(cor(rowMeans(ms[, g[[i-1]], drop=FALSE]), ms[, grp[[i]], drop=FALSE]));
    g[[i]]<-g[[i]][rev(order(corr))]; 
  }
  ms[, as.vector(unlist(g))]; 
}
d<-sortColumn(mtrx, grps); 
#rownames(d)<-paste(rownames(d), anno[rownames(d), 1], sep=':'); 
PlotColoredBlock(d, min=-1*max(abs(d), na.rm=TRUE), max=max(abs(d), na.rm=TRUE), num.breaks = 127, groups=grp, key='Expression level'); 
abline(v=c(0, cumsum(sapply(grp, length))), lwd=1); 
```

**Figure 2** Heatmap using all genes and samples

```{r series_1, include=TRUE, out.width='750px'}
ms<-colMeans(mtrx); 
m<-sapply(grps, function(g) mean(ms[g])); 
se<-sapply(grps, function(g) sd(ms[g])/sqrt(length(g))); 
m<-matrix(m, nr=1); 
se<-matrix(se, nr=1); 
colnames(m)<-names(grps); 
PlotSeries(m, se, labs=c('', 'Average expression'), title=name, draw.legend = FALSE); 
```

**Figure 3** Mean and standard errors across groups

## Gene set enrichment analysis

<div align="right">**_[Go back to parent page](`r home`)_**</div>

Find predefined gene sets enriched in gene cluster comparing to the background. 

```{r gse_1, include=FALSE}
gse<- TestGSE(rownames(mtrx), univ, gset[[2]])[[1]]; 
gse<-gse[, c(2, 4, 5, 6)]; 
colnames(gse)<-c('N', 'Enrichment', 'P', 'FDR'); 
gse<-cbind(gset[[1]][rownames(gse), ], gse); 
saveRDS(gse, file=paste(path, 'GSE.rds', sep='/')); 
write.csv(gse, file=paste(path, 'GSE.csv', sep='/')); 
gse$Name<- AddHref(gse$Name, gse$URL); 
gse<-gse[, colnames(gse)!='URL'];  

gse<-split(gse[, -1], gse[, 1]); 
fn<-sapply(names(gse), function(nm) CreateDatatable(gse[[nm]], paste(path, paste('GSE_', nm, '.html', sep=''), sep='/'), caption=nm)); 

lns<-paste('  - [', names(gse), '](', paste('GSE_', names(gse), '.html', sep=''), ')', sep=''); 
lns<-paste(lns, collapse='\n'); 
```

Click links to view full tabls of gene sets from different sources.

`r lns`

<div align="right">**_[Go back to parent page](`r home`)_**</div>

***
_END OF DOCUMENT_