PCA/Eigensoft/Eigenstrat

eigenstrat<-function(geno,maxMis=0,minMaf=0.01){
## geno: snp x ind matrix of genotypes \in 0,1,2
##maxMis maximum allowed missing genotypes for a site

nMis <- rowSums(is.na(geno))
freq <- rowMeans(geno,na.rm=T)/2               # get allele frequency
keep <- freq>minMaf&freq<(1-minMaf) & nMis<=maxMis         # remove sites with non-polymorphic data
freq<-freq[keep]
geno<-geno[keep,]
snp<-nrow(geno)                           #number of snps used in analysis
ind<-ncol(geno)                           #number of individuals used in analuysis
M <- (geno-freq*2)/sqrt(freq*(1-freq))       #normalize the genotype matrix
M[is.na(M)]<-0
X<-t(M)%*%M                               #get the (almost) covariance matrix
X<-X/(sum(diag(X))/(snp-1))
E<-eigen(X)

mu<-(sqrt(snp-1)+sqrt(ind))^2/snp         #for testing significance (assuming no LD!)
sigma<-(sqrt(snp-1)+sqrt(ind))/snp*(1/sqrt(snp-1)+1/sqrt(ind))^(1/3)
E$TW<-(E$values[1]*ind/sum(E$values)-mu)/sigma E$mu<-mu
E$sigma<-sigma E$nSNP <- nrow(geno)
E$nInd <- ncol(geno) class(E)<-"eigenstrat" E } plot.eigenstrat<-function(x,col=1,...) plot(x$vectors[,1:2],col=col,...)

print.eigenstrat<-function(x)


LD pruning in R

Install R package

wget http://www.popgen.dk/albrecht/misc_Rpackages/Rpakker/pruning_0.51.tar.gz
R CMD INSTALL pruning_0.51.tar.gz