Computerintensive Methoden - Hausübung Markus

From StatWiki
Jump to: navigation, search

Nehmen wir an, es gäbe 2 dominante Marker A und B, sodass wir erkennen können, wenn mindestens ein Allel von Linie 1 in den Nachkommen vorliegt. Wir sind interessiert die Distanz zwischen den Markern (\delta\,) oder äquivalent die Rekombinationsrate (\rho \,) zu schätzen (\rho\, ist Schätzer von r\,)

Komplette Info

Übergangswahrscheinlichkeiten
BB Bb bb bB
AA  (1-r)^2\,  (1-r) r\,  r^2\,  r (1-r)\,
Aa  (1-r) r\,  (1-r)^2\,  r (1-r)\,  r^2\,
aa  r^2\,  r (1-r)\,  (1-r)^2\,  (1-r) r\,
aA  r (1-r) \,  r^2\,  (1-r) r\,  (1-r)^2\,

Nach Angabe sind AA, Aa, aA (A_) voneinander nicht zu unterscheiden, selbiges für BB, Bb, bB (B_)


F_2\, Resign???

2 ingezüchtete Linien, daher gentisch ident innerhalb, werden gekreuzt, um eine F_1\, Generation zu erhalten. Diese ist genetisch wieder einheitlich, weil sie ja ein Chromosom vom Vater und eines von der Mutter erhält.

Von dieser F_1\, Generation wird eine F_2\, Generation durch Bruder-Schwester Paarung erhalten, die nun genetische Variabilität zeigen kann.

\Rightarrow \, 4 Gruppen unterscheidbar, daher

Freq(A_, B_) = 30
Freq(A_, bb) = 5
Freq(aa, B_) = 2
Freq(aa, bb) = 10

ComputerIntensiveMethoden.Homework.Angabe.png

Lösung

Verteilungen

 \rho \sim \mathrm{Beta}(n_r+1, N-n_r+1)\,

 N_{A\star,B\star} \sim \mathrm{Multinomial} \left(n=30, 
\vec p = \begin{pmatrix} 
(1-r)^2 \\
(1-r) r \\
r (1-r) \\
(1-r) r \\
(1-r)^2 \\
r^2 \\
r (1-r) \\
r^2 \\
(1-r)^2
\end{pmatrix} \cdot a \right)\,

<axiom echo="true"> solve(1 = a*((1-r)^2 + (1-r)*r + r*(1-r) + (1-r)*r + (1-r)^2 + r^2 + r*(1-r) + r^2 + (1-r)^2),a) </axiom>

 N_{A\star,bb} \sim \mathrm{Multinomial} \left(n=5, 
\vec p = \begin{pmatrix} 
r^2 \\
r (1-r) \\
(1-r) r
\end{pmatrix} \cdot a \right)\,

<axiom echo="true"> solve(1 = a*(r^2 + r*(1-r) + (1-r)*r),a) </axiom>

 N_{aa,B\star} \sim \mathrm{Multinomial} \left(n=2, 
\vec p = \begin{pmatrix} 
r^2 \\
r (1-r) \\
r (1-r) 
\end{pmatrix} \cdot a \right)\,

<axiom echo="true"> solve(1 = a*r^2 + a*r*(1-r) + a*(1-r)*r,a) </axiom>

Um die Anzahl der Rekombinationen zu zählen wird die Matrix der absoluten Häufigkeiten mit der folgenden Matrix positionsweise multipliziert (also nicht Matrizenmultiplikation)


\begin{pmatrix}
0&1&2&1 \\
1&0&1&2 \\
2&1&0&1 \\
1&2&1&0
\end{pmatrix}
\,

Program

> pdf(rpdf)
>
> library(coda)
Loading required package: lattice
>
> gibbs = function(r, n=5e3) {
+         #wir haben vergessen, dass jedes Individuum diploid ist, also
+         #2 Rekombinationen repraesentiert!
+         #deshalb der 2er
+         N = 2*47
+
+         # wird mit absoluten frequenzen multipliziert
+         # um nr zu berechnen (also die anzahl der rekombinationen)
+         # Wie entsteht diese? z.b. AA->BB = 0, keine Rekombination, AA->Bb = 1, eine Rekombination
+         nr.key=matrix(c(
+                  0,1,2,1,
+                  1,0,1,2,
+                  2,1,0,1,
+                  1,2,1,0), byrow=T, ncol=4)
+
+         replicate(n,
+                 {
+                         # print(r)
+
+                         AxBx=rmultinom(1, 30,
+                          c((1-r)^2, (1-r)*r,     r*(1-r),
+                            (1-r)*r, (1-r)^2,     r^2,
+
+                            r*(1-r), r^2,         (1-r)^2)*(1/(r^2 - 2*r + 3)))
+
+                         Axbb=rmultinom(1, 5,
+                          c(                  r^2,
+                                              r*(1-r),
+
+                                              (1-r)*r)*(-1/(r^2-2*r)))
+
+                         aaBx=rmultinom(1, 2,
+                          c(r^2,
+                            r*(1-r),
+                            r*(1-r)) * (-1/(r^2 - 2*r)))
+
+                         abs.freq=matrix(c(
+                             AxBx[1:2], Axbb[1], AxBx[3],
+                             AxBx[4:5], Axbb[2], AxBx[6],
+                             aaBx[1:2], 10,      aaBx[3],
+                             AxBx[7:8], Axbb[3], AxBx[9]), byrow=T, ncol=4)
+
+                         # print(abs.freq)
+
+                         # Anzahl der Rekombinations
+                         nr = sum(abs.freq * nr.key)
+
+                         # print(nr)
+
+                         # wichtig <<- ersetzt global variable r
+                         r <<- rbeta(1, nr+1, N-nr+1)
+                         r
+                 })
+ }
>
> hist(sapply(seq(0.001,1,length.out=10), gibbs, 5000))