学习笔记，仅供参考，有错必纠
参考资料：https://satijalab.org/seurat/index.html；
备注：这是我第一次学习的笔记，对此了解更深入之后，会重新整理
文章目录
- Seurat Package
- 加载包
- 读取数据
- 数据预处理
- 质量控制，选择需要分析的细胞
- 可视化
- 标准化数据
- 鉴定高度变化的基因
- Scaling the data
- 线性降维
- 决定数据集维度
- JackStrawPlot
- ElbowPlot
Seurat Package
在本教程中，我们将分析10X Genomics公司免费提供的外周血单核细胞(PBMC)数据集，该数据集包含2,700个单细胞在Illumina NextSeq 500上的测序信息. 原始数据指路:Peripheral Blood Mononuclear Cells;
加载包
加载Seurat包及其他需要使用的包.
library(Seurat)
library(dplyr)
library(patchwork)查看Seurat包的版本信息.
packageVersion("Seurat")[1] ‘4.1.1’读取数据
用Read10X()方法读取文件，该方法可以便捷的获取由10X genomics提供的稀疏数据矩阵.
Read10X(
data.dir,
gene.column = 2,
cell.column = 1,
unique.features = TRUE,
strip.suffix = FALSE
)- 参数
data.dir: 包含 10X 提供的 matrix.mtx、genes.tsv（或 features.tsv）和barcodes.tsv 文件的目录. 可以用一个向量或命名向量加载多个数据目录. 如果给定了一个命名向量，则cell barcodes的名称将以该名称为前缀.
gene.column: 指定genes.tsv或features.tsv的哪一列用于基因名称，默认为 2.
cell.column: 指定要用于单元格名称的barcodes.tsv 的哪一列，默认为 1.
unique.features: 使特征名称唯一(默认为TRUE)
strip.suffix: 如果所有cell barcodes中都存在尾随-1，请删除.- 返回值
如果features.tsv数据包含了多种数据类型，则将返回一个包含每种类型数据的稀疏矩阵的列表. 否则将返回一个包含表达数据的稀疏矩阵.
现在，我们读取下载好的数据文件.
data_dir <- "./data/filtered_gene_bc_matrices/hg19"
list.files(data_dir)
pbmc_data <- Read10X(data.dir = data_dir)[1] "barcodes.tsv" "genes.tsv"
[3] "Homo_sapiens.GRCh38.106.chr.gtf.gz" "matrix.mtx"读取文件后，我们利用CreateSeuratObject()函数从原始数据中建立Seurat对象.
CreateSeuratObject(
counts,
project = "CreateSeuratObject",
assay = "RNA",
names.field = 1,
names.delim = "_",
meta.data = NULL,
...
)- 参数
counts: 一个类似矩阵的对象，其为非标准化数据(细胞作为列，特征作为行).
project: Seurat对象的项目名称.
assay: 初始试验的名称.
min.cells: 基因过滤指标，包括至少在多个细胞中检测到的特征.
min.feature: 细胞过滤指标.
row.names：当counts是data.frame或data.frame派生对象时，其为可选的特征名称向量.- 返回值
一个Seurat对象.现在我们构造Seurat对象，并将其实例化为pbmc.
pbmc <- CreateSeuratObject(counts = pbmc.data, min.cells = 3,
min.features = 200, project = "10X_PBMC")
# min.cells = 3, 基因至少在3个细胞中表达
# min.features = 200, 一个细胞中表达基因至少有200
# project = "10X_PBMC", 置项目名称为10X_PBMC
pbmcAn object of class Seurat
13714 features across 2700 samples within 1 assay
Active assay: RNA (13714 features, 0 variable features)查看数据.
pbmc.data[c("CD3D", "TCL1A", "MS4A1"), 1:20]3 x 20 sparse Matrix of class "dgCMatrix"
[[ suppressing 20 column names ��AAACATACAACCAC-1��, ��AAACATTGAGCTAC-1��, ��AAACATTGATCAGC-1�� ... ]]
CD3D 4 . 10 . . 1 2 3 1 . . 2 7 1 . . 1 3 . 2
TCL1A . . . . . . . . 1 . . . . . . . . . . .
MS4A1 . 6 . . . . . . 1 1 1 . . . . . . . . .上表中列为细胞，行为基因，这里的点代表0，意思是这些基因在这个细胞里没有被检测到. 因为大部分基因在单细胞RNA-seq矩阵里都是0，所以Seurat尽可能使用稀疏矩阵（sparse-matrix)表示.
数据预处理
质量控制，选择需要分析的细胞
利用Seurat可以轻松探索QC指标并根据任何用户定义的标准过滤细胞. 常用的QC指标包括：
- 每个细胞中检测到的唯一基因数
- 低质量的细胞或 empty droplets 往往会有很少的基因
- doublets 或者 multiplets细胞可能表现出异常的高基因数
- 在一个细胞内检测到的分子数(与唯一基因密切相关)
- 读取到的线粒体基因组比例
- 低质量/濒临死亡的细胞往往有很高的线粒体污染
- 使用
PercentageFeatureSet()函数计算线粒体QC指标 - 设置以
MT-开头的基因作为线粒体基因
PercentageFeatureSet()能够轻松计算可能的特征子集的占所有特征集计数的百分比. 例如，计算映射到线粒体基因的转录本百分比时，这里的计算是属于特征子集的计数槽中存在的矩阵的列和除以所有列和再乘以100.
PercentageFeatureSet(
object,
pattern = NULL,
features = NULL,
col.name = NULL,
assay = NULL
)- 参数
object: Seurat对象
pattern: 用来匹配特性的正则表达式
features: 定义的特征集, 如果提供了特性，将忽略pattern参数
col.name: meta.data列中要分配的名称，如果这不是空的，则返回一个Seurat对象，其特征集的比例存储在metadata中.- 返回值
返回一个包含特征集比例的向量，如果md.name被设置，则返回一个Seurat对象，其特征集的比例存储在元数据中.现在计算映射到线粒体基因的转录本比例.
pbmc[["percent.mt"]] <- PercentageFeatureSet(object = pbmc, pattern = "^MT-")dim(pbmc@meta.data)
head(pbmc@meta.data, 5)[1] 2700 4
可视化
我们可以可视化QC指标，并用它们来过滤细胞.
- 将有唯一特征数的小于200 或 大于2500 的细胞过滤掉
- 将线粒体数占5%以上的细胞过滤
# Visualize QC metrics as a violin plot
VlnPlot(pbmc, features = c("nFeature_RNA", "nCount_RNA", "percent.mt"), ncol = 3)
FeatureScatter()可以创建两个特征之间的散点图. 细胞按其身份类别着色，两个特征之间的皮尔逊相关关系显示在图的上方.
我们用FeatureScatter()函数可视化特征与特征之间的关系.
plot1 <- FeatureScatter(object = pbmc,
feature1 = 'nCount_RNA', feature2 = 'percent.mt')
plot2 <- FeatureScatter(object = pbmc,
feature1 = 'nCount_RNA', feature2 = 'nFeature_RNA')
CombinePlots(plots = list(plot1, plot2))
现在，我们将有唯一特征数的小于200 或 大于2500 的细胞过滤掉.
pbmc <- subset(pbmc, subset = nFeature_RNA > 200 & nFeature_RNA < 2500 & percent.mt < 5)标准化数据
从数据集中去除不需要的细胞后，下一步是对数据进行归一化. 默认情况下，我们采用NormalizeData()方法的全局缩放的归一化 LogNormalize，该方法利用总的表达量对每个细胞里的基因表达值进行标准化，再乘以一个比例因子(默认为10000)，并将结果进行对数转换. 归一化的数值存储在pbmc[["RNA"]]@data中.
pbmc <- NormalizeData(pbmc,
normalization.method = "LogNormalize", # 默认
scale.factor = 10000)鉴定高度变化的基因
接下来我们通过直接对单细胞数据的均值-变异关系进行建模，计算在pbmc里细胞之间变化度很高的基因集(这些基因在一些细胞中高表达，在另一些细胞中低表达)，这一步用FindVariableFeatures()函数来执行.
默认情况下，我们为每个数据集返回2000个特征。这些将被用于下游分析，如PCA.
FindVariableFeatures()可识别mean variability图上的异常特征.
VariableFeatures()获取和设置一个实验对象的可变特征信息. HVFInfo和VariableFeatures一般利用可变特征，而SVFInfo和SpatiallyVariableFeatures则被限制于空间可变特征.
pbmc <- FindVariableFeatures(pbmc, selection.method = "vst", nfeatures = 2000)
# Identify the 10 most highly variable genes
top10 <- head(VariableFeatures(pbmc), 10)dim(pbmc)[1] 13714 2638top10[1] "PPBP" "LYZ" "S100A9" "IGLL5" "GNLY" "FTL" "PF4" "FTH1" "GNG11" "S100A8"以上为10个变异程度最大的特征.
现在我们使用VariableFeaturePlot()函数将变量特征可视化.
# plot variable features with and without labels
plot1 <- VariableFeaturePlot(pbmc)
plot2 <- LabelPoints(plot = plot1, points = top10, repel = TRUE)
plot2
Scaling the data
接下来，我们应用ScaleData()函数进行线性变换(“缩放”)， 这是应用在PCA等降维技术之前的一个标准预处理步骤. 该方法Shift以及Scales每个基因表达，使得细胞间的平均表达为0，变异为1. 结果将存储在pbmc[["RNA"]]@scale.data中.
all.genes <- rownames(pbmc)
pbmc <- ScaleData(pbmc, features = all.genes)线性降维
接下来，我们利用PCA对标准化后的数据降维. 在默认值下，只对前面选出的2000个基因降维，
但是也可以选择想要降维的基因集.
利用RunPCA()函数可进行PCA降维，PCA计算的细节请参见PrintPCAParams.
pbmc <- RunPCA(pbmc, features = VariableFeatures(object = pbmc))Seurat 提供几种方法可以对PCA分析后的细胞和基因进行可视化，即VizDimReduction(), DimPlot()和 DimHeatmap().
用几种不同的方法检查PCA的结果.
# Examine and visualize PCA results a few different ways
print(pbmc[["pca"]], dims = 1:5, nfeatures = 5)PC_ 1
Positive: CST3, TYROBP, LST1, AIF1, FTL
Negative: MALAT1, LTB, IL32, IL7R, CD2
PC_ 2
Positive: CD79A, MS4A1, TCL1A, HLA-DQA1, HLA-DQB1
Negative: NKG7, PRF1, CST7, GZMB, GZMA
PC_ 3
Positive: HLA-DQA1, CD79A, CD79B, HLA-DQB1, HLA-DPB1
Negative: PPBP, PF4, SDPR, SPARC, GNG11
PC_ 4
Positive: HLA-DQA1, CD79B, CD79A, MS4A1, HLA-DQB1
Negative: VIM, IL7R, S100A6, IL32, S100A8
PC_ 5
Positive: GZMB, NKG7, S100A8, FGFBP2, GNLY
Negative: LTB, IL7R, CKB, VIM, MS4A7VizDimLoadings(pbmc, dims = 1:2, reduction = "pca")
DimPlot(pbmc, reduction = "pca")
DimHeatmap()函数用于绘制聚焦于主成分的热图，该函数可以方便地探索数据集中异质性的主要来源，在决定将哪些主成分纳入下游分析时非常有用. 细胞和特征都是根据它们的PCA分数来排序的.
DimHeatmap(
object,
dims = 1,
nfeatures = 30,
cells = NULL,
reduction = "pca",
disp.min = -2.5,
disp.max = NULL,
balanced = TRUE,
projected = FALSE,
ncol = NULL,
fast = TRUE,
raster = TRUE,
slot = "scale.data",
assays = NULL,
combine = TRUE
)- 参数
object: Seurat对象
dims: 绘制的主成分维度
nfeatures: 要绘制的基因数量
cells: 要绘制的细胞列表. 如果是数字，则只绘制top cells.
reduction: 使用哪种降维方法
disp.min: 最小显示值(低于该值的所有值都被剪切)
disp.max: 最大显示值(高于该值的所有值都被剪掉)，如果`slot='scale.data'`，默认为2.5，否则为6
balanced: 用`+`和`-`标记相等数量的基因DimHeatmap(pbmc, dims = 1, cells = 500, balanced = TRUE)
绘制多个主成分的热力图.
DimHeatmap(pbmc, dims = 1:9, cells = 500, balanced = TRUE)
决定数据集维度
本部分我们用2种方法研究到底应该选取几个主成分进行后续分析.
JackStrawPlot
JackStraw()将随机排列一个数据子集，并计算这些 "随机"基因的预测PCA得分. 然后将 "随机 "基因的PCA得分与观察到的PCA得分进行比较，以确定统计意义. 最终结果是每个基因与每个主成分的关联的P值.
JackStraw(
object,
reduction = "pca",
assay = NULL,
dims = 20,
num.replicate = 100,
prop.freq = 0.01,
verbose = TRUE,
maxit = 1000
)- 参数
object: Seurat对象
reduction: DimReduc to use. ONLY PCA CURRENTLY SUPPORTED.
dims: Number of PCs to compute significance for
num.replicate: 要执行的replicate采样数
prop.freq: 在每次replicate中随机排列的数据的比例
verbose: 是否打印进度条
maxit: RunPCA的irlba函数的最大迭代数- 返回值
返回一个Seurat对象，其中JS(object = object[['pca']], slot = 'empirical') 代表PCA分析中每个基因的p值. 如果随后运行ProjectPCA，JS(object = object[['pca']], slot = 'full')则代表所有基因的p值.与null distribution相比，显著的主成分应该表现出强烈向左倾斜的P值分布. The p-value for each PC is based on a proportion test comparing the number of features with a p-value below a particular threshold (score.thresh), compared with the proportion of features expected under a uniform distribution of p-values.
# approximate techniques such as those implemented in ElbowPlot() can be used to reduce
# computation time
pbmc <- JackStraw(pbmc, num.replicate = 100)
pbmc <- ScoreJackStraw(pbmc, dims = 1:20)JackStrawPlot()函数提供了一个可视化工具，用于比较每个主成分的p值分布与均匀分布(虚线). 显著的 主成分将展示出更低的p值(虚线上方的实线曲线). 在这种情况下，似乎在前10-12个主成分之后，显著性会急剧下降.
JackStrawPlot(pbmc, dims = 1:15)
ElbowPlot
另一种方法使用"Elbow plot"，即根据每个主成分所解释的方差百分比对主成分进行排名. 在这个例子中，我们可以观察到主成分9-10周围的"Elbow"，这表明大部分的真实信息是在前10个主成分中获取的.
ElbowPlot(pbmc)
这里选择前10个主成分进行后续分析.
