4.8 实战八 | Smart-seq2 | 人胰腺细胞
刘小泽写于2020.7.20

1 前言

这是也是来自多个供体的人类胰腺细胞,使用Smart-seq2建库技术,数据来自Segerstolpe et al. (2016)

数据准备

1
library(scRNAseq)
2
sce.seger <- SegerstolpePancreasData()
3
sce.seger
4
# class: SingleCellExperiment
5
# dim: 26179 3514
6
# metadata(0):
7
# assays(1): counts
8
# rownames(26179): SGIP1 AZIN2 ... BIVM-ERCC5 eGFP
9
# rowData names(2): symbol refseq
10
# colnames(3514): HP1502401_N13 HP1502401_D14 ...
11
# HP1526901T2D_O11 HP1526901T2D_A8
12
# colData names(8): Source Name individual ... age
13
# body mass index
14
# reducedDimNames(0):
15
# altExpNames(1): ERCC
Copied!
看到3500多个细胞,包含ERCC,使用Symbol ID
看下样本信息:

ID转换

选择的方式是:将没有匹配的NA去掉,并且去掉重复的行
1
# 首先得到symbol ID和对应的Ensembl ID(其中会存在无对应的NA情况)
2
library(AnnotationHub)
3
edb <- AnnotationHub()[["AH73881"]]
4
symbols <- rowData(sce.seger)$symbol
5
ens.id <- mapIds(edb, keys=symbols, keytype="SYMBOL", column="GENEID")
6
7
# 之前见到的方法是:
8
# keep <- !is.na(gene.ids) & !duplicated(gene.ids)
9
10
# 这里使用了另一种方法(不是直接将NA去掉,而且替换成了symbol)
11
ens.id <- ifelse(is.na(ens.id), symbols, ens.id)
12
keep <- !duplicated(ens.id)
13
14
sce.seger <- sce.seger[keep,]
15
rownames(sce.seger) <- ens.id[keep]
Copied!
小结一下:至此见到了三种ID转换的方式,根据最后保留的基因数量,可以排个序:
保留基因最多(保留了NA和重复):uniquifyFeatureNames 中等(保留了NA,去掉重复):ifelse(is.na(ens.id), symbols, ens.id) 最少(去掉了NA以及重复):!is.na(gene.ids) & !duplicated(gene.ids)

编辑样本信息

之前有8列样本的信息,有点冗余了。这里只保留3列关心的,并重新命名
1
emtab.meta <- colData(sce.seger)[,c("cell type",
2
"individual", "single cell well quality")]
3
colnames(emtab.meta) <- c("CellType", "Donor", "Quality")
4
colData(sce.seger) <- emtab.meta
Copied!
另外把细胞类型这一列中的“cell”字符去掉,并把首字母大写
1
sce.seger$CellType <- gsub(" cell", "", sce.seger$CellType)
2
sce.seger$CellType <- paste0(
3
toupper(substr(sce.seger$CellType, 1, 1)),
4
substring(sce.seger$CellType, 2))
Copied!

2 质控

依然是备份一下,把unfiltered数据主要用在质控的探索上
1
unfiltered <- sce.seger
Copied!
之前作者在数据中已经标注了细胞质量,可以看到有问题的细胞还是很多的:
1
table(sce.seger$Quality)
2
#
3
# control, 2-cell well control, empty well low quality cell OK
4
# 32 96 1177 2209
Copied!
因此就要注意了,这里的数据会不会满足“大部分细胞都是高质量的”这个假设?
还是需要试一下,看看结果先
1
library(scater)
2
stats <- perCellQCMetrics(sce.seger)
3
qc1 <- quickPerCellQC(stats, percent_subsets="altexps_ERCC_percent",
4
batch=sce.seger$Donor)
5
6
7
colData(unfiltered) <- cbind(colData(unfiltered), stats)
8
unfiltered$discard <- qc1$discard
9
10
gridExtra::grid.arrange(
11
plotColData(unfiltered, x="Donor", y="sum", colour_by="discard") +
12
scale_y_log10() + ggtitle("Total count") +
13
theme(axis.text.x = element_text(angle = 90)),
14
plotColData(unfiltered, x="Donor", y="detected", colour_by="discard") +
15
scale_y_log10() + ggtitle("Detected features") +
16
theme(axis.text.x = element_text(angle = 90)),
17
plotColData(unfiltered, x="Donor", y="altexps_ERCC_percent",
18
colour_by="discard") + ggtitle("ERCC percent") +
19
theme(axis.text.x = element_text(angle = 90)),
20
ncol=3
21
)
Copied!
看到HP1509101在过滤时存在过滤不完全的情况,HP1504901过滤的ERCC数量太多,推测这两个批次效果可能并不是很好,可能存在大量的低质量细胞
因此,再次指定subset 参数,重新画图
1
library(scater)
2
stats <- perCellQCMetrics(sce.seger)
3
qc <- quickPerCellQC(stats, percent_subsets="altexps_ERCC_percent",
4
batch=sce.seger$Donor,
5
subset=!sce.seger$Donor %in% c("HP1504901", "HP1509101"))
Copied!
看看过滤掉多少
1
colSums(as.matrix(qc))
2
## low_lib_size low_n_features high_altexps_ERCC_percent
3
## 788 1056 1031
4
## discard
5
## 1246
Copied!
最后将qc过滤的与本来标注低质量的一同过滤
1
low.qual <- sce.seger$Quality == "low quality cell"
2
sce.seger <- sce.seger[,!(qc$discard | low.qual)]
3
4
# 过滤了大概1500个细胞
5
> dim(unfiltered);dim(sce.seger)
6
[1] 26179 3514
7
[1] 26179 2090
Copied!

