疑似1細胞発現データを作る
- 複数の細胞からのデータがある
- 細胞数が分かっている
- リード数を細胞数に分けて、疑似1細胞発現データを作りたい
- sample()関数を使って疑似作成できることを示し
# 検体数 N # ある遺伝子のNGSリード数が Mだった # m1 + ... + mN = M; mi >= 0となるような{mi}をランダムに生成したい N <- 100 M <- 50 # M本のリードを1:Nの検体にランダムに(等確率で)割り付ける s <- sample(1:N,M,replace=TRUE) # 検体IDごとに集計する t <- tabulate(s) # 関数にする my.divide.reads <- function(N,M){ s <- sample(1:N,M,replace=TRUE) t <- tabulate(s) ret <- rep(0,N) ret[1:length(t)] <- t return(ret) }
- そのうえで、多項分布でもできることを示す
# 別法 # 多項分布を使う t2 <- rmultinom(1,M,rep(1,N)) T2 <- rmultinom(10000,M,rep(1,N)) apply(T2,1,mean) S <- matrix(0,10000,N) for(i in 1:10000){ S[i,] <- my.divide.reads(N,M) }