.Rhistory

if(ls!="012" & ls!="0129"){print("Error: Your snp matrix does not have 0,1, and 2"); break}
### Checking that locusNames and popNames have the same lengths as the columns and rows of SNPmat
if(dim(SNPmat)[1]!=length(popname) ){
print("Error: your population names do not match your SNP matrix")
break}
if(dim(SNPmat)[2]!=length(locusname)){
print("Error:  your locus names do not match your SNP matrix")
break}
writeLines("Calculating FSTs, may take a few minutes...")
nloci <- length(locusname)
FSTmat2 <- matrix(NA, nrow=nloci, ncol=8)
return (SNPmat[1,])
# for (i in 1:nloci){
# 	   FSTmat[i,]=getFSTs_diploids(popname,SNPmat[i,])
# 	   if (i==10000){print(paste(i, "done of", nloci))}
# }
# outTemp=as.data.frame(FSTmat2)
# return (cbind(locusname,outTemp))
}
MakeDiploidFSTMat(tempdata,locusNamesTemp,popNamesTemp)
MakeDiploidFSTMat = function(SNPmat,locusNames,popNames){
# SNPmat is a matrix with individuals in rows and snps in columns
# 0, 1, or 2 represent the number of copies of the focal allele, and 9 is for missing data
# locusNames is a character vector of names of each SNP
# popNames is a character vector with the population identifier for each individual
locusname <- unlist(locusNames)
popname <- unlist(popNames)
### Check that SNPmat has appropriate values (0, 1, 2, or 9, only)
snplevs <- levels(as.factor(unlist(SNPmat)))
ls <- paste(snplevs, collapse="")
if(ls!="012" & ls!="0129"){print("Error: Your snp matrix does not have 0,1, and 2"); break}
### Checking that locusNames and popNames have the same lengths as the columns and rows of SNPmat
if(dim(SNPmat)[1]!=length(popname) ){
print("Error: your population names do not match your SNP matrix")
break}
if(dim(SNPmat)[2]!=length(locusname)){
print("Error:  your locus names do not match your SNP matrix")
break}
writeLines("Calculating FSTs, may take a few minutes...")
nloci <- length(locusname)
FSTmat2 <- matrix(NA, nrow=nloci, ncol=8)
return (SNPmat[,1])
# for (i in 1:nloci){
# 	   FSTmat[i,]=getFSTs_diploids(popname,SNPmat[,i])
# 	   if (i==10000){print(paste(i, "done of", nloci))}
# }
# outTemp=as.data.frame(FSTmat2)
# return (cbind(locusname,outTemp))
}
MakeDiploidFSTMat(tempdata,locusNamesTemp,popNamesTemp)
MakeDiploidFSTMat = function(SNPmat,locusNames,popNames){
# SNPmat is a matrix with individuals in rows and snps in columns
# 0, 1, or 2 represent the number of copies of the focal allele, and 9 is for missing data
# locusNames is a character vector of names of each SNP
# popNames is a character vector with the population identifier for each individual
locusname <- unlist(locusNames)
popname <- unlist(popNames)
### Check that SNPmat has appropriate values (0, 1, 2, or 9, only)
snplevs <- levels(as.factor(unlist(SNPmat)))
ls <- paste(snplevs, collapse="")
if(ls!="012" & ls!="0129"){print("Error: Your snp matrix does not have 0,1, and 2"); break}
### Checking that locusNames and popNames have the same lengths as the columns and rows of SNPmat
if(dim(SNPmat)[1]!=length(popname) ){
print("Error: your population names do not match your SNP matrix")
break}
if(dim(SNPmat)[2]!=length(locusname)){
print("Error:  your locus names do not match your SNP matrix")
break}
writeLines("Calculating FSTs, may take a few minutes...")
nloci <- length(locusname)
FSTmat2 <- matrix(NA, nrow=nloci, ncol=8)
for (i in 1:nloci){
FSTmat[i,]=getFSTs_diploids(popname,SNPmat[,i])
if (i==10000){print(paste(i, "done of", nloci))}
}
outTemp=as.data.frame(FSTmat2)
return (cbind(locusname,outTemp))
}
MakeDiploidFSTMat(tempdata,locusNamesTemp,popNamesTemp)
MakeDiploidFSTMat = function(SNPmat,locusNames,popNames){
# SNPmat is a matrix with individuals in rows and snps in columns
# 0, 1, or 2 represent the number of copies of the focal allele, and 9 is for missing data
# locusNames is a character vector of names of each SNP
# popNames is a character vector with the population identifier for each individual
locusname <- unlist(locusNames)
popname <- unlist(popNames)
### Check that SNPmat has appropriate values (0, 1, 2, or 9, only)
snplevs <- levels(as.factor(unlist(SNPmat)))
ls <- paste(snplevs, collapse="")
if(ls!="012" & ls!="0129"){print("Error: Your snp matrix does not have 0,1, and 2"); break}
### Checking that locusNames and popNames have the same lengths as the columns and rows of SNPmat
if(dim(SNPmat)[1]!=length(popname) ){
print("Error: your population names do not match your SNP matrix")
break}
if(dim(SNPmat)[2]!=length(locusname)){
print("Error:  your locus names do not match your SNP matrix")
break}
writeLines("Calculating FSTs, may take a few minutes...")
nloci <- length(locusname)
FSTmat <- matrix(NA, nrow=nloci, ncol=8)
for (i in 1:nloci){
FSTmat[i,]=getFSTs_diploids(popname,SNPmat[,i])
if (i==10000){print(paste(i, "done of", nloci))}
}
outTemp=as.data.frame(FSTmat2)
return (cbind(locusname,outTemp))
}
MakeDiploidFSTMat(tempdata,locusNamesTemp,popNamesTemp)
MakeDiploidFSTMat = function(SNPmat,locusNames,popNames){
# SNPmat is a matrix with individuals in rows and snps in columns
# 0, 1, or 2 represent the number of copies of the focal allele, and 9 is for missing data
# locusNames is a character vector of names of each SNP
# popNames is a character vector with the population identifier for each individual
locusname <- unlist(locusNames)
popname <- unlist(popNames)
### Check that SNPmat has appropriate values (0, 1, 2, or 9, only)
snplevs <- levels(as.factor(unlist(SNPmat)))
ls <- paste(snplevs, collapse="")
if(ls!="012" & ls!="0129"){print("Error: Your snp matrix does not have 0,1, and 2"); break}
### Checking that locusNames and popNames have the same lengths as the columns and rows of SNPmat
if(dim(SNPmat)[1]!=length(popname) ){
print("Error: your population names do not match your SNP matrix")
break}
if(dim(SNPmat)[2]!=length(locusname)){
print("Error:  your locus names do not match your SNP matrix")
break}
writeLines("Calculating FSTs, may take a few minutes...")
nloci <- length(locusname)
FSTmat <- matrix(NA, nrow=nloci, ncol=8)
for (i in 1:nloci){
FSTmat[i]=getFSTs_diploids(popname,SNPmat[,i])
if (i==10000){print(paste(i, "done of", nloci))}
}
outTemp=as.data.frame(FSTmat2)
return (cbind(locusname,outTemp))
}
MakeDiploidFSTMat(tempdata,locusNamesTemp,popNamesTemp)
MakeDiploidFSTMat = function(SNPmat,locusNames,popNames){
# SNPmat is a matrix with individuals in rows and snps in columns
# 0, 1, or 2 represent the number of copies of the focal allele, and 9 is for missing data
# locusNames is a character vector of names of each SNP
# popNames is a character vector with the population identifier for each individual
locusname <- unlist(locusNames)
popname <- unlist(popNames)
### Check that SNPmat has appropriate values (0, 1, 2, or 9, only)
snplevs <- levels(as.factor(unlist(SNPmat)))
ls <- paste(snplevs, collapse="")
if(ls!="012" & ls!="0129"){print("Error: Your snp matrix does not have 0,1, and 2"); break}
### Checking that locusNames and popNames have the same lengths as the columns and rows of SNPmat
if(dim(SNPmat)[1]!=length(popname) ){
print("Error: your population names do not match your SNP matrix")
break}
if(dim(SNPmat)[2]!=length(locusname)){
print("Error:  your locus names do not match your SNP matrix")
break}
writeLines("Calculating FSTs, may take a few minutes...")
nloci <- length(locusname)
FSTmat <- matrix(NA, nrow=nloci, ncol=8)
for (i in 1:nloci){
FSTmat[i,]=unlist(getFSTs_diploids(popname,SNPmat[,i]))
if (i==10000){print(paste(i, "done of", nloci))}
}
outTemp=as.data.frame(FSTmat2)
return (cbind(locusname,outTemp))
}
MakeDiploidFSTMat(tempdata,locusNamesTemp,popNamesTemp)
MakeDiploidFSTMat = function(SNPmat,locusNames,popNames){
# SNPmat is a matrix with individuals in rows and snps in columns
# 0, 1, or 2 represent the number of copies of the focal allele, and 9 is for missing data
# locusNames is a character vector of names of each SNP
# popNames is a character vector with the population identifier for each individual
locusname <- unlist(locusNames)
popname <- unlist(popNames)
### Check that SNPmat has appropriate values (0, 1, 2, or 9, only)
snplevs <- levels(as.factor(unlist(SNPmat)))
ls <- paste(snplevs, collapse="")
if(ls!="012" & ls!="0129"){print("Error: Your snp matrix does not have 0,1, and 2"); break}
### Checking that locusNames and popNames have the same lengths as the columns and rows of SNPmat
if(dim(SNPmat)[1]!=length(popname) ){
print("Error: your population names do not match your SNP matrix")
break}
if(dim(SNPmat)[2]!=length(locusname)){
print("Error:  your locus names do not match your SNP matrix")
break}
writeLines("Calculating FSTs, may take a few minutes...")
nloci <- length(locusname)
FSTmat <- matrix(NA, nrow=nloci, ncol=8)
for (i in 1:nloci){
FSTmat[i,]=unlist(getFSTs_diploids(popname,SNPmat[,i]))
if (i==10000){print(paste(i, "done of", nloci))}
}
outTemp=as.data.frame(FSTmat2)
return (cbind(locusname,outTemp))
}
MakeDiploidFSTMat(tempdata,locusNamesTemp,popNamesTemp)
MakeDiploidFSTMat = function(SNPmat,locusNames,popNames){
# SNPmat is a matrix with individuals in rows and snps in columns
# 0, 1, or 2 represent the number of copies of the focal allele, and 9 is for missing data
# locusNames is a character vector of names of each SNP
# popNames is a character vector with the population identifier for each individual
locusname <- unlist(locusNames)
popname <- unlist(popNames)
### Check that SNPmat has appropriate values (0, 1, 2, or 9, only)
snplevs <- levels(as.factor(unlist(SNPmat)))
ls <- paste(snplevs, collapse="")
if(ls!="012" & ls!="0129"){print("Error: Your snp matrix does not have 0,1, and 2"); break}
### Checking that locusNames and popNames have the same lengths as the columns and rows of SNPmat
if(dim(SNPmat)[1]!=length(popname) ){
print("Error: your population names do not match your SNP matrix")
break}
if(dim(SNPmat)[2]!=length(locusname)){
print("Error:  your locus names do not match your SNP matrix")
break}
writeLines("Calculating FSTs, may take a few minutes...")
nloci <- length(locusname)
FSTmat <- matrix(NA, nrow=nloci, ncol=8)
for (i in 1:nloci){
FSTmat[i,]=unlist(getFSTs_diploids(popname,SNPmat[,i]))
if (i==10000){print(paste(i, "done of", nloci))}
}
outTemp=as.data.frame(FSTmat)
return (cbind(locusname,outTemp))
}
MakeDiploidFSTMat(tempdata,locusNamesTemp,popNamesTemp)
MakeDiploidFSTMat = function(SNPmat,locusNames,popNames){
# SNPmat is a matrix with individuals in rows and snps in columns
# 0, 1, or 2 represent the number of copies of the focal allele, and 9 is for missing data
# locusNames is a character vector of names of each SNP
# popNames is a character vector with the population identifier for each individual
locusname <- unlist(locusNames)
popname <- unlist(popNames)
### Check that SNPmat has appropriate values (0, 1, 2, or 9, only)
snplevs <- levels(as.factor(unlist(SNPmat)))
ls <- paste(snplevs, collapse="")
if(ls!="012" & ls!="0129"){print("Error: Your snp matrix does not have 0,1, and 2"); break}
### Checking that locusNames and popNames have the same lengths as the columns and rows of SNPmat
if(dim(SNPmat)[1]!=length(popname) ){
print("Error: your population names do not match your SNP matrix")
break}
if(dim(SNPmat)[2]!=length(locusname)){
print("Error:  your locus names do not match your SNP matrix")
break}
writeLines("Calculating FSTs, may take a few minutes...")
nloci <- length(locusname)
FSTmat <- matrix(NA, nrow=nloci, ncol=8)
for (i in 1:nloci){
FSTmat[i,]=unlist(getFSTs_diploids(popname,SNPmat[,i]))
if (i==10000){print(paste(i, "done of", nloci))}
}
outTemp=as.data.frame(FSTmat)
outTemp = cbind(locusname,outTemp)
colnames(outTemp)= c(LocusName,He, FST, T1, T2, FSTNoCorr,T1NoCorr,T2NoCorr,meanAlleleFreq)
return (outTemp)
}
MakeDiploidFSTMat(tempdata,locusNamesTemp,popNamesTemp)
MakeDiploidFSTMat = function(SNPmat,locusNames,popNames){
# SNPmat is a matrix with individuals in rows and snps in columns
# 0, 1, or 2 represent the number of copies of the focal allele, and 9 is for missing data
# locusNames is a character vector of names of each SNP
# popNames is a character vector with the population identifier for each individual
locusname <- unlist(locusNames)
popname <- unlist(popNames)
### Check that SNPmat has appropriate values (0, 1, 2, or 9, only)
snplevs <- levels(as.factor(unlist(SNPmat)))
ls <- paste(snplevs, collapse="")
if(ls!="012" & ls!="0129"){print("Error: Your snp matrix does not have 0,1, and 2"); break}
### Checking that locusNames and popNames have the same lengths as the columns and rows of SNPmat
if(dim(SNPmat)[1]!=length(popname) ){
print("Error: your population names do not match your SNP matrix")
break}
if(dim(SNPmat)[2]!=length(locusname)){
print("Error:  your locus names do not match your SNP matrix")
break}
writeLines("Calculating FSTs, may take a few minutes...")
nloci <- length(locusname)
FSTmat <- matrix(NA, nrow=nloci, ncol=8)
for (i in 1:nloci){
FSTmat[i,]=unlist(getFSTs_diploids(popname,SNPmat[,i]))
if (i==10000){print(paste(i, "done of", nloci))}
}
outTemp=as.data.frame(FSTmat)
outTemp = cbind(locusname,outTemp)
colnames(outTemp)= c("LocusName","He", "FST", "T1", "T2", "FSTNoCorr", "T1NoCorr", "T2NoCorr", "meanAlleleFreq")
return (outTemp)
}
MakeDiploidFSTMat(tempdata,locusNamesTemp,popNamesTemp)
MakeDiploidFSTMat(tempdata3,locusNamesTemp,popNamesTemp)
tempdata3[[1,4]]=9
MakeDiploidFSTMat(tempdata3,locusNamesTemp,popNamesTemp)
tempdata3
tempdata2
MakeDiploidFSTMat(tempdata,locusNamesTemp,popNamesTemp)
MakeDiploidFSTMat(tempdata2,locusNamesTemp,popNamesTemp)
locusNamesTemp=c("Locus1","Locus2","Locus3","Locus4","Locus5","Locus6")
popNamesTemp = read.csv("/Users/whitlock/Desktop/diploidtestPopNames.csv",header=FALSE)
tempdata=read.csv("/Users/whitlock/Desktop/dipFsttestfile.csv",header=FALSE)
tempdata2=tempdata
tempdata2[[6,1]]=9
tempdata3=tempdata
tempdata3[[1,4]]=9
tempdata
source('~/Desktop/OutFLANK/R/Fst Diploids.R', echo=TRUE)
source('~/Desktop/OutFLANK/R/Fst Diploids.R', echo=TRUE)
ReturnDiploidFSTmat(tempdata,locusNamesTemp,popNamesTemp)
locus1testdata=c(4,2,0,2,4,2,1,6,6,0,2,2)
dim(locus1testdata)=c(3,4)
locus1testdata=t(locus1testdata)
locus1testdata
WC_FST_FiniteSample_Diploids_2Alleles(locus1testdata)
WC_FST_FiniteSample_Diploids_2Alleles_NoCorr(locus1testdata)
WC_FST_Diploids_2Alleles(locus1testdata)
library(devtools)
library(roxygen2)
document("~/Desktop/OutFLANK")
install_github("whitlock/OutFLANK")
library(OutFLANK)
?OutFLANK
document("~/Desktop/OutFLANK")
document("~/Desktop/OutFLANK")
datatemp=read.csv("~/Desktop/katietest.csv", header=TRUE,stringsAsFactors=FALSE)
summary(datatemp)
install_github("whitlock/OutFLANK")
library(OutFLANK)
?OutFLANK
inputtemp=MakeDiploidFSTMat(datatemp)
locusnamestemp=names(datatemp)
locusnamestemp
popnamestemp=c(rep("A",1000),rep("B",200))
popnamestemp
inputtemp=MakeDiploidFSTMat(datatemp,locusnamestemp,popnamestemp)
inputtemp
datatemp[[1]]
datatemp[[,1]]
datatemp[,1]
datatemp[1,]
filename= "sam_sims/katie_m0.005_s1_n1000_N500_a0.01_r0.5_u0.0001_h22222long1_50000_01.qfreqtot"
Gen.LFMM <- read.table(gsub(".qfreqtot", ".genmat", filename), header=TRUE)
Gen.LFMM <- read.table("katie_m0.005_s1_n1000_N500_a0.01_r0.5_u0.0001_h22222long1_50000_01.genmat", header=TRUE)
Gen.LFMM <- read.table("/Desktop/katie_m0.005_s1_n1000_N500_a0.01_r0.5_u0.0001_h22222long1_50000_01.genmat", header=TRUE)
Gen.LFMM <- read.table("~/Desktop/katie_m0.005_s1_n1000_N500_a0.01_r0.5_u0.0001_h22222long1_50000_01.genmat", header=TRUE)
head(Gen.LFMM[,1:10])
dim(Gen.LFMM)
Pop=gl(2,1000)
par(mfrow=c(1,1))
OutFLANKResultsPlotter(outFST)
OutFLANKResultsPlotter(outFST, Zoom=TRUE, RightZoomFraction = 0.2)
FSTmat <- MakeDiploidFSTMat(Gen.LFMM, as.character(1:1000), Pop)
head(FSTmat)
outFST <- OutFLANK(FSTmat, NumberOfSamples=2)
str(outFST)
outFST$FSTbar # mean FST
outFST$T1
sum(outFST$results$T1)
which(outFST$OutlierFlag==TRUE)
sum(outFST$results$T1, na.rm=TRUE)/sum(outFST$results$T2, na.rm=TRUE)
outFST$FSTbar # mean FST
which(outFST$results$qvalues<0.05)
which(outFST$results$He<0.05)
out$FST$results[272]
outFST$results[272]
outFST$results[[272]]
outFST$results
outFST$results[,272]
outFST$results[272,]
which(outFST$results$OutlierFlag==TRUE)
which(outFST$results$qvalues<0.05)
iters = 1000
sd = 2
precision = 20
results = rep(0,iters)
for(i in 1:iters) {
x = floor(rnorm(20,5,sd) %% 10)
results[i] = paste(c('.',x),sep="",collapse="")
}
results = as.numeric(results)
plot(density(results,bw=.01),col="blue",lwd=3,bty="n")
iters = 1000
sd = 2
precision = 20
results = rep(0,iters)
for(i in 1:iters) {
x = floor(rnorm(20,5,sd) %% 10)
results[i] = paste(c('.',x),sep="",collapse="")
}
results = as.numeric(results)
plot(density(results,bw=.01),col="blue",lwd=3,bty="n")
iters = 1000
sd = 2
precision = 20
results = rep(0,iters)
for(i in 1:iters) {
x = floor(rnorm(20,5,sd) %% 10)
results[i] = paste(c('.',x),sep="",collapse="")
}
results = as.numeric(results)
plot(density(results,bw=.01),col="blue",lwd=3,bty="n")
iters = 1000
sd = 2
precision = 20
results = rep(0,iters)
for(i in 1:iters) {
x = floor(rnorm(20,5,sd) %% 10)
results[i] = paste(c('.',x),sep="",collapse="")
}
results = as.numeric(results)
plot(density(results,bw=.01),col="blue",lwd=3,bty="n")
oor(rnorm(20,5,2) %% 10)
floor(rnorm(20,5,2) %% 10)
? %%
152%%10
?rnorm
8%%3
9%%3
10%%3
iters = 1000
sd = 2
precision = 20
results = rep(0,iters)
for(i in 1:iters) {
x = floor(rnorm(20,5,sd) %% 10)
results[i] = paste(c('.',x),sep="",collapse="")
}
results
floor(rnorm(20,5,2) %% 10)
floor(rnorm(20,5,2) %% 10)
floor(rnorm(20,5,2) %% 10)
xx=floor(rnorm(20,5,2) %% 10)
xx
paste(c('.',xx),sep="",collapse="")
iters = 1000
sd = 2
precision = 20
results = rep(0,iters)
for(i in 1:iters) {
x = floor(rnorm(20,5,sd) %% 10)
results[i] = paste(c('.',x),sep="",collapse="")
}
results = as.numeric(results)
plot(density(results,bw=.01),col="blue",lwd=3)
plot(density(results,bw=.01),col="blue",lwd=3,bty="n")
falsepos=read.csv("~/Desktop/ Get to Work!/Papers/Papers in progress/Trimming - OutFLANK/figures/comparative false positive.csv", header=TRUE,stringsAsFactors=FALSE)
falsepos=as.data.frame(falsepos)
barplot(falsepos$IM, names.arg = falsepos$Method, ylab = "False positive rate",las = 1, col = "firebrick", ylim = c(0,0.016), axes=TRUE,width=.25)
barplot(falsepos$IBD, names.arg = falsepos$Method, ylab = "False positive rate",las = 1, col = "firebrick", ylim = c(0,0.016), axes=TRUE,width=.25)
barplot(as.matrix(cbind(falsepos$IM,falsepos$IBD,falsepos$X1R,falsepos$X2R)), col = c("darkolivegreen","goldenrod1","darkkhaki","darkorchid4"),beside = TRUE,space = c(0.005, 0.8),xlab = "Island model                     IBD                         One refuge expansion                         Two refuge expansion",  ylab = "False positive rate", las=1,legend.text = falsepos$Method)
power=read.csv("~/Desktop/ Get to Work!/Papers/Papers in progress/Trimming - OutFLANK/figures/comparative power.csv", header=TRUE,stringsAsFactors=FALSE)
power=as.data.frame(power)
barplot(as.matrix(cbind(power$IM,power$IBD,power$X1R,power$X2R)), col = c("darkolivegreen","goldenrod1","darkkhaki","darkorchid4"),beside = TRUE,space = c(0.005, 0.8),xlab = "Island model                     IBD                         One refuge expansion                         Two refuge expansion",  ylab = "Power", las=1,legend.text = power$Method)
?OutFLANK
?OutFLANK
install_github("whitlock/OutFLANK")
library(OutFLANK)
?OutFLANK
?OutFLANK
library(OutFLANK)
?OutFLANK
falsepos=read.csv("~/Desktop/ Get to Work!/Papers/Papers in progress/Trimming - OutFLANK/figures/comparative false positive.csv", header=TRUE,stringsAsFactors=FALSE)
falsepos=as.data.frame(falsepos)
barplot(falsepos$IM, names.arg = falsepos$Method, ylab = "False positive rate",las = 1, col = "firebrick", ylim = c(0,0.016), axes=TRUE,width=.25)
barplot(falsepos$IBD, names.arg = falsepos$Method, ylab = "False positive rate",las = 1, col = "firebrick", ylim = c(0,0.016), axes=TRUE,width=.25)
barplot(as.matrix(cbind(falsepos$IM,falsepos$IBD,falsepos$X1R,falsepos$X2R)), col = c("darkolivegreen","goldenrod1","darkkhaki","darkorchid4"),beside = TRUE,space = c(0.005, 0.8),xlab = "Island model                     IBD                         One refuge expansion                         Two refuge expansion",  ylab = "False positive rate", las=1,legend.text = falsepos$Method)
###  power
power=read.csv("~/Desktop/ Get to Work!/Papers/Papers in progress/Trimming - OutFLANK/figures/comparative power.csv", header=TRUE,stringsAsFactors=FALSE)
power=as.data.frame(power)
barplot(as.matrix(cbind(power$IM,power$IBD,power$X1R,power$X2R)), col = c("darkolivegreen","goldenrod1","darkkhaki","darkorchid4"),beside = TRUE,space = c(0.005, 0.8),xlab = "Island model                     IBD                         One refuge expansion                         Two refuge expansion",  ylab = "Power", las=1,legend.text = power$Method)
falsepos
power
alleleFreqListMaker = function(pbar = 0.25, bottleneckN = 10 , numPops){
ptemp=rnorm(n = numPops,mean = pbar,sd = sqrt(pbar*(1-pbar)/(2*bottleneckN)))
ptemp=ifelse(ptemp<0, 0, ptemp)
ptemp=ifelse(ptemp>1, 1, ptemp)
ptemp
}
rowMaker = function(popNumber, p, nPerPop){
#p is the allele frequnecy in this populaiton; nPerPop is the sample size of diploid individuals
NumAA = rbinom(n=1, size =  nPerPop, prob = p*p)
NumAa = rbinom(n = 1, size = nPerPop-NumAA, prob = 2*p*(1-p)/(2*p*(1-p)+(1-p)*(1-p)))
c(NumAA,NumAa,nPerPop - NumAA - NumAa)
}
test=array(dim=c(10000,3))
for(i in 1:10000) test[i,]=rowMaker(i,.25,12)
#works.
countSampler = function(plist, nPerPop){
numPops=length(plist)
out = array(dim=c(numPops,3))
for(i in 1:numPops) out[i,] = rowMaker(i,plist[i],nPerPop)
out
}
FstFromplist = function(plist,nPerPop){
data=countSampler(plist,nPerPop)
WC_FST_FiniteSample_Diploids_2Alleles(data)$FST
}
FstFromplistVECT = function(plist,nPerPop){
data=countSampler(plist,nPerPop)
WC_FST_FiniteSample_Diploids_2Alleles(data)
}
######
FstFrompbar = function(pbar = 0.25, bottleneckN = 10 , numPops,nPerPop){
plist = alleleFreqListMaker(pbar = pbar, bottleneckN = bottleneckN , numPops)
FstFromplist(plist,nPerPop)
}
FstFrompbarVECT = function(pbar = 0.25, bottleneckN = 10 , numPops,nPerPop){
plist = alleleFreqListMaker(pbar = pbar, bottleneckN = bottleneckN , numPops)
FstFromplistVECT(plist,nPerPop)
}
###########
#tests
reps=1000
temp=rep(1:reps)
for(j in 1:reps) temp[j]=FstFrompbar(pbar=0.25, numPops = 20, nPerPop = 20)
mean(temp)
var(temp/mean(temp))