rna_seq单细胞测序单细胞测序

单细胞测序分析之scater包学习笔记

2019-11-23  本文已影响0人  生信start_site

这次学习的是单细胞测序的第三个比较常用的包:Scater。本文主要参考教程:单细胞转录组学习笔记-15-利用scRNAseq包学习scater来走分析流程。

scater包的官方网站:here
这篇文章写的些许混乱,因为官网里的很多函数与教程里的并不一样,是因为有些函数是最新版scater新更新的,所以我在学习的时候有的部分按照官网的来,有的部分则按照教程里的来单细胞转录组学习笔记-15-利用scRNAseq包学习scater。但是实际上大同小异,主要的流程都是一样的。

这个包像之前介绍的两个包一样,是为了分析单细胞测序分析而开发的,也可以进行质控、可视化、下游分析前的预处理。scater包的流程是这样的:

简单的来说,scater包可以做下面几件事:
(1)QC
(2)从read数据里进行转录本定量(根据伪时间排序)
(3)数据格式标准化
(4)可视化
(5)Seamless integration into the Bioconductor universe
(6)标准化数据

安装scater包

>BiocManager::install("scater")

创建SingleCellExperiment对象

和之前介绍的Monocle2和Seurat包一样,scater包也需要一个对象。
SingleCellExperiment这个对象,画出来大概是这么个意思(图片来自:单细胞转录组学习笔记-15-利用scRNAseq包学习scater):

在导入数据之前,请将表达数据存为矩阵。如果你的数据是10×Genomics的测序结果,可以用[DropletUtils]包的read10xCounts函数,它可以自动的生成SingleCellExperiment对象。这里我们用scater包自带的示例数据来走一遍流程。
首先,你要把数据导入进来:
第一个要导入的是基因的表达量矩阵,行是基因,列是细胞:

> rm(list = ls()) 
> options(stringsAsFactors = F) 
> data("sc_example_counts") #导入表达矩阵
> dim(sc_example_counts) #看一下表达矩阵里多少行多少列
[1] 2000   40
> sc_example_counts[1:3,1:3] # 简单的看一下表达矩阵长什么样
          Cell_001 Cell_002 Cell_003
Gene_0001        0      123        2
Gene_0002      575       65        3
Gene_0003        0        0        0

第二,导入细胞信息,注释信息可以是细胞名称,组织来源,收集样品的时间点,处理条件等等, 行是细胞,列是属性:

> data("sc_example_cell_info")
> sc_example_cell_info[1:3,1:3] # 样本注释信息
             Cell Mutation_Status Cell_Cycle
Cell_001 Cell_001        positive          S
Cell_002 Cell_002        positive         G0
Cell_003 Cell_003        negative         G1

有了上述的文件,就可以构建对象了:

> example_sce <- SingleCellExperiment(
   assays = list(counts = sc_example_counts), 
   colData = sc_example_cell_info)
#瞅一眼对象里有啥
> example_sce
class: SingleCellExperiment 
dim: 2000 40 
metadata(0):
assays(1): counts
rownames(2000): Gene_0001 Gene_0002 ... Gene_1999 Gene_2000
rowData names(0):
colnames(40): Cell_001 Cell_002 ... Cell_039 Cell_040
colData names(4): Cell Mutation_Status Cell_Cycle Treatment
reducedDimNames(0):
spikeNames(0):

上面在构建对象的时候,我们将表达矩阵起了一个名字是“counts”,是为了方便有的时候你如果需要把矩阵里的counts提取出来,用counts函数比较方便:

> str(counts(example_sce))
 int [1:2000, 1:40] 0 575 0 0 0 0 0 0 416 12 ...
 - attr(*, "dimnames")=List of 2
  ..$ : chr [1:2000] "Gene_0001" "Gene_0002" "Gene_0003" "Gene_0004" ...
  ..$ : chr [1:40] "Cell_001" "Cell_002" "Cell_003" "Cell_004" ...

你可以通过下面的方法来新增一列meatdata:

> example_sce$whee <- sample(LETTERS, ncol(example_sce), replace=TRUE)
> colData(example_sce)
DataFrame with 40 rows and 5 columns
                Cell Mutation_Status  Cell_Cycle   Treatment        whee
         <character>     <character> <character> <character> <character>
