In zhezhangsh/DEGandMore:

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`

**_[Go back to parent page](`r home`)_**

This procedure is usually run as a subroutine of 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 to view gene list and data of the cluster.

Visualization

**_[Go back to parent page](`r home`)_**

plot(hclust(as.dist(1-cor(mtrx))), main=name, xlab='Sample', sub='');

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

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

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

**_[Go back to parent page](`r home`)_**

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

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');