genetic.dist <- function(Data=Loci, Pop)
{
Result <- list()
# Calculate allele frequency
# 1) number of alleles in each population at each locus
Nr.ind.Pop <- tapply(Pop, Pop, length)
allele.frequency  <- list()
test <- list()
for(i in 1:ncol(Loci))
	{
	b <- list(Loci[,i], Pop)
	allele.frequency [[i]] <- tapply(Loci[,i], b, length)
	}
allele.frequency 
# 2) calculate allele frequency
for(i in 1:ncol(Loci))
	{
	for(j in 1:length(Nr.ind.Pop))
		{
		allele.frequency[[i]][,j] <-allele.frequency [[i]][,j]/Nr.ind.Pop[j]
		}
	}
allele.frequency	# allele frequency per locus and population 	
# 3) replace the NAs with zero
b <- list()
for(i in 1:ncol(Loci))
	{
	b <- allele.frequency[[i]]
	b[is.na(b)]=0
	allele.frequency[[i]] <- b
	}
test <-allele.frequency 
test # allele.frequency at each locus, list.
# allele.frequency: allele frequency for each locus
# for haploid data, for diploid data only if the two alleles are in rows


# Jxy (average homozygosity between pop x and pop y over r loci)
# Multipliziere Spalte 1:2 mit Spalte 2:3 von test1
# Summiere das Ergebnis auf 
# Teile jede Summe durch die Anzahl loci

C.nr <-(length(unique(Pop))^2 - length(unique(Pop)))/2;C.nr
Jx <- rep(0, length(unique(Pop)));Jx
Jxy <- rep(0, C.nr);Jxy
SQRT.xy <- rep(0, C.nr);SQRT.xy
for(x in 1:ncol(Loci))
{
test1 <-test[[x]]
#test1 <-test[[1]]
b <- list()
n <- 1
for(j in 1:(length(unique(Pop))-1))
	{
	for(k in (j+1):length(unique(Pop)))
		{
		b[[n]] <- (test1[,j]*test1[,k]) # Produkt von Pop x und y locus x 
		n <- n+1
		}
	}
a <- list()
for(i in 1:C.nr)
	{
	a[[i]] <- (sum(b[[i]])) # Summiere die Produkte b in die Liste a[[i]]
	}
a
# Schreibe a[[i]] in Liste d.1
n <- 1
d.1 <- c()
for(j in 1:C.nr)
	{
	d.1[j] <-c(a[[n]])
	n <- n+1
	}
Jxy <- Jxy + d.1 
#Jxy entspricht Jij

# Berechne Jx ("average homozygosity within pop")
n <- 1
b.1 <- list()
for(j in 1:length(unique(Pop)))
	{
	b.1[[n]] <- (test1[,j]^2) # Quadrat of the allele frquencies xij2 
	n <- n+1
	}
b.1
f <- list()
for(i in 1:length(unique(Pop)))
	{
	f[[i]] <- (sum(b.1[[i]])) # Sum products b into list a[[i]]
	}
f
# Write f[[i]] into vektor d.2
n <- 1
d.2 <- c()
for(j in 1:length(unique(Pop)))
	{
	d.2[j] <-c(f[[n]])
	n <- n+1
	}
Jx <- Jx + d.2

# Calculate sqrt(xij*yij)
b <- list()
n <- 1 
for(j in 1:(length(unique(Pop))-1))
	{
	for(k in (j+1):length(unique(Pop)))
		{
		b[[n]] <- sqrt(test1[,j]*test1[,k]) # Produkt von Pop x und y locus x 
		n <- n+1
		}
	}
b
a <- list()
for(i in 1:C.nr)
	{
	a[[i]] <- sum(b[[i]])
	}
a
# Write a[[i]] into vektor d.3
n <- 1
d.3 <- c()
for(j in 1:C.nr)
	{
	d.3[j] <-c(a[[n]])
	n <- n+1
	}
SQRT.xy <- SQRT.xy + d.3

}

Da <- 1-1/(ncol(Loci))* SQRT.xy; Da # Nei's (1983) Da distance in einem Vektor
Jxy.1 <- Jxy/ncol(Loci);Jxy.1 # entspricht Jxy in Nei 1972, gemittelt über Loci.
Jx.1 <- Jx/ncol(Loci);Jxy.1 # entspricht Jx und Jy in Nei 1972, gemittelt über Loci.

# Write Jxy in lower.tri of matrix J 
J <- matrix(NA, length(unique(Pop)), length(unique(Pop)), dimnames=list(unique(Pop),unique(Pop)))
diag(J) <- Jx.1
n <-1 
for(j in 1:(length(unique(Pop))-1))
	{
	for(k in (j+1):length(unique(Pop)))
		{
		J[k,j] <- Jxy.1[n]
		n <- n+1
		}
	}
J
# Write Jx into diagonal of matrix J


# Nei's 1973 gene diversity within pop i (=1-Jx)
# This is the probability of nonidentity within each deme.
Di <- c()
Di <- 1- (diag(J))
Di.Nei <- data.frame(Di)
Di.Nei


# Write Da into the matrix D.Nei.Da
D.Nei.Da <- matrix(NA, length(unique(Pop)), length(unique(Pop)), dimnames=list(unique(Pop),unique(Pop)))
n <-1 
for(j in 1:(length(unique(Pop))-1))
	{
	for(k in (j+1):length(unique(Pop)))
		{
		D.Nei.Da[k,j] <- Da[n]
		n <- n+1
		}
	}
D.Nei.Da


#D.Nei.Std is Nei's standard genetic distance (1972)
D.Nei.Std <- matrix(NA, length(unique(Pop)), length(unique(Pop)), dimnames=list(unique(Pop),unique(Pop)))
for(i in 1:(length(unique(Pop))-1))
	{
	for(j in (i+1):length(unique(Pop)))
		{
		D.Nei.Std[j,i] <- -log(J[j,i]/(sqrt(J[i,i]*J[j,j])))
		}
	}
D.Nei.Std

D.Nei.Min <- matrix(NA, length(unique(Pop)), length(unique(Pop)), dimnames=list(unique(Pop),unique(Pop)))
for(i in 1:(length(unique(Pop))-1))
	{
	for(j in (i+1):length(unique(Pop)))
		{
		D.Nei.Min[j,i] <- (J[i,i]+J[j,j])/2 - J[j,i]
		}
	}
D.Nei.Min

# Latter's 1972 Fst(=phi*)
Fst <- matrix(NA, length(unique(Pop)), length(unique(Pop)), dimnames=list(unique(Pop),unique(Pop)))
for(i in 1:(length(unique(Pop))-1))
	{
	for(j in (i+1):length(unique(Pop)))
		{
		Fst[j,i] <- ((J[i,i]+ J[j,j])/2 - J[j,i])/ (1-J[j,i])
		}
	}
Fst

# standardised Fst: Fst/(1-Fst)
#Fst.std <- matrix(NA, length(unique(Pop)), length(unique(Pop)), dimnames=list(unique(Pop),unique(Pop)))
#for(i in 1:(length(unique(Pop))-1))
#	{
#	for(j in (i+1):length(unique(Pop)))
#		{
#		Fst.std[j,i] <- Fst[j,i]/(1-Fst[j,i])
#		}
#	}
#Fst.std
Result <- list()
#Result$Fst.std <- Fst.std
Result$Fst <- Fst
#Result$D.Nei.Std <- D.Nei.Std
Result$D.Nei.Min <- D.Nei.Min
Result$D.Nei.Da <- D.Nei.Da
Result$GeneDiversity <- Di.Nei
Result$J <- J #average within (Jx) and between (Jxy) population homozygosity
Result
}


#Usage
#gendiff <- genetic.dist(Loci, Pop)
#gendiff
#gendiff$D.Nei.Min