Seurat是一个针对于单细胞RNA-seq数据处理,探索以及分析的R包。
目录
在Seurat4.0版本中,加入了综合多模态分析,采用加权近邻分析来定义细胞状态;还加入了数据集间的快速比对方法,引入Azimuth,即利用高质量参考数据集快速映射新的scRNA-seq数据集。
下面以官网指导教程(https://satijalab.org/seurat/articles/pbmc3k_tutorial.html)为例,对Seurat进行初步熟悉。
教程中分析2700例由二代测序产生出的单细胞外周血细胞数据。
#导入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指标:
- 每个细胞测得独特基因的数目(unique feature)--低质量与空液滴细胞中一般不含基因,细胞的二重体或多重体可能表现出异常高的基因计数。
- 细胞内检测到的分子总数
- 与线粒体基因组相匹配的读取率--低质量/死亡细胞常表现出线粒体污染,可用PercentageFeatureSet()函数计算线粒体QC度量
- 低质量/死亡的细胞常表现出广泛的线粒体污染
- 使用从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)