In zhezhangsh/DEGandMore:
cat("Setting up and loading data\n");
knitr::opts_chunk$set(eval=TRUE, dpi=300, fig.pos="H", fig.width=8, fig.height=6, echo=FALSE, warning=FALSE, message=FALSE, fig.path='FIGURE/');
library(MASS);
library(gplots);
library(knitr);
library(rmarkdown);
library(yaml);
library(awsomics);
library(rchive);
library(DEGandMore);
#fn.yaml<-'~/R/source/DEGandMore/examples/DdReport/DdReport.yml';
#yml<-yaml.load_file(fn.yaml);
path<-yml$output;
if (!file.exists(path)) dir.create(path, recursive = TRUE);
comp<-names(yml$input$comparison);
# Load results
cat("Loading results from pairwsie DEG analysis ...\n");
res<-lapply(paste(yml$input$comparison, 'result.rds', sep='/'), readRDS);
grp<-lapply(res, function(res) names(res$input$comparison));
smp<-lapply(res, function(x) x$input$comparison);
names(res)<-names(grp)<-names(smp)<-comp;
gset<-readRDS(yml$input$geneset$source);
smpl<-readRDS(yml$input$sample)[unlist(smp, use.names=FALSE), ];
grps<-c(res[[1]][[1]]$comparison, res[[2]][[1]]$comparison);
path.r<-paste(path, 'R', sep='/');
path.tbl<-paste(path, 'TABLE', sep='/');
if (!file.exists(path.r)) dir.create(path.r, recursive = TRUE);
if (!file.exists(path.tbl)) dir.create(path.tbl, recursive = TRUE);
**_[Go back to project home](`r yml$home`)_**
Introduction: This analysis is based on the outputs of pairwise comparison of differential gene expression generated by this template. It uses results from 2 pairwise comparisons of 2 sample groups vs. their corresponding control groups and compares how these 2 sample groups are different from each other in terms of their sample-control differences (delta-delta). An example of such analysis is the different responses of two cell types to the treatment of the same drug.
`r comp[1]` vs. `r comp[2]`
**_[Go back to project home](`r yml$home`)_**
Project
cat('Writing project background\n');
lns<-lapply(names(yml$project), function(nm) {
c(paste('##', nm), '\n', yml$project[[nm]], '\n');
});
lns<-paste(do.call('c', lns), collapse='\n');
r lns
Pairwise comparisons
This report compares the results of the 2 following pairwise comparisons. Click links to view full results of individual comparisons:
**_[Go back to project home](`r yml$home`)_**
Gene-level comparison
de<-lapply(res, function(res) res$de);
stat<-lapply(de, function(de) de$result$stat);
deg<-lapply(de, function(de) de$DEG);
anno<-lapply(res, function(res) res$input$anno);
gid<-lapply(stat, rownames);
if (setequal(rownames(stat[[1]]), rownames(stat[[2]]))) gid<-rownames(stat[[1]]) else
if (yml$input$gene$union) gid<-Reduce('union', lapply(stat, rownames)) else gid<-Reduce('intersect', lapply(stat, rownames))
names(stat)<-names(de)<-names(anno)<-comp;
anno[[2]]<-anno[[2]][!(rownames(anno[[2]]) %in% rownames(anno[[1]])), , drop=FALSE];
anno<-rbind(anno[[1]], anno[[2]])[gid, , drop=FALSE];
l2r<-sapply(stat, function(stat) stat[, 'LogFC'][gid]);
pvl<-sapply(stat, function(stat) stat[, 'Pvalue'][gid]);
fdr<-sapply(stat, function(stat) stat[, 'FDR'][gid]);
rnk<-apply(-sign(l2r)*log10(pvl), 2, rank);
l2r[is.na(l2r)]<-0;
pvl[is.na(pvl)]<-1;
fdr[is.na(fdr)]<-1;
corr<-round(cor(l2r[, 1], l2r[, 2]), 3);
p0<-yml$input$gene$pvalue;
tbl<-lapply(stat, function(s) sapply(colnames(s)[c(1,2,4,5,6)], function(cnm) s[, cnm][gid]));
tbl<-cbind(tbl[[1]], tbl[[2]]);
colnames(tbl)[c(3:5, 8:10)]<-paste(rep(colnames(tbl)[3:5], 2), rep(names(grp), each=3), sep='_');
tbl<-cbind(anno, FormatNumeric(tbl));
CreateDatatable(tbl, paste(path, 'gene_table.html', sep='/'), caption='Gene-level differential expression of 2 pairwise comparisons')->fn;
write.csv(tbl, paste(path.tbl, 'gene_stat.csv', sep='/'));
saveRDS(tbl, paste(path.r, 'gene_stat.rds', sep='/'));
**_[Go back to project home](`r yml$home`)_**
Global delta-delta correlation
Both comparisons reported the log ratio of 2 group means for each gene. The global agreement of log ratios of all genes indicates how much the results of these 2 comparisons are similar to or different from each other. Full table of gene-level statistics side-by-side is here.
Figure 1 This plot shows the global correlation (correlation coefficient = r corr) between the 2 pairwise comparisons: r comp[1] and r comp[2]. DEGs identified by both comparisons are highlighted.
wzxhzdk:3
**_[Go back to project home](`r yml$home`)_**
Differentially expression genes (DEGs)
Both comparisons identified DEGs of 2 compared groups. Overlapped DEGs identified by both comparisons are worthy of a closer look.
Table 1 Number of DEGs identified by both comparisons:
deg.n<-sapply(deg, function(deg) sapply(deg, nrow));
g<-paste(sapply(grp, function(x) x[2]), collapse='/');
rownames(deg.n)<-paste(c('Higher', 'Lower'), 'in', g);
pander::pander(deg.n, align=c('l', 'c', 'c'), );
Figure 2 Overlapping of DEGs. All combinations of differential expression towards opposite directions were plotted and Fisher's exact test was performed to evaluate the significance of overlapping or lack of overlapping.
wzxhzdk:5
Click links to view tables of overlapping DEGs:
r lns
**_[Go back to project home](`r yml$home`)_**
ANOVA
d<-lapply(res, function(res) {
input<-res$input;
s<-c(input$comparison[[1]], input$comparison[[2]]);
sapply(s, function(s) input$expr[, s][gid]);
});
d<-cbind(d[[1]], d[[2]]);
s<-smpl[colnames(d), ]
n<-do.call('c', lapply(de, function(d) c(length(d$group0), length(d$group1))));
for (i in 1:ncol(s)) s[[i]]<-as.factor(s[[i]]);
s$d<-d[1, ];
f<-formula(yml$input$gene$anova$formula);
d0<-lapply(1:nrow(d), function(i) d[i, ]);
paov<- parallel::mclapply(d0, function(x) {
s$d<-x;
smm<-summary(aov(f, data=s))[[1]];
if (length(smm)==1) smm[[1]] else smm;
}, mc.cores=yml$input$gene$anova$core);
cnm<-rownames(paov[[1]]);
cnm<-gsub(' ', '', cnm[-length(cnm)]);
paov<-t(sapply(paov, function(x) x[-nrow(x), 5]));
dimnames(paov)<-list(rownames(d), cnm);
cnm0<-c(yml$input$gene$anova$f1, yml$input$gene$anova$f2, paste(yml$input$gene$anova$f1, yml$input$gene$anova$f2, sep=':'));
paov<-paov[, c(cnm0, setdiff(cnm, cnm0))];
sig.aov<-paov[paov[, 3]<=yml$input$gene$anova$p, , drop=FALSE];
sig.aov<-sig.aov[order(sig.aov[, 3]), ]
#CreateDatatable(cbind(anno[rownames(sig.aov), ], FormatNumeric(sig.aov)), paste(path, 'anova_table.html', sep='/'), caption='ANOVA p value');
m<-sapply(grps, function(g) rowMeans(d[, g, drop=FALSE]));
colnames(paov)<-paste('P', colnames(paov), sep='_');
aov.stat<-cbind(m, paov);
aov.tbl<-data.frame(anno, FormatNumeric(aov.stat), stringsAsFactors = FALSE);
CreateDatatable(aov.tbl, paste(path, 'anova_table.html', sep='/'), caption='ANOVA results');
saveRDS(d, paste(path.r, 'expr.rds', sep='/'));
saveRDS(anno, paste(path.r, 'anno.rds', sep='/'));
saveRDS(aov.stat, paste(path.r, 'anova_stat.rds', sep='/'));
write.csv(d, paste(path.tbl, 'expr.csv', sep='/'));
write.csv(anno, paste(path.tbl, 'anno.csv', sep='/'));
write.csv(aov.stat, paste(path.tbl, 'anova_stat.csv', sep='/'));
2-way ANOVA analysis was performed to identify genes responding to r yml$input$gene$anova$f2 differently in different r yml$input$gene$anova$f1. The analysis reported 3 p values, corresponding to the effect of r yml$input$gene$anova$f2, r yml$input$gene$anova$f1, and their interaction. The analysis identified r nrow(sig.aov) significant genes with interaction p values less than r yml$input$gene$anova$p. The ANOVA results were summarized in a table here.
Figure 3 This is the top gene having the most significant interactive p value (r sig.aov[1, 3]).
wzxhzdk:7
**_[Go back to project home](`r yml$home`)_**
Gene set-level comparison
Genes are often grouped into pre-defined gene sets according to their function, interaction, location, etc. Analysis then can be performed on genes in the same gene set as a unit instead of individual genes.
Gene set average
# Gene set anno
lst<-gset$list;
gns<-unlist(lst, use.names=FALSE);
gst<-rep(names(lst), sapply(lst, length));
lst<-split(gns[gns %in% gid], gst[gns %in% gid]);
sz<-sapply(lst, length);
lst<-lst[sz>=yml$input$geneset$min & sz<=yml$input$geneset$max];
sz<-sapply(lst, length);
# split gene log-ratio into gene sets
gns<-unlist(lst, use.names=FALSE);
gst<-rep(names(lst), sapply(lst, length));
lr1<-l2r[gns, ];
splt1<-apply(lr1, 2, function(lr) split(lr, gst));
mean.gs1<-sapply(splt1, function(x) sapply(x, mean));
lr2<-abs(lr1);
splt2<-apply(lr2, 2, function(lr) split(lr, gst));
mean.gs2<-sapply(splt2, function(x) sapply(x, mean));
# statistics of log-ratio means
splt<-cbind(splt1[[1]], splt1[[2]]);
pprs<-apply(splt, 1, function(x) t.test(x[[1]], x[[2]], paired=TRUE)$p.value[[1]]);
dprs<-mean.gs1[, 2]-mean.gs1[, 1];
# ranked log-ratio and rank sum test
rk1<-rnk[gns, ];
splt3<-apply(rk1, 2, function(x) split(x, gst));
prnk<-apply(cbind(splt3[[1]], splt3[[2]]), 1, function(x) wilcox.test(x[[1]], x[[2]], paired = TRUE)$p.value[[1]]);
rnk.mn<-sapply(splt3, function(x) sapply(x, min));
rnk.mx<-sapply(splt3, function(x) sapply(x, max));
mrnk<-sapply(splt3, function(x) sapply(x, mean));
stat.rnk<-cbind(rnk.mn, rnk.mx);
mean.gs<-list("Mean of log-ratio"=mean.gs1, "Mean of absolute(log-ratio)"=mean.gs2);
stat.gs<-cbind(Size=sapply(splt[, 1], length), mean.gs1, Diff=mean.gs1[,2]-mean.gs1[,1], P_Ttest=pprs, FDR_Ttest=p.adjust(pprs, method='BH'), P_RST=prnk, FDR_RST=p.adjust(prnk, method='BH'));
stat.gs<-stat.gs[order(stat.gs[, 'P_RST']), ];
gs.tbl<-gset[[1]][rownames(stat.gs), ];
gs.tbl$Name<-AddHref(gs.tbl$Name, gs.tbl$URL);
gs.tbl<-cbind(gs.tbl[, 1:3], FormatNumeric(stat.gs));
CreateDatatable(gs.tbl, paste(path, 'geneset_table.html', sep='/'), rownames = FALSE, caption='Gene set average');
saveRDS(gs.tbl, paste(path.r, 'geneset_stat.rds', sep='/'));
write.csv(gs.tbl, paste(path.tbl, 'geneset_stat.csv', sep='/'));
Average differential expression of genes in the same gene set. The gene set-level statistics were fully summarized in this table here.
Figure 4 Each dot represents a gene set and the average log-ratio of all genes in this gene set. The averages were calculated with the log-ratio value of genes (left panel) and the absolute of the log-ratios (right panel). The correlation coefficients are r round(cor(mean.gs[[1]][, 1], mean.gs[[1]][, 2]), 4) and r round(cor(mean.gs[[2]][, 1], mean.gs[[2]][, 2]), 4) respectively.
wzxhzdk:9
**_[Go back to project home](`r yml$home`)_**
path.ora<-paste(path, 'ORA', sep='/');
if (!file.exists(path.ora)) dir.create(path.ora, recursive = TRUE);
ora<-lapply(res, function(x) x$ora);
ora.stat<-lapply(ora, function(x) lapply(x[1:2], function(x) x[[1]]));
ora.sets<-lapply(ora.stat, function(x) lapply(x, rownames));
ora.univ<-lapply(ora, function(x) lapply(x[1:2], function(x) rownames(x[[3]])));
ora.univ<-Reduce('union', c(ora.univ[[1]], ora.univ[[2]]));
# Combine stats
ora.n<-lapply(ora, function(o) sapply(o[1:2], function(o) o[[3]][, 'n11'][ora.univ]));
ora.n<-cbind(ora.n[[1]], ora.n[[2]]);
ora.n[is.na(ora.n)]<-0;
rownames(ora.n)<-ora.univ;
colnames(ora.n)<-paste('N', colnames(ora.n), sep='_');
ora.or<-lapply(ora, function(o) sapply(o[1:2], function(o) {
c<-o[[3]];
c<-(c[, 1]/c[, 2])*(c[, 4]/c[, 3]);
c<-c[ora.univ];
c;
}));
ora.or<-cbind(ora.or[[1]], ora.or[[2]]);
ora.or[is.na(ora.or)]<-1;
rownames(ora.or)<-ora.univ;
colnames(ora.or)<-paste('OR', colnames(ora.or), sep='_');
ora.or<-cbind(ora.or[, c(1, 3)], OR_High_Diff=ora.or[, 3]-ora.or[, 1], ora.or[, c(2, 4)], OR_Low_Diff=ora.or[, 4]-ora.or[, 2]);
ora.p<-lapply(ora.stat, function(o) sapply(o, function(o) o[, 'P_HyperGeo'][ora.univ]));
ora.p<-cbind(ora.p[[1]], ora.p[[2]]);
rownames(ora.p)<-ora.univ;
colnames(ora.p)<-paste('P', colnames(ora.p), sep='_');
ora.p[is.na(ora.p)]<-0.5;
ora.p<-cbind(ora.p[, c(1, 3)], P_High_Ratio=ora.p[, 3]/ora.p[, 1], ora.p[, c(2, 4)], P_Low_ratio=ora.p[, 4]/ora.p[, 2]);
ora.tbl<-FormatNumeric(cbind(ora.n, ora.or, ora.p));
ora.tbl<-cbind(gset[[1]][rownames(ora.tbl), ], ora.tbl);
CreateDatatable(ora.tbl, paste(path, 'ora_table.html', sep='/'), caption='Combined table of ORA statistics');
saveRDS(ora.tbl, paste(path.r, 'ora_stat.rds', sep='/'));
write.csv(ora.tbl, paste(path.tbl, 'ora_stat.csv', sep='/'));
gset.src<-sort(unique(gset[[1]][, 'Source']));
ora.s<-c(ora.stat[[1]], ora.stat[[2]]);
fn.tbl<-lapply(names(ora.s), function(nm) {
s<-ora.s[[nm]];
lapply(gset.src, function(src) {
a<-gset[[1]][rownames(s), ];
a$Name<-AddHref(a$Name, a$URL)
a<-a[, 1:3];
t<-cbind(a[a$Source==src, , drop=FALSE], FormatNumeric(s[a$Source==src, , drop=FALSE]))[, -1];
f<-paste('ORA/', src, '_', nm, '.html', sep='');
CreateDatatable(t, paste(path, f, sep='/'), rownames = FALSE, caption=paste(src, nm, sep=': '));
c<-nrow(t);
names(c)<-f;
c;
});
});
lnk<-sapply(fn.tbl, function(fn) {
c<-as.vector(unlist(fn));
f<-sapply(fn, names);
paste('[', c, '](', f, ')', sep='');
});
lnk<-lnk;
dimnames(lnk)<-list(gset.src, names(ora.s));
Gene set over-representation analysis (ORA)
Each 2-group comparison performs gene set over-representation analysis (ORA) that identifies gene sets over-represented with differentially expressed genes. The results of ORA of both 2-group comparisons are summarized and compared here. The ORA of each gene set reports an odds ratio and p value. These statistics from both comparisons were combined and listed side-by-side, as well as the difference of their odds ratios and ratio of their p values (p set to 0.5 when not available), in this table here
Table 2 Gene sets were broken down into subgroups by their sources. Click on the numbers of over-represented gene sets to see a full list.
r kable(lnk, align='c')
Figure 5 The overlapping of over-represented gene sets from both comparisons.
wzxhzdk:11
Click links to view tables of overlapping significant gene sets from ORA:
r lns
**_[Go back to project home](`r yml$home`)_**
Gene set enrichment analysis (GSEA)
path.gsea<-paste(path, 'GSEA', sep='/');
if (!file.exists(path.gsea)) dir.create(path.gsea, recursive = TRUE);
gsea<-lapply(res, function(x) x$gsea);
gsea.stat<-lapply(gsea, function(x) x$stat);
gsea.stat<-lapply(gsea.stat, function(x) {
rownames(x)<-paste(x[[1]], x[[2]], sep=':');
x;
})
gsea.univ<-Reduce('union', lapply(gsea.stat, rownames));
gsea.stat<-lapply(gsea.stat, function(x) x[gsea.univ, ]);
gsea.tbl<-cbind(gsea.stat[[1]][, 1:3], gsea.stat[[1]][gsea.univ, 4:7], gsea.stat[[2]][gsea.univ, 4:7]);
colnames(gsea.tbl)[c(4:6, 8:10)]<-paste(colnames(gsea.tbl)[c(4:6, 8:10)], rep(names(grp), each=3), sep='_');
gsea.tbl<-FormatNumeric(gsea.tbl);
rownames(gsea.tbl)<-1:nrow(gsea.tbl);
CreateDatatable(gsea.tbl, paste(path, 'gsea_table.html', sep='/'), caption='GSEA statistics');
saveRDS(gsea.tbl, paste(path.r, 'gsea_stat.rds', sep='/'));
write.csv(gsea.tbl, paste(path.tbl, 'gsea_stat.csv', sep='/'));
gsea.st<-lapply(gsea.stat, function(x) split(x[x[, 'PValue']<=0.05, -ncol(x)], x[x[, 'PValue']<=0.05, ncol(x)])[2:1]);
gsea.s<-c(gsea.st[[1]], gsea.st[[2]]);
names(gsea.s)<-names(ora.s);
gset.cll<-sort(unique(unlist(lapply(gsea.s, function(x) x$Collection), use.names=FALSE)));
fn.tbl<-lapply(names(gsea.s), function(nm) {
s<-gsea.s[[nm]];
lapply(gset.cll, function(cll) {
t<-s[s$Collection==cll, , drop=FALSE];
t<-FormatNumeric(t);
f<-paste('GSEA/', cll, '_', nm, '.html', sep='');
CreateDatatable(t, paste(path, f, sep='/'), rownames = FALSE, caption=paste(cll, nm, sep=': '));
c<-nrow(t);
names(c)<-f;
c;
});
});
lnk<-sapply(fn.tbl, function(fn) {
c<-as.vector(unlist(fn));
f<-sapply(fn, names);
paste('[', c, '](', f, ')', sep='');
});
dimnames(lnk)<-list(gset.cll, names(ora.s));
Each 2-group comparison performs gene set enrichment analysis (GSEA) on genes ranked by their differential expression. The results of GSEA of both 2-group comparisons are summarized and compared here. The GSEA of each gene set reports an enrichment score and p value. These statistics from both comparisons were combined and listed side-by-side in this table here
Table 3 Gene sets were broken down into subgroups by collections. Click on the numbers of enriched gene sets to see a full list.
r kable(lnk, align='c')
Figure 6 Nominal enrichment scores from both comparisons. Each dot represents a gene set. Gene sets with p values less than 0.01 from both comparisons are highlighted.
wzxhzdk:13
Figure 7 The overlapping of over-represented gene sets from both comparisons.
wzxhzdk:14
Click links to view tables of overlapping significant gene sets from GSEA:
r lns
**_[Go back to project home](`r yml$home`)_**
Gene clustering
path.cl<-paste(path, 'CLUSTER', sep='/');
if (!file.exists(path.cl)) dir.create(path.cl, recursive = TRUE);
# normalize data
d1<-lapply(res, function(x) {
y<-x$input;
z<-y$expr[, unlist(y$comparison)];
z-rowMeans(z[, y$comparison[[1]]])
});
gs<-Reduce('intersect', lapply(d1, rownames));
d1<-do.call('cbind', lapply(d1, function(x) x[gs, , drop=FALSE]));
d1<-t(apply(d1, 1, function(x) x/sd(x)));
x<-paov[order(paov[, 3]), ];
x<-x[x[,3]<=yml$input$geneset$cluster$panova, , drop=FALSE];
d0<-d1[rownames(x), ];
d0<-d0[1: min(nrow(d0), yml$input$geneset$cluster$top), , drop=FALSE];
# initiate clusters
hc<-hclust(as.dist(1-cor(t(d0))));
cl<-cutree(hc, k=yml$input$geneset$cluster$seed);
cl<-split(names(cl), cl);
# merge similar clusters
flag<-TRUE;
while(flag) {
cat('Number of clusters ', length(cl), '\n');
ms<-sapply(cl, function(cl) colMeans(d0[cl, ]));
tr<-cutree(hclust(as.dist(1-cor(ms))), k=length(cl)-1); # find the 2 most similar clusters
i<-tr[duplicated(tr)];
c<-ms[, tr==i];
r<-cor(c[, 1], c[, 2]);
p<-cor.test(c[, 1], c[, 2])$p.value[[1]];
if (r>yml$input$geneset$cluster$merge$corr & p<yml$input$geneset$cluster$merge$p) {
cl[tr==i][[1]]<-as.vector(unlist(cl[tr==i]));
cl<-cl[names(cl)!=names(i)];
} else flag<-FALSE;
}
# Sort clusters
m<-sapply(cl, function(cl) colMeans(d1[cl, ]));
ind<-apply(m, 2, function(x) which(x==max(x)));
cl<-cl[order(ind)];
names(cl)<-paste('Cluster', 1:length(cl), sep='_');
# re-cluster genes
reCl<-function(d, cl, r, dif) {
md<-sapply(cl, function(cl) apply(d[cl[cl %in% rownames(d)], , drop=FALSE], 2, median));
corr<-cor(t(d), md);
c<-lapply(1:ncol(corr), function(i) {
mx<-apply(corr[, -i, drop=FALSE], 1, max);
rownames(corr)[corr[, i]>=r & (corr[, i]-mx)>dif];
});
c;
}
cls<-reCl(d1, cl, yml$input$geneset$cluster$recluster$corr, yml$input$geneset$cluster$recluster$diff);
ms<-t(sapply(cls, function(x) colMeans(d1[x, ])));
rownames(ms)<-names(cls)<-paste('Cluster', 1:length(cls), sep='_');
fn.htmp<-sapply(names(cls), function(nm) {
x<-d1[cls[[nm]], ];
f<-paste(path.cl, '/Heatmap_', nm, '.pdf', sep='');
sz<-CalculateColoredBlockSize(x);
pdf(f, width = max(sz[2]/3, sz[1]), height = sz[2]);
PlotColoredBlock(x, num.breaks = 31, key = 'Normalized expression', groups = grps);
dev.off();
f;
});
saveRDS(list(cluster=cls, data=d1), file=paste(path.r, 'cluster.rds', sep='/'));
The top r nrow(d0) genes with significant ANOVA p values (p <= 'r yml$input$geneset$cluster$panova') were used as seeds to perform a gene-gene clustering analysis and r length(cls) clusters were identified. ORA was performed on the clusters to identify their functional association (see table below);
fn.tbl<-lapply(names(cls), function(nm) {
cl<-cls[[nm]];
s<-TestGSE(cl, gid, gset[[2]])[[1]];
t<-gset[[1]][rownames(s), ];
t$Name<-AddHref(t$Name, t$URL);
t<-cbind(t[, 1:3], FormatNumeric(s));
f<-paste('CLUSTER/ORA_', nm, '.html', sep='');
CreateDatatable(t, paste(path, f, sep='/'), rownames=FALSE, caption=paste('Over-represented gene set in', nm));
c<-nrow(t);
names(c)<-f;
t<-do.call('cbind', lapply(stat, function(s) s[cl, c(1, 2, 4, 5, 6)]));
t<-FormatNumeric(t);
colnames(t)[c(3,4,5,8,9,10)]<-paste(colnames(t)[c(3,4,5,8,9,10)], rep(names(grp), each=3), sep='_');
t<-data.frame(anno[rownames(t), ], t, stringsAsFactors = FALSE);
CreateDatatable(t, paste(path.cl, '/', nm, '.html', sep=''), caption=nm);
c;
});
mm<-round(sapply(grps, function(x) rowMeans(ms[, x, drop=FALSE])), 4);
n<-as.vector(unlist(fn.tbl));
sz<-sapply(cls, length);
lnk<-paste('[', n, '](', sapply(fn.tbl, names), ')', sep='');
cl.tbl<-data.frame(Size=sz, mm, Gene_set=lnk, stringsAsFactors = FALSE);
cl.tbl[, 1]<-paste('[', sz, '](CLUSTER/', names(cls), '.html)', sep='');
cl.tbl<-cbind(ID=rownames(cl.tbl), cl.tbl);
cl.tbl[, 1]<-paste('[', cl.tbl[, 1], '](CLUSTER/Heatmap_', names(cls), '.pdf)', sep='');
Table 4 This table lists the number of genes in each cluster (click the numbers to see gene lists), the average expression of all genes in a cluster of all sample groups, and then the gene sets over-represented in each cluster (click the numbers to see gene set lists). The gene expression levels were normalized so the mean of the control groups equals to 0 and the mean of the treatment groups is the number of standard deviations.
r kable(cl.tbl, row.names=FALSE, align='c')
Figure 8 This plot shows below the average expression levels of each cluster. Data was normalized before the analysis, so the mean of the control groups was zero and the standard deviation of all samples of each gene was 1.0. Values indicate number of standard deviation from mean of relative control group.
wzxhzdk:17
**_[Go back to project home](`r yml$home`)_**
Figure 9 This plot summarizes the group means and standard errors of all clusters.
wzxhzdk:18
**_[Go back to project home](`r yml$home`)_**
END OF DOCUMENT
zhezhangsh/DEGandMore documentation built on Sept. 22, 2022, 9:55 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
cat("Setting up and loading data\n"); knitr::opts_chunk$set(eval=TRUE, dpi=300, fig.pos="H", fig.width=8, fig.height=6, echo=FALSE, warning=FALSE, message=FALSE, fig.path='FIGURE/'); library(MASS); library(gplots); library(knitr); library(rmarkdown); library(yaml); library(awsomics); library(rchive); library(DEGandMore); #fn.yaml<-'~/R/source/DEGandMore/examples/DdReport/DdReport.yml'; #yml<-yaml.load_file(fn.yaml); path<-yml$output; if (!file.exists(path)) dir.create(path, recursive = TRUE); comp<-names(yml$input$comparison); # Load results cat("Loading results from pairwsie DEG analysis ...\n"); res<-lapply(paste(yml$input$comparison, 'result.rds', sep='/'), readRDS); grp<-lapply(res, function(res) names(res$input$comparison)); smp<-lapply(res, function(x) x$input$comparison); names(res)<-names(grp)<-names(smp)<-comp; gset<-readRDS(yml$input$geneset$source); smpl<-readRDS(yml$input$sample)[unlist(smp, use.names=FALSE), ]; grps<-c(res[[1]][[1]]$comparison, res[[2]][[1]]$comparison); path.r<-paste(path, 'R', sep='/'); path.tbl<-paste(path, 'TABLE', sep='/'); if (!file.exists(path.r)) dir.create(path.r, recursive = TRUE); if (!file.exists(path.tbl)) dir.create(path.tbl, recursive = TRUE);
**_[Go back to project home](`r yml$home`)_**
Introduction: This analysis is based on the outputs of pairwise comparison of differential gene expression generated by this template. It uses results from 2 pairwise comparisons of 2 sample groups vs. their corresponding control groups and compares how these 2 sample groups are different from each other in terms of their sample-control differences (delta-delta). An example of such analysis is the different responses of two cell types to the treatment of the same drug.
`r comp[1]` vs. `r comp[2]`
**_[Go back to project home](`r yml$home`)_**
Project
cat('Writing project background\n'); lns<-lapply(names(yml$project), function(nm) { c(paste('##', nm), '\n', yml$project[[nm]], '\n'); }); lns<-paste(do.call('c', lns), collapse='\n');
r lns
Pairwise comparisons
This report compares the results of the 2 following pairwise comparisons. Click links to view full results of individual comparisons:
**_[Go back to project home](`r yml$home`)_**
Gene-level comparison
de<-lapply(res, function(res) res$de); stat<-lapply(de, function(de) de$result$stat); deg<-lapply(de, function(de) de$DEG); anno<-lapply(res, function(res) res$input$anno); gid<-lapply(stat, rownames); if (setequal(rownames(stat[[1]]), rownames(stat[[2]]))) gid<-rownames(stat[[1]]) else if (yml$input$gene$union) gid<-Reduce('union', lapply(stat, rownames)) else gid<-Reduce('intersect', lapply(stat, rownames)) names(stat)<-names(de)<-names(anno)<-comp; anno[[2]]<-anno[[2]][!(rownames(anno[[2]]) %in% rownames(anno[[1]])), , drop=FALSE]; anno<-rbind(anno[[1]], anno[[2]])[gid, , drop=FALSE]; l2r<-sapply(stat, function(stat) stat[, 'LogFC'][gid]); pvl<-sapply(stat, function(stat) stat[, 'Pvalue'][gid]); fdr<-sapply(stat, function(stat) stat[, 'FDR'][gid]); rnk<-apply(-sign(l2r)*log10(pvl), 2, rank); l2r[is.na(l2r)]<-0; pvl[is.na(pvl)]<-1; fdr[is.na(fdr)]<-1; corr<-round(cor(l2r[, 1], l2r[, 2]), 3); p0<-yml$input$gene$pvalue; tbl<-lapply(stat, function(s) sapply(colnames(s)[c(1,2,4,5,6)], function(cnm) s[, cnm][gid])); tbl<-cbind(tbl[[1]], tbl[[2]]); colnames(tbl)[c(3:5, 8:10)]<-paste(rep(colnames(tbl)[3:5], 2), rep(names(grp), each=3), sep='_'); tbl<-cbind(anno, FormatNumeric(tbl)); CreateDatatable(tbl, paste(path, 'gene_table.html', sep='/'), caption='Gene-level differential expression of 2 pairwise comparisons')->fn; write.csv(tbl, paste(path.tbl, 'gene_stat.csv', sep='/')); saveRDS(tbl, paste(path.r, 'gene_stat.rds', sep='/'));
**_[Go back to project home](`r yml$home`)_**
Global delta-delta correlation
Both comparisons reported the log ratio of 2 group means for each gene. The global agreement of log ratios of all genes indicates how much the results of these 2 comparisons are similar to or different from each other. Full table of gene-level statistics side-by-side is here.
Figure 1 This plot shows the global correlation (correlation coefficient = r corr) between the 2 pairwise comparisons: r comp[1] and r comp[2]. DEGs identified by both comparisons are highlighted.
wzxhzdk:3
**_[Go back to project home](`r yml$home`)_**
Differentially expression genes (DEGs)
Both comparisons identified DEGs of 2 compared groups. Overlapped DEGs identified by both comparisons are worthy of a closer look.
Table 1 Number of DEGs identified by both comparisons:
deg.n<-sapply(deg, function(deg) sapply(deg, nrow)); g<-paste(sapply(grp, function(x) x[2]), collapse='/'); rownames(deg.n)<-paste(c('Higher', 'Lower'), 'in', g); pander::pander(deg.n, align=c('l', 'c', 'c'), );
Figure 2 Overlapping of DEGs. All combinations of differential expression towards opposite directions were plotted and Fisher's exact test was performed to evaluate the significance of overlapping or lack of overlapping.
wzxhzdk:5
Click links to view tables of overlapping DEGs:
r lns
**_[Go back to project home](`r yml$home`)_**
ANOVA
d<-lapply(res, function(res) { input<-res$input; s<-c(input$comparison[[1]], input$comparison[[2]]); sapply(s, function(s) input$expr[, s][gid]); }); d<-cbind(d[[1]], d[[2]]); s<-smpl[colnames(d), ] n<-do.call('c', lapply(de, function(d) c(length(d$group0), length(d$group1)))); for (i in 1:ncol(s)) s[[i]]<-as.factor(s[[i]]); s$d<-d[1, ]; f<-formula(yml$input$gene$anova$formula); d0<-lapply(1:nrow(d), function(i) d[i, ]); paov<- parallel::mclapply(d0, function(x) { s$d<-x; smm<-summary(aov(f, data=s))[[1]]; if (length(smm)==1) smm[[1]] else smm; }, mc.cores=yml$input$gene$anova$core); cnm<-rownames(paov[[1]]); cnm<-gsub(' ', '', cnm[-length(cnm)]); paov<-t(sapply(paov, function(x) x[-nrow(x), 5])); dimnames(paov)<-list(rownames(d), cnm); cnm0<-c(yml$input$gene$anova$f1, yml$input$gene$anova$f2, paste(yml$input$gene$anova$f1, yml$input$gene$anova$f2, sep=':')); paov<-paov[, c(cnm0, setdiff(cnm, cnm0))]; sig.aov<-paov[paov[, 3]<=yml$input$gene$anova$p, , drop=FALSE]; sig.aov<-sig.aov[order(sig.aov[, 3]), ] #CreateDatatable(cbind(anno[rownames(sig.aov), ], FormatNumeric(sig.aov)), paste(path, 'anova_table.html', sep='/'), caption='ANOVA p value'); m<-sapply(grps, function(g) rowMeans(d[, g, drop=FALSE])); colnames(paov)<-paste('P', colnames(paov), sep='_'); aov.stat<-cbind(m, paov); aov.tbl<-data.frame(anno, FormatNumeric(aov.stat), stringsAsFactors = FALSE); CreateDatatable(aov.tbl, paste(path, 'anova_table.html', sep='/'), caption='ANOVA results'); saveRDS(d, paste(path.r, 'expr.rds', sep='/')); saveRDS(anno, paste(path.r, 'anno.rds', sep='/')); saveRDS(aov.stat, paste(path.r, 'anova_stat.rds', sep='/')); write.csv(d, paste(path.tbl, 'expr.csv', sep='/')); write.csv(anno, paste(path.tbl, 'anno.csv', sep='/')); write.csv(aov.stat, paste(path.tbl, 'anova_stat.csv', sep='/'));
2-way ANOVA analysis was performed to identify genes responding to r yml$input$gene$anova$f2 differently in different r yml$input$gene$anova$f1. The analysis reported 3 p values, corresponding to the effect of r yml$input$gene$anova$f2, r yml$input$gene$anova$f1, and their interaction. The analysis identified r nrow(sig.aov) significant genes with interaction p values less than r yml$input$gene$anova$p. The ANOVA results were summarized in a table here.
Figure 3 This is the top gene having the most significant interactive p value (r sig.aov[1, 3]).
wzxhzdk:7
**_[Go back to project home](`r yml$home`)_**
Gene set-level comparison
Genes are often grouped into pre-defined gene sets according to their function, interaction, location, etc. Analysis then can be performed on genes in the same gene set as a unit instead of individual genes.
Gene set average
# Gene set anno lst<-gset$list; gns<-unlist(lst, use.names=FALSE); gst<-rep(names(lst), sapply(lst, length)); lst<-split(gns[gns %in% gid], gst[gns %in% gid]); sz<-sapply(lst, length); lst<-lst[sz>=yml$input$geneset$min & sz<=yml$input$geneset$max]; sz<-sapply(lst, length); # split gene log-ratio into gene sets gns<-unlist(lst, use.names=FALSE); gst<-rep(names(lst), sapply(lst, length)); lr1<-l2r[gns, ]; splt1<-apply(lr1, 2, function(lr) split(lr, gst)); mean.gs1<-sapply(splt1, function(x) sapply(x, mean)); lr2<-abs(lr1); splt2<-apply(lr2, 2, function(lr) split(lr, gst)); mean.gs2<-sapply(splt2, function(x) sapply(x, mean)); # statistics of log-ratio means splt<-cbind(splt1[[1]], splt1[[2]]); pprs<-apply(splt, 1, function(x) t.test(x[[1]], x[[2]], paired=TRUE)$p.value[[1]]); dprs<-mean.gs1[, 2]-mean.gs1[, 1]; # ranked log-ratio and rank sum test rk1<-rnk[gns, ]; splt3<-apply(rk1, 2, function(x) split(x, gst)); prnk<-apply(cbind(splt3[[1]], splt3[[2]]), 1, function(x) wilcox.test(x[[1]], x[[2]], paired = TRUE)$p.value[[1]]); rnk.mn<-sapply(splt3, function(x) sapply(x, min)); rnk.mx<-sapply(splt3, function(x) sapply(x, max)); mrnk<-sapply(splt3, function(x) sapply(x, mean)); stat.rnk<-cbind(rnk.mn, rnk.mx); mean.gs<-list("Mean of log-ratio"=mean.gs1, "Mean of absolute(log-ratio)"=mean.gs2); stat.gs<-cbind(Size=sapply(splt[, 1], length), mean.gs1, Diff=mean.gs1[,2]-mean.gs1[,1], P_Ttest=pprs, FDR_Ttest=p.adjust(pprs, method='BH'), P_RST=prnk, FDR_RST=p.adjust(prnk, method='BH')); stat.gs<-stat.gs[order(stat.gs[, 'P_RST']), ]; gs.tbl<-gset[[1]][rownames(stat.gs), ]; gs.tbl$Name<-AddHref(gs.tbl$Name, gs.tbl$URL); gs.tbl<-cbind(gs.tbl[, 1:3], FormatNumeric(stat.gs)); CreateDatatable(gs.tbl, paste(path, 'geneset_table.html', sep='/'), rownames = FALSE, caption='Gene set average'); saveRDS(gs.tbl, paste(path.r, 'geneset_stat.rds', sep='/')); write.csv(gs.tbl, paste(path.tbl, 'geneset_stat.csv', sep='/'));
Average differential expression of genes in the same gene set. The gene set-level statistics were fully summarized in this table here.
Figure 4 Each dot represents a gene set and the average log-ratio of all genes in this gene set. The averages were calculated with the log-ratio value of genes (left panel) and the absolute of the log-ratios (right panel). The correlation coefficients are r round(cor(mean.gs[[1]][, 1], mean.gs[[1]][, 2]), 4) and r round(cor(mean.gs[[2]][, 1], mean.gs[[2]][, 2]), 4) respectively.
wzxhzdk:9
**_[Go back to project home](`r yml$home`)_**
path.ora<-paste(path, 'ORA', sep='/'); if (!file.exists(path.ora)) dir.create(path.ora, recursive = TRUE); ora<-lapply(res, function(x) x$ora); ora.stat<-lapply(ora, function(x) lapply(x[1:2], function(x) x[[1]])); ora.sets<-lapply(ora.stat, function(x) lapply(x, rownames)); ora.univ<-lapply(ora, function(x) lapply(x[1:2], function(x) rownames(x[[3]]))); ora.univ<-Reduce('union', c(ora.univ[[1]], ora.univ[[2]])); # Combine stats ora.n<-lapply(ora, function(o) sapply(o[1:2], function(o) o[[3]][, 'n11'][ora.univ])); ora.n<-cbind(ora.n[[1]], ora.n[[2]]); ora.n[is.na(ora.n)]<-0; rownames(ora.n)<-ora.univ; colnames(ora.n)<-paste('N', colnames(ora.n), sep='_'); ora.or<-lapply(ora, function(o) sapply(o[1:2], function(o) { c<-o[[3]]; c<-(c[, 1]/c[, 2])*(c[, 4]/c[, 3]); c<-c[ora.univ]; c; })); ora.or<-cbind(ora.or[[1]], ora.or[[2]]); ora.or[is.na(ora.or)]<-1; rownames(ora.or)<-ora.univ; colnames(ora.or)<-paste('OR', colnames(ora.or), sep='_'); ora.or<-cbind(ora.or[, c(1, 3)], OR_High_Diff=ora.or[, 3]-ora.or[, 1], ora.or[, c(2, 4)], OR_Low_Diff=ora.or[, 4]-ora.or[, 2]); ora.p<-lapply(ora.stat, function(o) sapply(o, function(o) o[, 'P_HyperGeo'][ora.univ])); ora.p<-cbind(ora.p[[1]], ora.p[[2]]); rownames(ora.p)<-ora.univ; colnames(ora.p)<-paste('P', colnames(ora.p), sep='_'); ora.p[is.na(ora.p)]<-0.5; ora.p<-cbind(ora.p[, c(1, 3)], P_High_Ratio=ora.p[, 3]/ora.p[, 1], ora.p[, c(2, 4)], P_Low_ratio=ora.p[, 4]/ora.p[, 2]); ora.tbl<-FormatNumeric(cbind(ora.n, ora.or, ora.p)); ora.tbl<-cbind(gset[[1]][rownames(ora.tbl), ], ora.tbl); CreateDatatable(ora.tbl, paste(path, 'ora_table.html', sep='/'), caption='Combined table of ORA statistics'); saveRDS(ora.tbl, paste(path.r, 'ora_stat.rds', sep='/')); write.csv(ora.tbl, paste(path.tbl, 'ora_stat.csv', sep='/')); gset.src<-sort(unique(gset[[1]][, 'Source'])); ora.s<-c(ora.stat[[1]], ora.stat[[2]]); fn.tbl<-lapply(names(ora.s), function(nm) { s<-ora.s[[nm]]; lapply(gset.src, function(src) { a<-gset[[1]][rownames(s), ]; a$Name<-AddHref(a$Name, a$URL) a<-a[, 1:3]; t<-cbind(a[a$Source==src, , drop=FALSE], FormatNumeric(s[a$Source==src, , drop=FALSE]))[, -1]; f<-paste('ORA/', src, '_', nm, '.html', sep=''); CreateDatatable(t, paste(path, f, sep='/'), rownames = FALSE, caption=paste(src, nm, sep=': ')); c<-nrow(t); names(c)<-f; c; }); }); lnk<-sapply(fn.tbl, function(fn) { c<-as.vector(unlist(fn)); f<-sapply(fn, names); paste('[', c, '](', f, ')', sep=''); }); lnk<-lnk; dimnames(lnk)<-list(gset.src, names(ora.s));
Gene set over-representation analysis (ORA)
Each 2-group comparison performs gene set over-representation analysis (ORA) that identifies gene sets over-represented with differentially expressed genes. The results of ORA of both 2-group comparisons are summarized and compared here. The ORA of each gene set reports an odds ratio and p value. These statistics from both comparisons were combined and listed side-by-side, as well as the difference of their odds ratios and ratio of their p values (p set to 0.5 when not available), in this table here
Table 2 Gene sets were broken down into subgroups by their sources. Click on the numbers of over-represented gene sets to see a full list.
r kable(lnk, align='c')
Figure 5 The overlapping of over-represented gene sets from both comparisons.
wzxhzdk:11
Click links to view tables of overlapping significant gene sets from ORA:
r lns
**_[Go back to project home](`r yml$home`)_**
Gene set enrichment analysis (GSEA)
path.gsea<-paste(path, 'GSEA', sep='/'); if (!file.exists(path.gsea)) dir.create(path.gsea, recursive = TRUE); gsea<-lapply(res, function(x) x$gsea); gsea.stat<-lapply(gsea, function(x) x$stat); gsea.stat<-lapply(gsea.stat, function(x) { rownames(x)<-paste(x[[1]], x[[2]], sep=':'); x; }) gsea.univ<-Reduce('union', lapply(gsea.stat, rownames)); gsea.stat<-lapply(gsea.stat, function(x) x[gsea.univ, ]); gsea.tbl<-cbind(gsea.stat[[1]][, 1:3], gsea.stat[[1]][gsea.univ, 4:7], gsea.stat[[2]][gsea.univ, 4:7]); colnames(gsea.tbl)[c(4:6, 8:10)]<-paste(colnames(gsea.tbl)[c(4:6, 8:10)], rep(names(grp), each=3), sep='_'); gsea.tbl<-FormatNumeric(gsea.tbl); rownames(gsea.tbl)<-1:nrow(gsea.tbl); CreateDatatable(gsea.tbl, paste(path, 'gsea_table.html', sep='/'), caption='GSEA statistics'); saveRDS(gsea.tbl, paste(path.r, 'gsea_stat.rds', sep='/')); write.csv(gsea.tbl, paste(path.tbl, 'gsea_stat.csv', sep='/')); gsea.st<-lapply(gsea.stat, function(x) split(x[x[, 'PValue']<=0.05, -ncol(x)], x[x[, 'PValue']<=0.05, ncol(x)])[2:1]); gsea.s<-c(gsea.st[[1]], gsea.st[[2]]); names(gsea.s)<-names(ora.s); gset.cll<-sort(unique(unlist(lapply(gsea.s, function(x) x$Collection), use.names=FALSE))); fn.tbl<-lapply(names(gsea.s), function(nm) { s<-gsea.s[[nm]]; lapply(gset.cll, function(cll) { t<-s[s$Collection==cll, , drop=FALSE]; t<-FormatNumeric(t); f<-paste('GSEA/', cll, '_', nm, '.html', sep=''); CreateDatatable(t, paste(path, f, sep='/'), rownames = FALSE, caption=paste(cll, nm, sep=': ')); c<-nrow(t); names(c)<-f; c; }); }); lnk<-sapply(fn.tbl, function(fn) { c<-as.vector(unlist(fn)); f<-sapply(fn, names); paste('[', c, '](', f, ')', sep=''); }); dimnames(lnk)<-list(gset.cll, names(ora.s));
Each 2-group comparison performs gene set enrichment analysis (GSEA) on genes ranked by their differential expression. The results of GSEA of both 2-group comparisons are summarized and compared here. The GSEA of each gene set reports an enrichment score and p value. These statistics from both comparisons were combined and listed side-by-side in this table here
Table 3 Gene sets were broken down into subgroups by collections. Click on the numbers of enriched gene sets to see a full list.
r kable(lnk, align='c')
Figure 6 Nominal enrichment scores from both comparisons. Each dot represents a gene set. Gene sets with p values less than 0.01 from both comparisons are highlighted.
wzxhzdk:13
Figure 7 The overlapping of over-represented gene sets from both comparisons.
wzxhzdk:14
Click links to view tables of overlapping significant gene sets from GSEA:
r lns
**_[Go back to project home](`r yml$home`)_**
Gene clustering
path.cl<-paste(path, 'CLUSTER', sep='/'); if (!file.exists(path.cl)) dir.create(path.cl, recursive = TRUE); # normalize data d1<-lapply(res, function(x) { y<-x$input; z<-y$expr[, unlist(y$comparison)]; z-rowMeans(z[, y$comparison[[1]]]) }); gs<-Reduce('intersect', lapply(d1, rownames)); d1<-do.call('cbind', lapply(d1, function(x) x[gs, , drop=FALSE])); d1<-t(apply(d1, 1, function(x) x/sd(x))); x<-paov[order(paov[, 3]), ]; x<-x[x[,3]<=yml$input$geneset$cluster$panova, , drop=FALSE]; d0<-d1[rownames(x), ]; d0<-d0[1: min(nrow(d0), yml$input$geneset$cluster$top), , drop=FALSE]; # initiate clusters hc<-hclust(as.dist(1-cor(t(d0)))); cl<-cutree(hc, k=yml$input$geneset$cluster$seed); cl<-split(names(cl), cl); # merge similar clusters flag<-TRUE; while(flag) { cat('Number of clusters ', length(cl), '\n'); ms<-sapply(cl, function(cl) colMeans(d0[cl, ])); tr<-cutree(hclust(as.dist(1-cor(ms))), k=length(cl)-1); # find the 2 most similar clusters i<-tr[duplicated(tr)]; c<-ms[, tr==i]; r<-cor(c[, 1], c[, 2]); p<-cor.test(c[, 1], c[, 2])$p.value[[1]]; if (r>yml$input$geneset$cluster$merge$corr & p<yml$input$geneset$cluster$merge$p) { cl[tr==i][[1]]<-as.vector(unlist(cl[tr==i])); cl<-cl[names(cl)!=names(i)]; } else flag<-FALSE; } # Sort clusters m<-sapply(cl, function(cl) colMeans(d1[cl, ])); ind<-apply(m, 2, function(x) which(x==max(x))); cl<-cl[order(ind)]; names(cl)<-paste('Cluster', 1:length(cl), sep='_'); # re-cluster genes reCl<-function(d, cl, r, dif) { md<-sapply(cl, function(cl) apply(d[cl[cl %in% rownames(d)], , drop=FALSE], 2, median)); corr<-cor(t(d), md); c<-lapply(1:ncol(corr), function(i) { mx<-apply(corr[, -i, drop=FALSE], 1, max); rownames(corr)[corr[, i]>=r & (corr[, i]-mx)>dif]; }); c; } cls<-reCl(d1, cl, yml$input$geneset$cluster$recluster$corr, yml$input$geneset$cluster$recluster$diff); ms<-t(sapply(cls, function(x) colMeans(d1[x, ]))); rownames(ms)<-names(cls)<-paste('Cluster', 1:length(cls), sep='_'); fn.htmp<-sapply(names(cls), function(nm) { x<-d1[cls[[nm]], ]; f<-paste(path.cl, '/Heatmap_', nm, '.pdf', sep=''); sz<-CalculateColoredBlockSize(x); pdf(f, width = max(sz[2]/3, sz[1]), height = sz[2]); PlotColoredBlock(x, num.breaks = 31, key = 'Normalized expression', groups = grps); dev.off(); f; }); saveRDS(list(cluster=cls, data=d1), file=paste(path.r, 'cluster.rds', sep='/'));
The top r nrow(d0) genes with significant ANOVA p values (p <= 'r yml$input$geneset$cluster$panova') were used as seeds to perform a gene-gene clustering analysis and r length(cls) clusters were identified. ORA was performed on the clusters to identify their functional association (see table below);
fn.tbl<-lapply(names(cls), function(nm) { cl<-cls[[nm]]; s<-TestGSE(cl, gid, gset[[2]])[[1]]; t<-gset[[1]][rownames(s), ]; t$Name<-AddHref(t$Name, t$URL); t<-cbind(t[, 1:3], FormatNumeric(s)); f<-paste('CLUSTER/ORA_', nm, '.html', sep=''); CreateDatatable(t, paste(path, f, sep='/'), rownames=FALSE, caption=paste('Over-represented gene set in', nm)); c<-nrow(t); names(c)<-f; t<-do.call('cbind', lapply(stat, function(s) s[cl, c(1, 2, 4, 5, 6)])); t<-FormatNumeric(t); colnames(t)[c(3,4,5,8,9,10)]<-paste(colnames(t)[c(3,4,5,8,9,10)], rep(names(grp), each=3), sep='_'); t<-data.frame(anno[rownames(t), ], t, stringsAsFactors = FALSE); CreateDatatable(t, paste(path.cl, '/', nm, '.html', sep=''), caption=nm); c; }); mm<-round(sapply(grps, function(x) rowMeans(ms[, x, drop=FALSE])), 4); n<-as.vector(unlist(fn.tbl)); sz<-sapply(cls, length); lnk<-paste('[', n, '](', sapply(fn.tbl, names), ')', sep=''); cl.tbl<-data.frame(Size=sz, mm, Gene_set=lnk, stringsAsFactors = FALSE); cl.tbl[, 1]<-paste('[', sz, '](CLUSTER/', names(cls), '.html)', sep=''); cl.tbl<-cbind(ID=rownames(cl.tbl), cl.tbl); cl.tbl[, 1]<-paste('[', cl.tbl[, 1], '](CLUSTER/Heatmap_', names(cls), '.pdf)', sep='');
Table 4 This table lists the number of genes in each cluster (click the numbers to see gene lists), the average expression of all genes in a cluster of all sample groups, and then the gene sets over-represented in each cluster (click the numbers to see gene set lists). The gene expression levels were normalized so the mean of the control groups equals to 0 and the mean of the treatment groups is the number of standard deviations.
r kable(cl.tbl, row.names=FALSE, align='c')
Figure 8 This plot shows below the average expression levels of each cluster. Data was normalized before the analysis, so the mean of the control groups was zero and the standard deviation of all samples of each gene was 1.0. Values indicate number of standard deviation from mean of relative control group.
wzxhzdk:17
**_[Go back to project home](`r yml$home`)_**
Figure 9 This plot summarizes the group means and standard errors of all clusters.
wzxhzdk:18
**_[Go back to project home](`r yml$home`)_**
END OF DOCUMENT
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.