Cell_001    Cell_001        positive           S      treat1           Z
Cell_002    Cell_002        positive          G0      treat1           I
Cell_003    Cell_003        negative          G1      treat1           W
Cell_004    Cell_004        negative           S      treat1           O
Cell_005    Cell_005        negative          G1      treat2           I
...              ...             ...         ...         ...         ...
Cell_036    Cell_036        negative          G0      treat1           C
Cell_037    Cell_037        negative          G0      treat1           E
Cell_038    Cell_038        negative          G0      treat2           D
Cell_039    Cell_039        negative          G1      treat1           F
Cell_040    Cell_040        negative          G0      treat2           F

新增一行的方法:

#比如我想给每一行加上一个随机数,随机数生成的函数是runif()
> rowData(example_sce)$stuff <- runif(nrow(example_sce))
> rowData(example_sce)
DataFrame with 2000 rows and 1 column
                       stuff
                   <numeric>
Gene_0001 0.0491068968549371
Gene_0002  0.611523448256776
Gene_0003  0.282845608890057
Gene_0004 0.0250186347402632
Gene_0005  0.182038959115744
...                      ...
Gene_1996  0.633494263514876
Gene_1997  0.813403354492038
Gene_1998  0.987525492208079
Gene_1999  0.204807848902419
Gene_2000  0.369947757339105

还有一些其他的函数,以后如果有需要可以用的上,比如(摘自:单细胞转录组学习笔记-15-利用scRNAseq包学习scater):
(1)如果你的数据里有spike-in control,可以用isSpike 对spike-in操作
(2)sizeFactors 是进行标准化时对细胞文库大小计算的结果
(3)reducedDim 对降维结果(reduced dimensionality results)操作

其他导入数据的方法

如果你的count矩阵是存在一个CSV文件里的,可以用read.table() 来读取。
如果是大的数据集,可以用"Matrix"包里的readSparseCounts()函数读取。
或者,如果你的数据来自10X Genomics测序,你可以用 read10xCounts函数读取,它会自动产生一个sparse矩阵。
另外,kallisto和Salmon的转录本丰度数据可以使用tximeta包来进行读取。

质量控制

scater可以从三个水平上进行质量控制(QC):
(1)过滤细胞
(2)过滤基因
(3)实验变量的质控

细胞水平QC

细胞水平指标的计算是利用perCellQCMetrics()函数进行的,它包括:
(1)sum: 细胞里的总count数(文库大小)
(2)detected: 细胞里含有count的基因数
(3)subsets_X_percent: count的百分比(某个基因集)

NOTE这里需要安装最新版的scater(1.14.4),如果是老版的scater,这一步是不能完成的
由于我的破电脑用了N种方法也没安装上最新版的scater,所以质控这些步骤我只把代码放了上来,自己并没有亲自运行一遍,如果有同学可以安装,最好要试一遍质控的过程,因为像下一步perCellQCMetrics函数是最新版里新增的功能:

library(scater)
per.cell <- perCellQCMetrics(example_sce, 
    subsets=list(Mito=grep("mt-", rownames(example_sce))))
summary(per.cell$sum)
summary(per.cell$detected)
summary(per.cell$subsets_Mito_percent)
colData(example_sce) <- cbind(colData(example_sce), per.cell)

教程里并没有用这个函数,而是用了另外一种:calculateQCMetrics
这个函数对细胞和基因都计算了很多种的指标,分别存在colData和rowData里,默认情况下对count值进行计算,你也可以通过参数exprs_values进行修改:

> example_sce <- calculateQCMetrics(example_sce)
> colnames(colData(example_sce)) #列是样品
 [1] "Cell"                           "Mutation_Status"                "Cell_Cycle"                    
 [4] "Treatment"                      "whee"                           "is_cell_control"               
 [7] "total_features_by_counts"       "log10_total_features_by_counts" "total_counts"                  
[10] "log10_total_counts"             "pct_counts_in_top_50_features"  "pct_counts_in_top_100_features"
[13] "pct_counts_in_top_200_features" "pct_counts_in_top_500_features"
> colnames(rowData(example_sce)) #行是基因
[1] "stuff"                 "is_feature_control"    "mean_counts" #平均表达量          "log10_mean_counts"    
[5] "n_cells_by_counts"#多少个细胞表达了该基因     "pct_dropout_by_counts"#基因在细胞中表达量为0的细胞数占比 "total_counts"          "log10_total_counts" 

