本章将介绍使用Seurat（V3）的一些典型任务。尽管Seurat目前已经更新至V5版本，但仍不妨碍我们从此教程中学习一些基本操作及分析思想。
现在让我们加载本教程所需的所有库。

> library(Seurat)
> library(ggplot2)
> library(SingleR)
> library(dplyr)
> library(celldex)
> library(RColorBrewer)
> library(SingleCellExperiment)

7.1 基本质控和过滤
在开始分析之前要进行两步操作：1）使用soupX进行环境RNA校正；2）使用scrublet进行双胞检测。但本教程不对soupX和scrublet的使用做介绍。
首先读入SoupX校正矩阵。SoupX的输出仅包含可用的基因名称，所以不需要额外指定。如果从典型的Cell Ranger输出开始，则可以选择在表达矩阵中使用Ensemble ID或基因名称。这是通过gene.column选项实现的；默认值为“2”，即基因名称。

> adj.matrix <- Read10X("data/update/soupX_pbmc10k_filt")

在这之后，我们将创建一个Seurat对象。Seurat对象向我们展示了：1）细胞数量（samples）大致与每个数据集的描述（10,194）相匹配；2）有36,601个基因（features）。

> srat <- CreateSeuratObject(adj.matrix,project = "pbmc10k")
> srat
An object of class Seurat 
36601 features across 10194 samples within 1 assay 
Active assay: RNA (36601 features, 0 variable features)

让我们从内存中删除adj.matrix以节省RAM，并仔细查看Seurat对象。str命令允许我们查看该类的所有字段：

> adj.matrix <- NULL
> str(srat)
Formal class 'Seurat' [package "SeuratObject"] with 13 slots
  ..@ assays      :List of 1
  .. ..$ RNA:Formal class 'Assay' [package "SeuratObject"] with 8 slots
  .. .. .. ..@ counts       :Formal class 'dgCMatrix' [package "Matrix"] with 6 slots
  .. .. .. .. .. ..@ i       : int [1:24330253] 25 30 32 42 43 44 51 59 60 62 ...
  .. .. .. .. .. ..@ p       : int [1:10195] 0 4803 7036 11360 11703 15846 18178 20413 22584 27802 ...
  .. .. .. .. .. ..@ Dim     : int [1:2] 36601 10194
  .. .. .. .. .. ..@ Dimnames:List of 2
  .. .. .. .. .. .. ..$ : chr [1:36601] "MIR1302-2HG" "FAM138A" "OR4F5" "AL627309.1" ...
  .. .. .. .. .. .. ..$ : chr [1:10194] "AAACCCACATAACTCG-1" "AAACCCACATGTAACC-1" "AAACCCAGTGAGTCAG-1" "AAACCCAGTGCTTATG-1" ...
  .. .. .. .. .. ..@ x       : num [1:24330253] 1 2 1 1 1 3 1 1 1 1 ...
  .. .. .. .. .. ..@ factors : list()
  .. .. .. ..@ data         :Formal class 'dgCMatrix' [package "Matrix"] with 6 slots
  .. .. .. .. .. ..@ i       : int [1:24330253] 25 30 32 42 43 44 51 59 60 62 ...
  .. .. .. .. .. ..@ p       : int [1:10195] 0 4803 7036 11360 11703 15846 18178 20413 22584 27802 ...
  .. .. .. .. .. ..@ Dim     : int [1:2] 36601 10194
  .. .. .. .. .. ..@ Dimnames:List of 2
  .. .. .. .. .. .. ..$ : chr [1:36601] "MIR1302-2HG" "FAM138A" "OR4F5" "AL627309.1" ...
  .. .. .. .. .. .. ..$ : chr [1:10194] "AAACCCACATAACTCG-1" "AAACCCACATGTAACC-1" "AAACCCAGTGAGTCAG-1" "AAACCCAGTGCTTATG-1" ...
  .. .. .. .. .. ..@ x       : num [1:24330253] 1 2 1 1 1 3 1 1 1 1 ...
  .. .. .. .. .. ..@ factors : list()
  .. .. .. ..@ scale.data   : num[0 , 0 ] 
  .. .. .. ..@ key          : chr "rna_"
  .. .. .. ..@ assay.orig   : NULL
  .. .. .. ..@ var.features : logi(0) 
  .. .. .. ..@ meta.features:'data.frame':      36601 obs. of  0 variables
  .. .. .. ..@ misc         : list()
  ..@ meta.data   :'data.frame':        10194 obs. of  3 variables:
  .. ..$ orig.ident  : Factor w/ 1 level "pbmc10k": 1 1 1 1 1 1 1 1 1 1 ...
  .. ..$ nCount_RNA  : num [1:10194] 22196 7630 21358 857 15007 ...
  .. ..$ nFeature_RNA: int [1:10194] 4734 2191 4246 342 4075 2285 2167 2151 5134 3037 ...
  ..@ active.assay: chr "RNA"
  ..@ active.ident: Factor w/ 1 level "pbmc10k": 1 1 1 1 1 1 1 1 1 1 ...
  .. ..- attr(*, "names")= chr [1:10194] "AAACCCACATAACTCG-1" "AAACCCACATGTAACC-1" "AAACCCAGTGAGTCAG-1" "AAACCCAGTGCTTATG-1" ...
  ..@ graphs      : list()
  ..@ neighbors   : list()
  ..@ reductions  : list()
  ..@ images      : list()
  ..@ project.name: chr "pbmc10k"
  ..@ misc        : list()
  ..@ version     :Classes 'package_version', 'numeric_version'  hidden list of 1
  .. ..$ : int [1:3] 4 1 0
  ..@ commands    : list()
  ..@ tools       : list()

