my.dir <- "C:/Documents and Settings/Rich/My Documents/projects2/Zoo-MyWeb/Rdata"
#source(paste(my.dir, "/scripts_stat295F11.R", sep=""))
 source("http://www.uvm.edu/~rsingle/Rdata/scripts_stat295F14.R")

#ASSIGN THE EXAMPLE DATA TO THE OBJECT x 
 #z <- otherdata("example.dat",web=F, Rdata.dir=my.dir)
  z <- otherdata("example.dat") 

#-------------------------------------------------------------------
 dim(z)
 
 length(z$locus)
 
 z$locus

 levels(z$locus)

 unique(z$locus)

 unique(as.character(z$locus))

 locus.list <- unique(as.character(z$locus))

 locus.list[1]

 locus.list[2:3]

 locus.list[-2]

 z[z$locus==locus.list[1],]

 tab <- xtabs(z$allele.count ~ z$allele + z$locus)
 tab
 
 barplot(tab, col=rainbow(21), legend=T, beside=T)
 title("Barplot of Allele Frequencies")
 
 barplot(tab, col=rainbow(21), beside=F)
 legend("topright",rownames(tab),fill=rainbow(21))
 
 tmp<-barplot(tab, col=rainbow(21), beside=F)
 tmp #[1] 0.7 1.9 3.1 4.3
 legend(4.3,34,rownames(tab),fill=rainbow(21))
 
 
#A function for computing heterozygosity 
 hetz.function <- function(x){1-sum(x^2)}
 
 locus.name <- "D12S1638"
 sum( z$allele.freq[z$locus==locus.name] )           #Sum should be 1.0
 hetz.function( z$allele.freq[z$locus==locus.name] ) 
 
 for (locus.name in locus.list)
 {
   cat(locus.name, "\n")
   print( hetz.function( z$allele.freq[z$locus==locus.name] ) )
 }
 
#---------- 
 freqs.tab <- xtabs(z$allele.freq ~ z$allele + z$locus )
 barplot(freqs.tab)     

 freqs.tab <- xtabs(z$allele.freq ~ z$allele + z$locus + z$pop) #if different populations too
 barplot(freqs.tab)     #Error -> barplot() requires a matrix argument
 freqs.tab
 freqs.tab[,,1]
 barplot(freqs.tab[,,1])

 
#Hetz for each locus (column) & population
 apply(freqs.tab, 2, sum           ) #second arg=2 --> column operation, 1-->row operation
 apply(freqs.tab, 2, hetz.function ) #second arg=2 --> column operation, 1-->row operation

#----------
 z$locus
 locus
 
 search()
 attach(z)
 search()
 locus
 xtabs(allele.count ~ allele + locus)
 detach(z)
 locus
 
 attach(z) 
 attach(z) #attach a 2nd time --> masking