1 前言
使用的是来自永生化小鼠骨髓祖细胞系的2个96孔板的416B cells（Lun et al. 2017 ），并且在细胞裂解后文库制备前，在每个细胞中加入一定量的外源RNA（ERCC）
ERCC就是外源RNA对照联盟开发的人工设计好的已知序列和数量的mRNA，高ERCC占比与低质量数据相关，并且通常是排除的标准 。假如看到ERCC在很多样本中占比超过了20%（假如认为20%是一个较高的比例），那就是说这些样本中有超过20%的reads都”浪费“在了这种不相关的外源序列上，而减少了自身内源转录本的量，这对下游分析是不利的
之后再进行高通量测序，得到每个基因的表达量（这个是通过计算比对到外显子区域的reads数得到的，就像featureCounts的操作）；同样，spike-in的含量也是计算有多少reads比对到了spike-in的参考序列
先简单了解一下scRNAseq数据集
早期版本中只是内置了三个数据集Fluidig、Th2、Allen：当时也写了这一篇单细胞转录组学习笔记-14-学习scRNAseq这个R包 ，但后来这个包进行了更新，原来的数据获取方法也发生了变化
复制 # 原来使用
data(fluidigm)
# 现在使用
ReprocessedFluidigmData() #provides 65 cells
ReprocessedTh2Data() #provides 96 T helper cells
ReprocessedAllenData() #provides 379 cells from the mouse visual cortex 另外新版本的包还加入了其他很多数据集，在官网帮助文档 可以看到全部
2 数据加载
我们会使用416b这个数据
复制 library(scRNAseq)
sce.416b <- LunSpikeInData(which="416b")
> sce.416b
# class: SingleCellExperiment
# dim: 46604 192
# metadata(0):
# assays(1): counts
# rownames(46604): ENSMUSG00000102693
# ENSMUSG00000064842 ... ENSMUSG00000095742
# CBFB-MYH11-mcherry
# rowData names(1): Length
# colnames(192):
# SLX-9555.N701_S502.C89V9ANXX.s_1.r_1
# SLX-9555.N701_S503.C89V9ANXX.s_1.r_1 ...
# SLX-11312.N712_S508.H5H5YBBXX.s_8.r_1
# SLX-11312.N712_S517.H5H5YBBXX.s_8.r_1
# colData names(9): Source Name cell line ...
# spike-in addition block
# reducedDimNames(0):
# altExpNames(2): ERCC SIRV
> table(sce.416b$block)
# 20160113 20160325
# 96 96
# 设置分组信息（因为来自2个96孔板），也就是为后面的批次处理做铺垫
sce.416b$block <- factor(sce.416b$block)
# 当前关于行（基因）的信息只有基因长度
> rowData(sce.416b)
# DataFrame with 46604 rows and 1 column
# Length
# <integer>
# ENSMUSG00000102693 1070
# ENSMUSG00000064842 110
# ENSMUSG00000051951 6094
# ENSMUSG00000102851 480
# ENSMUSG00000103377 2819
# ... ...
# ENSMUSG00000094621 121
# ENSMUSG00000098647 99
# ENSMUSG00000096730 3077
# ENSMUSG00000095742 243
# CBFB-MYH11-mcherry 2998 增加基因的信息（symbol ID、染色体位置）
首先是基因ID转换
这个函数是我非常喜欢的，方便了基因ID转换
将Ensembl ID转为Symbol，使用uniquifyFeatureNames，参数很简单，第一个是ID，第二是names
This function will attempt to use names if it is unique. If not, it will append the _ID to any non-unique value of names. Missing names will be replaced entirely by ID.
含义就是：
如果有symbol name的，它会将Ensembl ID替换为symbol name，没有name的仍然使用ID；
如果name相同、ID不同（即两个Ensembl ID对应同一个Symbol name），它会将ID与name组合（name_id）保证特异性
复制 library(AnnotationHub)
ens.mm.v97 <- AnnotationHub()[["AH73905"]]
columns(ens.mm.v97)
# [1] "DESCRIPTION" "ENTREZID" "EXONID" "EXONIDX" "EXONSEQEND"
# [6] "EXONSEQSTART" "GENEBIOTYPE" "GENEID" "GENEIDVERSION" "GENENAME"
# [11] "GENESEQEND" "GENESEQSTART" "INTERPROACCESSION" "ISCIRCULAR" "PROTDOMEND"
# [16] "PROTDOMSTART" "PROTEINDOMAINID" "PROTEINDOMAINSOURCE" "PROTEINID" "PROTEINSEQUENCE"
# [21] "SEQCOORDSYSTEM" "SEQLENGTH" "SEQNAME" "SEQSTRAND" "SYMBOL"
# [26] "TXBIOTYPE" "TXCDSSEQEND" "TXCDSSEQSTART" "TXID" "TXIDVERSION"
# [31] "TXNAME" "TXSEQEND" "TXSEQSTART" "TXSUPPORTLEVEL" "UNIPROTDB"
# [36] "UNIPROTID" "UNIPROTMAPPINGTYPE"
rowData(sce.416b)$ENSEMBL <- rownames(sce.416b)
rowData(sce.416b)$SYMBOL <- mapIds(ens.mm.v97, keys=rownames(sce.416b),
keytype="GENEID", column="SYMBOL") 如果需要把这个函数应用在日常使用中，也很简单，它只需要两个参数：一个ID，一个ID对应的名字
复制 # 一个ID
head(rowData(sce.416b)$ENSEMBL)
# [1] "ENSMUSG00000102693" "ENSMUSG00000064842"
# [3] "ENSMUSG00000051951" "ENSMUSG00000102851"
# [5] "ENSMUSG00000103377" "ENSMUSG00000104017"
# 一个ID对应的名字(这一步是问题的根源：存在无对应NA或者一对多、多对一的问题)
head(rowData(sce.416b)$SYMBOL)
# ENSMUSG00000102693 ENSMUSG00000064842
# "4933401J01Rik" "Gm26206"
# ENSMUSG00000051951 ENSMUSG00000102851
# "Xkr4" "Gm18956"
# ENSMUSG00000103377 ENSMUSG00000104017
# "Gm37180" "Gm37363"
uniquifyFeatureNames(head(rowData(sce.416b)$ENSEMBL),
head(rowData(sce.416b)$SYMBOL)) 然后增加基因的染色体编号（方便后面识别线粒体基因）
复制 rowData(sce.416b)$SEQNAME <- mapIds(ens.mm.v97, keys=rownames(sce.416b),
keytype="GENEID", column="SEQNAME")
> table(rowData(sce.416b)$SEQNAME=='MT')
# FALSE TRUE
# 46004 37 3 质控
备份数据
一般来说，有spike-in其实可以不需要线粒体过滤（因为它们的目的一致），但是这里还是加上了，双重保险
复制 mito <- which(rowData(sce.416b)$SEQNAME=="MT")
# 先计算各个细胞的QC指标
stats <- perCellQCMetrics(sce.416b, subsets=list(Mt=mito))
colnames(stats)
# [1] "sum" "detected"
# [3] "percent_top_50" "percent_top_100"
# [5] "percent_top_200" "percent_top_500"
# [7] "subsets_Mt_sum" "subsets_Mt_detected"
# [9] "subsets_Mt_percent" "altexps_ERCC_sum"
# [11] "altexps_ERCC_detected" "altexps_ERCC_percent"
# [13] "altexps_SIRV_sum" "altexps_SIRV_detected"
# [15] "altexps_SIRV_percent" "total"
# 再根据QC指标进行判断，哪些该去除
qc <- quickPerCellQC(stats, percent_subsets=c("subsets_Mt_percent",
"altexps_ERCC_percent"), batch=sce.416b$block)
sce.416b <- sce.416b[,!qc$discard]
> dim(sce.416b)
[1] 46604 185
> dim(unfiltered)
[1] 46604 192 根据原来的数据，加上质控标准作图
复制 colData(unfiltered) <- cbind(colData(unfiltered), stats)
unfiltered$discard <- qc$discard
gridExtra::grid.arrange(
plotColData(unfiltered, x="block", y="sum",
colour_by="discard") + scale_y_log10() + ggtitle("Total count"),
plotColData(unfiltered, x="block", y="detected",
colour_by="discard") + scale_y_log10() + ggtitle("Detected features"),
plotColData(unfiltered, x="block", y="subsets_Mt_percent",
colour_by="discard") + ggtitle("Mito percent"),
plotColData(unfiltered, x="block", y="altexps_ERCC_percent",
colour_by="discard") + ggtitle("ERCC percent"),
nrow=2,
ncol=2
) 再看下文库大小和ERCC分别和线粒体含量的关系
复制 gridExtra::grid.arrange(
plotColData(unfiltered, x="sum", y="subsets_Mt_percent",
colour_by="discard") + scale_x_log10(),
plotColData(unfiltered, x="altexps_ERCC_percent", y="subsets_Mt_percent",
colour_by="discard"),
ncol=2
) 然后检查一下被过滤的原因
复制 colSums(as.matrix(qc))
## low_lib_size low_n_features high_subsets_Mt_percent
## 5 0 2
## high_altexps_ERCC_percent discard
## 2 7 4 归一化
之前在 3.2 归一化 中介绍：去卷积方法计算size factor时，操作是这样的
先混合：Pool counts from many cells to increase the size of the counts
后去卷积：Pool-based size factors are then “deconvolved” into cell-based factors for normalization of each cell’s expression profile.
复制 set.seed(100)
# 这里一看有pre-clustering就知道是利用的去卷积的size factor
clust.zeisel <- quickCluster(sce.zeisel)
sce.zeisel <- computeSumFactors(sce.zeisel, cluster=clust.zeisel, min.mean=0.1)
sce.zeisel <- logNormCounts(sce.zeisel) 但是这里的数据集由于细胞数量较少，并且所有的细胞都是源自一个细胞系，所以不需要提前进行聚类quickCluster操作，直接进行去卷积即可。目的只有一个，将技术因素（如测序）导致的文库大小差异降到最低
复制 library(scran)
sce.416b <- computeSumFactors(sce.416b)
sce.416b <- logNormCounts(sce.416b)
summary(sizeFactors(sce.416b))
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.347 0.711 0.921 1.000 1.152 3.604 看看两种归一化方法的差异
复制 # 常规：最简单的只考虑文库大小
summary(librarySizeFactors(sce.416b))
# Min. 1st Qu. Median Mean 3rd Qu. Max.
# 0.3996 0.7601 0.9406 1.0000 1.1207 3.6720
plot(librarySizeFactors(sce.416b), sizeFactors(sce.416b), pch=16,
xlab="Library size factors", ylab="Deconvolution factors",
# 这是一个不错的主意，可以学习，直接将不同类型赋予不同颜色
# +1操作是为了将原来的False（0）、True（1）变成1、2
# 因此得到的black就是uninduced，red就是induced
col=c("black", "red")[grepl("induced", sce.416b$phenotype)+1],
log="xy") 看到去卷积的方法的确能暴露出细胞类型对文库矫正的影响，因为如果只考虑文库大小，那么常规的方法和去卷积方法应该结果一致，也不会出现红点和黑点分离的情况，但现在既然由于细胞类型不同而导致size factor出现了偏差，那么就是去卷积方法还考虑了其他因素，做的更加精细。
当然，使用常规的文库大小归一化方法对细胞分群和鉴定每群的marker基因是足够的。这里使用去卷积得到更准确的结果，虽然对分群改善不大，但对于每个基因的表达量数据估计与解释还是很重要的，因为它考虑了细胞类型组成差异的影响。
5 找表达量高变化基因
这里会使用到之前设置的批次信息，即：用spike-in来拟合更准确的技术偏差，同时考虑上批次信息【在之前 3.3 挑选高变化基因 的2.4 考虑批次效应部分提到】
复制 dec.416b <- modelGeneVarWithSpikes(sce.416b, "ERCC", block=sce.416b$block)
chosen.hvgs <- getTopHVGs(dec.416b, prop=0.1) 分别作图
复制 par(mfrow=c(1,2))
blocked.stats <- dec.416b$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)
} 黑点是基因，红点是spike-in，蓝线是根据红点画的拟合线。看到这里的两个批次之间差异很小，表示重复效果不错
6 批次矫正
从上面的效果看，这两个细胞板之间的细胞细胞组成应该是差不多的，因此使用简单的批次处理即可，不需要动用太复杂的方法；如果数据集比较大，需要考虑的批次效应更多，可能要动用 batchelor包的regressBatches()
复制 library(limma)
assay(sce.416b, "corrected") <- removeBatchEffect(logcounts(sce.416b),
design=model.matrix(~sce.416b$phenotype), batch=sce.416b$block) 7 降维
这个小数据集我们估计也不会有太大的异质性，因此设置10个主成分即可
然后使用精确的 SVD方法，来避免 irlba包对于处理小数据可能会出现的一些warning信息
复制 sce.416b <- runPCA(sce.416b, ncomponents=10, subset_row=chosen.hvgs,
exprs_values="corrected", BSPARAM=BiocSingular::ExactParam())
set.seed(1010)
sce.416b <- runTSNE(sce.416b, dimred="PCA", perplexity=10) 8 聚类
复制 # 看到这里使用了层次聚类的方法
my.dist <- dist(reducedDim(sce.416b, "PCA"))
my.tree <- hclust(my.dist, method="ward.D2")
library(dynamicTreeCut)
my.clusters <- unname(cutreeDynamic(my.tree, distM=as.matrix(my.dist),
minClusterSize=10, verbose=0))
colLabels(sce.416b) <- factor(my.clusters) 看一下现在的分群和批次之间的联系
复制 table(Cluster=colLabels(sce.416b), Plate=sce.416b$block)
## Plate
## Cluster 20160113 20160325
## 1 40 38
## 2 37 32
## 3 10 14
## 4 6 8 再比较一下现在的分群和表型之间的联系
复制 ## Oncogene
## Cluster induced CBFB-MYH11 oncogene expression wild type phenotype
## 1 78 0
## 2 0 69
## 3 1 23
## 4 14 0 作图
复制 plotTSNE(sce.416b, colour_by="label") 可以借助”轮廓图“（silhouette width）检查分群的质量
会对每个细胞都计算一个silhouette width值，如果一个细胞的width值为正并且越大 ，表示相对于其他亚群的细胞，这个细胞和它所在亚群中的细胞更接近，分群效果越好 ；如果width为负，就表示这个亚群的这个细胞和其他亚群的细胞更接近，即分群效果不太理想。
复制 library(cluster)
sil <- silhouette(my.clusters, dist = my.dist)
> colnames(sil)
[1] "cluster" "neighbor" "sil_width"
# 设置每个cluster的颜色
clust.col <- scater:::.get_palette("tableau10medium")
# 下面这一句的意思是：如果sil_width>0，就属于和自己接近，否则属于和其他亚群接近
sil.cols <- clust.col[ifelse(sil[,3] > 0, sil[,1], sil[,2])]
sil.cols <- sil.cols[order(-sil[,1], sil[,3])]
plot(sil, main = paste(length(unique(my.clusters)), "clusters"),
border=sil.cols, col=sil.cols, do.col.sort=FALSE) 图中可以反映的问题：
这个图是对上面得到的各个cluster中的细胞做的barplot，每个cluster是一种颜色。
这里看到cluster2都是正值并且横坐标width数值还是最大的，因此说明了分群没有问题，之前t-SNE的结果的确是它本身的问题
如果发现width值较小，表示分群结果一般，还有可能是分群过度，本来属于一个群的，又被拆分成小群
利用这个图，我们就能调整之前的参数，来调整分群效果，不过也不需要太过纠结完美的分群结果。因为即使图上看似合理的分群，可能实际上也不会得到更多的生物信息
9 找marker基因并解释结果
复制 markers <- findMarkers(sce.416b, my.clusters, block=sce.416b$block)
# 找cluster1的marker基因
marker.set <- markers[["1"]]
# cluster1 的Top10
head(marker.set, 10) 使用cluster1的Top10基因（但不一定只是10个）画热图
看到这些基因在cluster1的分布确实和其他几个clusters不同，表示差异基因找的比较可靠；另外看到，cluster 1含有致癌细胞，而且和DNA复制和细胞周期相关基因有强烈的下调，道理上也说得通，毕竟细胞衰老也是诱导癌细胞的一个因素
复制 top.markers <- rownames(marker.set)[marker.set$Top <= 10]
> length(top.markers)
[1] 29
plotHeatmap(sce.416b, features=top.markers, order_columns_by="label",
colour_columns_by=c("label", "block", "phenotype"),
center=TRUE, symmetric=TRUE, zlim=c(-5, 5)) 看着是不是很像pheatmap的结果，其实看一下这个函数，就是基于pheatmap做的
