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

背景

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

例

例1

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

例２

6つのリードを２つにダウンサンプルする。このうち geneA の数を3とする
$A_1,A_2,A_3,B,C,D$から２つとる組み合わせは
$$\begin{array}{lllll}{A}_1{A}_2& A_1{A}_3& A_{1}B& A_{1}C& A_{1}D\\ A_2{A}_3& A_{2}B& A_{2}C& A_{2}D\\ A_{3}B& A_{3}C& A_{3}D\\ BC& BD\\ CD\end{array}$$
このうち
$$\begin{array}{l}{A}^2:A_1{A}_2,A_1{A}_3,A_2{A}_3:3\\ A^1:A_{1}B,A_{2}B,A_{3}B,A_{1}C,A_{2}C,A_{3}C,A_{1}D,A_{2}D,A_{3}D: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+\\ ABC{x}^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+\\ ABCD{x}^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{array}{rcl}(1+px{)}^2(1+x{)}^2& =& (1+2px+p^2{x}^2)(1+2x+x^2)\\ & =& 1+2(p+1)x+(p^2+4p+1)x^2+2p(p+1)x^3+p^2{x}^4\end{array}$$
この式の$x^2$の係数は$p^2+4p+1$これの$p$の係数が例1を表している。

例２

$A=B=C=p,D=E=F=1$とおく
$$\begin{array}{r}(1+px{)}^3(1+x{)}^3=(1+3px+3p^2{x}^2+p^3{x}^3)(1+3x+3x^2+x^3)\end{array}$$
$$\begin{array}{cccc}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^3{x}^6\end{array}$$
これの$x^2$の項は$3p^2+9p+3$この式は例２を表している。

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

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

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

$$\begin{array}{rcl}f(x)& =& a+bx+c{x}^2+d{x}^3+e{x}^4+..\\ f^{\prime }(x)& =& b+2cx+3d{x}^2+4e{x}^3+..\\ f^{″}(x)& =& 2c+6dx+12e{x}^2+..\\ f^{‴}(x)& =& 6d+24ex+..\end{array}$$

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

べき乗式の微分

$P=(1+px{)}^A$とすると
$$\begin{array}{rcl}P& =& (1+px{)}^A\\ P^{\prime }& =& Aa(1+px{)}^{A-1}\\ P^{″}& =& A(A-1)p^2(1+px{)}^{A-2}\\ P^{‴}& =& A(A-1)(A-2)p^3(1+px{)}^{A-3}\\ ..\\ P^{(s)}& =& \frac{A!p^s}{(A-s)!}(1+px{)}^{A-s}\end{array}$$
また
$$\begin{array}{r}{P}^{(s)}(0)=\frac{A!p^s}{(A-s)!}\end{array}$$

$Q=(1+x{)}^{N-A}$とすると
$$\begin{array}{rcl}Q& =& (1+x{)}^{N-A}\\ Q^{\prime }& =& (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{array}$$
また
$$\begin{array}{r}{Q}^{(t)}(0)=\frac{(N-A)!}{(N-A-t)!}\end{array}$$

積の微分

$$\begin{array}{rcl}(PQ{)}^{\prime }& =& P^{\prime }Q+P{Q}^{\prime }\\ (PQ{)}^{″}& =& P^{″}Q+2P^{\prime }{Q}^{\prime }+P{Q}^{″}\\ (PQ{)}^{‴}& =& P^{‴}Q+3P^{″}Q+3P{Q}^{″}+P{Q}^{‴}\\ ..\\ (PQ{)}^{(n)}& =& P^{(n)}Q+n{P}^{(n-1)}{Q}^{\prime }+{\displaystyle (\genfrac{}{}{0ex}{}{n}{2})}{P}^{(n-2)}{Q}^{″}+..+{\displaystyle (\genfrac{}{}{0ex}{}{n}{s})}{P}^{(s)}{Q}^{(n-s)}+..\end{array}$$
$x^n$の係数は
$$\begin{array}{r}(PQ{)}^{(n)}(0)/n!\end{array}$$
これまで結果より
$$\begin{array}{r}\frac{A!}{n!}(\frac{p^n}{(A-n)!}+n\frac{(N-A)p^{n-1}}{(A-n+1)!}+{\displaystyle (\genfrac{}{}{0ex}{}{n}{2})}\frac{(N-A)(N-A-1)p^{n-2}}{(A-n+2)!}+..)\end{array}$$
これの$p^a$の係数は
$$\begin{array}{rcl}& & {\displaystyle (\genfrac{}{}{0ex}{}{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)!}\\ & =& {\displaystyle (\genfrac{}{}{0ex}{}{A}{a})}{\displaystyle (\genfrac{}{}{0ex}{}{N-A}{n-a})}\end{array}$$

例1

$N=4,n=2,A=2$
$$\begin{array}{ll}a=0:& {\displaystyle (\genfrac{}{}{0ex}{}{2}{0})}{\displaystyle (\genfrac{}{}{0ex}{}{2}{2})}=1\\ a=1:& {\displaystyle (\genfrac{}{}{0ex}{}{2}{1})}{\displaystyle (\genfrac{}{}{0ex}{}{2}{1})}=4\\ a=2:& {\displaystyle (\genfrac{}{}{0ex}{}{2}{2})}{\displaystyle (\genfrac{}{}{0ex}{}{2}{0})}=1\end{array}$$

例2

$N=6,n=2,A=3$
$$\begin{array}{ll}a=0:& {\displaystyle (\genfrac{}{}{0ex}{}{3}{0})}{\displaystyle (\genfrac{}{}{0ex}{}{3}{2})}=3\\ a=1:& {\displaystyle (\genfrac{}{}{0ex}{}{3}{1})}{\displaystyle (\genfrac{}{}{0ex}{}{3}{1})}=9\\ a=2:& {\displaystyle (\genfrac{}{}{0ex}{}{3}{2})}{\displaystyle (\genfrac{}{}{0ex}{}{3}{0})}=3\end{array}$$

組み合わせ論的考察

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

一般化

総リード$N$にgeneAが$A$,geneBが$B$,geneCが$C$,...含まれている。これから$n$リードとってくる。そのリードにそれぞれのgeneが$a,b,c,..$含まれる場合の数を考える。母関数は
$$\begin{array}{r}(1+\alpha x{)}^A(1+\beta x{)}^B(1+\gamma x{)}^C...(1+\zeta x{)}^Z\end{array}$$
ここで$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{array}{rcl}& & (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{array}$$
$x^4$の係数は
$$\begin{array}{r}{\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{array}$$
なので$a,b,c$の場合の数は以下のようになる
$$\overline{)\begin{array}{ccc}a& b& c\\ 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{array}{rcl}P(x)& =& (1+\alpha x{)}^A(1+\beta x{)}^B..(1+\zeta x{)}^Z\\ & =& P_A{P}_B..P_Z\end{array}$$

積の微分

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

例

$A=2,B=4,C=2$,$a=1,b=2,c=1$のとき
$$\begin{array}{r}{\displaystyle (\genfrac{}{}{0ex}{}{2}{1})}{\displaystyle (\genfrac{}{}{0ex}{}{4}{2})}{\displaystyle (\genfrac{}{}{0ex}{}{2}{1})}=24\end{array}$$
正しく例を表している。