今天在造dropout模拟数据集的时候,发现了一个问题,Rtsne的结果和skelarn中的tsne画图结果非常不同,让我感到很好奇,到底是什么原因导致这个的
%load_ext rpy2.ipython
%%R -o counts -o truecounts -o geneinfo -o cellinfo
# make sure that splatter is installed: https://github.com/Oshlack/splatter
suppressPackageStartupMessages({
library(splatter)
library(Rtsne)
library(ggplot2)
library(repr)
})
simulate <- function(nGroups=2, nGenes=200, batchCells=2000, dropout=3)
{
if (nGroups > 1) method <- 'groups'
else method <- 'single'
group.prob <- rep(1, nGroups) / nGroups
# new splatter requires dropout.type
if ('dropout.type' %in% slotNames(newSplatParams())) {
if (dropout)
dropout.type <- 'experiment'
else
dropout.type <- 'none'
sim <- splatSimulate(group.prob=group.prob, nGenes=nGenes, batchCells=batchCells,
dropout.type=dropout.type, method=method,
seed=42, dropout.shape=-1, dropout.mid=dropout)
} else {
sim <- splatSimulate(group.prob=group.prob, nGenes=nGenes, batchCells=batchCells,
dropout.present=!dropout, method=method,
seed=42, dropout.shape=-1, dropout.mid=dropout)
}
counts <- as.data.frame(t(counts(sim)))
truecounts <- as.data.frame(t(assays(sim)$TrueCounts))
dropout <- assays(sim)$Dropout
mode(dropout) <- 'integer'
cellinfo <- as.data.frame(colData(sim))
geneinfo <- as.data.frame(rowData(sim))
list(counts=counts,
cellinfo=cellinfo,
geneinfo=geneinfo,
truecounts=truecounts)
}
sim <- simulate(nGroups=6, dropout=1) ## 因为这个地方的dropout从1变成了3,导致两者不一致
counts <- sim$counts
geneinfo <- sim$geneinfo
cellinfo <- sim$cellinfo
truecounts <- sim$truecounts
X <- t(counts) ## counts with dropout
Y <- as.integer(substring(cellinfo$Group,6))
Y <- Y-1
X.normalized <- apply(X, 2, function(z) z/sum(z)) ## 我大概知道这个图的原因了,因为sc.ppnor
X.normalized <- log(X.normalized + 1)
tsne.X <- Rtsne(t(X.normalized), max_iter = 1000)
tsne_plot.X <- data.frame(`x-tsne` = tsne.X$Y[,1], `y-tsne` = tsne.X$Y[,2],
truelabel = Y, check.names = F)
tsne_plot.X$truelabel <- factor(tsne_plot.X$truelabel, levels = c(0:max(Y)))
ggplot(tsne_plot.X) + geom_point(aes(x=`x-tsne`, y=`y-tsne`, color=truelabel), size=0.5) +
labs(color='true label') +
ggtitle("Simulated data") +
theme_bw() +
theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
legend.key = element_rect(fill = 'white', colour = 'white'), legend.position="none",
axis.title.y=element_blank(), axis.title.x=element_blank())
# sim <- simulate()
# counts <- sim$counts
# geneinfo <- sim$geneinfo
# cellinfo <- sim$cellinfo
# truecounts <- sim$truecounts
## 我也可以在R中显示的
结果如下
visulize in python
import scanpy as sc
adata = sc.AnnData(counts.values, obs=cellinfo, var=geneinfo)
adata.obs_names = cellinfo.Cell
adata.var_names = geneinfo.Gene
adata.obs.index = list(adata.obs.index) ## add this line to avoid error
adata.var.index = list(adata.var.index) ## add this line to avoid error
#sc.pp.filter_genes(adata, min_counts=1)
adata_true = sc.AnnData(truecounts.values, obs=cellinfo, var=geneinfo)
adata_true.obs_names = cellinfo.Cell
adata_true.var_names = geneinfo.Gene
adata_true.obs.index = list(adata_true.obs.index) ## add this line to avoid error
adata_true.var.index = list(adata_true.var.index) ## add this line to avoid error
adata_true = adata_true[:, adata.var_names].copy()
sc.pp.normalize_per_cell(adata)
sc.pp.normalize_per_cell(adata_true)
sc.pp.log1p(adata)
sc.pp.log1p(adata_true)
sc.pp.pca(adata)
sc.pp.pca(adata_true)
sc.tl.tsne(adata)
sc.tl.tsne(adata_true)
sc.pp.neighbors(adata)
sc.pp.neighbors(adata_true)
sc.pl.tsne(adata,color=["Group"],title="dropout_tsne")
sc.pl.tsne(adata_true,color=["Group"],title="true_tsne")
sc.tl.umap(adata)
sc.tl.umap(adata_true)
sc.pl.umap(adata,color=["Group"],title="dropout_umap")
sc.pl.umap(adata_true,color=["Group"],title="true_out")
可以看到这结果的tsne看起来是不一样的,这个结果很奇怪,Rtsne有两个cluster竟然显示出来了,而python中的tsne直接混合了数据
我首先排除了可能是Rtsne的原因
# X_raw=adata.X.copy()
# Y=adata.obs["Group"].cat.codes.values
X_raw=adata_true.X.copy()
Y=adata_true.obs["Group"].cat.codes.values
%%R -i X_raw -i Y
tsne.X <- Rtsne(X_raw, max_iter = 1000)
tsne_plot.X <- data.frame(`x-tsne` = tsne.X$Y[,1], `y-tsne` = tsne.X$Y[,2],
truelabel = Y, check.names = F)
tsne_plot.X$truelabel <- factor(tsne_plot.X$truelabel, levels = c(0:max(Y)))
ggplot(tsne_plot.X) + geom_point(aes(x=`x-tsne`, y=`y-tsne`, color=truelabel), size=0.5) +
labs(color='true label') +
ggtitle("Simulated data") +
theme_bw() +
theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
legend.key = element_rect(fill = 'white', colour = 'white'), legend.position="none",
axis.title.y=element_blank(), axis.title.x=element_blank())
结果如下
所以唯一的问题,就是这个预处理的过程
R的预处理流程如下
而scanpy的处理流程如下
根据我前几天刚看的sc.pp.normal_per_cell的实现原理,里面涉及到中位数,但是R中的标准化并没有,也就是说scanpy得到的结果可能会更稳定,这也就是结果的原因
我验证了这个结果,我在R中实现了中位数的标准化,
%%R -o counts -o truecounts -o geneinfo -o cellinfo
# make sure that splatter is installed: https://github.com/Oshlack/splatter
suppressPackageStartupMessages({
library(splatter)
library(Rtsne)
library(ggplot2)
library(repr)
})
simulate <- function(nGroups=2, nGenes=200, batchCells=2000, dropout=3)
{
if (nGroups > 1) method <- 'groups'
else method <- 'single'
group.prob <- rep(1, nGroups) / nGroups
# new splatter requires dropout.type
if ('dropout.type' %in% slotNames(newSplatParams())) {
if (dropout)
dropout.type <- 'experiment'
else
dropout.type <- 'none'
sim <- splatSimulate(group.prob=group.prob, nGenes=nGenes, batchCells=batchCells,
dropout.type=dropout.type, method=method,
seed=42, dropout.shape=-1, dropout.mid=dropout)
} else {
sim <- splatSimulate(group.prob=group.prob, nGenes=nGenes, batchCells=batchCells,
dropout.present=!dropout, method=method,
seed=42, dropout.shape=-1, dropout.mid=dropout)
}
counts <- as.data.frame(t(counts(sim)))
truecounts <- as.data.frame(t(assays(sim)$TrueCounts))
dropout <- assays(sim)$Dropout
mode(dropout) <- 'integer'
cellinfo <- as.data.frame(colData(sim))
geneinfo <- as.data.frame(rowData(sim))
list(counts=counts,
cellinfo=cellinfo,
geneinfo=geneinfo,
truecounts=truecounts)
}
sim <- simulate(nGroups=6, dropout=1) ## 因为这个地方的dropout从1变成了3,导致两者不一致
counts <- sim$counts
geneinfo <- sim$geneinfo
cellinfo <- sim$cellinfo
truecounts <- sim$truecounts
X <- t(counts) ## counts with dropout
Y <- as.integer(substring(cellinfo$Group,6))
Y <- Y-1
#X.normalized <- apply(X, 2, function(z) z/sum(z))
X_t=t(X)
me=median(apply(X_t,1,sum))
sf=me/apply(X_t,1,sum)
X.normalized=diag(sf)%*%X_t
X.normalized <- log(X.normalized + 1)
tsne.X <- Rtsne(X.normalized, max_iter = 1000)
tsne_plot.X <- data.frame(`x-tsne` = tsne.X$Y[,1], `y-tsne` = tsne.X$Y[,2],
truelabel = Y, check.names = F)
tsne_plot.X$truelabel <- factor(tsne_plot.X$truelabel, levels = c(0:max(Y)))
ggplot(tsne_plot.X) + geom_point(aes(x=`x-tsne`, y=`y-tsne`, color=truelabel), size=0.5) +
labs(color='true label') +
ggtitle("Simulated data") +
theme_bw() +
theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
legend.key = element_rect(fill = 'white', colour = 'white'), legend.position="none",
axis.title.y=element_blank(), axis.title.x=element_blank())
# sim <- simulate()
# counts <- sim$counts
# geneinfo <- sim$geneinfo
# cellinfo <- sim$cellinfo
# truecounts <- sim$truecounts
## 我也可以在R中显示的
结果如下