【R code】SeuratObjectに格納したsingle-cell 10x dataを用いたbulk qPCRとの比較のためのコード

10x dataを解析しているときに、ふと、10x dataの平均とってうまく解析したら、bulkでcDNA作成し、qPCRした時の結果と同じ様になるのではないか？と考えました。

また、それはある種のquality checkにもなるのではないか、と考えたので、codeを紹介します。

codeの紹介

# your.data; nFeature, nCountなどで上限、下限を設定しsubset.また、必要ならmultiplet除去後のseuratobject
# "geneA","geneB", "geneC", "geneD"の発言量をみたい


> your.data <- NormalizeData(your.data, normalization.method = "LogNormalize", scale.factor=10000)

> avg.data <- as.data.frame(log(AverageExpression(your.data, verbose = FALSE)$RNA+1),base=2) 

> avg.data$gene <- rownames(avg.data)

> genes.to.label = c("geneA","geneB", "geneC", "geneD")

> 10x.gene <- subset(avg.data, gene == genes.to.label[1])
> for (i in genes.to.label){
  10x.add <- subset(avg.data, gene == i)
  10x.gene <- rbind(10x.gene, 10x.add)
}

> 10x.gene$gene <- as.factor(10x.gene$gene)
> 10x.gene <- transform(10x.gene, gene= factor(gene, levels=c("geneA","geneB", "geneC", "geneD")))

> ggplot(10x.gene, aes(x= gene, y= all)) +geom_col(width=0.5)

何をしていうのか解説

まずは、normalizeを行います。single cellだと各細胞で読まれるUMI数が異なり、UMIの多少で著明な遺伝子の発現量に違いが生じてしまいます。ですのでnormalizationとして、全ての細胞においてUMIを10Kと仮定します。そして、その値にlog(1+p)とします。こうすることで遺伝子発現の分布が正規分布となり、解析手法が使える様になる（のだと理解しています。）https://medium.com/@kyawsawhtoon/log-transformation-purpose-and-interpretation-9444b4b049c9

> your.data <- NormalizeData(your.data, normalization.method = “LogNormalize”, scale.factor=10000)

次に、your.data内にある全ての細胞に対してのRNA発現量を検討します。AverageExpressionはまず指数変換してから平均値を取る関数だそうで、それに対してlog2(x+1)を行い、各遺伝子のlog2ベースの発現量をavg.dataに格納

avg.data <- as.data.frame(log(AverageExpression(your.data, verbose = FALSE)$RNA+1),base=2)

今、all列に遺伝子発現量が格納され、列名に遺伝子名が入っているので、新しくgene列を作成し、格納

avg.data$gene <- rownames(avg.data)

10x.geneにしりたい遺伝子のみを抽出する。このとき、”geneA”が二回入ることになるが、outputは特に問題にならないので気にしない。

> 10x.gene <- subset(avg.data, gene == genes.to.label[1])

> for (i in genes.to.label){ 10x.add <- subset(avg.data, gene == i) 10x.gene <- rbind(10x.gene, 10x.add) }

plotの時にgeneの順番がバラバラになったので、あらかじめ順番を指定しておく。

> 10x.gene$gene <- as.factor(10x.gene$gene)

> 10x.gene <- transform(10x.gene, gene= factor(gene, levels=c(“geneA”,”geneB”, “geneC”, “geneD”)))

棒グラフで出力

ggplot(10x.gene, aes(x= gene, y= all)) +geom_col(width=0.5)

まとめ

これでqPCRのサイクルのように、log2をベースとした遺伝子発現量の違いを抽出することができるはずです。

解釈等間違っていましたら、教えてください。

通常のNormalizeDataのlogの底がeなのは統計解析上の都合なのだと解釈していますが、あってるのでしょうか、、