Skip to content

Commit 2d14ac1

Browse files
cjbatteycjbattey
committed
removed replicates, fixed pos for graphs, added multiple pop sizes and allele freqs
1 parent eab1960 commit 2d14ac1

File tree

3 files changed

+102
-62
lines changed

3 files changed

+102
-62
lines changed

‎dev.R

+26-6
Original file line numberDiff line numberDiff line change
@@ -9,14 +9,34 @@ binomialDraw <- function(n,p){
99
}
1010

1111
#main simulation function
12-
runPopSim <- function(gen=100,p=0.5,Waa=1,Wab=1,Wbb=1,n=100,nPop=2,m=0,stats=c("p","Fst"),Uab=0,Uba=0,infinitePop=F){
12+
runPopSim <- function(gen=100,p=0.5,Waa=1,Wab=1,Wbb=1,n=100,nPop=2,m=0,stats=c("p","Fst"),Uab=0,Uba=0,
13+
infinitePop=F,continue=F){
14+
if(continue==T){ #continue function currently broken... will update when possible.
15+
#allele.freq <- endState
16+
} else {
17+
allele.freq <- data.frame(matrix(ncol=3*nPop)) #initialize summary stat matrix
18+
if(length(p)==nPop){
19+
allele.freq[1,(1:nPop)] <- p
20+
} else {
21+
allele.freq[1,(1:nPop)] <- rep(p,nPop) #starting allele freqs
22+
}
23+
}
24+
1325
withProgress(message="simulating populations...",value=0,{
14-
allele.freq <- data.frame(matrix(ncol=3*nPop))
15-
allele.freq[1,(1:nPop)] <- rep(p,nPop) #starting allele freqs
1626
for(i in 1:gen){
17-
mean.p <- as.numeric(rowMeans(allele.freq[i,(1:nPop)]))
27+
if(length(n)==nPop){ #weight mean allele freq by relative population size if different
28+
ps <- allele.freq[i,(1:nPop)] %>% unlist()
29+
mean.p <- mean(ps*(n/sum(n)))
30+
} else {
31+
mean.p <- as.numeric(rowMeans(allele.freq[i,(1:nPop)]))
32+
}
1833
for(j in 1:nPop){
1934
p <- allele.freq[i,j]
35+
if(length(n)==nPop){ #get population-specific pop size if multiple listed
36+
n <- n[j]
37+
} else {
38+
n <- n
39+
}
2040
p <- (1 - Uab) * p + Uba * (1 - p) #mutation
2141
p <- p*(1-m)+m*mean.p # migration
2242
q <- 1-p
@@ -60,6 +80,7 @@ runPopSim <- function(gen=100,p=0.5,Waa=1,Wab=1,Wbb=1,n=100,nPop=2,m=0,stats=c("
6080
} #end populations loop
6181
incProgress(1/gen)
6282
} #end generations loop
83+
}) #end progress bar
6384
#summary stats
6485
names <- c()
6586
for(i in 1:nPop){names[i]<-paste0("p",i)}
@@ -68,15 +89,14 @@ runPopSim <- function(gen=100,p=0.5,Waa=1,Wab=1,Wbb=1,n=100,nPop=2,m=0,stats=c("
6889
colnames(allele.freq) <- names
6990
allele.freq$meanHo <- rowMeans(allele.freq[(nPop+1):(nPop*2)])
7091
allele.freq$meanHe <- rowMeans(allele.freq[(nPop*2+1):(nPop*3)])
71-
allele.freq$Fis <- abs(1-(allele.freq$meanHo/allele.freq$meanHe))
92+
allele.freq$Fis <- 1-(allele.freq$meanHo/allele.freq$meanHe)
7293
allele.freq$mean.p <- rowMeans(allele.freq[1:nPop])
7394
allele.freq$Hs <- rowMeans(allele.freq[(nPop*2+1):(nPop*3)])
7495
allele.freq$Ht <- 2*allele.freq$mean.p*(1 - allele.freq$mean.p)
7596
allele.freq$Fst <- (allele.freq$Ht-allele.freq$Hs)/allele.freq$Ht
7697
allele.freq$Fst[allele.freq$Fst<0] <- 0
7798
allele.freq$gen <- 0:gen
7899
return(allele.freq)
79-
})
80100
}
81101

82102
#format for plotting

‎server.R

+62-44
Original file line numberDiff line numberDiff line change
@@ -4,17 +4,31 @@ shinyServer(function(input,output,session){
44

55
source("./dev.R")
66

7+
p <- reactive({as.numeric(unlist(strsplit(input$p,",")))})
8+
n <- reactive({as.numeric(unlist(strsplit(input$n,",")))})
9+
10+
#run single sim
711
sim.data <- eventReactive(input$go,ignoreNULL = F,{
812
validate(
13+
need((input$nPop==length(p())|length(p())==1),"number of populations must equal number of starting allele frequencies."),
914
need(input$gen<5001,"Please select < 5000 generations."),
1015
need(input$nPop<101,"Please select < 100 populations"),
1116
need(input$n<1000001,"Please select n < 1,000,000"),
1217
need(input$plotStats!="","Select a variable to plot.")
1318
)
14-
runPopSim(gen=input$gen,p=input$p,Waa=input$Waa,Wab=input$Wab,Wbb=input$Wbb,n=input$n,
15-
nPop=input$nPop,m=input$m,stats=input$plotStats,infinitePop=input$infinitePop,Uab=input$Uab,Uba=input$Uab)
19+
tmp <- runPopSim(gen=input$gen,p=p(),Waa=input$Waa,Wab=input$Wab,Wbb=input$Wbb,n=n(),
20+
nPop=input$nPop,m=input$m,stats=input$plotStats,infinitePop=input$infinitePop,Uab=input$Uab,Uba=input$Uab,
21+
continue=F)
22+
tmp
1623
})
1724

25+
#continue sim
26+
# sim.data <- eventReactive(input$continue,ignoreNULL = F,{
27+
# runPopSim(gen=input$gen,p=input$p,Waa=input$Waa,Wab=input$Wab,Wbb=input$Wbb,n=input$n,
28+
# nPop=input$nPop,m=input$m,stats=input$plotStats,infinitePop=input$infinitePop,Uab=input$Uab,Uba=input$Uab,
29+
# continue=T)
30+
# })
31+
1832
plot.data <- eventReactive(sim.data(),{
1933
meltPlotData(allele.freq.df = sim.data(),gen=input$gen,nPop=input$nPop,stats=input$plotStats)
2034
})
@@ -34,50 +48,54 @@ shinyServer(function(input,output,session){
3448
nLost.text()
3549
})
3650

37-
output$endStateTable <- renderTable({
38-
endState <- sim.data()[input$gen,c("Fis","Hs","Ht","Fst")]
51+
endStateTable <- eventReactive(sim.data(),{
52+
#pNames <- c()
53+
#for(i in 1:input$nPop){pNames[i] <- paste0("p",i)}
54+
sim.data()[(input$gen+1),c("Fis","Hs","Ht","Fst")]
3955
})
4056

41-
sumTable <- eventReactive(input$runSim,{
42-
validate(
43-
need(input$n<=100000,"Please select n <= 100,000")
44-
)
45-
sumTable <- data.frame(matrix(ncol=14))
46-
withProgress(message="simulating populations...",value=0,{
47-
for(i in 1:100){
48-
df <- runPopSim2(gen=100,p=input$p,Waa=input$Waa,Wab=input$Wab,Wbb=input$Wbb,n=input$n,
49-
nPop=2,m=input$m,infinitePop=input$infinitePop,Uab=input$Uab,Uba=input$Uab)
50-
names(sumTable) <- names(df)
51-
sumTable[i,] <- df[nrow(df),]
52-
incProgress(1/100)
53-
}
54-
})
55-
sumTable
56-
})
57-
58-
meanTable <- reactive({
59-
tbl <- colMeans(sumTable(),na.rm=T)
60-
tbl <- tbl[c("Fis","Hs","Ht","Fst")]
61-
t(tbl)
62-
})
63-
64-
varTable <- reactive({
65-
tbl <- apply(sumTable(),2,function(e) var(e,na.rm=T))
66-
tbl <- tbl[c("Fis","Hs","Ht","Fst")]
67-
t(tbl)
57+
output$endStateTable <- renderTable({
58+
endStateTable()
6859
})
6960

70-
71-
72-
output$meanTable <- renderTable(meanTable(),colnames = T,digits=4,caption = "Mean state at final generation:",
73-
caption.placement = getOption("xtable.caption.placement", "top"),
74-
caption.width = getOption("xtable.caption.width", NULL))
75-
76-
output$varTable <- renderTable(varTable(),colnames = T,digits=4,caption = "Variance:",
77-
caption.placement = getOption("xtable.caption.placement", "top"),
78-
caption.width = getOption("xtable.caption.width", NULL))
79-
80-
output$sumTable <- renderTable(sumTable(),caption = "Final Generation States:",
81-
caption.placement = getOption("xtable.caption.placement", "top"),
82-
caption.width = getOption("xtable.caption.width", NULL))
61+
# sumTable <- eventReactive(input$runSim,{
62+
# validate(
63+
# need(input$n<=100000,"Please select n <= 100,000")
64+
# )
65+
# sumTable <- data.frame(matrix(ncol=14))
66+
# withProgress(message="simulating populations...",value=0,{
67+
# for(i in 1:100){
68+
# df <- runPopSim2(gen=100,p=input$p,Waa=input$Waa,Wab=input$Wab,Wbb=input$Wbb,n=input$n,
69+
# nPop=2,m=input$m,infinitePop=input$infinitePop,Uab=input$Uab,Uba=input$Uab)
70+
# names(sumTable) <- names(df)
71+
# sumTable[i,] <- df[nrow(df),]
72+
# incProgress(1/100)
73+
# }
74+
# })
75+
# sumTable
76+
# })
77+
#
78+
# meanTable <- reactive({
79+
# tbl <- colMeans(sumTable(),na.rm=T)
80+
# tbl <- tbl[c("Fis","Hs","Ht","Fst")]
81+
# t(tbl)
82+
# })
83+
#
84+
# varTable <- reactive({
85+
# tbl <- apply(sumTable(),2,function(e) var(e,na.rm=T))
86+
# tbl <- tbl[c("Fis","Hs","Ht","Fst")]
87+
# t(tbl)
88+
# })
89+
#
90+
# output$meanTable <- renderTable(meanTable(),colnames = T,digits=4,caption = "Mean state at final generation:",
91+
# caption.placement = getOption("xtable.caption.placement", "top"),
92+
# caption.width = getOption("xtable.caption.width", NULL))
93+
#
94+
# output$varTable <- renderTable(varTable(),colnames = T,digits=4,caption = "Variance:",
95+
# caption.placement = getOption("xtable.caption.placement", "top"),
96+
# caption.width = getOption("xtable.caption.width", NULL))
97+
#
98+
# output$sumTable <- renderTable(sumTable(),caption = "Final Generation States:",
99+
# caption.placement = getOption("xtable.caption.placement", "top"),
100+
# caption.width = getOption("xtable.caption.width", NULL))
83101
})

‎ui.R

+14-12
Original file line numberDiff line numberDiff line change
@@ -2,14 +2,14 @@ shinyUI(fluidPage(
22
titlePanel("driftR: Population Genetic Simulations in R"),
33
sidebarLayout(
44
sidebarPanel(
5-
numericInput("p","Starting allele frequency",value=0.5,min=0,max=1),
5+
textInput("p","Starting allele frequency",value="0.5"),
66
sliderInput("Uab","Mutation Rate",value=0,min=0,max=0.1),
77
sliderInput("Waa","Fitness of genotype AA",value=1,min=0,max=1),
88
sliderInput("Wab","Fitness of genotype AB",value=1,min=0,max=1),
99
sliderInput("Wbb","Fitness of genotype BB",value=1,min=0,max=1),
1010
sliderInput("m","Migration Rate",0,min=0,max=0.35),
11-
numericInput("n","Population Size",100,min=1,max=1e5),
1211
numericInput("nPop","Number of Populations",2,min=1,max=100),
12+
textInput("n","Population Size",value=100),
1313
numericInput("gen","Number of Generations",100,min=1,max=5000),
1414
checkboxInput("infinitePop","Infinite Population (no drift)",value = F),
1515
checkboxGroupInput(inputId="plotStats",label="plot:",choices=c("p","He","Hs","Ht","Fst"),inline=T,selected="p"),
@@ -20,18 +20,20 @@ shinyUI(fluidPage(
2020
Full code available on github: https://github.com/cjbattey/driftR"),style="font-size:75%")
2121
),
2222

23-
mainPanel(
23+
mainPanel(style="position:fixed;margin-left:32vw;",
2424
plotOutput("plot"),
2525
textOutput("nLost"),
26-
tableOutput("endStateTable"),
27-
div(style="border-top:1px solid black;"),
28-
helpText("Click the button below to run 100, 100-generation simulations of 2 populations using the current
29-
parameters."),
30-
actionButton("runSim","Run Replicate Simulations"),
31-
tableOutput("meanTable"),
32-
tableOutput("varTable"),
33-
div(tableOutput("sumTable"), style = "font-size: 75%; width: 75%;")
26+
div(style="height:5vh;"),
27+
div("Final Generation State:"),
28+
tableOutput("endStateTable")
29+
# actionButton("continue","Continue Simulation"),
30+
# div(style="border-top:1px solid black;"),
31+
#helpText("Click the button below to run 100, 100-generation simulations of 2 populations using the current
32+
# parameters.")
33+
#actionButton("runSim","Run Replicate Simulations"),
34+
#tableOutput("meanTable"),
35+
#tableOutput("varTable"),
36+
#div(tableOutput("sumTable"), style = "font-size: 75%; width: 75%;")
3437
)
3538
)
36-
3739
))

0 commit comments

Comments
 (0)