data <- read.table("O://Genetik//MA//Silke//RCode//TestData//2pop.txt", header=T)
data[1,]
Loci <- data[,3:8]
Pop <- data[,9]

############################
# Calculate Nei's unbiased #
# gene diversity (1978)    #
# without resampling 	   #
############################

##### _WITHOUT_ resampling:
# Number of individuals per population = number of thalli
Ind.Nr.Pop.Thalli <- rep(NA,length(unique(Pop)))
data.split <- split(data,data$Pop)
for(i in 1:length(unique(Pop)))
	{
	Ind.Nr.Pop.Thalli[i] <-nrow(data.split[[i]])
	}
Ind.Nr.Pop.Thalli

# Alternative option of calculating Ind.Nr.Pop.Thalli 
# Ind.Nr.Pop.Thalli <- tapply(Pop,Pop,length)

##### Calculate Nei's unbiased gene diversity H
source("O://Genetik//MA//Silke//RCode//GeneDiversity//unbiased.gene.diversity.Nei1978_source-v2.R")
# change path accordingly
gendiversity <- gene.div(Loci, Pop, Ind.Nr.Pop=Ind.Nr.Pop.Thalli)
gendiversity


##############################################################################
##### Diesen Teil hier brauchen Leute, die schlau gesampelt haben, nicht! ####
##### The following code is only useful if you have unbalanced sampling   ####
##### design. Then, you would characteristically calculate a mean value   ####
##### of gene diversity over several resamplings of the data. 		  ####
##############################################################################

############################
# Calculate Nei's unbiased #
# gene diversity (1978)    #
# with resampling 	   #
############################

##### _WITH_ resampling (=randomly draw 1 thallus per tree)
# Number of individuals per population = number of trees in this case!!
Ind.Nr.Pop.Tree <- rep(NA,length(unique(Pop)))
data.split <- split(data,data$Pop)
for(i in 1:length(unique(Pop)))
	{
	Ind.Nr.Pop.Tree[i] <-length(unique(data.split[[i]]$TreeID)) 
	}
Ind.Nr.Pop.Tree

## use Resampling
#source("O://Genetik//MA//Silke//RCode//Resampling//Code_resampling.R")

resample <- function(d=data)
{
    Trees <- split(1:nrow(d), data$TreeID)
    Trees[[1]] # first data point

    a <- rep(NA,length(Trees))

    for(i in 1:length(Trees))
    {
    	a[i] <- Trees[[i]][1]
    	if(length(Trees[[i]]) > 1) 
    	{
	    a[i] <- sample(Trees[[i]], 1, replace=FALSE)
    	}
    }
    d[a,]
}

nr.res <- 3 # enter the number of resamplings 
data.res <- list()
for(i in 1:nr.res)
	{
	data.res[[i]] <- resample(data)
	}



#Assign "Loci.1" and "Pop.1" to each resampled data set
Loci.1 <- list()
Pop.1 <- list()
for(n in 1:nr.res)
	{
	Loci.1[[n]] <- data.res[[n]][,3:8] #change the columns if your data file differs
	Pop.1[[n]] <- data.res[[n]][,9] #change the columns if your data file differs
	}
#Calculate gene diversity
rm(Loci)
rm(Pop)
rm(data)
source("O://Genetik//MA//Silke//RCode//GeneDiversity//unbiased.gene.diversity.Nei1978_source-v2.R")
# change pathname according to the location where you saved this file

test1 <- list()
n <- 1
for(i in 1:nr.res)
	{
	Loci <- Loci.1[[i]]
	Pop <- Pop.1[[i]]
	test1[[n]] <- gene.div(Loci, Pop, Ind.Nr.Pop=Ind.Nr.Pop.Tree)
	n <- n+1
	}
test1[[1]]$GeneDiversity.unbiased 
# returns gene diversity of the first resampled dataset, [[2]] of the second aso.









