# 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)