Seurat官网教程复现 以及相关问题解决

Seurat是一个针对于单细胞RNA-seq数据处理，探索以及分析的R包。

目录

#创建Seurat项目

#数据预处理

#质控（QC）

#标准化数据

#识别高度可变特征的基因

#规范化数据

#线性降维并可视化

#确定数据维度

#细胞聚类

#非线性降维

#识别簇中的细胞类型

---

在Seurat4.0版本中，加入了综合多模态分析，采用加权近邻分析来定义细胞状态；还加入了数据集间的快速比对方法，引入Azimuth，即利用高质量参考数据集快速映射新的scRNA-seq数据集。

下面以官网指导教程（https://satijalab.org/seurat/articles/pbmc3k_tutorial.html）为例，对Seurat进行初步熟悉。

教程中分析2700例由二代测序产生出的单细胞外周血细胞数据。

原始数据下载地址https://cf.10xgenomics.com/samples/cell/pbmc3k/pbmc3k_filtered_gene_bc_matrices.tar.gzhttps://cf.10xgenomics.com/samples/cell/pbmc3k/pbmc3k_filtered_gene_bc_matrices.tar.gz

#导入R包

dplyr包包含管道函数%>%，patchwork包可以将图形之间进行组合，二者辅助Seurat包。

install.packages("Seurat")
library(Seurat)
install.packages("dplyr")
library(dplyr)
install.packages("patchwork")
library(patchwork)

#创建Seurat项目

#导入数据

pbmc.data <- Read10X(data.dir = “F:/bio/AU/AU1/Seurat/filtered_gene_bc_matrices/hg19”)

#对原始数据构建Seurat项目--过滤数据集中测到的少于200个基因的细胞（min.features = 200）和少于3个细胞覆盖的基因（min.cells = 3）

pbmc <- CreateSeuratObject(counts = pbmc.data, project = "pbmc3k", min.cells = 3, min.features = 200)

#数据预处理

#质控（QC）

#QC指标：

1. 每个细胞测得独特基因的数目（unique feature）--低质量与空液滴细胞中一般不含基因，细胞的二重体或多重体可能表现出异常高的基因计数。
2. 细胞内检测到的分子总数
3. 与线粒体基因组相匹配的读取率--低质量/死亡细胞常表现出线粒体污染，可用PercentageFeatureSet（）函数计算线粒体QC度量
4. 低质量/死亡的细胞常表现出广泛的线粒体污染
5. 使用从MT-开始的所有基因集合作为一组线粒体基因

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

#绘制小提琴图--可视化QC指标，用来过滤细胞（过滤独特基因超过2500或低于200，线粒体计数>5%的细胞），nFeature_RNA代表每个细胞测到的基因数目，nCount代表每个细胞测到所有基因的表达量之和，percent.mt代表测到的线粒体基因的比例

VlnPlot(pbmc, features = c("nFeature_RNA", "nCount_RNA", "percent.mt"), ncol = 3)

#可视化特征间的关系--FeatureScatter创建两个特征间的散点图

plot1 <- FeatureScatter(pbmc, feature1 = "nCount_RNA", feature2 = "percent.mt")
plot2 <- FeatureScatter(pbmc, feature1 = "nCount_RNA", feature2 = "nFeature_RNA")
plot1 + plot2 

#保留独特基因低于2500或高于200，线粒体计数<5%的细胞并储存

pbmc <- subset(pbmc, subset = nFeature_RNA > 200 & nFeature_RNA < 2500 & percent.mt < 5)

#标准化数据

采用全局尺度归一化的方法将每个单元格的特征表达式度量归一化，并将其乘以一个尺度因子(默认情况下为10,000)，然后对结果进行对数转换。

pbmc <- NormalizeData(pbmc, normalization.method = "LogNormalize", scale.factor = 10000)

#识别高度可变特征的基因

鉴别出细胞间表达量差异较大的基因。通过FindVariableFeatures ()直接建模单细胞数据中固有的均值-方差关系。

pbmc <- FindVariableFeatures(pbmc, selection.method = "vst", nfeatures = 2000)

#识别前10个变异基因并保存
top10 <- head(VariableFeatures(pbmc), 10)

#可视化高度可变细胞
plot1 <- VariableFeaturePlot(pbmc)
plot2 <- LabelPoints(plot = plot1, points = top10, repel = TRUE)
plot1 + plot2 

以我为例，出图时可能会报错：

