doing some pca on snps from Axiom_aegypti1 chip

genlight_PCA.R

library("adegenet")
library("rgl")
library("phangorn")
library("ape")

setwd("C:/Users/ben/Desktop/temp/pca")
#Read in the data
X <- read.PLINK("pca6.raw", parallel=FALSE)

col <- rainbow(length(levels(pop(X))), start=.3, end=.1)
icol <- rainbow(length(levels(pop(X))), start=.3, end=.1)[as.integer(X@pop)]

#look at missingness, blocks of genotypes
# glPlot(X, posi="topright")

#allele Freqs
# myFreq <- glMean(X)
# hist(myFreq, proba=TRUE, col="lightblue", xlab="Allele frequencies",
# main="Distribution of (second) allele frequencies")
# temp <- density(myFreq)
# lines(temp$x, temp$y*1.8,lwd=3)

#NJ
# tre <- nj(dist(as.matrix(X)))
# plot(tre, typ="fan", cex=0.7)

##PCA
pca1 <- glPca(X, cent=TRUE, scale = FALSE, parallel=FALSE, nf = 3)

# myCol <- colorplot(pca1$scores,pca1$scores, transp=TRUE, cex=2)
# abline(h=0,v=0, col="grey")
# add.scatter.eig(pca1$eig[1:20],2,1,2, posi="bottomright", inset=.12, ratio=.2)

s.class(pca1$scores, pop(X), xax=1, yax=2, col=transp(col,.6), axesell=FALSE, cstar=0, cpoint=1, grid=FALSE)

plot3d(pca1$scores, size = 1, type="s", col=icol, xlab="PC1", ylab="PC2", zlab="PC3", box=FALSE)
#rgl.postscript("persp3dd.pdf","pdf")

##NJ Trees
# D <- dist(X)
# treeNJ <- NJ(D)
# plot.phylo(tre, type = "unrooted", use.edge.length=FALSE, show.tip = FALSE)
# legend("topright", X$pop, fill = col, ncol=2, cex = .8)
# tiplabels(X$ind.names, col = icol, cex = .6, frame = "none", adj = -.1)
# title(paste("Neighbor-Joining Tree (Euclidean Distance)\n", length(X$ind.names),"Individuals,", length(X$pop.names), "Sample Groups,",length(X$loc.names), " Loci"))

# Xpop <- genind2genpop(X)
# popD = dist.genpop(Xpop, method=2)
# popNJ <-NJ(popD)
# plot.phylo(popNJ, type = "unrooted",use.edge.length=FALSE, show.tip = FALSE)
# legend("topright", X$pop.names, fill = col, ncol=2, cex = .8)
# tiplabels(X$pop.names, col = col, cex = .6, frame = "none", adj = -.1)
# title("Neighbor-Joining Network (Euclidean Distance)")

# ##DAPC
# grps <- find.clusters(X, n.clust=length(X$pop.names), scale=FALSE, n.pca=100)
# dapc1 <- dapc(X, grps$grp)
# temp <- optim.a.score(dapc1)
# dapc2 <- dapc(X, grps$grp)

# compoplot(dapc2, cleg=.6, posi=list(x=0,y=1.13),lab=lab, col=col)