如果你的测序里添加了ERCC,那么这里还可以设置对照:

> example_sce <- calculateQCMetrics(example_sce, 
                                   feature_controls = list(ERCC = 1:20, mito = 500:1000),
                                   cell_controls = list(empty = 1:5, damaged = 31:40))
> 
> all_col_qc <- colnames(colData(example_sce))
> all_col_qc <- all_col_qc[grep("ERCC", all_col_qc)]
> all_col_qc #这里就是ERCC相关的信息了
[1] "total_features_by_counts_ERCC"       "log10_total_features_by_counts_ERCC"
[3] "total_counts_ERCC"                   "log10_total_counts_ERCC"            
[5] "pct_counts_ERCC"                     "pct_counts_in_top_50_features_ERCC" 
[7] "pct_counts_in_top_100_features_ERCC" "pct_counts_in_top_200_features_ERCC"
[9] "pct_counts_in_top_500_features_ERCC"

可视化

(一)细胞水平上的质控可以画一个总基因表达量 比上feature control的百分比,如果一个基因在总基因表达量上的比例多,在ERCC里少,就是正常的细胞,反之则不正常。

> plotColData(example_sce, x = "total_features_by_counts",
             y = "pct_counts_feature_control", colour = "Mutation_Status") +
             theme(legend.position = "top") +
             stat_smooth(method = "lm", se = FALSE, size = 2, fullrange = TRUE)

(二)plotScater(单细胞转录组学习笔记-15-利用scRNAseq包学习scater
先在表达量最高的基因中选一部分(默认是500个),然后从高到低累加,看看它们对每个细胞文库的贡献值如何。它将不同细胞的不同基因表达分布绘制出来,就像芯片数据或bulk转录组数据每个样本做一个箱线图,来看基因表达量分布。但是由于单细胞数据样本太多,没办法全画出箱线图。这个图在处理不同批次细胞时就可以很清楚地看到它们之间的差异。

为了看不同细胞的表达量差异,可以利用colData中的数据将细胞分类;默认使用counts值,但是只要在assays存在的表达矩阵,就可以在exprs_values参数中自定义:

> plotScater(example_sce, block1 = "Mutation_Status", block2 = "Treatment",
            colour_by = "Cell_Cycle", nfeatures = 300, exprs_values = "counts")
image.png

(三)QC结果两两比较(单细胞转录组学习笔记-15-利用scRNAseq包学习scater

# 比较feature的属性
> plotRowData(example_sce, x = "n_cells_by_counts", y = "mean_counts")
 # 比较样本的属性
> p1 <- plotColData(example_sce, x = "total_counts", 
                   y = "total_features_by_counts")
> p2 <- plotColData(example_sce, x = "pct_counts_feature_control",
                   y = "total_features_by_counts")
> p3 <- plotColData(example_sce, x = "pct_counts_feature_control",
                   y = "pct_counts_in_top_50_features")
 # 最后组合
> multiplot(p1, p2, p3, cols = 3)

过滤低质量细胞

举个栗子,我们只保留有10万个counts和500个基因表达的细胞:

> keep.total <- example_sce$total_counts > 1e5
> keep.n <- example_sce$total_features_by_counts > 500
> filtered <- example_sce[,keep.total & keep.n]
> dim(filtered)
[1] 1973   37

isOutlier 函数提供了一个更为合适的方法来选择阈值。它定义距离中位数一定数量的中位数绝对偏差(MADs)的阈值。超过这个值的细胞认为是离群的,从而被过滤掉:

> keep.total <- isOutlier(example_sce$total_counts, nmads=3, 
                         type="lower", log=TRUE)
> filtered <- example_sce[,keep.total]

使用quickPerCellQC()函数还可以更方便地检测几个常见指标的异常值。该方法使用总count数、检测到的基因数和某个基因集的count百分比(如线粒体基因、spik -in转录本)来确定要舍弃哪些细胞:

qc.stats <- quickPerCellQC(per.cell, percent_subsets="subsets_Mito_percent")
colSums(as.matrix(qc.stats))
filtered <- example_sce[,!qc.stats$discard]

isOutlier的方法需要根据你的实验需求进行调整,如测序深度、spike-in的加入、细胞类型。

你也可以对基因进行过滤:

> keep_feature <- nexprs(example_sce, byrow=TRUE) >= 4
> example_sce <- example_sce[keep_feature,]
> dim(example_sce)
[1] 1753   40

基因水平的质控

指标的定义

使用的perFeatureQCMetrics()函数包括下面几个指标:
(1) mean: 所有的细胞之间基因的平均count
(2)detected: 对于每一个基因都是非零的细胞百分比
(3)subsets_Y_ratio: 某一个子集的细胞和总细胞的平均count的比值

# 这里举例:前10个孔是空孔
per.feat <- perFeatureQCMetrics(example_sce, subsets=list(Empty=1:10))
summary(per.feat$mean)
summary(per.feat$detected)
summary(per.feat$subsets_Empty_ratio)

calculateAverage()函数提供了更精确的平均值计算,它根据取平均值之前的相对文库大小来调整计数:

#这里我用的是过滤前的对象来进行演示的
> ave <- calculateAverage(example_sce)
> summary(ave)
     Min.   1st Qu.    Median      Mean   3rd Qu.      Max. 
    0.009     5.529    67.730   197.135   190.286 12428.798 

我们还可以直接计算出表达基因的细胞数量:

> summary(nexprs(example_sce, byrow=TRUE))
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
   1.00    6.00   12.00   14.66   21.00   40.00 

可视化:
在默认值下,只显示前50个基因。行是基因,橙色的线(bar)代表该基因在每一个细胞里的表达量。圆圈代表这个基因表达量的中位数。

> plotHighestExprs(example_sce, exprs_values = "counts")

下面这个图是对于每一个基因都是非零的细胞百分比:

> plotExprsFreqVsMean(example_sce)

根据行来取子集

从SingleCellExperiment对象里移除一些基因,比如我们把在所有细胞中都不表达的基因删除:

> keep_feature <- rowSums(counts(example_sce) > 0) > 0
> example_sce <- example_sce[keep_feature,]
> dim(example_sce)
[1] 1973   40

也可以用现有的基因注释进行过滤,比如线粒体基因,通常我们对这些基因不感兴趣,对我们了解细胞异质性也没有多大帮助。

变量水平的质控

用getVarianceExplained()函数(经过标准化之后,请参阅下面的内容)。

example_sce <- logNormCounts(example_sce) # see below.
vars <- getVarianceExplained(example_sce, 
    variables=c("tissue", "total mRNA mol", "sex", "age"))
head(vars)
plotExplanatoryVariables(vars)

计算表达值

标准化文库大小差异

第一种方法:
最常用的函数是logNormCounts(),结果储存在"logcounts"里。它是用每一个count除以自身的size factor,加一个 pseudo-count 和log-transforming。

example_sce <- logNormCounts(example_sce)
assayNames(example_sce)
summary(librarySizeFactors(example_sce))

第二种方法:
用calculateCPM()函数,结果储存在"cpm"中。相关函数还包括:calculateTPM()和calculateFPKM()。

cpm(example_sce) <- calculateCPM(example_sce)
assay(example_sce, "normed") <- normalizeCounts(example_sce, 
    size_factors=runif(ncol(example_sce)), pseudo_count=1.5)

上面两种方法是官网的代码,我用的是下面这个,是教程里用的代码:

> cpm(example_sce) <- calculateCPM(example_sce)
> example_sce <- normalize(example_sce)
> assayNames(example_sce)
[1] "counts"    "cpm"       "logcounts"
> plotExplanatoryVariables(example_sce)

也可以选择一部分进行画图:

> plotExplanatoryVariables(example_sce,
                          variables = c("total_features_by_counts", "total_counts",
                                        "Mutation_Status", "Treatment", "Cell_Cycle"))
结果看到:total_counts 和 total_features_by_counts 对基因表达变化差异贡献很高(接近10%),真实数据中这两个占比应该在1-5%之间

可视化表达值

plotExpression()函数可以画某一个基因子集的表达值,默认使用"logcounts"里的数据,也可以通过添加exprs_values参数来修改。

> plotExpression(example_sce, rownames(example_sce)[1:6])
这是前6个基因的logcount值的图

你也可以按照分组的情况画基因表达情况,比如在这个包自带的数据info里,有Mutation_Status,Cell_cycle的分组情况,你也可以根据这些来画图:

> plotExpression(example_sce, rownames(example_sce)[1:6],
                x = "Mutation_Status", exprs_values = "logcounts")
> plotExpression(example_sce, rownames(example_sce)[1:6],
                x = "Cell_Cycle", exprs_values = "logcounts")

你还可以同时将两个变量一起画出来:

> plotExpression(example_sce, rownames(example_sce)[1:6],
                colour_by = "Cell_Cycle", shape_by = "Mutation_Status", 
                size_by = "Gene_0002")

降维

(一)PCA

最常用的是PCA降维的方法,PCA结果将存放在SingleCellExperiment对象的reducedDims里:

> example_sce <- runPCA(example_sce)
> str(reducedDim(example_sce, "PCA"))
 num [1:40, 1:2] -11.35 -1.41 2.01 -9.8 6.35 ...
 - attr(*, "dimnames")=List of 2
  ..$ : chr [1:40] "Cell_001" "Cell_002" "Cell_003" "Cell_004" ...
  ..$ : chr [1:2] "PC1" "PC2"
 - attr(*, "percentVar")= num [1:2] 0.0679 0.0422

默认情况下,runPCA()函数用第一个PCs(聚类)里变化量最大的500个基因来进行计算。你也可以用subset_row参数来调整你想要进行计算的基因数,或者基因子集。用name参数可以修改在reducedDims slot里的名称:

> example_sce <- runPCA(example_sce, name="PCA2",
    subset_row=rownames(example_sce)[1:1000],
    ncomponents=25)

(二)其他降维方法(tSNE, UMAP)

> set.seed(1000)
> example_sce <- runTSNE(example_sce, perplexity=10)
> head(reducedDim(example_sce, "TSNE"))
               [,1]       [,2]
Cell_001  3.7408876 -0.8718345
Cell_002  0.9961570 -1.0279169
Cell_003  0.3328361  1.1045595
Cell_004  4.9745275  0.2320463
Cell_005 -4.1067009  0.1909082
Cell_006  0.5936069 -2.1781319

这里官方网站推荐多用几个不同的随机种子和perplexity值(表示近邻细胞的数量)来跑tSNE,以确保你的结果是比较稳定的。

一个更常见的模式是用你之前跑的PCA结果,作为tSNE的输入,进行分析。这有助于提高运行速度,并且降低噪音,只集中在几个主要的变量上。下面这个代码是用的PCA前10个dimension的结果来进行tSNE计算:

> set.seed(1000)
> example_sce <- runTSNE(example_sce, perplexity=50, 
                       dimred="PCA", n_dimred=10)
> head(reducedDim(example_sce, "TSNE"))

你也可以用UMAP来计算:

> example_sce <- runUMAP(example_sce)
> head(reducedDim(example_sce, "UMAP"))
               [,1]       [,2]
Cell_001  2.1841452 -0.4125833
Cell_002  0.2940562  1.0461869
Cell_003 -0.5211948  1.2969395
Cell_004  2.1369329  0.1383032
Cell_005 -2.0110685 -0.4748272
Cell_006  1.0558627  0.7567186

你还可以用diffusion maps来降维:

> example_sce <- runDiffusionMap(example_sce)
Registered S3 method overwritten by 'xts':
  method     from
  as.zoo.xts zoo 
> plotDiffusionMap(example_sce, colour_by = "Gene_0001", size_by = "Gene_1000")

降维的可视化

#根据Mutation_Status来画
> plotReducedDim(example_sce, use_dimred = "PCA", colour_by = "Mutation_Status")

还可以利用表达量添加颜色、大小分组(参考单细胞转录组学习笔记-15-利用scRNAseq包学习scater):

> # 看特定基因在PCA过程中起到的作用
> plotReducedDim(example_sce, use_dimred = "PCA", 
                colour_by = "Gene_1000", size_by = "Gene_0500")

如果你想同时画2个以上成分的可视化图的时候,可以用下面的代码:

> example_sce <- runPCA(example_sce, ncomponents=20)
> plotPCA(example_sce, ncomponents = 4, colour_by = "Mutation_Status")
图里的对角线上的方块图,显示的是每一个component的密度
上一篇下一篇

猜你喜欢

热点阅读