3 归一化

此处会有一点小问题,值得注意!
本来有ERCC,操作应该是:
1
library(scran)
2
sce.seger = computeSpikeFactors(sce.seger, "ERCC")
3
sce.seger <- logNormCounts(sce.seger)
4
# Error in .local(x, ...) : size factors should be positive
Copied!
但由于存在几个细胞中一个ERCC都没有,所以会报错
此时面临两个选择:要么把这几个细胞去掉;要么就不借助ERCC,用另一种去卷积方法
1
> table(colSums(counts(altExp(sce.seger)))==0)
2
3
FALSE TRUE
4
2087 3
Copied!
如果要去掉这几个细胞:
1
test=sce.seger[,!colSums(counts(altExp(sce.seger)))==0]
2
sce.test = computeSpikeFactors(test, "ERCC")
3
sce.test <- logNormCounts(test)
Copied!
我们这里选择保守的方法,不去掉细胞,使用另一种去卷积方法:
1
clusters <- quickCluster(sce.seger)
2
sce.seger <- computeSumFactors(sce.seger, clusters=clusters)
3
sce.seger <- logNormCounts(sce.seger)
4
5
summary(sizeFactors(sce.seger))
6
# Min. 1st Qu. Median Mean 3rd Qu. Max.
7
# 0.0000 0.1832 0.4016 1.0000 1.0996 12.9607
Copied!

4 找表达量高变化基因

下面构建模型想使用modelGeneVarWithSpikes,于是首先应该把那几个没有ERCC的细胞去掉;另外由于AZ这个批次相对其他批次的细胞数量过于少,因此在模型构建中也把它去掉吧
1
for.hvg <- sce.seger[,librarySizeFactors(altExp(sce.seger)) > 0
2
& sce.seger$Donor!="AZ"]
3
dec.seger <- modelGeneVarWithSpikes(for.hvg, "ERCC", block=for.hvg$Donor)
4
chosen.hvgs <- getTopHVGs(dec.seger, n=2000)
Copied!
如果要批量作图检查的话
1
# 批次数量较多,因此设置多行多列显示
2
par(mfrow=c(3,3))
3
blocked.stats <- dec.seger$per.block
4
for (i in colnames(blocked.stats)) {
5
current <- blocked.stats[[i]]
6
plot(current$mean, current$total, main=i, pch=16, cex=0.5,
7
xlab="Mean of log-expression", ylab="Variance of log-expression")
8
curfit <- metadata(current)
9
points(curfit$mean, curfit$var, col="red", pch=16)
10
curve(curfit$trend(x), col='dodgerblue', add=TRUE, lwd=2)
11
}
Copied!
注意,这里在找完HVGs后,没有进行批次矫正,如果继续向下做,会发现什么?

5 降维聚类

降维

1
library(BiocSingular)
2
set.seed(101011001)
3
sce.seger <- runPCA(sce.seger, subset_row=chosen.hvgs, ncomponents=25)
4
sce.seger <- runTSNE(sce.seger, dimred="PCA")
Copied!

聚类

1
snn.gr <- buildSNNGraph(sce.seger, use.dimred="PCA")
2
colLabels(sce.seger) <- factor(igraph::cluster_walktrap(snn.gr)$membership)
Copied!

检查聚类分群与批次

1
tab <- table(Cluster=colLabels(sce.seger), Donor=sce.seger$Donor)
2
library(pheatmap)
3
pheatmap(log10(tab+10), color=viridis::viridis(100))
Copied!
结果真的是:批次效应影响了分群,因此最好还是做一遍fastMNN操作
tSNE图中也是显示出了强烈的批次效应
1
gridExtra::grid.arrange(
2
plotTSNE(sce.seger, colour_by="label"),
3
plotTSNE(sce.seger, colour_by="Donor"),
4
ncol=2
5
)
Copied!

6 补充矫正批次效应

上图看到很明显的批次效应,那么如果处理后,会有什么不同吗?

利用fastMNN矫正

1
library(batchelor)
2
set.seed(1001010)
3
merged.seger <- fastMNN(sce.seger, subset.row=chosen.hvgs,
4
batch=sce.seger$Donor)
5
merged.seger
6
# class: SingleCellExperiment
7
# dim: 2000 2090
8
# metadata(2): merge.info pca.info
9
# assays(1): reconstructed
10
# rownames(2000): GCG TTR ... MAP6 LCP1
11
# rowData names(1): rotation
12
# colnames(2090): HP1502401_H13 HP1502401_J14 ...
13
# HP1526901T2D_N8 HP1526901T2D_A8
14
# colData names(2): batch label
15
# reducedDimNames(2): corrected TSNE
16
# altExpNames(0):
17
18
# metadata(merged.seger)$merge.info$lost.var
19
# lost.var :值越大表示丢失的真实生物异质性越多
Copied!
因为fastMNN会包含PCA降维,所以下面继续进行tSNE即可

降维聚类

1
library(BiocSingular)
2
set.seed(101011001)
3
merged.seger <- runTSNE(merged.seger, dimred="corrected")
4
5
snn.gr <- buildSNNGraph(merged.seger, use.dimred="corrected")
6
colLabels(merged.seger) <- factor(igraph::cluster_walktrap(snn.gr)$membership)
Copied!
再次作图,是不是明显比之前好很多?
1
gridExtra::grid.arrange(
2
plotTSNE(merged.seger, colour_by="label"),
3
plotTSNE(merged.seger, colour_by="batch"),
4
ncol=2
5
)
Copied!
最近更新 1yr ago