# 4.1 实战一 | Smart-seq2 | 小鼠骨髓

## 1 前言

使用的是来自永生化小鼠骨髓祖细胞系的2个96孔板的416B cells（[Lun et al. 2017](https://pubmed.ncbi.nlm.nih.gov/29030468/)），并且在细胞裂解后文库制备前，在每个细胞中加入一定量的外源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包](https://www.jianshu.com/p/cc621292d29e)，但后来这个包进行了更新，原来的数据获取方法也发生了变化

```r
# 原来使用
data(fluidigm)
# 现在使用
ReprocessedFluidigmData() #provides 65 cells 
ReprocessedTh2Data() #provides 96 T helper cells 
ReprocessedAllenData() #provides 379 cells from the mouse visual cortex
```

另外新版本的包还加入了其他很多数据集，在[官网帮助文档](https://bioconductor.org/packages/release/data/experiment/vignettes/scRNAseq/inst/doc/scRNAseq.html#available-data-sets)可以看到全部

![](https://jieandze1314-1255603621.cos.ap-guangzhou.myqcloud.com/blog/2020-07-19-020516.png)

## 2 数据加载

我们会使用416b这个数据

```r
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`）保证特异性
>
>   ![image-20200719110335976](https://jieandze1314-1255603621.cos.ap-guangzhou.myqcloud.com/blog/2020-07-19-030336.png)

```r
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对应的名字**

```r
# 一个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))
```

**然后增加基因的染色体编号（方便后面识别线粒体基因）**

```r
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 质控

### 备份数据

```r
unfiltered <- sce.416b
```

一般来说，有spike-in其实可以不需要线粒体过滤（因为它们的目的一致），但是这里还是加上了，双重保险

```r
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
```

### **根据原来的数据，加上质控标准作图**

```r
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
)
```

![](https://jieandze1314-1255603621.cos.ap-guangzhou.myqcloud.com/blog/2020-07-19-032210.png)

### **再看下文库大小和ERCC分别和线粒体含量的关系**

```r
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
)
```

![](https://jieandze1314-1255603621.cos.ap-guangzhou.myqcloud.com/blog/2020-07-19-032650.png)

### **然后检查一下被过滤的原因**

```r
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 归一化](https://jieandze1314.osca.top/03/03-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.

```r
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`操作，直接进行去卷积即可。目的只有一个，将技术因素（如测序）导致的文库大小差异降到最低

```r
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
```

### **看看两种归一化方法的差异**

```r
# 常规：最简单的只考虑文库大小
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")
```

![](https://jieandze1314-1255603621.cos.ap-guangzhou.myqcloud.com/blog/2020-07-19-040144.png)

看到去卷积的方法的确能暴露出细胞类型对文库矫正的影响，因为如果只考虑文库大小，那么常规的方法和去卷积方法应该结果一致，也不会出现红点和黑点分离的情况，但现在既然由于细胞类型不同而导致size factor出现了偏差，那么就是去卷积方法还考虑了其他因素，做的更加精细。

当然，使用常规的文库大小归一化方法对细胞分群和鉴定每群的marker基因是足够的。这里使用去卷积得到更准确的结果，虽然对分群改善不大，但对于每个基因的表达量数据估计与解释还是很重要的，因为它考虑了细胞类型组成差异的影响。

## 5 找表达量高变化基因

这里会使用到之前设置的批次信息，即：用spike-in来拟合更准确的技术偏差，同时考虑上批次信息【在之前 [3.3 挑选高变化基因](https://jieandze1314.osca.top/03/03-3) 的2.4 考虑批次效应部分提到】

```r
dec.416b <- modelGeneVarWithSpikes(sce.416b, "ERCC", block=sce.416b$block)
chosen.hvgs <- getTopHVGs(dec.416b, prop=0.1)
```

分别作图

```r
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，蓝线是根据红点画的拟合线。看到这里的两个批次之间差异很小，表示重复效果不错

![](https://jieandze1314-1255603621.cos.ap-guangzhou.myqcloud.com/blog/2020-07-19-041126.png)

#### 6 批次矫正

从上面的效果看，这两个细胞板之间的细胞细胞组成应该是差不多的，因此使用简单的批次处理即可，不需要动用太复杂的方法；如果数据集比较大，需要考虑的批次效应更多，可能要动用 batchelor包的`regressBatches()`

```r
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信息

```r
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 聚类

```r
# 看到这里使用了层次聚类的方法
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)
```

**看一下现在的分群和批次之间的联系**

```r
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
```

**再比较一下现在的分群和表型之间的联系**

```r
##        Oncogene
## Cluster induced CBFB-MYH11 oncogene expression wild type phenotype
##       1                                     78                   0
##       2                                      0                  69
##       3                                      1                  23
##       4                                     14                   0
```

**作图**

```r
plotTSNE(sce.416b, colour_by="label")
```

![](https://jieandze1314-1255603621.cos.ap-guangzhou.myqcloud.com/blog/2020-07-19-062356.png)

### **可以借助”轮廓图“（silhouette width）检查分群的质量**

会对每个细胞都计算一个silhouette width值，如果一个细胞的**width值为正并且越大**，表示相对于其他亚群的细胞，这个细胞和它所在亚群中的细胞更接近，**分群效果越好**；如果width为负，就表示这个亚群的这个细胞和其他亚群的细胞更接近，即分群效果不太理想。

```r
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)
```

![](https://jieandze1314-1255603621.cos.ap-guangzhou.myqcloud.com/blog/2020-07-02-092004.png)

图中可以反映的问题：

* 这个图是对上面得到的各个cluster中的细胞做的barplot，每个cluster是一种颜色。
* 这里看到cluster2都是正值并且横坐标width数值还是最大的，因此说明了分群没有问题，之前t-SNE的结果的确是它本身的问题
* 如果发现width值较小，表示分群结果一般，还有可能是分群过度，本来属于一个群的，又被拆分成小群
* 利用这个图，我们就能调整之前的参数，来调整分群效果，不过也不需要太过纠结完美的分群结果。因为即使图上看似合理的分群，可能实际上也不会得到更多的生物信息

## 9 找marker基因并解释结果

```r
markers <- findMarkers(sce.416b, my.clusters, block=sce.416b$block)
# 找cluster1的marker基因
marker.set <- markers[["1"]]
# cluster1 的Top10
head(marker.set, 10)
```

![](https://jieandze1314-1255603621.cos.ap-guangzhou.myqcloud.com/blog/2020-07-19-063547.png)

**使用cluster1的Top10基因（但不一定只是10个）画热图**

看到这些基因在cluster1的分布确实和其他几个clusters不同，表示差异基因找的比较可靠；另外看到，cluster 1含有致癌细胞，而且和DNA复制和细胞周期相关基因有强烈的下调，道理上也说得通，毕竟细胞衰老也是诱导癌细胞的一个因素

```r
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))
```

![](https://jieandze1314-1255603621.cos.ap-guangzhou.myqcloud.com/blog/2020-07-19-063741.png)

看着是不是很像pheatmap的结果，其实看一下这个函数，就是基于pheatmap做的

![](https://jieandze1314-1255603621.cos.ap-guangzhou.myqcloud.com/blog/2020-07-19-063840.png)


---

# Agent Instructions: Querying This Documentation

If you need additional information that is not directly available in this page, you can query the documentation dynamically by asking a question.

Perform an HTTP GET request on the current page URL with the `ask` query parameter:

```
GET https://jieandze1314.osca.top/04/04-1.md?ask=<question>
```

The question should be specific, self-contained, and written in natural language.
The response will contain a direct answer to the question and relevant excerpts and sources from the documentation.

Use this mechanism when the answer is not explicitly present in the current page, you need clarification or additional context, or you want to retrieve related documentation sections.