Created
September 12, 2020 10:52
-
-
Save romunov/d78d27a0852e23ba80325994c4e1c495 to your computer and use it in GitHub Desktop.
modified version of scatter.dapc for plotting points as individual points
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
library(adegenet) | |
data(H3N2) | |
dapc1 <- dapc(H3N2, pop=H3N2$other$epid, n.pca=30,n.da=6) | |
custom.cols <- rep(c("red", "blue"), times = c(1000,903)) | |
scatter(dapc1, col = custom.cols) |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
scatter.dapc <- function(x, xax = 1, yax = 2, grp = x$grp, | |
col = seasun(length(levels(grp))), | |
pch = 20, bg = "white", solid = .7, | |
scree.da = TRUE, scree.pca = FALSE, | |
posi.da = "bottomright", | |
posi.pca = "bottomleft", | |
bg.inset = "white", | |
ratio.da = .25, ratio.pca = .25, | |
inset.da = 0.02, inset.pca = 0.02, inset.solid = .5, | |
onedim.filled = TRUE, mstree = FALSE, | |
lwd = 1, lty = 1, segcol = "black", | |
legend = FALSE, posi.leg = "topright", | |
cleg = 1, txt.leg = levels(grp), | |
cstar = 1, cellipse = 1.5, axesell = FALSE, | |
label = levels(grp), clabel = 1, xlim = NULL, ylim = NULL, | |
grid = FALSE, addaxes = TRUE, origin = c(0,0), | |
include.origin = TRUE, sub = "", csub = 1, possub = "bottomleft", | |
cgrid = 1, pixmap = NULL, contour = NULL, | |
area = NULL, label.inds = NULL, ...){ | |
ONEDIM <- xax==yax | ncol(x$ind.coord)==1 | |
## recycle color and pch | |
pch <- rep(pch, length(levels(grp))) | |
col <- transp(col, solid) | |
bg.inset <- transp(bg.inset, inset.solid) | |
## handle grp | |
if(is.null(grp)){ | |
grp <- x$grp | |
} | |
## handle xax or yax NULL | |
if(is.null(xax)||is.null(yax)){ | |
xax <- 1L | |
yax <- ifelse(ncol(x$ind.coord)==1L, 1L, 2L) | |
ONEDIM <- TRUE | |
} | |
## handle 1 dimensional plot | |
if(!ONEDIM){ | |
## set par | |
opar <- par(mar = par("mar")) | |
par(mar = c(0.1, 0.1, 0.1, 0.1), bg=bg) | |
on.exit(par(opar)) | |
axes <- c(xax,yax) | |
## basic empty plot | |
s.class(x$ind.coord[,axes], fac = grp, col = col, cpoint = 0, | |
cstar = cstar, cellipse = cellipse, axesell = axesell, | |
label = label, clabel = clabel, xlim = xlim, ylim = ylim, | |
grid = grid, addaxes = addaxes, origin = origin, | |
include.origin = include.origin, | |
sub = sub, csub = csub, possub = possub, | |
cgrid = cgrid, pixmap = pixmap, | |
contour = contour, area = area) | |
## add points | |
colfac <- pchfac <- grp | |
levels(colfac) <- col | |
levels(pchfac) <- pch | |
colfac <- as.character(colfac) | |
pchfac <- as.character(pchfac) | |
if(is.numeric(col)) colfac <- as.numeric(colfac) | |
if(is.numeric(pch)) pchfac <- as.numeric(pchfac) | |
points(x$ind.coord[,xax], | |
x$ind.coord[,yax], | |
col = col, pch = pchfac, ...) | |
s.class(x$ind.coord[,axes], fac = grp, col = col, cpoint = 0, | |
add.plot=TRUE, cstar = cstar, cellipse = cellipse, | |
axesell = axesell, label = label,clabel = clabel, | |
xlim = xlim, ylim = ylim, grid = grid, | |
addaxes = addaxes, origin = origin, | |
include.origin = include.origin, | |
sub = sub, csub = csub, possub = possub, | |
cgrid = cgrid, pixmap = pixmap, | |
contour = contour, area = area) | |
## Add labels of individuals if specified. Play around with "air" to get | |
## a satisfactory result. | |
if (!is.null(label.inds) & is.list(label.inds)) { | |
appendList <- function (x, val) { | |
# recursevly "bind" a list into a longer list, | |
# from http://stackoverflow.com/a/9519964/322912 | |
stopifnot(is.list(x), is.list(val)) | |
xnames <- names(x) | |
for (v in names(val)) { | |
x[[v]] <- if (v %in% xnames && is.list(x[[v]]) && is.list(val[[v]])) | |
appendList(x[[v]], val[[v]]) | |
else c(x[[v]], val[[v]]) | |
} | |
x | |
} | |
do.call("orditorp", | |
c(appendList(list(x = x$ind.coord[, c(xax, yax)], | |
display = "species"), | |
label.inds))) | |
} | |
## add minimum spanning tree if needed | |
if(mstree){ | |
meanposi <- apply(x$tab,2, tapply, grp, mean) | |
D <- dist(meanposi)^2 | |
tre <- ade4::mstree(D) | |
x0 <- x$grp.coord[tre[,1], axes[1]] | |
y0 <- x$grp.coord[tre[,1], axes[2]] | |
x1 <- x$grp.coord[tre[,2], axes[1]] | |
y1 <- x$grp.coord[tre[,2], axes[2]] | |
segments(x0, y0, x1, y1, lwd = lwd, lty = lty, col = segcol) | |
} | |
} else { | |
## set screeplot of DA to FALSE (just 1 bar) | |
scree.da <- FALSE | |
## get plotted axis | |
if(ncol(x$ind.coord)==1) { | |
pcLab <- 1 | |
} else{ | |
pcLab <- xax | |
} | |
## get densities | |
ldens <- tapply(x$ind.coord[,pcLab], grp, density) | |
allx <- unlist(lapply(ldens, function(e) e$x)) | |
ally <- unlist(lapply(ldens, function(e) e$y)) | |
par(bg=bg) | |
plot(allx, ally, type = "n", | |
xlab = paste("Discriminant function", pcLab), | |
ylab = "Density") | |
for(i in 1:length(ldens)){ | |
if(!onedim.filled) { | |
lines(ldens[[i]]$x, ldens[[i]]$y, col = col[i], lwd = 2) # add lines | |
} else { | |
polygon(c(ldens[[i]]$x, rev(ldens[[i]]$x)), | |
c(ldens[[i]]$y, rep(0,length(ldens[[i]]$x))), | |
col = col[i], lwd = 2, border = col[i]) # add lines | |
} | |
points(x = x$ind.coord[grp==levels(grp)[i], pcLab], | |
y = rep(0, sum(grp==levels(grp)[i])), | |
pch = "|", col = col[i]) # add points for indiv | |
} | |
} | |
## ADD INSETS ## | |
## group legend | |
if(legend){ | |
## add a legend | |
temp <- list(...)$cex | |
if(is.null(temp)) temp <- 1 | |
if(ONEDIM | temp<0.5 | all(pch == "")) { | |
legend(posi.leg, fill = col, legend = txt.leg, | |
cex = cleg, bg = bg.inset) | |
} else { | |
legend(posi.leg, col = col, legend = txt.leg, cex = cleg, | |
bg = bg.inset, pch = pch, pt.cex = temp) | |
} | |
} | |
## eigenvalues discriminant analysis | |
if(scree.da && ratio.da>.01) { | |
inset <- function(){ | |
myCol <- rep("white", length(x$eig)) | |
myCol[1:x$n.da] <- "grey" | |
myCol[c(xax, yax)] <- "black" | |
myCol <- transp(myCol, inset.solid) | |
barplot(x$eig, col=myCol, xaxt="n", yaxt="n", ylim=c(0, x$eig[1]*1.1)) | |
mtext(side=3, "DA eigenvalues", line=-1.2, adj=.8) | |
box() | |
} | |
add.scatter(inset(), posi = posi.da, ratio = ratio.da, | |
bg.col = bg.inset, inset = inset.da) | |
} | |
## eigenvalues PCA | |
if(scree.pca && !is.null(x$pca.eig) && ratio.pca>.01) { | |
inset <- function(){ | |
temp <- 100* cumsum(x$pca.eig) / sum(x$pca.eig) | |
myCol <- rep(c("black","grey"), c(x$n.pca, length(x$pca.eig))) | |
myCol <- transp(myCol, inset.solid) | |
plot(temp, col=myCol, ylim=c(0,115), | |
type="h", xaxt="n", yaxt="n", xlab="", ylab="", lwd=2) | |
mtext(side=3, "PCA eigenvalues", line=-1.2, adj=.1) | |
} | |
add.scatter(inset(), posi = posi.pca, ratio = ratio.pca, | |
bg.col = bg.inset, inset = inset.pca) | |
} | |
return(invisible(match.call())) | |
} |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment