population genetics

Population Genetics Chpt 1 Box

Chapter 1 Box Problems

Box B

A 9:3:3:1 test ratio is a total of 16 possibilities

1
2
3
(factorial(16)/(factorial(9)*factorial(3)*factorial(3)*factorial(1)))*(9/16)^9*(3/16)^3*(3/16)^3*(1/16)^1
## [1] 0.02452135

Ab/aB with r= 0.2 will produce nonrecombinants 80% of the time and recombinants 20% of the time.

1
2
3
4
5
6
7
8
9
10
11
12
13

Genotype <- c("Ab","aB","AB","ab")
Proportion <- c(4,4,1,1)
tbl1 <- cbind( Genotype, Proportion)

cat(hwrite(tbl1,
           border=NA,
           table.class="t1",
           row.class=list(c("header col","header col"),
                          c("col","col"),
                          c("col","col"),
                          c("col","col"),
                          c("col","col"))))
GenotypeProportion
Ab4
aB4
AB1
ab1
1
2
3
4

(factorial(10)/(factorial(4)*factorial(4)*factorial(1)*factorial(1)))*(4/10)^4*(4/10)^4*(1/10)^1*(1/10)^1

[1] 0.04128768
## Box C a 
1
2
3
4
5
N0 <-100
r <- 0.02
t <- c(50, 100, 150, 200, 250)
Nt <- N0*((1+r)^t)
plot(t, Nt)
## Box C b 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
N0 <-100
r <- 0.02
t <- c(50, 100, 150, 200, 250)
ra <- log(1+r)
rb <- r
rc <- r -( (r^2)/2 )
 
Na <- N0*(exp(ra*t))
Nb <- N0*(exp(rb*t))
Nc <- N0*(exp(rc*t))

Na

[1]   269.1588   724.4646  1949.9603  5248.4897 14126.7721

Nb

[1]   271.8282   738.9056  2008.5537  5459.8150 14841.3159

Nc

[1]   269.1234   724.2743  1949.1920  5245.7326 14117.4964

plot(t, Na, col="red", ylab="Nt")
points( t, Nb, col="black")
points( t, Nc, col="blue")

$$\frac{dN}{[N(1-\frac{N}{K})]}=rdt$$

## Box C c

Given that:

$$\int \frac{1}{x(1+ax)} dx = ln \frac{x}{1+ax}$$

then:

$$\int_{N_0}^{N_t} \frac{1}{N(1- \frac{1}{K}N)} dN = ln \int_{t=0}^{t}r dt$$

Solve to get:

$$ln \frac{N_t}{1- \frac{1}{K} N_t} - ln \frac{N_0}{1- \frac{1}{K}N_0} = rt$$

ln(x) - ln(y) = ln(x/y) = ln(x*1/y)

$$ln \frac{N_t}{1- \frac{1}{K} N_t} . ln \frac{1- \frac{1}{K}N_0}{N_0} = rt$$

Inverse natural log of both sides:

$$\frac{N_t}{N_0} . \frac{1- \frac{1}{K}N_0}{1-\frac{1}{K}N_t} = e^{rt} $$

$$\frac{N_t - \frac{N_tN_0}{K}}{N_0 - \frac{N_tN_0}{K}} = e^{rt}$$

Multiply the left side by K/K

$$\frac{N_tK - N_tN_0}{N_0K - N_tN_0} = e^{rt}$$

$$\frac{N_t(K-N_0)}{N_0(K-N_t)} = e^{rt} $$

$$N_0e^{rt}(K-N_t) = N_t(K-N_0)$$

$$\frac{N_0e^{rt}}{K-N_0} = \frac{N_t}{K-N_t}$$

## Box C d 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33

N.exact <- vector(mode="numeric", length=6)
N.exact[1] <- 1000
r <- 0.1
K <- 2000

for( t in 2:6){
     N.exact[t] <- r*N.exact[t-1]*((K - N.exact[t-1])/K) + N.exact[t-1]
}

N.exact

[1] 1000.000 1050.000 1099.875 1149.376 1198.261 1246.295

N.approx <- vector(mode="numeric", length=6)
N.approx[1] <- 1000
C <- (K-N.approx[1])/N.approx[1]

for( t in 2:6){
     N.approx[t] <- K/(1 + C*exp(-r*t))
}

N.approx

[1] 1000.000 1099.668 1148.885 1197.375 1244.919 1291.313



N.diff <- N.approx - N.exact

N.diff

[1]  0.00000 49.66799 49.01003 47.99907 46.65808 45.01739
## Box D preamble

With statistical software such as R, the manipulations described in Box D are unnecessary except to educate you as to what is going on behind the scenes. Here is one method to perform the manipulations as described in the box:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
x <- c(1,2,3,4)
y <- c(3,7,6,9)

x.bar <- mean(x) # the mean of the X values
x.bar

[1] 2.5


y.bar <- mean(y) # the mean of the Y values
y.bar

[1] 6.25

x2.bar <- mean( x^2) # the mean of the X squared values
x2.bar

[1] 7.5

xy.bar <- mean( x*y ) #mean of the sum of the xy products
xy.bar

[1] 17.75

xy.cov <- xy.bar - (x.bar*y.bar) #covariance of x and y
xy.cov

[1] 2.125

x.var <- x2.bar - x.bar^2 # variance of x
x.var

[1] 1.25

m <- xy.cov/x.var # slope
m

[1] 1.7

b <- y.bar - m*x.bar # y intercept
b 

[1] 2


# So the best fit line is y = mx + b or
# y = 1.7x + 2

# here is what you will normally do using R
# using the original x, y vectors from above

model <- lm(y~x)
summary(model)

Call:
lm(formula = y ~ x)

Residuals:
   1    2    3    4 
-0.7  1.6 -1.1  0.2 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)
(Intercept)   2.0000     1.7958   1.114    0.381
x             1.7000     0.6557   2.592    0.122

Residual standard error: 1.466 on 2 degrees of freedom
Multiple R-squared:  0.7707,	Adjusted R-squared:  0.656 
F-statistic: 6.721 on 1 and 2 DF,  p-value: 0.1221

##Now let's incorporate weights

x <- c(1,2,3,4)
y <- c(3,7,6,9)
w <- c(10,10,1,1)

w.sum <- sum(w) # sum of the weights
w.sum

[1] 22

xw.bar <- sum(x*w)/w.sum # the mean of the weighted X values
xw.bar

[1] 1.681818

yw.bar <- sum(y*w)/w.sum # the mean of the weighted Y values
yw.bar

[1] 5.227273

x2w.bar <- sum( x^2*w)/w.sum #mean of the sum of the weighted x squared values
x2w.bar

[1] 3.409091

xyw.bar <- sum( x*y*w )/w.sum #mean of the sum of the weighted xy products
xyw.bar

[1] 10.18182

xyw.cov <- xyw.bar - (xw.bar*yw.bar) #covariance of weighted x and y
xyw.cov

[1] 1.390496

xw.var <- x2w.bar - xw.bar^2 # variance of weighted x
xw.var

[1] 0.5805785

mw <- xyw.cov/xw.var # slope
mw

[1] 2.395018

bw <- yw.bar - mw*xw.bar # y intercept
bw 

[1] 1.199288

# So the best fit line to the weighted data is y = mx + b or
# y = 2.4x + 1.2

# lm has a weight option

model <- lm(y~x, weight= w)
summary(model)

Call:
lm(formula = y ~ x, weights = w)

Weighted Residuals:
     1      2      3      4 
-1.879  3.196 -2.384 -1.779 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)
(Intercept)   1.1993     1.7366   0.691    0.561
x             2.3950     0.9405   2.546    0.126

Residual standard error: 3.361 on 2 degrees of freedom
Multiple R-squared:  0.7643,	Adjusted R-squared:  0.6464 
F-statistic: 6.484 on 1 and 2 DF,  p-value: 0.1258

## Box D a & b