meta.data是接下来步骤中最重要的字段。可以使用@和[[]]运算符访问它。
目前每个细胞有3个字段：数据集ID、每个细胞检测到的UMI数（nCount_RNA）以及每个相同细胞表达（检测到）的基因数（nFeature_RNA）。

> meta <- srat@meta.data
> dim(meta)
[1] 10194     3
> head(meta)
                   orig.ident nCount_RNA nFeature_RNA
AAACCCACATAACTCG-1    pbmc10k      22196         4734
AAACCCACATGTAACC-1    pbmc10k       7630         2191
AAACCCAGTGAGTCAG-1    pbmc10k      21358         4246
AAACCCAGTGCTTATG-1    pbmc10k        857          342
AAACGAACAGTCAGTT-1    pbmc10k      15007         4075
AAACGAACATTCGGGC-1    pbmc10k       9855         2285
> summary(meta$nCount_RNA)
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
    499    5549    7574    8902   10730   90732 
> summary(meta$nFeature_RNA)
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
     47    1725    2113    2348    2948    7265

我们来添加更多可用于评估细胞质量的值。线粒体基因是细胞状态的有用指标。对于小鼠数据集，将匹配模式更改为“Mt-”，或使用features = ...选项明确列出基因ID。

> srat[["percent.mt"]] <- PercentageFeatureSet(srat, pattern = "^MT-") 

类似地，我们可以定义核糖体蛋白（名称以RPS或RPL开头），它们通常占据reads的大部分：

> srat[["percent.rb"]] <- PercentageFeatureSet(srat, pattern = "^RP[SL]") 

现在，让我们将scrublet生成的双胞注释添加到Seurat对象元数据中。

> doublets <- read.table("data/update/scrublet_calls.tsv",header = F,row.names = 1)
> colnames(doublets) <- c("Doublet_score","Is_doublet")
> srat <- AddMetaData(srat,doublets)
> head(srat[[]])
                   orig.ident nCount_RNA nFeature_RNA percent.mt percent.rb
AAACCCACATAACTCG-1    pbmc10k      22196         4734   5.275725  25.067580
AAACCCACATGTAACC-1    pbmc10k       7630         2191   8.833552  33.853211
AAACCCAGTGAGTCAG-1    pbmc10k      21358         4246   6.283360  19.276149
AAACCCAGTGCTTATG-1    pbmc10k        857          342  31.388565   1.750292
AAACGAACAGTCAGTT-1    pbmc10k      15007         4075   7.916306  14.986340
AAACGAACATTCGGGC-1    pbmc10k       9855         2285   7.762557  41.024860
                   Doublet_score Is_doublet
