Skip to content

Instantly share code, notes, and snippets.

Embed
What would you like to do?
# PCA
set.seed(100)
pca <- snpgdsPCA(genofile, sample.id = geno.sample.ids,
snp.id = snpset.ibd, num.thread = 1)
pctab <- data.frame(sample.id = pca$sample.id,
PC1 = pca$eigenvect[,1],
PC2 = pca$eigenvect[,2],
stringsAsFactors = F)
# Subset and/or reorder origin accordingly
origin <- origin[match(pca$sample.id, origin$sample.id),]
pcaCol <- rep(rgb(0,0,0,.3), length(pca$sample.id)) # Set black for chinese
pcaCol[origin$Country == "I"] <- rgb(1,0,0,.3) # red for indian
pcaCol[origin$Country == "M"] <- rgb(0,.7,0,.3) # green for malay
png("PCApopulation.png", width = 500, height = 500)
plot(pctab$PC1, pctab$PC2, xlab = "PC1", ylab = "PC2", col = pcaCol, pch = 16)
abline(h = 0, v = 0, lty = 2, col = "grey")
legend("top", legend = c("Chinese", "Indian", "Malay"), col = 1:3, pch = 16, bty = "n")
dev.off()
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment