サンプリングによる分布の変化の考察

背景

遺伝子シーケンサーからとれたリードデータをサンプリングしたとき、どのように元々の遺伝子カウントの分布が変わるかを考える。

例

例1

４つのリードを２つにダウンサンプルする。このうち geneA の数を 2とする
$A_1,A_2,B,C$から２つとる組み合わせは
\begin{array}{l} A_1 A_2 & A_1B & A_1 C\\ A_2B & A_2C\\ B C \end{array}
このうち$A$が2つ含まれるのは$A_1A_2$の1通り
$A$が1つ含まれるのは$A_1B,A_1C,A_2B,A_2C$の4通り
$A$が0つ含まれるのは$BC$の1通り

例２

6つのリードを２つにダウンサンプルする。このうち geneA の数を3とする
$A_1,A_2,A_3,B,C,D$から２つとる組み合わせは
\begin{array}{l} A_1 A_2&A_1 A_3&A_1 B&A_1 C&A_1 D\\ A_2A_3&A_2B&A_2C&A_2D\\ A_3B& A_3 C & A_3 D \\ B C & B D \\ C D \end{array}
このうち
\begin{array}{l} A^2 : A_1A_2,A_1A_3,A_2 A_3 : 3 \\ A^1 : A_1B,A_2B,A_3B,A_1C,A_2C,A_3C,A_1D,A_2D,A_3D :9\\ A^0 : BC,BD,CD : 3 \end{array}

母関数

\begin{array}{l} (1+Ax)(1+Bx)(1+Cx) =\\ 1+\\ (A+B+C)x+\\ (AB+AC+CB)x^2+\\ ABCx^3\\ \\\ (1+Ax)(1+Bx)(1+Cx)(1+Dx)=\\ 1+\\ (A+B+C+D)x+\\ (AB+AC+AD+BC+BD+CD)x^2+\\ (ABC+ABD+ACD+BCD)x^3+\\ ABCDx^4 \end{array}
$A,B,C,D$から２つとる組み合わせは$(1+Ax)(1+Bx)(1+Cx)(1+Dx)$の$x^2$の係数に現れる:$AB,AC,AD,BC,BD,CD$
$A=B=C=D=1$とすると$(1+x)^4=1+4x+6x^2+4x^3+1$
$x^2$の係数は6となる。

母関数による表現

例1

$A=B=p,C=D=1$とおく
\begin{eqnarray} (1+px)^2(1+x)^2&=&(1+2px+p^2x^2)(1+2x+x^2)\\ &=&1+2(p+1)x+(p^2+4p+1)x^2+2p(p+1)x^3+p^2x^4 \end{eqnarray}
この式の$x^2$の係数は$p^2+4p+1$これの$p$の係数が例1を表している。

例２

$A=B=C=p,D=E=F=1$とおく
\begin{eqnarray} (1+px)^3(1+x)^3=(1+3px+3p^2x^2+p^3x^3)(1+3x+3x^2+x^3)\\ \end{eqnarray}
\begin{array}{ccccccc} 1&3p&3p^2&p^3\\ &3&9p&9p^2&3p^3\\ &&3&9p&9p^2&3p^3\\ &&&1&3p&3p^2&p^3\\ \end{array}
なので
\begin{array}{l} (1+px)^3(1+x)^3=\\ 1+\\ 3(p+1)x+\\ 3(p^2+3p+1)x^2+\\ (p^3+9p^2+9p+1)x^3+\\ 3p(p^2+3p+1)x^4+\\ 3p^2(p+1)x^5+\\ p^3x^6 \end{array}
これの$x^2$の項は$3p^2+9p+3$この式は例２を表している。

母関数によるサンプリングの表現

総リード$N$にはgeneAが$A$含まれている。これからダウンサンプルしてリード数$n$とってくる。このうちgeneAが$a$個含まれる場合の数を考える。
これは$(1+p x)^A(1+x)^{N-A}$の$x^n$の係数から分かる。この係数は$p$の多項式になり、$p^0$の係数がgeneAが0つ含まれる数、$p^1$の係数がgeneAが1つ含まれる数、等々である。

微分による$x^n$の係数の求め方

\begin{eqnarray} f(x)&=&a+bx+cx^2+dx^3+ex^4+..\\ f'(x)&=&b+2cx+3dx^2+4ex^3+..\\ f''(x)&=&2c+6dx+12ex^2+..\\ f'''(x)&=&6d+24ex+.. \end{eqnarray}

\begin{eqnarray} a&=&f(0)\\ b&=&f'(0)\\ c&=&f''(0)/2\\ d&=&f'''(0)/6\\ .. \end{eqnarray}
なので$x^n$の係数は$f^{(n)}(0)/n!$

べき乗式の微分

$P=(1+p x)^A$とすると
\begin{eqnarray} P&=&(1+p x)^A\\ P'&=&Aa(1+p x)^{A-1}\\ P''&=&A(A-1)p^2(1+p x)^{A-2}\\ P'''&=&A(A-1)(A-2)p^3(1+p x)^{A-3}\\ ..\\ P^{(s)}&=&\frac{A!p ^s}{(A-s)!}(1+p x)^{A-s}\\ \end{eqnarray}
また
\begin{eqnarray} P^{(s)}(0) = \frac{A!p ^s}{(A-s)!} \end{eqnarray}

$Q=(1+x)^{N-A}$とすると
\begin{eqnarray} Q&=&(1+x)^{N-A}\\ Q'&=&(N-A)(1+x)^{N-A-1}\\ Q''&=&(N-A)(N-A-1)(1+x)^{N-A-2}\\ Q'''&=&(N-A)(N-A-1)(N-A-2)(1+x)^{N-A-3}\\ ..\\ Q^{(t)}&=&\frac{(N-A)!}{(N-A-t)!}(1+x)^{N-A-t}\\ \end{eqnarray}
また
\begin{eqnarray} Q^{(t)}(0)=\frac{(N-A)!}{(N-A-t)!} \end{eqnarray}

積の微分

\begin{eqnarray} (PQ)'&=&P'Q+PQ'\\ (PQ)''&=&P''Q+2P'Q'+PQ''\\ (PQ)'''&=&P'''Q+3P''Q+3PQ''+PQ'''\\ ..\\ (PQ)^{(n)}&=&P^{(n)}Q+nP^{(n-1)}Q'+\dbinom{n}{2}P^{(n-2)}Q''+..+\dbinom{n}{s}P^{(s)}Q^{(n-s)}+.. \end{eqnarray}
$x^n$の係数は
\begin{eqnarray} (PQ)^{(n)}(0)/n! \end{eqnarray}
これまで結果より
\begin{eqnarray} \frac{A!}{n!}\left(\frac{p^n}{(A-n)!}+n\frac{(N-A)p^{n-1}}{(A-n+1)!}+\dbinom{n}{2}\frac{(N-A)(N-A-1)p^{n-2}}{(A-n+2)!}+..\right) \end{eqnarray}
これの$p^a$の係数は
\begin{eqnarray} &&\dbinom{n}{a}\frac{A!(N-A)!}{n!(A-a)!(N-A-n+a)!}\\ &=&\frac{n!}{(n-a)!a!}\frac{A!}{n!(A-a)!}\frac{(N-A)!}{(N-A-n+a)!}\\ &=&\frac{A!}{(A-a)!a!}\frac{(N-A)!}{(N-A-n+a)!(n-a)!}\\ &=&\dbinom{A}{a}\dbinom{N-A}{n-a} \end{eqnarray}

例1

$N=4,~n=2,~A=2$
\begin{array}{l} a=0:&\dbinom{2}{0}\dbinom{2}{2}=1\\ a=1:&\dbinom{2}{1}\dbinom{2}{1}=4\\ a=2:&\dbinom{2}{2}\dbinom{2}{0}=1 \end{array}

例2

$N=6,~n=2,~A=3$
\begin{array}{l} a=0:&\dbinom{3}{0}\dbinom{3}{2}=3\\ a=1:&\dbinom{3}{1}\dbinom{3}{1}=9\\ a=2:&\dbinom{3}{2}\dbinom{3}{0}=3 \end{array}

組み合わせ論的考察

$N$から$n$を選ぶ。その$n$のうち$a$個を$A$から選び、その残りの$n-a$を$N-A$から選ぶ。なので
\begin{eqnarray} \dbinom{A}{a}\dbinom{N-A}{n-a} \end{eqnarray}
これは超幾何分布に関係するらしい

一般化

総リード$N$にgeneAが$A$,geneBが$B$,geneCが$C$,...含まれている。これから$n$リードとってくる。そのリードにそれぞれのgeneが$a,b,c,..$含まれる場合の数を考える。母関数は
\begin{eqnarray} (1+\alpha x)^A(1+\beta x)^B(1+\gamma x)^C ...(1+\zeta x)^Z \end{eqnarray}
ここで$A+B+C+..+Z=N$
これの$x^n$の係数の$\alpha^a \beta^b \gamma^c ...\zeta^z$の係数が求めるものである。

例

