Seurat是一个r包，被设计用于单细胞rna-seq数据的细胞质控和分析
标准分析的流程：
image.png
https://www.jianshu.com/p/0727cabbc9c0
1：创建Seurat对象
library(dplyr)
library(Seurat)
# 使用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
2: 细胞选择
#使用从MT-作为线粒体基因集
pbmc[["percent.mt"]] <- PercentageFeatureSet(pbmc, pattern = "^MT-")
#qc指标可视化
VlnPlot(pbmc, features = c("nFeature_RNA", "nCount_RNA", "percent.mt"), ncol = 3)
# 过滤具有超过2500或少于200个基因计数的细胞，同时过滤掉线粒体比例超过5%的细胞。
pbmc <- subset(pbmc, subset = nFeature_RNA > 200 & nFeature_RNA < 2500 & percent.mt < 5)
3: 数据标准化
# 默认情况下，我们使用“lognormalize”全局缩放的归一化方法，通过总表达值对每个细胞的基因表达值归一化，并将其乘以缩放因子(默认为10,000)，最后对结果进行对数变换，标准化数据存储在pbmc@assays$RNA@data中。
pbmc <- NormalizeData(pbmc, normalization.method = "LogNormalize", scale.factor = 10000)
4: 高变异基因选择
# 在下游分析中关注这些基因有助于突出单细胞数据集中的生物信号。默认情况下，每个数据集选择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))
5: 缩放数据
# 使用线性变换(“scaling”)处理数据，过程：改变每个基因的表达，使细胞间的平均表达为0；缩放每个基因的表达，使细胞间的差异为1，这样高表达基因就不会影响后续分析；这步是pca等降维技术之前的标准预处理步骤，缩放的数据储存在pbmc@assays$RNA@scale.data中。
all.genes <- rownames(pbmc)
pbmc <- ScaleData(pbmc, features = all.genes)
6: 线性降维
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)
7: 确定数据集的维度
每个pc本质上代表一个“元特征”，它将相关特征集中的信息组合在一起。因此，越在顶部的主成分越可能代表数据集。可以使用JackStraw函数来选择PC数目：
#将数据的1%打乱重新运行pca，构建特征得分的“零分布”，这个过程重复100次
pbmc <- JackStraw(pbmc, num.replicate = 100)
pbmc <- ScoreJackStraw(pbmc, dims = 1:20)
plot1<-JackStrawPlot(pbmc, dims = 1:15)
plot2<-ElbowPlot(pbmc)
CombinePlots(plots = list(plot1, plot2))
- 细胞聚类
#计算最邻近距离
pbmc <- FindNeighbors(pbmc, dims = 1:10)
#聚类，包含设置下游聚类的“间隔尺度”的分辨率参数resolution ，增加值会导致更多的聚类。
pbmc <- FindClusters(pbmc, resolution = 0.5)
- 运行非线性降维
几种非线性降维技术，例如tsne和umap，
pca是线性降维，能很好处理高维度数据，
tsne和umap是非线性降维方式，对低维度数据处理速度快，但难以处理高维度数据。相比tsne，umap能更好地反映原始数据结构，有着更好的连续性。
pbmc <- RunUMAP(pbmc, dims = 1:10)
pbmc <- RunTSNE(pbmc, dims = 1:10)
DimPlot(pbmc, reduction = "umap")
DimPlot(pbmc, reduction = "tsne")
Run开头的函数降维，Find开头的函数聚类
PCA将原来2000维的数据降到50维，dims参数表示使用多少个主成分（一般20左右就可以了，多几个少几个对结果影响不大），resolution参数表达聚类的分辨率，这个值大于0，一般都是在0-1范围里面调整，越大得到的cluster越多，这个值可以反复调整，并不会改变降维的结果（也就是tsne、umap图的二维坐标）
链接：https://www.jianshu.com/p/ce6205fe1ff9
10： 差异分析
seurat可以通过差异表达找到每个聚类的markgene，差异分析可以有多种形式，如找到所有聚类的markene（如cluster1中所有的markgene是指cluster1相对于其余所有cluster是差异的）、两个cluster之间的差异分析、某个cluster中两个样品之间差异分析等
#找到cluster1中的markgene
cluster1.markers <- FindMarkers(pbmc, ident.1 = 1, min.pct = 0.25)
#找到cluster5和cluster0、3之间的markgene
cluster5.markers <- FindMarkers(pbmc, ident.1 = 5, ident.2 = c(0, 3), min.pct = 0.25)
#找到每个cluster相比于其余cluster的markgene，只报道阳性的markgene
pbmc.markers <- FindAllMarkers(pbmc, only.pos = TRUE, min.pct = 0.25, logfc.threshold = 0.25)
网友评论