Seurat标准流程 / 开普饭

作者 | 单细胞天地小编刘小泽

课程链接在：http://jm.grazy.cn/index/mulitcourse/detail.html?cid=55
第二单元第7讲：走Seurat标准流程【文章结构总-分-总，结尾有完整的代码，熟悉者前面可以跳过，去看后面8min完成的代码】

前言

前面介绍了自己利用cellranger count的结果进行seurat分析，但是整合数据方面做得还是不如原作者优秀，虽然我们不知道他们是如何处理的，但还是可以继续向下进行，而且这一次将会使用他们合并好的数据

从下载GEO开始

https://www.ncbi.nlm.nih.gov//geo/query/acc.cgi?acc=GSE117988

下载：GSE117988_raw.expMatrix_PBMC.csv.gz
(ftp://ftp.ncbi.nlm.nih.gov/geo/series/GSE117nnn/GSE117988/suppl/GSE117988_raw.expMatrix_PBMC.csv.gz)

1rm(list = ls()) 2options(warn=-1) 3suppressMessages(library(Seurat)) 4 5# 读取表达矩阵 6start_time <- Sys.time() 7raw_dataPBMC <- read.csv('./GSE117988_raw.expMatrix_PBMC.csv.gz', header = TRUE, row.names = 1) 8end_time <- Sys.time() 9end_time - start_time 10#Time difference of 1.486653 mins 11 12> dim(raw_dataPBMC) 13[1] 17712 12874 14# 得到17712基因，12874细胞

复现作者的分群结果

下载后怎么进行归一化、时间点划分、创建对象、分群其实作者也已经给了答案：https://www.researchgate.net/publication/328016998_Supplementary_Material_6/data/5bb2eac9299bf13e605a0a74/41467-2018-6300-MOESM6-ESM.txt

首先对文库大小进行一个归一化

1dataPBMC <- log2(1 + sweep(raw_dataPBMC, 2, median(colSums(raw_dataPBMC))/colSums(raw_dataPBMC), '*'))

这个sweep函数很有趣

它可以在apply的统一操作基础上，增加个性化服务(比如对每行/列使用不同的数值进行计算，而且可以指定计算方法)

例如，想要将1-3行分别减去对应的行号1-3，就可以这样：

1> test <- matrix(1:12,ncol = 4,byrow = T) 2> test 3 [,1] [,2] [,3] [,4] 4[1,] 1 2 3 4 5[2,] 5 6 7 8 6[3,] 9 10 11 12 7> sweep(test,1,c(1,2,3),"-") 8 [,1] [,2] [,3] [,4] 9[1,] 0 1 2 3 10[2,] 3 4 5 6 11[3,] 6 7 8 9

第一个位置test这里需要是矩阵或数据框；
第二个位置1和2选一个，原理和apply一样
第三个位置是要操作的向量，如果要对行操作，那么这个向量长度就要和行数一样
第四个位置是计算符，比如：+ - * / < > 等

于是就能懂了这里的操作：先求每个细胞文库的总大小，然后用它的中位数除以总大小得到一个小数，然后按列乘以这个小数就相当于对文库进行了归一化，将文库本身的大小差异置之度外

自定义划分时间点

看看现在的细胞命名，其实我们看完文章是知道作者用了四个时间点的，当然，作者在做实验过程中，也会有自己的一套方法对不同时间点的细胞进行命名。比如作者这里的想法就是：

1> head(colnames(dataPBMC)) 2[1] "AAACCTGAGCGAAGGG.1" "AAACCTGAGGTCATCT.1" "AAACCTGAGTCCTCCT.1" 3[4] "AAACCTGCACCAGCAC.1" "AAACCTGGTAACGTTC.1" "AAACCTGGTAAGGATT.1" 4> tail(colnames(dataPBMC)) 5[1] "TTTGTCAAGCGAGAAA.4" "TTTGTCAAGGAATTAC.4" "TTTGTCAAGTGCGTGA.4" 6[4] "TTTGTCACACGAGGTA.4" "TTTGTCATCATTGCGA.4" "TTTGTCATCCACGCAG.4"

看到每个细胞barcode后面都有一个点号，然后接着数字1-4来区分四个时间点，并且是按照时间顺序来的，它们的对应顺序就是：

11 => PBMC_Pre 22 => PBMC_EarlyD27 33 => PBMC_RespD376 44 => PBMC_ARD614

作者的想法是：将数字提取出来，然后分别对应到具体时期名称（利用if else结构）

1# 怎么提取？ 2# 作者利用的是Seurat V2的ExtractField函数 3# for Seurat V2 4timePoints <- sapply(colnames(dataPBMC), 5 function(x) ExtractField(x, 2, '[.]')) 6# 但如果使用V3，就要用常规方法strsplit，并且注意点号是正则匹配符，需要用\\来转义 7timePoints <- sapply(colnames(dataPBMC), function(x) unlist(strsplit(x, "\\."))[2]) 8 9> table(timePoints) 10timePoints 11 1 2 3 4 122082 1592 4684 4516

完成对应

1timePoints <-ifelse(timePoints == '1', 'PBMC_Pre', 2 ifelse(timePoints == '2', 'PBMC_EarlyD27', 3 ifelse(timePoints == '3', 'PBMC_RespD376', 'PBMC_ARD614'))) 4 5> table(timePoints) 6timePoints 7 PBMC_ARD614 PBMC_EarlyD27 PBMC_Pre PBMC_RespD376 8 4516 1592 2082 4684

对表达矩阵进行质控

主要看两点：

1# 第一点：基因在多少细胞表达 2fivenum(apply(dataPBMC,1,function(x) sum(x>0) )) 3# RP4-669L17.10 LAMB3 NAT10 AC093673.5 RPL21 4# 1 6 50 207 12102 5boxplot(apply(dataPBMC,1,function(x) sum(x>0) )) 6 7# 第二点：细胞中有多少表达的基因 8fivenum(apply(dataPBMC,2,function(x) sum(x>0) )) 9# CAAGAAATCGATCCCT.2 GGCCGATTCCAGAGGA.3 TCAACGAAGAGCTGGT.3 10# 10 321 395 11# TGCCAAAGTCTGCGGT.4 TCTGGAATCTATCGCC.3 12# 481 2865 13hist(apply(dataPBMC,2,function(x) sum(x>0) ))

看到：在一万多个基因中，75%的基因只在200多个细胞中有表达；而在一万多个细胞中，75%的细胞也只表达不到500个基因。按照常理，10X的数据应该能做到平均表达800个基因

创建Seurat对象

1# Seurat V2 2PBMC <- CreateSeuratObject(raw.data = dataPBMC, 3 min.cells = 1, min.genes = 0, project = '10x_PBMC') 4# Seurat V3 5PBMC <- CreateSeuratObject(dataPBMC, 6 min.cells = 1, min.features = 0, project = '10x_PBMC') 7 8> PBMC 9An object of class Seurat 1017712 features across 12874 samples within 1 assay 11Active assay: RNA (17712 features)

添加metadata

1# Seurat V2 使用AddMetaData 2# 3.0版本可以直接使用 object$name <- vector，当然也可以用AddMetaData 3PBMC <- AddMetaData(object = PBMC, 4 metadata = apply(raw_dataPBMC, 2, sum), 5 col.name = 'nUMI_raw') 6PBMC <- AddMetaData(object = PBMC, metadata = timePoints, col.name = 'TimePoints')

聚类标准流程

标准化=》找高变异基因=》根据这些基因进行PCA降维=》根据PCA结果找分群=》TSNE降维=》可视化

1# Seurat V2 2PBMC <- ScaleData(object = PBMC, vars.to.regress = c('nUMI_raw'), model.use = 'linear', use.umi = FALSE) 3 4PBMC <- FindVariableGenes(object = PBMC, mean.function = ExpMean, dispersion.function = LogVMR, x.low.cutoff = 0.0125, x.high.cutoff = 3, y.cutoff = 0.5) 5 6PBMC <- RunPCA(object = PBMC, pc.genes = PBMC@var.genes) 7 8PBMC <- FindClusters(object = PBMC, reduction.type = "pca", dims.use = 1:10, resolution = 1, k.param = 35, save.SNN = TRUE) # 13 clusters 9 10PBMC <- RunTSNE(object = PBMC, dims.use = 1:10) 11 12TSNEPlot(PBMC, colors.use = c('green4', 'pink', '#FF7F00', 'orchid', '#99c9fb', 'dodgerblue2', 'grey30', 'yellow', 'grey60', 'grey', 'red', '#FB9A99', 'black'))

如果使用Seurat V3，会发生一些变化

3.0版本将FindVariableGenes换为FindVariableFeatures，另外将原来的cutoff进行整合，x轴统一归到mean.cutoff中，y轴归到dispersion.cutoff中

1PBMC <- FindVariableFeatures(object = PBMC, mean.function = ExpMean, dispersion.function = LogVMR, mean.cutoff = c(0.0125,3), dispersion.cutoff = c(0.5,Inf))
3.0版本将FindClusters拆分为FindNeighbors和FindClusters，而版本2只有一个函数FindClusters

1PBMC <- FindNeighbors(PBMC, reduction = "pca", dims = 1:10, 2 k.param = 35) 3PBMC <- FindClusters(object = PBMC, 4 resolution = 0.9, verbose=F)
3.0版本在最后的可视化上可以使用DimPlot、TSNEPlot，颜色参数变成了cols

1DimPlot(PBMC, cols = c('green4', 'pink', '#FF7F00', 'orchid', '#99c9fb', 'dodgerblue2', 'grey30', 'yellow', 'grey60', 'grey', 'red', '#FB9A99', 'black'))

3.0的结果

【注意】如果参数使用不正确，DimPlot或TSNEPlot会调用默认颜色设置，例如使用参数colors或colors.use都是V3不识别的，因此它不会按照我们的颜色操作，而是生成类似这种：

保存对象

1save(PBMC,file = 'patient1.PBMC.output.Rdata') 2# 结果有1.9G

附加一个问题：Seurat2、3得到的结果差别大吗？

就看最后的TSNE聚类图，上面用Seurat3操作了一遍，虽然使用了和作者一样的数据，但结果和原文还是差别很大。

原文聚类结果

那么原因真的是由于版本引起的吗？使用作者的V2会不会好一些？

下面进行V2测试

8min跑一遍

1rm(list = ls()) 2options(warn=-1) 3suppressMessages(library(Seurat)) 4 5############################# 6# 读取表达矩阵 7############################# 8start_time <- Sys.time() 9raw_dataPBMC <- read.csv('./GSE117988_raw.expMatrix_PBMC.csv.gz', header = TRUE, row.names = 1) 10 11## step1: 归一化 12dataPBMC <- log2(1 + sweep(raw_dataPBMC, 2, 13 median(colSums(raw_dataPBMC))/colSums(raw_dataPBMC), '*')) # Normalization 14 15## step2: 自定义划分时间点 16# 作者利用的是Seurat V2的ExtractField函数 17# for Seurat V2 18timePoints <- sapply(colnames(dataPBMC), 19 function(x) ExtractField(x, 2, '[.]')) 20 21timePoints <-ifelse(timePoints == '1', 'PBMC_Pre', 22 ifelse(timePoints == '2', 'PBMC_EarlyD27', 23 ifelse(timePoints == '3', 'PBMC_RespD376', 'PBMC_ARD614'))) 24 25## step3: 表达矩阵质控 26# 第一点：基因在多少细胞表达 27fivenum(apply(dataPBMC,1,function(x) sum(x>0) )) 28# 第二点：细胞中有多少表达的基因 29fivenum(apply(dataPBMC,2,function(x) sum(x>0) )) 30 31## step4: 创建Seurat对象 32PBMC <- CreateSeuratObject(dataPBMC, 33 min.cells = 1, min.features = 0, project = '10x_PBMC') 34PBMC # 17,712 genes and 12,874 cells 35 36## step5: 添加metadata (nUMI and timePoints) 37PBMC <- AddMetaData(object = PBMC, metadata = apply(raw_dataPBMC, 2, sum), col.name = 'nUMI_raw') 38PBMC <- AddMetaData(object = PBMC, metadata = timePoints, col.name = 'TimePoints') 39 40## step6: 聚类标准流程 41# Seurat V2 42PBMC <- ScaleData(object = PBMC, vars.to.regress = c('nUMI_raw'), model.use = 'linear', use.umi = FALSE) 43PBMC <- FindVariableGenes(object = PBMC, mean.function = ExpMean, dispersion.function = LogVMR, x.low.cutoff = 0.0125, x.high.cutoff = 3, y.cutoff = 0.5) 44PBMC <- RunPCA(object = PBMC, pc.genes = PBMC@var.genes) 45PBMC <- FindClusters(object = PBMC, reduction.type = "pca", dims.use = 1:10, resolution = 1, k.param = 35, save.SNN = TRUE) # 13 clusters 46PBMC <- RunTSNE(object = PBMC, dims.use = 1:10) 47TSNEPlot(PBMC, colors.use = c('green4', 'pink', '#FF7F00', 'orchid', '#99c9fb', 'dodgerblue2', 'grey30', 'yellow', 'grey60', 'grey', 'red', '#FB9A99', 'black')) 48 49end_time <- Sys.time() 50end_time - start_time 51# Time difference of 8.085459 mins

Seurat V2结果

探究一下Seurat2和3的分群结果

起初我认为这两个版本的差异蛮大的，因为看tsne图明显感觉V2更好一些，导致我得出了错误的结论，认为两个版本的包处理结果千差万别。实际上，站在作者的角度思考一下，他也不会允许自己的“孩子”在成长过程中出现太大的偏差

学到了一点：图片不能说明问题，需要用数据来表现它们之间的分群结果到底差异大不大

下面👇就用数据来证明：

V2的Seurat得到了13群，V3利用resolution = 0.9得到了13群，而使用和V2一样的resolution = 1 会得到14群

主要使用table看一下比较结果

如果直接使用table(PBMC@meta.data$res.1)，那么就会出现以下情况：

1> table(PBMC@meta.data$res.1) 2 3 0 1 10 11 12 2 3 4 5 6 7 8 9 41877 1859 219 212 191 1569 1414 1322 1321 1123 772 595 400 5# 发现虽然是13群，但10、11、12出现在了1以后，很明显这是由于字符串ascii排序的原因，在linux中常见 6> class(PBMC@meta.data$res.1) 7[1] "character" 8# 使用数值型就会和平常一样 9> table(as.numeric(PBMC@meta.data$res.1)) 10 11 0 1 2 3 4 5 6 7 8 9 10 11 12 121877 1859 1569 1414 1322 1321 1123 772 595 400 219 212 191

横着的0-12是V2得到的PBMC，竖着的0-12是V3得到的PBMC_V3。看到V2的第1群在V3中分成了第4和第2群【在图中显示就是原文图中粉色群 = V3得到的橙色和浅蓝色（如下图）】；V2的第6群在V3中分成了第6和9群

1> table(PBMC_V3$RNA_snn_res.0.9,as.numeric(PBMC@meta.data$res.1)) 2 3 0 1 2 3 4 5 6 7 8 9 10 11 12 4 0 1873 1 0 0 0 0 0 0 3 7 0 0 0 5 1 0 0 224 0 1094 277 0 0 10 1 3 0 1 6 2 0 420 1101 0 37 2 0 0 4 2 40 0 0 7 3 0 0 0 1379 0 0 0 82 60 17 0 5 0 8 4 0 1418 53 0 0 0 0 0 0 7 15 0 0 9 5 0 2 47 0 93 1029 0 0 6 0 1 1 0 10 6 0 0 1 0 1 0 887 0 0 0 0 0 32 11 7 0 0 0 13 0 0 0 648 0 6 0 2 0 12 8 1 17 3 3 21 11 0 12 500 5 0 0 0 13 9 1 1 0 1 0 0 236 0 2 0 0 0 158 14 10 2 0 1 11 1 0 0 12 1 354 2 1 0 15 11 0 0 139 0 75 2 0 0 9 1 158 0 0 16 12 0 0 0 7 0 0 0 18 0 0 0 203 0

但是，这个结果说明了：V2、V3处理同一组数据，结果其实相差并不大。并不是像我们想象的，图片差距大，结果差距就大。真实情况可能是本来的某一群在另一种处理中被分得离散一些。但是tSNE本来就是这样，图中的距离并不代表真实的差异，它的运行次数会直接导致最后的图片形态不同

关于tsne这个流行的算法，有必要了解一下：

tsne的作者Laurens强调，可以通过t-SNE的可视化图提出一些假设，但是不要用t-SNE来得出一些结论，想要验证你的想法，最好用一些其他的办法。
t-SNE中集群之间的距离并不表示相似度 ，同一个数据上运行t-SNE算法多次，很有可能得到多个不同“形态”的集群。但话说回来，真正有差异的群体之间，不管怎么变换形态，它们还是有差别
关于perplexity的使用：(默认值是30) 如果忽视了perplexity带来的影响，有的时候遇到t-SNE可视化效果不好时，对于问题无从下手。perplexity表示了近邻的数量，例如设perplexity为2，那么就很有可能得到很多两个一对的小集群。
有的时候会出现同一集群被分为两半的情况，但群间的距离并不能说明什么，解决这个问题，只需要跑多次找出效果最好的就可以了

引用自：https://bindog.github.io/blog/2018/07/31/t-sne-tips/
很好的tsne可视化：https://distill.pub/2016/misread-tsne/c

结尾

Seurat两个版本的结果的确存在不同，但不至于差异太大。它们的tSNE聚类结果“看似差异大”，其实是我们误认为tSNE图中的点之间距离代表了相似性。就像上面V3图中的两群分开的红点，它们其实还是一群，只不过此时此刻的映射坐标是这样的，而V2的映射坐标正好把它们映射到了相近的位置而已。也许多运行几次，使用不同的perplexity参数，V3的结果又会发生变化
数据胜于图片，使用table比较是相对于肉眼观察图片更有效的方式

“

Seurat标准流程

前言