本笔记主要参考
https://www.jianshu.com/p/75c6d55a63bb?from=timeline
主要是为了理清思路，记录一下，方便以后查找
Seurat软件介绍
导入分析需要的包
library(Seurat)
library(dplyr)
library(magrittr)
读入矩阵
#matrix.mtx, genes.tsv, and barcodes.tsv
#使用Read10X函数读取矩阵数据
pbmc.data <- Read10X(data.dir = "../data/pbmc3k/filtered_feature_bc_matrix/")
#使用CreateSeuratObject函数创建seurat对象，这里要求细胞中基因数目不能小于200，且基因至少在三个细胞中有表达
pbmc <- CreateSeuratObject(counts = pbmc.data, project = "pbmc3k", min.cells = 3, min.features = 200)
pbmc
默认情况下，我们的seurat对象中是一个叫RNA的Assay。
在我们处理数据的过程中，做整合（integration），或者做变换（SCTransform），或者做去除污染（SoupX），或者是融合velocity的数据等，我们可能会生成新的相关的Assay，用于存放这些处理之后的矩阵。
在之后的处理中，我们可以根据情况使用指定Assay下的数据。不指定Assay使用数据的时候， Seurat给我们调用的是Default Assay下的内容。
可以通过对象名@active.assay查看当前Default Assay，通过DefaultAssay函数更改当前Default Assay。
Assay数据中，counts为raw，data为normalized，scale为scaled
调用Assay中的数据的方式为，以调取一个名为PBMC的Seurat对象中Assay integrate中的nomalized数据为例：
PBMC@assays$RNA@data
细胞选择
由于低质量细胞或空液滴通常只有很少的基因、双细胞可能表现出异常高的基因数目、低质量/濒死细胞通常表现出广泛的线粒体污染三个原因，我们需要进行细胞数目的选择。
使用percentagefeatureset函数计算线粒体qc指标，该函数计算源自一组基因的计数百分比。
#使用从MT-作为线粒体基因集
pbmc[["percent.mt"]] <- PercentageFeatureSet(pbmc, pattern = "^MT-")
qc指标可视化
VlnPlot(pbmc, features = c("nFeature_RNA", "nCount_RNA", "percent.mt"), ncol = 3)
从可视化结果进行细胞选择，过滤超过2500或少于200基因的细胞，同时过滤线粒体比例超过5%的细胞。
数据标准化
过滤细胞后，之后进行标准化数据，默认情况下，我们使用“lognormalize”全局缩放的归一化方法，通过总表达值对每个细胞的基因表达值归一化，并将其乘以缩放因子(默认为10,000)，最后对结果进行对数变换，标准化数据存储在pbmc@assays$RNA@data中。
pbmc <- NormalizeData(pbmc, normalization.method = "LogNormalize", scale.factor = 10000)
高变异基因选择
计算数据集中表现出高细胞间差异的基因子集(即，它们在一些细胞中高表达，而在另一些细胞中低表达)
在下游分析中关注这些基因有助于突出单细胞数据集中的生物信号。默认情况下，每个数据集选择2000个高变异基因用于下游分析。
pbmc <- FindVariableFeatures(pbmc, selection.method = "vst", nfeatures = 2000)
#前10的高变基因
top10 <- head(VariableFeatures(pbmc), 10)
#将前10的高变基因在图中标注出来
plot1 <- VariableFeaturePlot(pbmc)
plot2 <- LabelPoints(plot = plot1, points = top10, repel = TRUE)
CombinePlots(plots = list(plot1, plot2))
缩放基因
使用线性变换(“scaling”)处理数据
过程：改变每个基因的表达，使细胞间的平均表达为0；缩放每个基因的表达，使细胞间的差异为1，这样高表达基因就不会影响后续分析；这步是pca等降维技术之前的标准预处理步骤，缩放的数据储存在pbmc@assays$RNA@scale.data中。
all.genes <- rownames(pbmc)
pbmc <- ScaleData(pbmc, features = all.genes)
线性降维
对缩放后的数据进行降维，默认降至50维。
PCA(Principal Component Analysis)，即主成分分析方法，是一种使用最广泛的数据降维算法。PCA的主要思想是将n维特征映射到k维上，在尽可能保留原始数据信息的同时降低数据维度来加速数据分析。
过程就是从原始高维的空间中按顺序地找一组相互正交的坐标轴系统，新的坐标轴的选择与数据本身是密切相关的。
其中，第一个新坐标轴选择是原始数据中方差最大的方向，第二个新坐标轴选取是与第一个坐标轴正交的平面中使得方差最大的，第三个轴是与第1,2个轴正交的平面中方差最大的。依次类推，可以得到n个这样的坐标轴。大部分方差都包含在前面k个坐标轴中，后面的坐标轴所含的方差几乎为0。于是，我们可以忽略余下的坐标轴，只保留前面k个含有绝大部分方差的坐标轴，实现对数据降维。
pbmc <- RunPCA(pbmc, features = VariableFeatures(object = pbmc))
#选择前5个维度进行查看
print(pbmc[["pca"]], dims = 1:5, nfeatures = 5)
使用DimPlot可视化函数查看样品中的细胞在pca中的分布情况：
DimPlot(pbmc, reduction = "pca")
使用DimHeatmap可视化函数查看样品中前500个细胞在前6个PCA中的热图：
DimHeatmap(pbmc, dims = 1:6, cells = 500, balanced = TRUE)
确定数据集的维度
每个pc本质上代表一个“元特征”，它将相关特征集中的信息组合在一起。因此，越在顶部的主成分越可能代表数据集。然而，我们应该选择多少个主成分才认为我们选择的数据包含了绝大部分的原始数据信息呢？我们可以使用JackStraw函数来选择PC数目
在JackStraw()函数中, 使用基于零分布的置换检验方法。随机抽取一部分基因(默认1%)然后进行pca分析得到pca分数，将这部分基因的pca分数与先前计算的pca分数进行比较得到显著性p-Value，。根据主成分（pc）所包含基因的p-value进行判断选择主成分。最终的结果是每个基因与每个主成分的关联的p-Value。保留下来的主成分是那些富集小的p-Value基因的主成分。
处理大数据时会花费大量时间；ElbowPlot()内置了一些其它的方法可以减少运行时间。
#将数据的1%打乱重新运行pca，构建特征得分的“零分布”，这个过程重复100次
pbmc <- JackStraw(pbmc, num.replicate = 100)
pbmc <- ScoreJackStraw(pbmc, dims = 1:20)
JackStrawPlot()函数提供可视化方法，用于比较每一个主成分的p-value的分布，虚线是均匀分布；显著的主成分富集有小p-Value基因，实线位于虚线左上方。下图表明保留10个pca主成分用于后续分析是比较合理的。
JackStrawPlot(pbmc, dims = 1:15)
ElbowPlot
ElbowPlot(pbmc)
启发式评估方法生成一个Elbow plot图。在图中展示了每个主成分对数据方差的解释情况（百分比表示），并进行排序。根据自己需要选择主成分，图中发现第9个主成分是一个拐点，后续的主成分(PC)变化都不大了。
注意：鉴别数据的真实维度不是件容易的事情；除了上面两种方法，Serat官当文档还建议将主成分(数据异质性的相关来源有关)与GSEA分析相结合。Dendritic cell 和 NK aficionados可能识别的基因与主成分 12 和 13相关，定义了罕见的免疫亚群 (i.e. MZB1 is a marker for plasmacytoid DCs)。如果不是事先知道的情况下，很难发现这些问题。
Serat文档因此鼓励用户使用不同数量的PC（10、15，甚至50）重复下游分析。其实也将观察到的，结果通常没有显著差异。因此，在选择此参数时，可以尽量选大一点的维度，维度太小的话对结果会产生不好的影响。
细胞聚类
Seurat v3应用基于图形的聚类方法，例如KNN方法。具有相似基因表达模式的细胞之间绘制边缘，然后将他们划分为一个内联群体。
在PhenoGraph中，首先基于pca维度中（先前计算的pca数据）计算欧式距离（the euclidean distance），然后根据两个细胞在局部的重合情况（Jaccard 相似系数）优化两个细胞之间的边缘权值。此步骤内置于FindNeighbors函数，输入时先前确定的pc数据。
为了聚类细胞，接下来应用模块化优化技术迭代将细胞聚集在一起。（the Louvain algorithm (default) or SLM [SLM, Blondel et al., Journal of Statistical Mechanics]），FindClusters函数实现这一功能，其中需要注意resolution参数，该参数设置下游聚类分析的“granularity”，更大的resolution会导致跟多的细胞类群。3000左右的细胞量，设置resolution为0.4-1.2是比较合适的。细胞数据集越大，需要更大的resolution参数, 会获得更多的细胞聚类。
查看细胞属于那个类群可以使用函数Idents。
#计算最邻近距离
pbmc <- FindNeighbors(pbmc, dims = 1:10)
#聚类，包含设置下游聚类的“间隔尺度”的分辨率参数resolution ，增加值会导致更多的聚类。
pbmc <- FindClusters(pbmc, resolution = 0.5)
#可以使用idents函数找到聚类情况：
#查看前5个细胞的聚类id
head(Idents(pbmc), 5)
非线性降维分析
Seurat提供了一些非线性降维度分析的方法，如tSNE和UMAP，以可视化和探索这些数据集；将上一步PC维度中的数据进一步降至二维，实现数据的可视化。之前pca是线性降维，能很好处理高维度数据，tsne和umap是非线性降维方式，对低维度数据处理速度快，但难以处理高维度数据。
相比tsne，umap能更好地反映原始数据结构，有着更好的连续性。
# install UMAP： reticulate::py_install(packages ='umap-learn')
pbmc <- RunUMAP(pbmc, dims = 1:10)
pbmc <- RunTSNE(pbmc, dims = 1:10)
#画图
DimPlot(pbmc, reduction = "umap")
DimPlot(pbmc, reduction = "tsne")
#添加细胞类群xiba的标签
DimPlot(pbmc, reduction = "umap",label = TRUE)
LabelClusters(DimPlot(pbmc, reduction = "umap"),id = 'ident')
saveRDS(pbmc, file = "../output/pbmc_tutorial.rds")
差异分析
seurat可以通过差异表达找到每个聚类的markgene，差异分析可以有多种形式，如找到所有聚类的markene（如cluster1中所有的markgene是指cluster1相对于其余所有cluster是差异的）、两个cluster之间的差异分析、某个cluster中两个样品之间差异分析等。
seurat可以通过差异表达分析寻找不同细胞类群的标记基因。FindMarkers函数可以进行此操作，但是默认寻找单个类群(参数ident.1)与其他所有类群阳性和阴性标记基因。FindAllMarkers函数会自动寻找每个类群和其他每个类群之间的标记基因。
min.pct参数:设定在两个细胞群中任何一个被检测到的百分比，通过此设定不检测很少表达基因来缩短程序运行时间。默认0.1
thresh.test参数：设定在两个细胞群中基因差异表达量。可以设置为0 ，程序运行时间会更长
max.cells.per.ident参数：每个类群细胞抽样设置；也可以缩短程序运行时间。
# find all markers of cluster 1
cluster1.markers <- FindMarkers(pbmc, ident.1 = 1, min.pct = 0.25)
head(cluster1.markers, n = 5)
# find all markers distinguishing cluster 5 from clusters 0 and 3(找到cluster5和cluster0、3之间的markgene)
cluster5.markers <- FindMarkers(pbmc, ident.1 = 5, ident.2 = c(0, 3), min.pct = 0.25)
head(cluster5.markers, n = 5)
# find markers for every cluster compared to all remaining cells, report only the positive ones(找到每个cluster相比于其余cluster的markgene，只报道阳性的markgene)
pbmc.markers <- FindAllMarkers(pbmc, only.pos = TRUE, min.pct = 0.25, logfc.threshold = 0.25)
pbmc.markers %>% group_by(cluster) %>% top_n(n = 2, wt = avg_logFC)
Seurat可以通过参数test.use设定检验差异表达的方法(详情见 DE vignett/ https://satijalab.org/seurat/archive/v3.0/de_vignette.html)。
cluster1.markers <- FindMarkers(pbmc, ident.1 = 0, logfc.threshold = 0.25, test.use = "roc", only.pos = TRUE)
head(cluster1.markers, n = 5)
添加细胞注释信息
如果你对每个cluster已经鉴定出细胞类型，就可以在可视化中将细胞类型标注在图中。
new.cluster.ids <- c("Naive CD4 T", "Memory CD4 T", "CD14+ Mono", "B", "CD8 T", "FCGR3A+ Mono",
"NK", "DC", "Platelet")
names(new.cluster.ids) <- levels(pbmc)
pbmc <- RenameIdents(pbmc, new.cluster.ids)
DimPlot(pbmc, reduction = "umap", label = TRUE, pt.size = 0.5) + NoLegend()
可视化Seurat有多种方法可视化标记基因的方法
VlnPlot: 基于细胞类群的基因表达概率分布
FeaturePlot：在tSNE 或 PCA图中画出基因表达情况
RidgePlot，CellScatter，DotPlot
VlnPlot(pbmc, features = c("MS4A1", "CD79A"))
# you can plot raw counts as well
VlnPlot(pbmc, features = c("NKG7", "PF4"), slot = "counts", log = TRUE)
#查看特定基因在聚类图中的表达量分布情况
FeaturePlot(pbmc, features = c("MS4A1", "GNLY", "CD3E", "CD14", "FCER1A", "FCGR3A", "LYZ", "PPBP",
"CD8A"))
#查看特定基因在cluster中的表达山脊图
features <- c("LYZ", "CCL5")
RidgePlot(pbmc, features = features, ncol = 2)
#查看特定基因在cluster中的分布点图，点的大小代表表达该基因的细胞比例，颜色代表平均表达水平
DotPlot(pbmc, features = features) + RotatedAxis()
DoHeatmap为指定的细胞和基因画表达热图。每个类群默认展示top 20标记基因
top10 <- pbmc.markers %>% group_by(cluster) %>% top_n(n = 10, wt = avg_logFC)
DoHeatmap(pbmc, features = top10$gene) + NoLegend()
Assigning cell type identity to clusters
根据已知细胞类型与基因标记的对应关系，可以为细胞类群找到对应的细胞类型。
new.cluster.ids <- c("Naive CD4 T", "Memory CD4 T", "CD14+ Mono", "B", "CD8 T", "FCGR3A+ Mono",
"NK", "DC", "Platelet")
names(new.cluster.ids) <- levels(pbmc)
pbmc <- RenameIdents(pbmc, new.cluster.ids)
DimPlot(pbmc, reduction = "umap", label = TRUE, pt.size = 0.5) + NoLegend()
saveRDS(pbmc, file = "../output/pbmc3k_final.rds")
参考文章：
https://www.jianshu.com/p/75c6d55a63bb?from=timeline
https://www.jianshu.com/p/44c1f526e450
https://www.bilibili.com/read/cv6411199/
https://www.jianshu.com/p/fef17a1babc2