疑似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)
}