AAACCCACATAACTCG-1    0.30542385       True
AAACCCACATGTAACC-1    0.01976120      False
AAACCCAGTGAGTCAG-1    0.03139876      False
AAACCCAGTGCTTATG-1    0.02755040      False
AAACGAACAGTCAGTT-1    0.36867997       True
AAACGAACATTCGGGC-1    0.09723644      False

绘制所选元数据特征的小提琴图。请注意，这些图按名为identity的类分组。Identity类可以在srat@active.ident中看到，或者使用Idents()函数，SetIdents()可以更改处于活动状态的identity。

> VlnPlot(srat, features = c("nFeature_RNA","nCount_RNA","percent.mt","percent.rb"), ncol = 4,pt.size = 0.1) & 
    theme(plot.title = element_text(size=10))

再来绘制一些元数据特征，看看它们之间的相关性。每个图上方的数字是皮尔逊相关系数。

> FeatureScatter(srat, feature1 = "nCount_RNA", feature2 = "percent.mt")

> FeatureScatter(srat, feature1 = "nCount_RNA", feature2 = "nFeature_RNA")

> FeatureScatter(srat, feature1 = "nCount_RNA", feature2 = "percent.rb")

> FeatureScatter(srat, feature1 = "percent.rb", feature2 = "percent.mt")

> FeatureScatter(srat, feature1 = "nFeature_RNA", feature2 = "Doublet_score")

上图清楚地表明，高线粒体百分比与低UMI密切相关，这通常被解释为死细胞。但高核糖体蛋白含量与线粒体强烈负相关，并且似乎包含生物信号。双胞得分和表达的基因数之间也存在很强的相关性。让我们在元数据中增加QC列，并充分的定义它。

> srat[['QC']] <- ifelse(srat@meta.data$Is_doublet == 'True','Doublet','Pass')
> srat[['QC']] <- ifelse(srat@meta.data$nFeature_RNA < 500 & srat@meta.data$QC == 'Pass','Low_nFeature',srat@meta.data$QC)
> srat[['QC']] <- ifelse(srat@meta.data$nFeature_RNA < 500 & srat@meta.data$QC != 'Pass' & srat@meta.data$QC != 'Low_nFeature',paste('Low_nFeature',srat@meta.data$QC,sep = ','),srat@meta.data$QC)
> srat[['QC']] <- ifelse(srat@meta.data$percent.mt > 15 & srat@meta.data$QC == 'Pass','High_MT',srat@meta.data$QC)
> srat[['QC']] <- ifelse(srat@meta.data$nFeature_RNA < 500 & srat@meta.data$QC != 'Pass' & srat@meta.data$QC != 'High_MT',paste('High_MT',srat@meta.data$QC,sep = ','),srat@meta.data$QC)
> table(srat[['QC']])
QC
                     Doublet                      High_MT 
                         546                          548 
        High_MT,Low_nFeature High_MT,Low_nFeature,Doublet 
                         275                            1 
                        Pass 
                        8824

我们可以看到，双胞通常不会与基因数较少的细胞重叠；同时，后者通常与线粒体含量高的细胞重叠。让我们仅绘制通QC列为pass的细胞的元数据：

> VlnPlot(subset(srat, subset = QC == 'Pass'), 
          features = c("nFeature_RNA", "nCount_RNA", "percent.mt","percent.rb"), ncol = 4, pt.size = 0.1) & 
    theme(plot.title = element_text(size=10))

往期内容：
重生之我在剑桥大学学习单细胞RNA-seq分析——6. 生物学分析（1）
重生之我在剑桥大学学习单细胞RNA-seq分析——6. 生物学分析（2）
重生之我在剑桥大学学习单细胞RNA-seq分析——6. 生物学分析（3）
重生之我在剑桥大学学习单细胞RNA-seq分析——6. 生物学分析（4）