$N=8$ から$n=4$ サンプリングする。$A=2,B=4,C=2$のとき $a=1,b=2,c=1$となる場合の数を求める
\begin{eqnarray} &&(1+\alpha x)^2(1+\beta x)^4 (1+ \gamma x)^2 \\ &=& 1 +\\ &&(2\alpha+4\beta+2\gamma)x +\\ &&(\alpha^2+6\beta^2+\gamma^2+8\alpha\beta+4\alpha\gamma+8\beta\gamma)x^2+\\ &&(4\beta^3+4\alpha^2\beta+2\alpha^2\gamma+12\alpha\beta^2+12\beta^2\gamma+2\alpha\gamma^2+4\beta\gamma^2+16\alpha\beta\gamma)x^3+\\ &&(\beta^4+8\alpha\beta^3+8\gamma\beta^3+6\alpha^2\beta^2+\alpha^2\gamma^2+6\beta^2\gamma^2+8\alpha^2\beta\gamma+24\alpha\beta^2\gamma+8\alpha\beta\gamma^2)x^4+\\ &&(2\alpha\beta^4+2\beta^4\gamma+4\alpha^2\beta^3+4\beta^3\gamma^2+16\alpha\beta^3\gamma+12\alpha^2\beta^2\gamma+12\alpha\beta^2\gamma^2+4\alpha^2\beta\gamma^2)x^5+\\ && (\alpha^2\beta^4+\beta^4\gamma^2+4\alpha\beta^4\gamma+8\alpha^2\beta^3\gamma+8\alpha\beta^3\gamma^2+6\alpha^2\beta^2\gamma^2)x^6+\\ && (2\alpha^2\beta^4\gamma+2\alpha\beta^4\gamma^2+4\alpha^2\beta^3\gamma^2)x^7+\\ && \alpha^2\beta^4\gamma^2 \end{eqnarray}
$x^4$の係数は
\begin{eqnarray} \beta^4+8\alpha\beta^3+8\gamma\beta^3+6\alpha^2\beta^2+\alpha^2\gamma^2+6\beta^2\gamma^2+8\alpha^2\beta\gamma+24\alpha\beta^2\gamma+8\alpha\beta\gamma^2 \end{eqnarray}
なので$a,b,c$の場合の数は以下のようになる
\begin{array}{|ccc|c} a&b&c\\ \hline 0&4&0&1\\ 1&3&0&8\\ 0&3&1&8\\ 2&2&0&6\\ 2&0&2&1\\ 0&2&2&6\\ 2&1&1&8\\ 1&2&1&24\\ 1&1&2&8\\ \end{array}
よって元の分布を反映する$a=1,b=2,c=1$になる場合の数が最も大きいことが分かる

母関数

\begin{eqnarray} P(x)&=&(1+\alpha x)^A(1+\beta x)^B .. (1+\zeta x)^Z\\ &=& P_AP_B..P_Z \end{eqnarray}

積の微分

微分作用素を$D$とする。積の微分は次のように考えられる。
\begin{eqnarray} D(P_AP_B..P_Z)=(D_A+D_B+~..~D_Z)(P_AP_B..P_Z) \end{eqnarray}
ここで$D_X$は$P_X$にだけ作用する微分作用素である。
\begin{eqnarray} (D_A+D_B+~..~D_Z)^n \end{eqnarray}
の$D_A^{a}D_B^b..D_Z^z$の係数は多項係数
\begin{eqnarray} \dbinom{n!}{a!,b!,..,z!} = \frac{n!}{a!b!,..z!} \end{eqnarray}
である。また
\begin{eqnarray} P_A^{(n)}&=&\frac{A!\alpha^n}{(A-n)!}(1+\alpha x)^{A-n}\\ P_A^{(n)}(0)&=&\frac{A!\alpha^n}{(A-n)!}\\ \end{eqnarray}
よって$D_A^n$から$\alpha^n$の項が生じる。$P$の$x^n$の係数は$P^{(n)}(0)/n!$なので、$P$の$x^n$の係数の$\alpha^a\beta^n..\zeta^z$の係数は
\begin{eqnarray} &&\frac{1}{n!}\frac{n!}{a!b!..z!} \frac{A!}{(A-a)!}\frac{B!}{(B-b)!}..\frac{Z!}{(Z-z)!}\\ &=&\dbinom{A}{a}\dbinom{B}{b}..\dbinom{C}{c} \end{eqnarray}
これは多変量超幾何分布に関係するらしい。

例

$A=2,B=4,C=2$,$a=1,b=2,c=1$のとき
\begin{eqnarray} \dbinom{2}{1}\dbinom{4}{2}\dbinom{2}{1}=24 \end{eqnarray}
正しく例を表している。