刘小泽写于19.8.20-第三单元第十讲：使用Seurat包(2、3版本)
笔记目的：根据生信技能树的单细胞转录组课程探索smart-seq2技术相关的分析技术
课程链接在：http://jm.grazy.cn/index/mulitcourse/detail.html?cid=53
载入数据，创建对象
rm(list = ls())
Sys.setenv(R_MAX_NUM_DLLS=999)
## 首先载入文章的数据
load(file='../input.Rdata')
# 原始表达矩阵
counts=a
counts[1:4,1:4];dim(counts) # 2万多个基因，768个细胞(需要下一步进行过滤)
library(stringr)
# 样本信息
meta=df
> head(meta)
g plate n_g all
SS2_15_0048_A3 1 0048 2624 all
SS2_15_0048_A6 1 0048 2664 all
SS2_15_0048_A5 2 0048 3319 all
SS2_15_0048_A4 3 0048 4447 all
SS2_15_0048_A1 2 0048 4725 all
SS2_15_0048_A2 3 0048 5263 all
# 按行/基因检查：每个基因在多少细胞中有表达量
fivenum(apply(counts,1,function(x) sum(x>0) ))
boxplot(apply(counts,1,function(x) sum(x>0) ))
# 按列/样本检查：每个细胞中存在多少表达的基因
fivenum(apply(counts,2,function(x) sum(x>0) ))
hist(apply(counts,2,function(x) sum(x>0) ))
按行+按列检查结果
左图显示了：有50%的基因只在低于20个细胞中有表达量，还有许多没有表达量的：
> table(apply(counts,1,function(x) sum(x>0) )==0)
FALSE TRUE
17429 7153
# 存在7000多个基因在任何一个细胞中都没表达
右图显示了：大部分细胞都包含2000个以上的有表达的基因
开始使用CreateSeuratObject构建Seurat对象：
因为Seurat存在两个版本，并且应用都比较多，所以这里会列出两种版本的代码
# 使用Seurat V3版本（以下简写"V3"）(版本3有大量的参数从之前的genes改成了features)
sce <- CreateSeuratObject(counts = counts,
meta.data = meta,
min.cells = 5,
min.features = 2000,
project = "sce")
# 使用V2版本(以下简写"V2")
sce <- CreateSeuratObject(raw.data = counts,
meta.data =meta,
min.cells = 5,
min.genes = 2000,
project = "sce")
> sce
An object of class Seurat
14406 features across 693 samples within 1 assay
Active assay: RNA (14406 features)
# 其中min.cells和min.features两个参数实际上做了下面这个事：
table(apply(counts,2,function(x) sum(x>0) )>2000)
table(apply(counts,1,function(x) sum(x>0) )>4)
有了对象，怎么从中得到矩阵和样本信息呢？
# 得到表达矩阵
# V3
dim(sce@assays$RNA)
[1] 14406 693
# 或者更简单的：
> dim(sce[['RNA']])
[1] 14406 693
# V2
dim(sce@data)
# 得到样本信息
> head(sce@meta.data)
orig.ident nCount_RNA nFeature_RNA g plate n_g all
SS2_15_0048_A3 SS2 128784 3067 1 0048 2624 all
SS2_15_0048_A6 SS2 131208 3037 1 0048 2664 all
SS2_15_0048_A5 SS2 207004 3739 2 0048 3319 all
SS2_15_0048_A4 SS2 140323 5004 3 0048 4447 all
SS2_15_0048_A1 SS2 301161 5117 2 0048 4725 all
SS2_15_0048_A2 SS2 311213 5673 3 0048 5263 all
预处理之可视化
对feature信息可视化
# V2
VlnPlot(object = sce,
features.plot = c("nGene", "nUMI" ),
group.by = 'plate',
nCol = 2)
# V3
plot1 <- VlnPlot(object = sce,
features = c("nFeature_RNA", "nCount_RNA" ),
group.by = 'plate',
ncol = 2)
plot2 <- VlnPlot(object = sce,
features = c("nFeature_RNA", "nCount_RNA" ),
group.by = 'g',
ncol = 2)
CombinePlots(plots = list(plot1, plot2))
可以看到：
- 左边两张图中不管是基因数量还是reads count数量，两个细胞板都没有批次效应；
- 右边两张图是根据hclust的层次聚类结果分类，可以看到这个分类很有可能是根据基因数量来的，有的样本表达的基因数量多就分到一类，表达基因数量少的就分为另一类，不过总体数据量没有太大区别
image
找其他几个基因再看看批次：
# V3
VlnPlot(sce,group.by = 'plate',c("Gapdh","Bmp3","Brca1","Brca2","nFeature_RNA"))
# V2
VlnPlot(sce,group.by = 'plate',c("Gapdh","Bmp3","Brca1","Brca2","nGene"))
image
然后增加一列样本信息—ERCC，再作图
# V2
ercc.genes <- grep(pattern = "^ERCC-", x = rownames(x = sce@raw.data), value = TRUE) # 注意这类的value = TRUE设置，如果不设置这个，结果就返回匹配位置；设置了可以直接返回匹配结果
percent.ercc <- Matrix::colSums(sce@raw.data[ercc.genes, ]) / Matrix::colSums(sce@raw.data)
sce <- AddMetaData(object = sce, metadata = percent.ercc,
col.name = "percent.ercc")
VlnPlot(object = sce,
features.plot = c("nGene", "nUMI", "percent.ercc" ),
group.by = 'g',
nCol = 3)
# V3(将多个过程重新整合)，可以直接利用PercentageFeatureSet函数计算来自某一个feature的reads count数量占比，它直接得到的就是百分比数值(比如V2得到的是0.234，那么V3结果就直接是23.4)
sce[["percent.ercc"]] <- PercentageFeatureSet(sce, pattern = "^ERCC-")
VlnPlot(object = sce,
features = c("nGene", "nUMI", "percent.ercc" ),
group.by = 'g',
ncol = 3)
下面这个图第一张和第三张就明显是负相关：也说明了检测的ERCC占比越多，检测到的基因数目就越少（因为第二张图中总共的reads数是一定的嘛，大约就0.5M）
image
看两组基因的相关性
# V2
GenePlot(object = sce, gene1 = "nUMI", gene2 = "nGene")
GenePlot(object = sce, gene1 = "Brca1", gene2 = "Brca2")
# V3
plot1 <- FeatureScatter(sce, feature1 = "nFeature_RNA", feature2 = "nCount_RNA")
plot2 <- FeatureScatter(sce, feature1 = "Brca1", feature2 = "Brca2")
CombinePlots(plots = list(plot1, plot2))
image
看两个细胞的相关性
# V2
CellPlot(sce,sce@cell.names[3],
sce@cell.names[4],
do.ident = FALSE)
# V3
CellScatter(sce,
colnames(sce[['RNA']])[3],
colnames(sce[['RNA']])[4])
image
预处理之归一化
关于Seurat归一化原理，可以看这一篇：https://www.biorxiv.org/content/biorxiv/early/2019/03/18/576827.full.pdf
默认是进行一个全局的LogNormalize操作：log1p(value/colSums[cell-idx] *scale_factor) ，其中log1p指的是log(x + 1) ，得到的结果存储在：pbmc[["RNA"]]@data
sce <- NormalizeData(object = sce,
normalization.method = "LogNormalize",
scale.factor = 10000)
# 当然其中都是默认参数，于是直接写这种也是可以的
sce <- NormalizeData(sce)
# 检查一下归一化后的数据
# V2
sce@data[1:3,1:3]
# V3
> sce[['RNA']][1:3,1:3]
3 x 3 sparse Matrix of class "dgCMatrix"
SS2_15_0048_A3 SS2_15_0048_A6 SS2_15_0048_A5
0610005C13Rik . . .
0610007P14Rik . . 0.6256969
0610009B22Rik . . .
# 结果将0存储为.(为了节省空间，这是松散矩阵sparse Matrix的一个特点)
鉴定差异基因HVGs（highly variable features）
在细胞与细胞间进行比较，选择表达量差别最大的，FindVariableFeatures函数，会计算一个mean-variance结果
# V2
sce <- FindVariableGenes(object = sce,
mean.function = ExpMean,
dispersion.function = LogVMR )
length( sce@var.genes) # 876(自动寻找了876个差异基因)
# 默认值是：x.low.cutoff = 0.1, x.high.cutoff = 8, y.cutoff = 1，就是说取log后的平均表达量(x轴)介于0.1-8之间的；分散程度(y轴，即标准差)至少为1的
# V3
sce <- FindVariableFeatures(sce, selection.method = "vst")
可见V3相比V2做了较大的改动：
- V3计算
mean.function和FastLogVMR均采用了加速的FastExpMean、FastLogVMR模式 - V3横坐标范围设定参数改成：
mean.cutoff；纵坐标改成：dispersion.cutoff - 鉴定差异基因的算法包含三种：
vst(默认)、mean.var.plot、dispersion- vst：首先利用loess对 log(variance) 和log(mean) 拟合一条直线，然后利用观测均值和期望方差对基因表达量进行标准化。最后根据保留最大的标准化的表达量计算方差
- mean.var.plot: 首先利用mean.function和 dispersion.function分别计算每个基因的平均表达量和离散情况；然后根据平均表达量将基因们分散到一定数量（默认是20个）的小区间（bin）中，并且计算每个bin中z-score
- dispersion：挑选最高离散值的基因
- V3默认选择2000个差异基因，检查方法也不同（V3用
VariableFeatures(sce)检查，V2用sce@var.genes检查）
标准化
这里去除一些技术误差，例如UMI、ERCC，之后就不需要考虑ERCC的影响了，下面的代码中vars.to.regress就是对一些技术因素的排除
关于scale的操作过程：
ScaleData这个函数有两个默认参数：do.scale = TRUE和do.center = TRUE，然后需要输入进行scale/center的基因名（默认是取前面FindVariableFeatures的结果）。scale和center这两个默认参数应该不陌生了，center就是对每个基因在不同细胞的表达量都减去均值；scale就是对每个进行center后的值再除以标准差（就是进行了一个z-score的操作）
# V2
sce <- ScaleData(object = sce,
vars.to.regress = c("nUMI","percent.ercc" ))
# 结果存储在sce@scale.data中
# V3
all.genes <- rownames(sce)
length(rownames(sce))
sce <- ScaleData(sce, features = all.genes,
vars.to.regress = c("nCount_RNA","percent.ercc" ))
# 结果存储在sce[["RNA"]]@scale.data中
实际上，scale函数默认是只针对之前鉴定的HVGs进行标准化（版本3中默认得到2000个HVGs），这样操作的结果中降维和聚类不会被影响，但是只对HVGs基因进行标准化，下面画的热图可能会受到全部基因中某些极值的影响。所以为了热图的结果，还是对所有基因进行归一化比较好。
很多时候，我们会混淆normalization和scale：那么看一句英文解释：
Normalization “normalizes” within the cell for the difference in sequenicng depth / mRNA throughput. 主要着眼于样本的文库大小差异
Scaling “normalizes” across the sample for differences in range of variation of expression of genes . 主要着眼于基因的表达分布差异
降维之PCA
最常用的PCA方法是一种线性降维，它使用标准化的数据
# V2 使用自动检测得到的sce@var.genes(876个)；它默认分析20个主成分
sce <- RunPCA(object = sce, pc.genes = sce@var.genes,
do.print = TRUE, pcs.print = 1:5,
genes.print = 5)
# 结果在 sce@dr$pca@gene.loadings
# V3 默认选择之前鉴定的差异基因（2000个）作为input，但可以使用features进行设置；默认分析的主成分数量是50个；然后输出到屏幕上5个主成分，每个主成分
sce <- RunPCA(sce, features = VariableFeatures(object = sce),
ndims.print = 1:5, nfeatures.print = 5)
# 这50个主成分的全部结果在 sce@reductions$pca@feature.loadings
对print的主成分结果可视化（由于V2、V3默认分析的主成分数量不同，它们的结果也有略微差别）：
# V2 它利用的也就是sce@dr$pca@gene.loadings结果
VizPCA( sce, pcs.use = 1:2)
# V3，它利用的也就是sce@reductions$pca@feature.loadings结果
VizDimLoadings(pbmc, dims = 1:2, reduction = "pca")
V3得到的前两个主成分结果
按批次检查前两个主成分
# V2
PCAPlot(sce, dim.1 = 1, dim.2 = 2,group.by = 'plate')
# V3
DimPlot(sce, reduction = "pca",group.by = 'plate')
可以看到，即使前面没有对批次进行校正，这里也没有批次效应
V3结果
再看下层次聚类cluster结果在PCA的可视化
# V2
PCAPlot(sce, dim.1 = 1, dim.2 = 2,group.by = 'g')
# V3 (DimPlot可以通过降维算法选择，默认首选umap，其次tsne，最后pca，来替代PCAPlot、TSNEPlot、UMAPPlot)
DimPlot(sce, reduction = "pca",group.by = 'g')
# 或者直接 PCAPlot(sce, group.by = 'g')
可以看到几个层次聚类的cluster混在了一起，个人推测可能是：
- PCA默认针对前2000个HVGs进行分析，而层次聚类是针对全部基因做的而分成4群，PCA又是一个线性降维模式，因此分的不是特别开
V3结果
然后又探索了使用全部基因做PCA：
# V3
sce <- RunPCA(sce, features = all.genes,
ndims.print = 1:5, nfeatures.print = 5)
PCAPlot(sce, group.by = 'g')
image
好像结果没有太大变化，没关系先继续向下探索
探索异质性来源—DimHeatmap
每个细胞和基因都根据PCA结果得分进行了排序，默认画前30个基因（nfeatures设置），1个主成分(dims设置)，细胞数没有默认值（cells设置）
# V2
PCHeatmap(object = sce, pc.use = 1:15, cells.use = 100,
do.balanced = TRUE, label.columns = FALSE)
# V3
DimHeatmap(sce, dims = 1:15, cells = 100, balanced = TRUE)
# 其中balanced表示正负得分的基因各占一半
那么降维后选多少个主成分来代表整个数据集？
最简单使用ElbowPlot(sce)看一看，主要关注”肘部“的PC，它是一个转折点
ElbowPlot(sce)
细胞聚类
# V2
sce1 <- FindClusters(object = sce, reduction.type = "pca",
dims.use = 1:15, force.recalc = T,
resolution = 0.4, print.output = 0,
save.SNN = TRUE)
# V3
# 先根据ElbowPlot挑选了15个PCs，所以这里dims定义为15个
sce <- FindNeighbors(sce, dims = 1:15)
# 然后使用FindClusters() 进行聚类。其中包含了一个resolution的选项，它会设置一个”间隔“值，该值越大，间隔越大，得到的cluster数量越多
sce <- FindClusters(sce, resolution = 0.4)
# # 发现700细胞，设resolution=0.4时，得到Number of communities: 4
> table(sce@meta.data$RNA_snn_res.0.4)
0 1 2 3
434 169 52 38
关于这个resolution参数：分辨率越高，看的越清楚，看到的细节越丰富（cluster越多）；反之，如果分辨率调的很低，结果就看的模模糊糊一大坨
另一种降维方法—非线性降维（tSNE）
sce <- RunTSNE(object = sce, dims.use = 1:15, do.fast = TRUE)
# V3
DimPlot(sce,reduction = "tsne",label=T)
# V2
TSNEPlot(object = sce, do.label=T)
tsne结果
看一下tSNE分类结果和之前层次聚类结果比较：
> table(sce$RNA_snn_res.0.4,sce$g)
1 2 3 4
0 184 249 0 1
1 40 3 108 18
2 4 48 0 0
3 9 0 13 16
# 左侧0-3是tSNE结果，顶部1-4是hclust结果：hclust的第1组对应tsne的第0组；hclust的第2组还是对应tsne的第0组；hclust的最后一组没有完全分开，说明它们差异还是很大的
找marker基因
# 记住这里tsne的分组是0、1、2、3这四组，因此下面取得ident.1 = 1实际上是第2组
# min.pct的意思是：一个基因在任何两群细胞中的占比最低不能低于多少
cluster2.markers <- FindMarkers(sce, ident.1 = 1, min.pct = 0.25)
> head(cluster2.markers, n = 5)
p_val avg_logFC pct.1 pct.2 p_val_adj
Lsp1 3.359064e-89 1.228134 0.828 0.084 4.839067e-85
Pdgfra 2.317021e-84 1.842903 0.864 0.153 3.337901e-80
Mfap5 3.442206e-84 2.782743 0.917 0.240 4.958843e-80
Svep1 7.297568e-84 1.876058 0.834 0.132 1.051288e-79
Tspan11 3.860428e-83 0.876791 0.757 0.065 5.561333e-79
对找到的marker基因进行可视化：
-
VlnPlot：用小提琴图对某些基因在全部cluster中表达量进行绘制markers_genes = rownames(head(cluster2.markers, n = 5)) # V2 VlnPlot(object = sce, features.plot =markers_genes, use.raw = TRUE, y.log = TRUE) # V3 VlnPlot(pbmc, features = markers_genes)第二组的marker基因 -
FeaturePlot：最常用的可视化=》将基因表达量投射到降维聚类结果中# V2 FeaturePlot(object = sce, features.plot =markers_genes, cols.use = c("grey", "blue"), reduction.use = "tsne") # V3 FeaturePlot(sce, features = markers_genes)第2群marker基因映射结果
可视化文献作者给出的基因
会了基本操作以后，可以将文章中的4个细胞亚群的marker基因拿过来，看看它们分别在我们自己结果中的那一组
就是根据这样图：
原文的72个marker基因
如何快速提取其中的基因名呢？自己一个一个手动敲入是个办法，但难免会眼花：
# 第一步：鼠标复制(前提是pdf中的图片可以选择文字)，
# 因为基因数太多，这里就不全展示了，只展示一部分
tmp <- c('Atp1b2 Notch3 Ano1 Gm13889 Des
Aoc3
Ndu fa4l2 Gucy1a3
Esam
Gdpd3
Mcam
Higd1b
Cpe
...')
# 第二步：这样粘贴过来的一定存在换行符，就像这样：
> tmp
[1] "Atp1b2 Notch3 Ano1 Gm13889 Des\nAoc3\nNdu fa4l2 Gucy1a3\nEsam\nGdpd3\nMcam\nHigd1b\nCpe\nKcnj8\nAbcc9\nRgs4\nSparcl1\nRgs5\nSmoc2\nItgbl1\nFbln1\nCdh11\nCrabp1\nPdgf ra\nSvep1\nPdpn\nLsp1\nCpxm1\nLrrc15\nCilp\nDcn\nLum\nMfap5\nFbln2\nOlfml3\nRnase4\nMki67\nAnln\nCdca3 2810417H13Rik Tpx2\nCdca8\nFam64a\nNuf2\nBirc5\nCep55\nSka1\nKif15\nTtk\nMelk\nTop2a\nPbk\nCcna2\nSpc25\nTfap2b\nCol4a6\nTfap2c\nEps8l2\nExtl1\nAim1\nIrx1\nGata3\nCol9a1\nCol4a5\nChad\nSmoc1\nCol9a3\nScrg1\nMia\nCd24a\nPerp\nTrf"
# 那么就需要先将换行符换成空格
tmp <- gsub("\n"," ",tmp)
> tmp
[1] "Atp1b2 Notch3 Ano1 Gm13889 Des Aoc3 Ndu fa4l2 Gucy1a3 Esam Gdpd3 Mcam Higd1b Cpe Kcnj8 Abcc9 Rgs4 Sparcl1 Rgs5 Smoc2 Itgbl1 Fbln1 Cdh11 Crabp1 Pdgf ra Svep1 Pdpn Lsp1 Cpxm1 Lrrc15 Cilp Dcn Lum Mfap5 Fbln2 Olfml3 Rnase4 Mki67 Anln Cdca3 2810417H13Rik Tpx2 Cdca8 Fam64a Nuf2 Birc5 Cep55 Ska1 Kif15 Ttk Melk Top2a Pbk Ccna2 Spc25 Tfap2b Col4a6 Tfap2c Eps8l2 Extl1 Aim1 Irx1 Gata3 Col9a1 Col4a5 Chad Smoc1 Col9a3 Scrg1 Mia Cd24a Perp Trf"
# 第三步：分割字符串
library(stringr)
paper_marker <- as.character(str_split(tmp,' ',simplify = T))
# 第四步：检查错误
> length(paper_marker)
[1] 74
# 这个不对，原文只有4组*18=72个基因，多了两个，应该是复制粘贴过来出的错误(因为这里小鼠基因名都是首字母大写，于是先找到基因名首字母不是大写的，再将它们替换掉)
# 其中一个就是‘Ndufa4l2’本来是一个，但是粘贴过来成了两个：'Ndu','fa4l2'
paper_marker_1 <- str_replace(paper_marker,c('Ndu','fa4l2'),"Ndufa4l2")
# 每次都要检查
> setdiff(paper_marker_1,paper_marker)
[1] "Ndufa4l2"
# 另一个就是"Pdgf" 和"ra"
paper_marker_2 <- str_replace(paper_marker_1,c('Pdgf','ra'),"Pdgfra")
# 每次修改都要检查：
> setdiff(paper_marker_2,paper_marker_1)
[1] "CPdgfrabp1" "Pdgfra"
# 发现paper_marker_2中由于替换了一个简单的字符"ra"，所以原来的"Crabp1"也被替换了，修改回来即可
paper_marker_2[paper_marker_2=='CPdgfrabp1'] <- 'Crabp1'
paper_marker <- unique(paper_marker_2)
> length(paper_marker)
[1] 72
使用上面得到的基因进行FeaturePlot
# 原文中的第一群(18个基因中，我们只取前6个基因看一下)
FeaturePlot(sce, features = paper_marker[1:6])
结合我们上面tsne的结果看：
tsne结果
发现原文中的第一群细胞在我们这里也是大体分到了第一群，当然不是很完美的再现，因为发现这些基因还有少量映射到了旁边的第3群、第4群
原文中的第一群细胞在我们这里也是第一群
# 同理对其他三群
FeaturePlot(sce, features = paper_marker[19:24])
FeaturePlot(sce, features = paper_marker[37:42])
FeaturePlot(sce, features = paper_marker[55:60])
再回到自己分析的结果，找全部的marker
# V2
sce.markers <- FindAllMarkers(object = sce, only.pos = TRUE,
min.pct = 0.25,
thresh.use = 0.25)
# V3
sce.markers <- FindAllMarkers(sce, only.pos = TRUE,
min.pct = 0.25, logfc.threshold = 0.25)
# 这一步过滤好好理解(进行了分类、排序、挑每个cluster的前20个)
top20 <- sce.markers %>% group_by(cluster) %>% top_n(n = 20, wt = avg_logFC)
# 看看和原文挑选的基因重合度(约50%重合度，还不错)
> intersect(top20$gene,paper_marker)
[1] "Rgs5" "Sparcl1" "Rgs4" "Atp1b2" "Abcc9" "Kcnj8"
[7] "Cpe" "Des" "Ano1" "Esam" "Gucy1a3" "Mfap5"
[13] "Smoc2" "Itgbl1" "Cpxm1" "Dcn" "Fbln1" "Rnase4"
[19] "Cdh11" "Lum" "Olfml3" "Lrrc15" "Fbln2" "Crabp1"
[25] "Cilp" "Top2a" "Pbk" "Mki67" "Spc25" "Anln"
[31] "Ccna2" "Col9a1" "Col4a5" "Chad" "Col9a3" "Trf"
画自己数据的top20基因热图：
# V3
DoHeatmap(sce, features = top20$gene) + NoLegend()
# V2
DoHeatmap(object = sce, genes.use = top20$gene,
slim.col.label = TRUE, remove.key = TRUE)
自己数据的top20基因热图
