ExampleStats.R

# plotfchheatmap function is a wrap around heatmap.2 from gplots package to plot fold changes between groups
library(gplots)
plotfchheatmap<-function(M, CNM, col=colorRampPalette(c("green","white","orange"))(21),Colv = FALSE,dendrogram="row",margins=c(7,15),...)
{
  heatmap.2(as.matrix(M), cellnote=CNM, notecol="black", col=col,trace="none", margins = margins,Colv = Colv,dendrogram=dendrogram, ...)
}#end plotfchheatmap

# findmetext function facilitates search for text string
findmetext<-function(x,y)
{
  return(which(regexpr(x,y)>-1))
}

# statistical testing of differences between groups
data(iris)

#making one variable not significant
iris[,1]<-rnorm(n = nrow(iris),mean = 5,sd = 0.1)
tmp<-iris[,1:4]
grp<-iris$Species
fdrthr<-0.1


kruskalres<-apply(tmp,2,function(x) kruskal.test(x,g=grp))

padjres<-p.adjust(sapply(kruskalres,function(x) x$p.value), method="BH")

library(FSA)
#dunn tests only on significant
dunnTestres<-apply(rbind(tmp[,names(which(padjres<fdrthr))]),2,function(x) dunnTest(x~grp))
names(dunnTestres)<-names(which(padjres<fdrthr))


#stats table creating

d<-data.table(iris)
names(d)
dd<-melt(d)


library(reshape2)
dmed<-reshape2::dcast(dd,variable~Species,median)
dmed2<-dmed[,-1]
rownames(dmed2)<-dmed[,1]

coms<-combn(names(dmed2),2,simplify = FALSE) # pairwise combinations of groups to compare
meddiffs<-sapply(coms,function(x) dmed2[,x[1]]-dmed2[,x[2]])
colnames(meddiffs)<-sapply(coms, function(x) paste0(x,collapse = " - "))
rownames(meddiffs)<-dmed$variable

adjpvalsdunn<-t(sapply(dunnTestres,function(x) x$res$P.adj))
colnames(adjpvalsdunn)<-as.character(dunnTestres[[1]]$res$Comparison)


CNM<-round(meddiffs,3)
rownames(CNM)<-rownames(dmed2)
CNM[,1:ncol(CNM)]<-apply(round(meddiffs,3),2,function(x) as.character(x))

# this for cycle marks * in all significant comparisons
for (i in rownames(adjpvalsdunn))
{
  CNM[i,which(adjpvalsdunn[i,colnames(CNM)]<0.1)]<-paste0(CNM[i,which(adjpvalsdunn[i,colnames(CNM)]<0.1)],"*")
}

#colmap<-colorRampPalette(c("green","white","orange"))(101) #for green to orange
colmap<-colorRampPalette(c("dark blue","white","dark red"))(101) #for blue to red

pdf(file=paste0(path,"Analyses/Results/Results-September-2020/median_diffs_NMRstool.pdf"), height = 7+nrow(CNM)/7, width=10)
print(plotfchheatmap(meddiffs, CNM=CNM, col=colmap,margins=c(20,12)))
dev.off()

#only on significant

sel<-which(apply(CNM, 1, function(x) length(findmetext("\\*",x)))>0) # selects only significant variables

pdf(file=paste0(path,"Analyses/Results/Results-September-2020/median_diffs_NMRstool-onlySignif.pdf"), height = 7+nrow(CNM)/7, width=10)
print(plotfchheatmap(meddiffs[sel,], CNM=CNM[sel,], col=colmap,margins=c(20,12)))
dev.off()