查看原文
其他

Seurat标准流程

刘小泽 单细胞天地 2022-06-06

课程笔记




粉丝:有单细胞线上课程吗?

小编:什么? 我们的单细胞转录组分析线上课程已经上线好久了,你们竟然都不知道吗,每篇推文后面的课程推荐没人看的吗,小编已哭晕在厕所

好了,戏演完了,下面郑重介绍下我们的单细胞线上课程:(详情戳下方链接) 

全网第二个单细胞视频课程预售


这个课程笔记栏目记录了学员们学习单细胞转录组课程的学习笔记

希望大家能有所收获!



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


课程链接在: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[117712 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这里需要是矩阵或数据框;

  • 第二个位置12选一个,原理和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拆分为FindNeighborsFindClusters,而版本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版本在最后的可视化上可以使用DimPlotTSNEPlot,颜色参数变成了cols

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

    3.0的结果

    【注意】如果参数使用不正确,DimPlot或TSNEPlot会调用默认颜色设置,例如使用参数colorscolors.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的R语言环境
10X scRNA免疫治疗学习笔记

学习scRNAseq这个R包

利用scRNAseq包学习scater

用Scater包分析文章数据

用Seurat包分析文章数据(二)

scRNA包学习Monocle2

使用monocle2分析文章数据

使用作者代码重复结果

基因在任意癌症表达量相关性

评估任意基因集在癌症的表现

多个基因集相关性热图


如果你对单细胞转录组研究感兴趣,但又不知道如何入门,也许你可以关注一下下面的课程

生信技能树(爆款入门培训课)全国巡讲约你

(福州、上海见!)全国巡讲第19-20站(生信入门课加量不加价)



您可能也对以下帖子感兴趣

文章有问题?点此查看未经处理的缓存