PCA in R

plot_pca.R

first_time <- function() {
    t1 <- read.table("../ALL.chrY.phase3_integrated_v1b.20130502.genotypes.vcf.gz")
    print(nrow(t1))
    nucs <- c('A', 'C', 'G', 'T')
    t1 <- t1[t1$V4 %in% nucs & t1$V5 %in% nucs, ]
    print(nrow(t1))
    t2 <- t1[, 10:ncol(t1)]
    t2 <- sapply(t2, as.character)
    save(t2, file = "t2")
}

load("t2")
t2.pca <- prcomp(t2)

info <- read.csv("sample_pop.tsv")
population <- info$V2
print(population)

pdf("pca-y.pdf")
plot(t2.pca$x[, 1:2], xlab = "PC1", ylab = "PC2", col = population)
dev.off()