## MDS plot (cluster algorithm)  E.Korsching  2024

mds.eR <- function(expr, coldata, col, pch, main="", p.n=T, norm="none", cex=1){
	# MDS from limma via edgeR
	# expr: matrix: columns = samples
	# coldata: create df: one column: condition=group labels (char), row.names=dimnames(expr)[[2]]
	# col: color vector,  pch: pch number vector -- both vectors corresponding to coldata and group structure
	# p.n: T: points F: names,  main: title
	# norm: "none": do not scale columns [sometimes errors] , "TMMwsp","TMM","RLE","upperquartile": scaling [sometimes errors]
	require(ggplot2)
	require(edgeR)
	require(scales)
	
	rn <- rownames(coldata)
	# create colors & join with coldata
	coldata <- cbind(condition=coldata, col=col, pch=pch)
	rownames(coldata) <- rn
	
	# DGEList object -- library sizes should be finite and non-negative ... and non 0 + not too sparse ... try
	y <- DGEList(counts=expr, group=coldata$condition)
	# scaling
	if(norm!="none"){ y <- normLibSizes(y, method=norm, Acutoff=-1e10, p=0.75) }
	
	# generate the layout matrix
	lmat <- matrix(c(1,2), nrow=1, byrow=T)
	layout(lmat, width=c(0.8,0.2), height=1)	
#	layout.show(1); return();
	
	# plot similarity
	pm <- par()$mar
	par(mar=c(5.1,4.1,4.1,0.3))
	if(p.n){	# points
		plotMDS(y, pch=coldata$pch, col=coldata$col, cex=cex)	# bg > 21 fill
	}else{		# names
		plotMDS(y, col=coldata$col, cex=cex)
	}
	title(paste("MDS,",main))
	
	par(mar=c(0.1,0.1,0.1,0.1))
	plot(c(0,1),c(0,1), type="n", axes=F, xlab="", ylab="")
	legend(x=0, y=0.8,		# tune in case
			title=NULL,
			legend=unique(coldata$condition),
			col=unique(coldata$col),
			pch=unique(coldata$pch),
			bty="n", xjust=0, cex=0.8)
	par(mar=pm)
	return()
}

# errors: https://support.bioconductor.org/p/70347/

#aa <- as.matrix(x[,c(k.p.thal.first,k.p.thal.second,k.p.thal.third, k.Npn.thal$gray, k.n.thal.first,k.n.thal.second,k.n.thal.third)])
#ab <- data.frame(condition=
#				c(rep("th1a",length(k.p.thal.first)),rep("th2a",length(k.p.thal.second)),rep("th3a",length(k.p.thal.third)),rep("other",length(k.Npn.thal$gray)),
#						rep("th1n",length(k.n.thal.first)),rep("th2n",length(k.n.thal.second)),rep("th3n",length(k.n.thal.third)) ),
#		row.names=dimnames(aa)[[2]] )
#ac <- color.groups(c("deepskyblue1","chartreuse1","tan1","gray","blue1","chartreuse4","tan3"),
#		k.p.thal.first,k.p.thal.second,k.p.thal.third, k.Npn.thal$gray, k.n.thal.first,k.n.thal.second,k.n.thal.third)
#ad <- rep(20,length(ac))
#mds.eR(expr=aa, coldata=ab, col=ac, pch=ad, main="Thalassemia, raw", p.n=T, norm="none", cex=0.9)