单细胞交响乐25-实战八 Smart-seq2 胰腺细胞
刘小泽写于2020.7.20
为何取名叫“交响乐”？因为单细胞分析就像一个大乐团，需要各个流程的协同配合
单细胞交响乐1-常用的数据结构SingleCellExperiment
单细胞交响乐2-scRNAseq从实验到下游简介
单细胞交响乐3-细胞质控
单细胞交响乐4-归一化
单细胞交响乐5-挑选高变化基因
单细胞交响乐6-降维
单细胞交响乐7-聚类分群
单细胞交响乐8-marker基因检测
单细胞交响乐9-细胞类型注释
单细胞交响乐9-细胞类型注释
单细胞交响乐10-数据集整合后的批次矫正
单细胞交响乐11-多样本间差异分析
单细胞交响乐12-检测Doublet
单细胞交响乐13-细胞周期推断
单细胞交响乐14-细胞轨迹推断
单细胞交响乐15-scRNA与蛋白丰度信息结合
单细胞交响乐16-处理大型数据
单细胞交响乐17-不同单细胞R包的数据格式相互转换
单细胞交响乐18-实战一 Smart-seq2
单细胞交响乐19-实战二 STRT-Seq
单细胞交响乐20-实战三 10X 未过滤的PBMC数据
单细胞交响乐21-实战三 批量处理并整合多个10X PBMC数据
单细胞交响乐22-实战五 CEL-seq2
单细胞交响乐23-实战六 CEL-seq
单细胞交响乐24-实战七 SMARTer
1 前言
前面的种种都是作为知识储备，但是不实战还是记不住前面的知识
这是第八个实战练习
这是也是来自多个供体的人类胰腺细胞，使用Smart-seq2建库技术，数据来自Segerstolpe et al. (2016)
数据准备
library(scRNAseq)
sce.seger <- SegerstolpePancreasData()
sce.seger
# class: SingleCellExperiment
# dim: 26179 3514
# metadata(0):
# assays(1): counts
# rownames(26179): SGIP1 AZIN2 ... BIVM-ERCC5 eGFP
# rowData names(2): symbol refseq
# colnames(3514): HP1502401_N13 HP1502401_D14 ...
# HP1526901T2D_O11 HP1526901T2D_A8
# colData names(8): Source Name individual ... age
# body mass index
# reducedDimNames(0):
# altExpNames(1): ERCC
看到3500多个细胞，包含ERCC，使用Symbol ID
看下样本信息：
ID转换
选择的方式是：将没有匹配的NA去掉，并且去掉重复的行
# 首先得到symbol ID和对应的Ensembl ID（其中会存在无对应的NA情况）
library(AnnotationHub)
edb <- AnnotationHub()[["AH73881"]]
symbols <- rowData(sce.seger)$symbol
ens.id <- mapIds(edb, keys=symbols, keytype="SYMBOL", column="GENEID")
# 之前见到的方法是：
# keep <- !is.na(gene.ids) & !duplicated(gene.ids)
# 这里使用了另一种方法（不是直接将NA去掉，而且替换成了symbol）
ens.id <- ifelse(is.na(ens.id), symbols, ens.id)
keep <- !duplicated(ens.id)
sce.seger <- sce.seger[keep,]
rownames(sce.seger) <- ens.id[keep]
小结一下：至此见到了三种ID转换的方式，根据最后保留的基因数量，可以排个序：
保留基因最多（保留了NA和重复）：
uniquifyFeatureNames
中等（保留了NA，去掉重复）：ifelse(is.na(ens.id), symbols, ens.id)
最少（去掉了NA以及重复）：!is.na(gene.ids) & !duplicated(gene.ids)
编辑样本信息
之前有8列样本的信息，有点冗余了。这里只保留3列关心的，并重新命名
emtab.meta <- colData(sce.seger)[,c("cell type",
"individual", "single cell well quality")]
colnames(emtab.meta) <- c("CellType", "Donor", "Quality")
colData(sce.seger) <- emtab.meta
另外把细胞类型这一列中的“cell”字符去掉，并把首字母大写
sce.seger$CellType <- gsub(" cell", "", sce.seger$CellType)
sce.seger$CellType <- paste0(
toupper(substr(sce.seger$CellType, 1, 1)),
substring(sce.seger$CellType, 2))
2 质控
依然是备份一下，把unfiltered数据主要用在质控的探索上
unfiltered <- sce.seger
之前作者在数据中已经标注了细胞质量，可以看到有问题的细胞还是很多的：
table(sce.seger$Quality)
#
# control, 2-cell well control, empty well low quality cell OK
# 32 96 1177 2209
因此就要注意了，这里的数据会不会满足“大部分细胞都是高质量的”这个假设？
还是需要试一下，看看结果先
library(scater)
stats <- perCellQCMetrics(sce.seger)
qc1 <- quickPerCellQC(stats, percent_subsets="altexps_ERCC_percent",
batch=sce.seger$Donor)
colData(unfiltered) <- cbind(colData(unfiltered), stats)
unfiltered$discard <- qc1$discard
gridExtra::grid.arrange(
plotColData(unfiltered, x="Donor", y="sum", colour_by="discard") +
scale_y_log10() + ggtitle("Total count") +
theme(axis.text.x = element_text(angle = 90)),
plotColData(unfiltered, x="Donor", y="detected", colour_by="discard") +
scale_y_log10() + ggtitle("Detected features") +
theme(axis.text.x = element_text(angle = 90)),
plotColData(unfiltered, x="Donor", y="altexps_ERCC_percent",
colour_by="discard") + ggtitle("ERCC percent") +
theme(axis.text.x = element_text(angle = 90)),
ncol=3
)
看到HP1509101在过滤时存在过滤不完全的情况，HP1504901过滤的ERCC数量太多，推测这两个批次效果可能并不是很好，可能存在大量的低质量细胞
因此，再次指定subset 参数，重新画图
library(scater)
stats <- perCellQCMetrics(sce.seger)
qc <- quickPerCellQC(stats, percent_subsets="altexps_ERCC_percent",
batch=sce.seger$Donor,
subset=!sce.seger$Donor %in% c("HP1504901", "HP1509101"))
看看过滤掉多少
colSums(as.matrix(qc))
## low_lib_size low_n_features high_altexps_ERCC_percent
## 788 1056 1031
## discard
## 1246
最后将qc过滤的与本来标注低质量的一同过滤
low.qual <- sce.seger$Quality == "low quality cell"
sce.seger <- sce.seger[,!(qc$discard | low.qual)]
# 过滤了大概1500个细胞
> dim(unfiltered);dim(sce.seger)
[1] 26179 3514
[1] 26179 2090
3 归一化
此处会有一点小问题，值得注意！
本来有ERCC，操作应该是：
library(scran)
sce.seger = computeSpikeFactors(sce.seger, "ERCC")
sce.seger <- logNormCounts(sce.seger)
# Error in .local(x, ...) : size factors should be positive
但由于存在几个细胞中一个ERCC都没有，所以会报错
此时面临两个选择：要么把这几个细胞去掉；要么就不借助ERCC，用另一种去卷积方法
> table(colSums(counts(altExp(sce.seger)))==0)
FALSE TRUE
2087 3
如果要去掉这几个细胞：
test=sce.seger[,!colSums(counts(altExp(sce.seger)))==0]
sce.test = computeSpikeFactors(test, "ERCC")
sce.test <- logNormCounts(test)
我们这里选择保守的方法，不去掉细胞，使用另一种去卷积方法：
clusters <- quickCluster(sce.seger)
sce.seger <- computeSumFactors(sce.seger, clusters=clusters)
sce.seger <- logNormCounts(sce.seger)
summary(sizeFactors(sce.seger))
# Min. 1st Qu. Median Mean 3rd Qu. Max.
# 0.0000 0.1832 0.4016 1.0000 1.0996 12.9607
4 找高变异基因
下面构建模型想使用modelGeneVarWithSpikes，于是首先应该把那几个没有ERCC的细胞去掉；另外由于AZ这个批次相对其他批次的细胞数量过于少，因此在模型构建中也把它去掉吧
for.hvg <- sce.seger[,librarySizeFactors(altExp(sce.seger)) > 0
& sce.seger$Donor!="AZ"]
dec.seger <- modelGeneVarWithSpikes(for.hvg, "ERCC", block=for.hvg$Donor)
chosen.hvgs <- getTopHVGs(dec.seger, n=2000)
如果要批量作图检查的话
# 批次数量较多，因此设置多行多列显示
par(mfrow=c(3,3))
blocked.stats <- dec.seger$per.block
for (i in colnames(blocked.stats)) {
current <- blocked.stats[[i]]
plot(current$mean, current$total, main=i, pch=16, cex=0.5,
xlab="Mean of log-expression", ylab="Variance of log-expression")
curfit <- metadata(current)
points(curfit$mean, curfit$var, col="red", pch=16)
curve(curfit$trend(x), col='dodgerblue', add=TRUE, lwd=2)
}
注意，这里在找完HVGs后，没有进行批次矫正，如果继续向下做，会发现什么？
5 降维聚类
降维
library(BiocSingular)
set.seed(101011001)
sce.seger <- runPCA(sce.seger, subset_row=chosen.hvgs, ncomponents=25)
sce.seger <- runTSNE(sce.seger, dimred="PCA")
聚类
snn.gr <- buildSNNGraph(sce.seger, use.dimred="PCA")
colLabels(sce.seger) <- factor(igraph::cluster_walktrap(snn.gr)$membership)
检查聚类分群与批次
tab <- table(Cluster=colLabels(sce.seger), Donor=sce.seger$Donor)
library(pheatmap)
pheatmap(log10(tab+10), color=viridis::viridis(100))
结果真的是：批次效应影响了分群，因此最好还是做一遍fastMNN操作
tSNE图中也是显示出了强烈的批次效应
gridExtra::grid.arrange(
plotTSNE(sce.seger, colour_by="label"),
plotTSNE(sce.seger, colour_by="Donor"),
ncol=2
)
6 补充矫正批次效应
上图看到很明显的批次效应，那么如果处理后，会有什么不同吗？
利用fastMNN矫正
library(batchelor)
set.seed(1001010)
merged.seger <- fastMNN(sce.seger, subset.row=chosen.hvgs,
batch=sce.seger$Donor)
merged.seger
# class: SingleCellExperiment
# dim: 2000 2090
# metadata(2): merge.info pca.info
# assays(1): reconstructed
# rownames(2000): GCG TTR ... MAP6 LCP1
# rowData names(1): rotation
# colnames(2090): HP1502401_H13 HP1502401_J14 ...
# HP1526901T2D_N8 HP1526901T2D_A8
# colData names(2): batch label
# reducedDimNames(2): corrected TSNE
# altExpNames(0):
# metadata(merged.seger)$merge.info$lost.var
# lost.var ：值越大表示丢失的真实生物异质性越多
因为fastMNN会包含PCA降维，所以下面继续进行tSNE即可
降维聚类
library(BiocSingular)
set.seed(101011001)
merged.seger <- runTSNE(merged.seger, dimred="corrected")
snn.gr <- buildSNNGraph(merged.seger, use.dimred="corrected")
colLabels(merged.seger) <- factor(igraph::cluster_walktrap(snn.gr)$membership)
再次作图，是不是明显比之前好很多？
gridExtra::grid.arrange(
plotTSNE(merged.seger, colour_by="label"),
plotTSNE(merged.seger, colour_by="batch"),
ncol=2
)
欢迎关注我们的公众号～_～
我们是两个农转生信的小硕，打造生信星球，想让它成为一个不拽术语、通俗易懂的生信知识平台。需要帮助或提出意见请后台留言或发送邮件到jieandze1314@gmail.com
Welcome to our bioinfoplanet!
