4.11 实战十一 | Smart-seq2 | 小鼠造血干细胞
刘小泽写于2020.7.21
1 前言
数据来自 (Nestorowa et al. 2016) 的小鼠造血干细胞 haematopoietic stem cell (HSC) ，使用的技术是Smart-seq2
准备数据
library(scRNAseq)
sce.nest <- NestorowaHSCData()
sce.nest
# class: SingleCellExperiment
# dim: 46078 1920
# metadata(0):
# assays(1): counts
# rownames(46078): ENSMUSG00000000001
# ENSMUSG00000000003 ... ENSMUSG00000107391
# ENSMUSG00000107392
# rowData names(0):
# colnames(1920): HSPC_007 HSPC_013 ... Prog_852
# Prog_810
# colData names(2): cell.type FACS
# reducedDimNames(1): diffusion
# altExpNames(1): ERCC
counts(sce.nest)[1:3,1:3]
# 3 x 3 sparse Matrix of class "dgCMatrix"
# HSPC_007 HSPC_013 HSPC_019
# ENSMUSG00000000001 . 7 1
# ENSMUSG00000000003 . . .
# ENSMUSG00000000028 4 1 2看到使用了ERCC、Ensembl ID
ID转换
library(AnnotationHub)
ens.mm.v97 <- AnnotationHub()[["AH73905"]]
anno <- select(ens.mm.v97, keys=rownames(sce.nest),
keytype="GENEID", columns=c("SYMBOL", "SEQNAME"))
# 这里全部对应
> sum(is.na(anno$SYMBOL))
[1] 0
> sum(is.na(anno$SEQNAME))
[1] 0
# 接下来只需要匹配顺序即可
rowData(sce.nest) <- anno[match(rownames(sce.nest), anno$GENEID),]
sce.nest
# class: SingleCellExperiment
# dim: 46078 1920
# metadata(0):
# assays(1): counts
# rownames(46078): ENSMUSG00000000001
# ENSMUSG00000000003 ... ENSMUSG00000107391
# ENSMUSG00000107392
# rowData names(3): GENEID SYMBOL SEQNAME
# colnames(1920): HSPC_007 HSPC_013 ... Prog_852
# Prog_810
# colData names(2): cell.type FACS
# reducedDimNames(1): diffusion
# altExpNames(1): ERCC2 质控
依然是备份一下，把unfiltered数据主要用在质控的探索上
unfiltered <- sce.nest这里没有线粒体基因，因此只能用ERCC计算过滤条件
library(scater)
stats <- perCellQCMetrics(sce.nest)
qc <- quickPerCellQC(stats, percent_subsets="altexps_ERCC_percent")
sce.nest <- sce.nest[,!qc$discard]
# 看下过滤的细胞
colSums(as.matrix(qc))
# low_lib_size low_n_features
# 146 28
# high_altexps_ERCC_percent discard
# 241 264做个图
colData(unfiltered) <- cbind(colData(unfiltered), stats)
unfiltered$discard <- qc$discard
gridExtra::grid.arrange(
plotColData(unfiltered, y="sum", colour_by="discard") +
scale_y_log10() + ggtitle("Total count"),
plotColData(unfiltered, y="detected", colour_by="discard") +
scale_y_log10() + ggtitle("Detected features"),
plotColData(unfiltered, y="altexps_ERCC_percent",
colour_by="discard") + ggtitle("ERCC percent"),
ncol=2
)最后对数据进行过滤
sce.nest <- sce.nest[,!qc$discard]
# 过滤前后
> dim(unfiltered);dim(sce.nest)
[1] 46078 1920
[1] 46078 16563 归一化
library(scran)
set.seed(101000110)
clusters <- quickCluster(sce.nest)
sce.nest <- computeSumFactors(sce.nest, clusters=clusters)
sce.nest <- logNormCounts(sce.nest)4 找表达量高变化基因
使用基于ERCC的构建模型方法
set.seed(00010101)
dec.nest <- modelGeneVarWithSpikes(sce.nest, "ERCC")
top.nest <- getTopHVGs(dec.nest, prop=0.1)
class(dec.nest)
# [1] "DFrame"
# attr(,"package")
# [1] "S4Vectors"
# 其中ERCC的信息就存储在dec.nest的metadata中
curfit <- metadata(dec.nest)
class(curfit)
# [1] "list"
names(curfit)
# [1] "mean" "var" "trend" "std.dev"
length(unique(names(curfit$mean))) # 一共92个ERCC spike-in
# [1] 92
# 其中的mean、var就定义了横纵坐标
head(curfit$mean)
# ERCC-00002 ERCC-00003 ERCC-00004 ERCC-00009 ERCC-00012 ERCC-00013
# 14.91183375 11.27060119 13.31197197 11.94866319 0.02211546 0.21249156
head(curfit$var)
# ERCC-00002 ERCC-00003 ERCC-00004 ERCC-00009 ERCC-00012 ERCC-00013
# 0.02375131 0.29308411 0.05376959 0.41814635 0.14928826 1.08599155然后把基因（黑点）、ERCC（红点）、根据ERCC拟合的线（蓝线）画出来
plot(dec.nest$mean, dec.nest$total, pch=16, cex=0.5,
xlab="Mean of log-expression", ylab="Variance of log-expression")
curve(curfit$trend(x), col='dodgerblue', add=TRUE, lwd=2)
points(curfit$mean, curfit$var, col="red")5 降维聚类
降维
set.seed(101010011)
sce.nest <- denoisePCA(sce.nest, technical=dec.nest, subset.row=top.nest)
sce.nest <- runTSNE(sce.nest, dimred="PCA")
# 检查PC的数量
ncol(reducedDim(sce.nest, "PCA"))
## [1] 9聚类
snn.gr <- buildSNNGraph(sce.nest, use.dimred="PCA")
colLabels(sce.nest) <- factor(igraph::cluster_walktrap(snn.gr)$membership)
table(colLabels(sce.nest))
##
## 1 2 3 4 5 6 7 8 9
## 203 472 258 175 142 229 20 83 74作图
plotTSNE(sce.nest, colour_by="label")6 marker基因检测
markers <- findMarkers(sce.nest, colLabels(sce.nest),
test.type="wilcox", direction="up", lfc=0.5,
row.data=rowData(sce.nest)[,"SYMBOL",drop=FALSE])比如检测一下cluster8：
chosen <- markers[['8']]
best <- chosen[chosen$Top <= 10,]
length(best)
# [1] 13
# 将cluster8与其他clusters对比的AUC结果提取出来
aucs <- getMarkerEffects(best, prefix="AUC")
rownames(aucs) <- best$SYMBOL
library(pheatmap)
pheatmap(aucs, color=viridis::plasma(100))看到其中血红蛋白相关基因（Hba1、Hba2、Hbb）、Car2、Hebp1基因上调，说明clsuter8可能包含红细胞前体细胞
7 细胞类型注释
将会使用内置的参考注释数据，
SingleR中就包含了一些内置数据集，大部分是bulk RNA-Seq或芯片数据中经过筛选的细胞类型。
准备参考数据
library(SingleR)
mm.ref <- MouseRNAseqData()
mm.ref
# class: SummarizedExperiment
# dim: 21214 358
# metadata(0):
# assays(1): logcounts
# rownames(21214): Xkr4 Rp1 ... LOC100039574
# LOC100039753
# rowData names(0):
# colnames(358): ERR525589Aligned
# ERR525592Aligned ... SRR1044043Aligned
# SRR1044044Aligned
# colData names(3): label.main label.fine
# label.ont进行转换
renamed <- sce.nest
# 参考数据集中使用的是symbol name，这里也转换一下
rownames(renamed) <- uniquifyFeatureNames(rownames(renamed),
rowData(sce.nest)$SYMBOL)
# 然后把我们的细胞在参考数据集中找对应的细胞类型
# 返回的pred结果是一个数据框，每行是我们自己数据的一个细胞
pred <- SingleR(test=renamed, ref=mm.ref, labels=mm.ref$label.fine)
table(pred$labels)
#
# B cells Endothelial cells Erythrocytes Granulocytes Macrophages Monocytes NK cells T cells
# 61 1 1005 1 2 500 1 85这里也看到cluster8与红细胞更相近
最后更新于
