RNAseqのためのDownsampling - ryamadaの遺伝学・遺伝統計学メモ

Downsamplingはそれなりに数のある標本から、標本数を減らしてリサンプリングする作業のこと
RNAseqのデータにおいて、複数の細胞からの１細胞RNAseqのデータを相互比較したりするときに、ダウンサンプリングした方がよい、という話が、こちらの論文にあったので、Rでシミュレーションしてどんな具合かを眺めてみる
１細胞RNAseq では、１細胞内にmRNAのプールがあって、Ng個(たとえば25000個とか)の遺伝子のそれぞれのmRNA分子が、n_1,n_2,...,n_{Ng}個というように整数値で存在していると考える。そして１細胞シークエンスという条件では、このn_1+n_2+...+n_{Ng}の総mRNA分子のプールから、それなりの数が多項分布的に選び出されてそれが「読み本数」として報告される、という仕組みになっている
こんなときに、実験を複数の細胞で実施すると、これまた実験条件ブレが結構入るのは少量系の宿命なのだが、「多項分布的な選び出し」の部分にかなりのばらつきが入ってしまう
また、この論文では、典型的な細胞状態、というものを複数想定して、その想定細胞から、RNAseq結果が得られたとして、どの想定細胞である確率がどれくらい高いか、という評価をしているのだが、それをするときに、細胞ごとに、読み本数がばらつくと、そのばらつきが、「比率評価」の曖昧度に影響することもあり、「読み本数」をそろえたい
「読み本数」をそろえるためにダウンサンプリングしましょう、となっている
リプレイスなしで指定本数までダウンサンプリングしている
書き散らしソース

# 遺伝子数
Ng <- 1000
# １細胞シークエンスする細胞数
Nc <- 1000
# 遺伝子別のmRNA割合の設定
p <- dexp(1:Ng,0.01)
p <- p/sum(p)

plot(p)
# RNAseqでの読み本数
X <- matrix(0,Nc,Ng)
Y <- list()
Z <- rep(0,Ng)
for(i in 1:Nc){
	# 細胞ごとに総読み本数をばらつかせる
	N.mRNA <- round(10^(runif(1)*2+3))
	#print(N.mRNA)
	s <- sample(1:Ng,N.mRNA,replace=TRUE,prob=p)
	X[i,] <- tabulate(s,Ng)
	Y[[i]] <- s
	Z <- Z + X[i,]
}

# ダウンサンプリングの本数
thres <- 10^4
P1 <- P2 <- matrix(0,Nc,Ng)
oks <- rep(0,Nc)
n.rna <- rep(0,Nc)
for(i in 1:Nc){
	P1[i,] <- (X[i,]+1)/sum(X[i,]+1)
	n.rna[i] <- sum(X[i,])
	# リプレースなしでダウンサンプリングの目標数に達しないときは「使わない」
	if(sum(X[i,]) > thres){
		tmp <- sample(Y[[i]],thres)
		tmp.tab <- tabulate(tmp,Ng)
		P2[i,] <- (tmp.tab+1)/sum(tmp.tab+1)
		oks[i] <- 1
	}
}

mean1 <- apply(P1,2,mean)
mean1.1 <- apply(P1[which(oks==1),],2,mean)
var1 <- apply(P1,2,var)
mean2 <- apply(P2[which(oks==1),],2,mean)
var2 <- apply(P2[which(oks==1),],2,var)



plot(mean1,p)
plot(mean2,p)

matplot(log(cbind(p,mean1,mean1.1,mean2,Z/sum(Z))),type="l")
legend(0, -10, c("theory","simple.mean","mean.of.manys","down.sampled","total.mean"), col = 1:5, lty = 1:5)

plot(n.rna,P1[,1])