## This script calculates ratio between detected number of alleles and effetcive number of allele ## Pojskic et Kalamujic message("ALRATIO - Alleles-ratio, by Naris Pojskic (naris.pojskic@ingeb.unsa.ba)") message("University of Sarajevo, Institute for Genetic Engineering and Biotechnology") message("Laboratory for Bioinformatics and Biostatistics") message("Bioinfo Research Group - www.bioinfo.ba") message("---------------------------------------------------------------------------") #File manipulation message("R script and CSV input data file must be in the same directory") message("Select CSV input data file") cfi <- file.choose(new = FALSE) cfd <-dirname(cfi) setwd(cfd) getwd() #Import data mydata <- read.csv2(cfi) a <- mydata[[1]] loc <- as.character(a) b <- mydata[[2]] bb <- factor(b) He <- as.numeric(sub(",", ".", paste(bb), fixed = TRUE)) c <- mydata[[3]] cc <- factor(c) An <- as.numeric(sub(",", ".", paste(cc), fixed = TRUE)) #Ae and R calculation" Ae <- round(1/(1-He),2) R <- round(Ae/An,3) #Z value Ro<-rep(c(1),times=ncol(t(R))) z<-(R-Ro)/sqrt((R*(1-R)/2)) #P value pval <-pnorm(-abs(z)) round(pval,4) #Graph and data tabela <- cbind(loc, He, An, Ae, R, round(z,3), round(pval,5)) colnames(tabela) <- c("Loc","He","An","Ae","R", "Z", "P") print(tabela) write.table(format(tabela), file = "output.txt", dec = ",", sep = "\t", row.names = FALSE, quote=FALSE) tabgraf <- cbind(He, R) colnames(tabgraf) <- c("He","R") jpeg("plot1.jpg") barplot(as.matrix(tabgraf), beside=TRUE, legend.text=loc, args.legend=list(bty="n",horiz=FALSE, "topright"),border="white", ylim=c(0,1),ylab="Values", main="Expected heterozygosity, Ae/An ratio") dev.off() jpeg('plot2.jpg') plot(An, Ae, col='blue',main='An vs Ae', pch=1, cex=0.9, type='p',xlab='Observed number of alleles',ylab='Effective number of alleles') dev.off() message("Statistical significance level at P<0.01") h <- "File (output.txt) and graphs (plot1.jpg; plot2.jpg) are saved in" print(paste(h, cfd))