4.11 实战十一 | Smart-seq2 | 小鼠造血干细胞
刘小泽写于2020.7.21

1 前言

数据来自 (Nestorowa et al. 2016) 的小鼠造血干细胞 haematopoietic stem cell (HSC) ,使用的技术是Smart-seq2

准备数据

1
library(scRNAseq)
2
sce.nest <- NestorowaHSCData()
3
4
sce.nest
5
# class: SingleCellExperiment
6
# dim: 46078 1920
7
# metadata(0):
8
# assays(1): counts
9
# rownames(46078): ENSMUSG00000000001
10
# ENSMUSG00000000003 ... ENSMUSG00000107391
11
# ENSMUSG00000107392
12
# rowData names(0):
13
# colnames(1920): HSPC_007 HSPC_013 ... Prog_852
14
# Prog_810
15
# colData names(2): cell.type FACS
16
# reducedDimNames(1): diffusion
17
# altExpNames(1): ERCC
18
19
counts(sce.nest)[1:3,1:3]
20
# 3 x 3 sparse Matrix of class "dgCMatrix"
21
# HSPC_007 HSPC_013 HSPC_019
22
# ENSMUSG00000000001 . 7 1
23
# ENSMUSG00000000003 . . .
24
# ENSMUSG00000000028 4 1 2
Copied!
看到使用了ERCC、Ensembl ID

ID转换

1
library(AnnotationHub)
2
ens.mm.v97 <- AnnotationHub()[["AH73905"]]
3
anno <- select(ens.mm.v97, keys=rownames(sce.nest),
4
keytype="GENEID", columns=c("SYMBOL", "SEQNAME"))
5
# 这里全部对应
6
> sum(is.na(anno$SYMBOL))
7
[1] 0
8
> sum(is.na(anno$SEQNAME))
9
[1] 0
10
# 接下来只需要匹配顺序即可
11
rowData(sce.nest) <- anno[match(rownames(sce.nest), anno$GENEID),]
12
13
sce.nest
14
# class: SingleCellExperiment
15
# dim: 46078 1920
16
# metadata(0):
17
# assays(1): counts
18
# rownames(46078): ENSMUSG00000000001
19
# ENSMUSG00000000003 ... ENSMUSG00000107391
20
# ENSMUSG00000107392
21
# rowData names(3): GENEID SYMBOL SEQNAME
22
# colnames(1920): HSPC_007 HSPC_013 ... Prog_852
23
# Prog_810
24
# colData names(2): cell.type FACS
25
# reducedDimNames(1): diffusion
26
# altExpNames(1): ERCC
Copied!

2 质控

依然是备份一下,把unfiltered数据主要用在质控的探索上
1
unfiltered <- sce.nest
Copied!
这里没有线粒体基因,因此只能用ERCC计算过滤条件
1
library(scater)
2
stats <- perCellQCMetrics(sce.nest)
3
qc <- quickPerCellQC(stats, percent_subsets="altexps_ERCC_percent")
4
sce.nest <- sce.nest[,!qc$discard]
5
6
# 看下过滤的细胞
7
colSums(as.matrix(qc))
8
# low_lib_size low_n_features
9
# 146 28
10
# high_altexps_ERCC_percent discard
11
# 241 264
Copied!
做个图
1
colData(unfiltered) <- cbind(colData(unfiltered), stats)
2
unfiltered$discard <- qc$discard
3
4
gridExtra::grid.arrange(
5
plotColData(unfiltered, y="sum", colour_by="discard") +
6
scale_y_log10() + ggtitle("Total count"),
7
plotColData(unfiltered, y="detected", colour_by="discard") +
8
scale_y_log10() + ggtitle("Detected features"),
9
plotColData(unfiltered, y="altexps_ERCC_percent",
10
colour_by="discard") + ggtitle("ERCC percent"),
11
ncol=2
12
)
Copied!
最后对数据进行过滤
1
sce.nest <- sce.nest[,!qc$discard]
2
3
# 过滤前后
4
> dim(unfiltered);dim(sce.nest)
5
[1] 46078 1920
6
[1] 46078 1656
Copied!

3 归一化

1
library(scran)
2
set.seed(101000110)
3
clusters <- quickCluster(sce.nest)
4
sce.nest <- computeSumFactors(sce.nest, clusters=clusters)
5
sce.nest <- logNormCounts(sce.nest)
Copied!

4 找表达量高变化基因