Obtaining quantities 8 and 9 (slope and y intercept) for problems a and b is now trivial using the built in R function lm

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
x <- c(2,4,6,8)
y <- c(11,21,22,32)

model <- lm(y~x)
summary(model)

Call:
lm(formula = y ~ x)

Residuals:
   1    2    3    4 
-0.9  2.7 -2.7  0.9 


Coefficients:
            Estimate Std. Error t value Pr(>|t|)  
(Intercept)   5.5000     3.4857   1.578   0.2553  
x             3.2000     0.6364   5.028   0.0373 *
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 2.846 on 2 degrees of freedom
Multiple R-squared:  0.9267,	Adjusted R-squared:   0.89 
F-statistic: 25.28 on 1 and 2 DF,  p-value: 0.03735


w <- c(1,1,20,20)

model <- lm(y~x, weight= w)
summary(model)


Call:
lm(formula = y ~ x, weights = w)

Weighted Residuals:
     1      2      3      4 
 4.021  5.913 -5.342  3.121 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)  
(Intercept)  -1.1288     5.4466  -0.207   0.8550  
x             4.0539     0.7854   5.162   0.0355 *
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 6.686 on 2 degrees of freedom
Multiple R-squared:  0.9302,	Adjusted R-squared:  0.8953 
F-statistic: 26.64 on 1 and 2 DF,  p-value: 0.03554
# Box E a $$0 = r_1N_1[ 1 - \frac{N_1 + \alpha_{21} N_2}{K_1} ]$$ $$0 = r_1N_1 - r_1N_1[\frac{N_1 + \alpha_{21} N_2}{K_1} ]$$ $$0 = r_1N_1 - \frac{r_1N_1N_1}{K_1} - \frac{r_1N_1 \alpha_{21} N_2}{K_1}$$ Subtract $r_1N_1$ from both sides and multiply by -1 $$r_1N_1= \frac{r_1N_1N_1}{K_1} + \frac{r_1N_1 \alpha_{21} N_2}{K_1}$$ Multiply by $\frac{K_1}{r_1N_1}$ $$K_1= N_1 + \alpha_{21} N_2$$ $$N_1= K_1 - \alpha_{21} N_2$$ Substitute for $N_1$ in $$0 = r_2N_2[ 1 - \frac{N_2 + \alpha_{12} N_1}{K_2} ]$$ $$0 = r_2N_2[ 1 - \frac{N_2 + \alpha_{12} (K_1 - \alpha_{21} N_2) }{K_2} ]$$ $$0 = r_2N_2 - r_2N_2[\frac{N_2 + \alpha_{12}K_1 - \alpha_{12}\alpha_{21} N_2 }{K_2} ]$$ $$0 = r_2N_2 + \frac{-r_2N_2N_2 - r_2N_2\alpha_{12}K_1 + r_2N_2\alpha_{12}\alpha_{21} N_2 }{K_2} $$ $$-r_2N_2K_2= -r_2N_2N_2 - r_2N_2\alpha_{12}K_1 + r_2N_2\alpha_{12}\alpha_{21} N_2$$ Divide by $-r_2N_2$
$$K_2= N_2 +\alpha_{12}K_1 - \alpha_{12}\alpha_{21} N_2$$ $$K_2 -\alpha_{12}K_1= N_2 - \alpha_{12}\alpha_{21} N_2$$ $$K_2 -\alpha_{12}K_1= N_2(1 - \alpha_{12}\alpha_{21})$$ $$\frac{K_2 -\alpha_{12}K_1}{1 - \alpha_{12}\alpha_{21}}= N_2$$
## Box E b 
1
2
3
4
5
6
7
8
9
K1 <-1500
K2 <- 2500
a21 <- 0.25
a12 <- 0.5

N1 <- (K1-a21*K2)/(1-a21*a12)
N2 <- (K2-a12*K1)/(1-a21*a12)

cat( "The value of N1 is:", N1)
1
## The value of N1 is: 1000
1
cat( "The value of N2 is:", N2)
1
## The value of N2 is: 2000

$r_1$ and $r_2$ affect the rate of achieving equilibrium, not the final population count.

Share