# open an example dataset (HapMap)
genofile <- snpgdsOpen(snpgdsExampleFileName())
# get population information
# or pop_code <- scan("pop.txt", what=character())
# if it is stored in a text file "pop.txt"
pop_code <- read.gdsn(index.gdsn(genofile, "sample.annot/pop.group"))
# get sample id
samp.id <- read.gdsn(index.gdsn(genofile, "sample.id"))
# run eigen-analysis
RV <- snpgdsEIGMIX(genofile)
# eigenvalues
RV$eigenval
# make a data.frame
tab <- data.frame(sample.id = samp.id, pop = factor(pop_code),
EV1 = RV$eigenvect[,1], # the first eigenvector
EV2 = RV$eigenvect[,2], # the second eigenvector
stringsAsFactors = FALSE)
head(tab)
# draw
plot(tab$EV2, tab$EV1, col=as.integer(tab$pop),
xlab="eigenvector 2", ylab="eigenvector 1")
legend("topleft", legend=levels(tab$pop), pch="o", col=1:4)
# define groups
groups <- list(CEU = samp.id[pop_code == "CEU"],
YRI = samp.id[pop_code == "YRI"],
CHB = samp.id[is.element(pop_code, c("HCB", "JPT"))])
prop <- snpgdsAdmixProp(RV, groups=groups)
# draw
plot(prop[, "YRI"], prop[, "CEU"], col=as.integer(tab$pop),
xlab = "Admixture Proportion from YRI",
ylab = "Admixture Proportion from CEU")
abline(v=0, col="gray25", lty=2)
abline(h=0, col="gray25", lty=2)
abline(a=1, b=-1, col="gray25", lty=2)
legend("topright", legend=levels(tab$pop), pch="o", col=1:4)
# run eigen-analysis
RV <- snpgdsEIGMIX(genofile, sample.id=samp.id[pop_code=="JPT"],
need.ibdmat=TRUE)
z <- RV$ibdmat
mean(c(z))
mean(diag(z))
# close the genotype file
snpgdsClose(genofile)
Run the code above in your browser using DataLab