使用基于ERCC的构建模型方法
1
set.seed(00010101)
2
dec.nest <- modelGeneVarWithSpikes(sce.nest, "ERCC")
3
top.nest <- getTopHVGs(dec.nest, prop=0.1)
4
5
class(dec.nest)
6
# [1] "DFrame"
7
# attr(,"package")
8
# [1] "S4Vectors"
9
10
# 其中ERCC的信息就存储在dec.nest的metadata中
11
curfit <- metadata(dec.nest)
12
class(curfit)
13
# [1] "list"
14
names(curfit)
15
# [1] "mean" "var" "trend" "std.dev"
16
length(unique(names(curfit$mean))) # 一共92个ERCC spike-in
17
# [1] 92
18
19
# 其中的mean、var就定义了横纵坐标
20
head(curfit$mean)
21
# ERCC-00002 ERCC-00003 ERCC-00004 ERCC-00009 ERCC-00012 ERCC-00013
22
# 14.91183375 11.27060119 13.31197197 11.94866319 0.02211546 0.21249156
23
head(curfit$var)
24
# ERCC-00002 ERCC-00003 ERCC-00004 ERCC-00009 ERCC-00012 ERCC-00013
25
# 0.02375131 0.29308411 0.05376959 0.41814635 0.14928826 1.08599155
Copied!
然后把基因(黑点)、ERCC(红点)、根据ERCC拟合的线(蓝线)画出来
1
plot(dec.nest$mean, dec.nest$total, pch=16, cex=0.5,
2
xlab="Mean of log-expression", ylab="Variance of log-expression")
3
curve(curfit$trend(x), col='dodgerblue', add=TRUE, lwd=2)
4
points(curfit$mean, curfit$var, col="red")
Copied!

5 降维聚类

降维

1
set.seed(101010011)
2
sce.nest <- denoisePCA(sce.nest, technical=dec.nest, subset.row=top.nest)
3
sce.nest <- runTSNE(sce.nest, dimred="PCA")
4
5
# 检查PC的数量
6
ncol(reducedDim(sce.nest, "PCA"))
7
## [1] 9
Copied!

聚类

1
snn.gr <- buildSNNGraph(sce.nest, use.dimred="PCA")
2
colLabels(sce.nest) <- factor(igraph::cluster_walktrap(snn.gr)$membership)
3
4
table(colLabels(sce.nest))
5
##
6
## 1 2 3 4 5 6 7 8 9
7
## 203 472 258 175 142 229 20 83 74
Copied!

作图

1
plotTSNE(sce.nest, colour_by="label")
Copied!

6 marker基因检测

1
markers <- findMarkers(sce.nest, colLabels(sce.nest),
2
test.type="wilcox", direction="up", lfc=0.5,
3
row.data=rowData(sce.nest)[,"SYMBOL",drop=FALSE])
Copied!
比如检测一下cluster8:
1
chosen <- markers[['8']]
2
best <- chosen[chosen$Top <= 10,]
3
length(best)
4
# [1] 13
5
6
# 将cluster8与其他clusters对比的AUC结果提取出来
7
aucs <- getMarkerEffects(best, prefix="AUC")
8
rownames(aucs) <- best$SYMBOL
9
10
library(pheatmap)
11
pheatmap(aucs, color=viridis::plasma(100))
Copied!
看到其中血红蛋白相关基因(Hba1、Hba2、Hbb)、Car2、Hebp1基因上调,说明clsuter8可能包含红细胞前体细胞

7 细胞类型注释

将会使用内置的参考注释数据,SingleR中就包含了一些内置数据集,大部分是bulk RNA-Seq或芯片数据中经过筛选的细胞类型。

准备参考数据

1
library(SingleR)
2
mm.ref <- MouseRNAseqData()
3
mm.ref
4
# class: SummarizedExperiment
5
# dim: 21214 358
6
# metadata(0):
7
# assays(1): logcounts
8
# rownames(21214): Xkr4 Rp1 ... LOC100039574
9
# LOC100039753
10
# rowData names(0):
11
# colnames(358): ERR525589Aligned
12
# ERR525592Aligned ... SRR1044043Aligned
13
# SRR1044044Aligned
14
# colData names(3): label.main label.fine
15
# label.ont
Copied!

进行转换

1
renamed <- sce.nest
2
# 参考数据集中使用的是symbol name,这里也转换一下
3
rownames(renamed) <- uniquifyFeatureNames(rownames(renamed),
4
rowData(sce.nest)$SYMBOL)
5
# 然后把我们的细胞在参考数据集中找对应的细胞类型
6
# 返回的pred结果是一个数据框,每行是我们自己数据的一个细胞
7
pred <- SingleR(test=renamed, ref=mm.ref, labels=mm.ref$label.fine)
8
table(pred$labels)
9
#
10
# B cells Endothelial cells Erythrocytes Granulocytes Macrophages Monocytes NK cells T cells
11
# 61 1 1005 1 2 500 1 85
Copied!
这里也看到cluster8与红细胞更相近
最近更新 1yr ago