Error in grid.Call(C_convert, x, as.integer(whatfrom),as.integer(whatto), : Viewport has zero dimension(s) ---可能是视图窗口太小。调整视图窗口就ok

#规范化数据

使细胞间平均表达为0，方差为1。

all.genes <- rownames(pbmc)
pbmc <- ScaleData(pbmc, features = all.genes)

#简化方法（可只执行这一步）
pbmc <- ScaleData(pbmc)

#线性降维并可视化

pbmc <- RunPCA(pbmc, features = VariableFeatures(object = pbmc))

#用不同的方式对PCA结果进行检查与可视化--Seurat包中有VizDimReduction(), DimPlot(),和DimHeatmap()，三种可视化PCA结果的方法。

print(pbmc[["pca"]], dims = 1:5, nfeatures = 5)
VizDimLoadings(pbmc, dims = 1:2, reduction = "pca")

DimPlot(pbmc, reduction = "pca")

 

 DimHeatmap()单元格和特征都根据它们的PCA分数进行排序。将细胞设置为一个数字可以绘制出光谱两端的“极端”细胞

DimHeatmap(pbmc, dims = 1, cells = 500, balanced = TRUE)

DimHeatmap(pbmc, dims = 1:15, cells = 500, balanced = TRUE)

 

#确定数据维度

 研究pc来确定异构性的相关来源。

JackStraw() 统计PCA评分的统计意义，设置重复取样100次。

ScoreJackStraw()计算Jackstraw的得分差异。

pbmc <- JackStraw(pbmc, num.replicate = 100)
pbmc <- ScoreJackStraw(pbmc, dims = 1:20)

•       #法二--JackStrawPlot()--比较每台PC的p值分布是否均匀(虚线)，基于随机空模型的统计测试，但对于大型数据集来说很耗时，而且可能不会返回一个明确的PC截断值

JackStrawPlot(pbmc, dims = 1:15)

• 法三“肘部图”--ElbowPlot()基于每个成分所解释的方差百分比对主要成分进行排序，属于启发式算法，可以即时计算

ElbowPlot(pbmc)

#细胞聚类

FindNeighbors() 在主成分空间中基于欧氏距离构建KNN图，然后基于相邻区域的共享重叠，定义任意两个单元格之间的边权值。输入值：前10个PCs

pbmc <- FindNeighbors(pbmc, dims = 1:10)

基于共享近邻（SNN）法来识别细胞簇。参数设置在0.4-1.2之间，对3k细胞的数据集有良好的结果，数据集越大，最佳分辨率越高

pbmc <- FindClusters(pbmc, resolution = 0.5) 

 Idents()可查看聚簇结果

head(Idents(pbmc), 5)

#非线性降维

tSNE和UMAP是Seurat进行非线性降维的方法。

pbmc <- RunUMAP(pbmc, dims = 1:10)
head(pbmc@reductions$umap@cell.embeddings)  #提取UMAP坐标值
DimPlot(pbmc, reduction = "umap")  #降维方法首选umap，再tsne，最后再pca

 #在聚簇中找差异表达的特征

识别单个集群中阳性及阴性的标记

ident.1参数设置待分析的细胞类别，min.pct表示该基因表达数目占该类细胞总数的比例，找出簇2中的所有标记物

cluster2.markers <- FindMarkers(pbmc, ident.1 = 2, min.pct = 0.25)
head(cluster2.markers, n = 5)

找出簇5与簇0/3之间的差异表达标志物

cluster5.markers <- FindMarkers(pbmc, ident.1 = 5, ident.2 = c(0, 3), min.pct = 0.25)
head(cluster5.markers, n = 5)

在所有剩余细胞中找到每簇的标记物，并只报道阳性位点

pbmc.markers <- FindAllMarkers(pbmc, only.pos = TRUE, min.pct = 0.25, logfc.threshold = 0.25)
pbmc.markers %>%
    group_by(cluster) %>%
    slice_max(n = 2, order_by = avg_log2FC)
#这里的管道函数%>%需要调用dplyr包，将左边的数据导入右边，按照cluster进行排序，并返回avg_log2FC中变量最高的2列

ROC测试返回任何单个标记的“分类能力”(范围从0 -随机，到1 -完美)。

cluster0.markers <- FindMarkers(pbmc, ident.1 = 0, logfc.threshold = 0.25, test.use = "roc", only.pos = TRUE)

VlnPlot()和FeaturePlot()用来可视化tSNE和PCA的位点

VlnPlot(pbmc, features = c("MS4A1", "CD79A"))
VlnPlot(pbmc, features = c("NKG7", "PF4"), slot = "counts", log = TRUE)

FeaturePlot(pbmc, features = c("MS4A1", "GNLY", "CD3E", "CD14", "FCER1A", "FCGR3A", "LYZ", "PPBP", "CD8A"))

DoHeatmap()为给定的细胞和特征生成表达热图

pbmc.markers %>%
    group_by(cluster) %>%
    top_n(n = 10, wt = avg_log2FC) -> top10
DoHeatmap(pbmc, features = top10$gene) + NoLegend()
#产生top10基因的表达热图，并去掉图例

#识别簇中的细胞类型

new.cluster.ids <- c("Naive CD4 T", "CD14+ Mono", "Memory CD4 T", "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")

 Seurat的基本用法如以上所示。

敬请批评与指点。

相关问题

 解决方法：

根据自身数据量大小进行设置

options(future.globals.maxSize = 1000000 * 1024^2)