单细胞Seurat V5分析流程

作者: KS科研分享与服务 | 来源:发表于2024-03-17 23:59 被阅读0次

单细胞实战2: Seurat+SingleR--从矩阵到细胞类型
Harmony整合单细胞数
跟着Cell学单细胞转录组分析(二):单细胞转录组测序文件的读入
Seurat使用教程（v3.0）
实用干货 | 运用Seurat包做单细胞数据的合并分析
Seurat包学习笔记（三）：Analysis of spati
单细胞转录组Seurat
答读者问(6)：单细胞TPM矩阵如何分析？
sc-RAN-seq 数据分析||Seurat新版教程:Guid
实用干货|Seurat-空间转录组数据基本分析教学

自从seurat V5更新之后呢，很多小伙伴，初学者居多吧，都有点不适应，再加上网上有些人的“煽风点火”，导致大家望而却步，好像这次更新非常可怕一样。其实不然，seurat的更新在我看来并没有多大的变化，不必望而生畏。此外，他的更新也是非常好的，首先第一点如他官网上所述，数据结构发生了很大的改变，这样在运行的时候不会耗费太多的内容。其次我认为最好的地方就是数据整合这里，将目前比较优秀的方法通过一句代码实现，非常方便。这次我们演示一下它的基本分析，其实很简单，也不会有太大的问题。

首先我们下载安装相关的软件，读入数据，数据读取没有什么变化！

#install packages
install.packages('Seurat')
library(Seurat)
#安装一些额外的包
setRepositories(ind = 1:3, addURLs = c('https://satijalab.r-universe.dev', 'https://bnprks.r-universe.dev/'))
install.packages(c("BPCells", "presto", "glmGamPoi"))
remotes::install_github("satijalab/seurat-data", quiet = TRUE)
remotes::install_github("satijalab/azimuth", quiet = TRUE)
remotes::install_github("satijalab/seurat-wrappers", quiet = TRUE)
# If users encounter any errors related to the Matrix package, please resolve by re-installing the TFBSTools package using the command below and opening a fresh R session:
BiocManager::install("TFBSTools", type = "source", force = TRUE)

setwd("/home/ks_ts/data_analysis/seuratV5_test/scRNA_analysis/")

#read data and creat seurat obj
WT <- Read10X("./scRNA_data/WT_E18/")
WT <- WT$`Gene Expression`
WT <- CreateSeuratObject(counts = WT, project = "WT", min.cells = 3, min.features = 200)

GO <- Read10X("./scRNA_data/GO_E18/")
GO <- GO$`Gene Expression`
GO <- CreateSeuratObject(counts = GO, project = "GO", min.cells = 3, min.features = 200)

数据质控什么的和V4一样：


#线粒体比例
WT[["percent.mt"]] <- PercentageFeatureSet(WT, pattern = "^mt-")
GO[["percent.mt"]] <- PercentageFeatureSet(GO, pattern = "^mt-")
#血红蛋白基因
WT[["percent.hb"]] <- PercentageFeatureSet(WT, pattern = "^Hb[^(p)]")
GO[["percent.hb"]] <- PercentageFeatureSet(GO, pattern = "^Hb[^(p)]")

#核糖体基因
WT[["percent.rb"]] <- PercentageFeatureSet(WT, pattern = "^Rbs|Rpl")
GO[["percent.rb"]] <- PercentageFeatureSet(GO, pattern = "^Rbs|Rpl")

p1= VlnPlot(WT, features = c("nFeature_RNA", "nCount_RNA", "percent.mt","percent.hb","percent.rb"),pt.size = 0.1, ncol = 5)

p2 = VlnPlot(GO, features = c("nFeature_RNA", "nCount_RNA", "percent.mt","percent.hb","percent.rb"),pt.size = 0.1, ncol = 5)
p1/p2

#QC质控
WT <- subset(WT, subset = nFeature_RNA > 200 & nFeature_RNA < 6000 & nCount_RNA < 30000 &  percent.mt < 5)
GO <- subset(GO, subset = nFeature_RNA > 200 & nFeature_RNA < 6000 & nCount_RNA < 30000 &  percent.mt < 5)

p3= VlnPlot(WT, features = c("nFeature_RNA", "nCount_RNA", "percent.mt","percent.hb","percent.rb"),pt.size = 0.1, ncol = 5)
p4 = VlnPlot(GO, features = c("nFeature_RNA", "nCount_RNA", "percent.mt","percent.hb","percent.rb"),pt.size = 0.1, ncol = 5)
p3/p4

接下来就是数据整合了，也是他更新的地方。提供了很多方法，例如CCA、Harmony、scVI、RPCA等等，选择适合自己的即可！


#merge data
sce <- merge(WT, y=GO,add.cell.ids = c("WT", "GO"))


#NormalizeData & ScaleData
sce <- NormalizeData(sce)
sce <- FindVariableFeatures(sce)
sce <- ScaleData(sce, vars.to.regress = c("percent.mt"))
sce <- RunPCA(sce, verbose=F)

# methods of Integration
# CCA integration (method=CCAIntegration)
# RPCA integration (method=RPCAIntegration)
# Harmony (method=HarmonyIntegration)
# JointPCA (method= JointPCAIntegration)

# FastMNN (method= FastMNNIntegration)
# scVI (method=scVIIntegration)

###Integrated with CCA
sce_cca <- IntegrateLayers(object = sce, 
                           method = CCAIntegration, 
                           orig.reduction = "pca", 
                           new.reduction = "integrated.cca",verbose = FALSE)

# re-join layers after integration
sce_cca[["RNA"]] <- JoinLayers(sce_cca[["RNA"]])


sce_scvi <- IntegrateLayers(object = sce, 
                            method = scVIIntegration, 
                            orig.reduction = "pca", 
                            new.reduction = "integrated.scvi",
                            conda_env="/home/ks_ts/miniconda3/envs/scvi-env",
                            verbose = FALSE)

# re-join layers after integration
sce_scvi[["RNA"]] <- JoinLayers(sce_scvi[["RNA"]])

#================================================================================
#Perform cca reduction
Seurat::ElbowPlot(sce_cca, ndims = 50)
sce_cca <- FindNeighbors(sce_cca, reduction = "integrated.cca", dims = 1:20)
sce_cca <- FindClusters(sce_cca, resolution = seq(from = 0.1, to = 1.0, by = 0.1))
sce_cca <- RunUMAP(sce_cca, dims = 1:20, reduction = "integrated.cca")
# clustree(sce_cca)

#Perform scVI reduction
Seurat::ElbowPlot(sce_scvi, ndims = 50)
sce_scvi <- FindNeighbors(sce_scvi, reduction = "integrated.scvi", dims = 1:20)
sce_scvi <- FindClusters(sce_scvi, resolution = seq(from = 0.1, to = 1.0, by = 0.1))
sce_scvi <- RunUMAP(sce_scvi, dims = 1:20, reduction = "integrated.scvi")
# clustree(sce_scvi)


#cluster plot
DimPlot(sce_cca, reduction = "umap", group.by = "orig.ident")+
  ggtitle("CCA")
DimPlot(sce_scvi, reduction = "umap", group.by = "orig.ident")+
  ggtitle("scvi")


DimPlot(sce_cca, reduction = "umap", label = T)+
  ggtitle("CCA")
DimPlot(sce_scvi, reduction = "umap", label = T)+
  ggtitle("scvi")

然后就是细胞注释了：我的建议还是手动！


#================================================================================
library(Seurat)
library(ggplot2)

Allmarkers <- FindAllMarkers(sce_cca, logfc.threshold = 0.3, min.pct = 0.3, only.pos = T)
write.csv(Allmarkers, file = 'Allmarkers.csv')
#================================================================================
#Manual annotation, reference to published articles
markers <- c("Pparg", "Myh11", "Mrc1", "Flt1", "Col11a1", "Mymk", "Pax7", "Pdgfra","Ttn","Sox2")
DotPlot(sce_cca, features = markers, col.min = 0)+coord_flip()
FeaturePlot(sce_cca, features = )

#20 Adipocytes
#19 SMC
#13 Macrophages
#14 Endothelial
#7,10,21 Tenocytes
#9 Myoblasts
#11,12 MuSCs
#0,1,3,22 Mesenchymal
#2,4,5,6,8,15,16,18 Myonuclei
#17 NPCs

差异基因的分析和V4一样：

#所有细胞类型中两组的差异
celltypes <- unique(sce_cca$celltype)
DEGs_celltype <- list()

for (i in 1:length(celltypes)) {
  
  data = subset(sce_cca, celltype==celltypes[i])
  deg = FindMarkers(data,
                    group.by="orig.ident",
                    ident.1 = 'GO',
                    ident.2 = "WT",
                    logfc.threshold=0.25,
                    min.pct = 0.25)
  
  DEGs_celltype[[i]] <- deg
  names(DEGs_celltype)[i] <- celltypes[i]
  
}

这就是Seurat V5的基本分析了，就这么简单，没有什么难的地方。其他详细内容请在官网观看，给出的步骤很详细了！详细请参考：https://mp.weixin.qq.com/s/s0FlOruxzPEYcXfwfJd33Q

单细胞实战2: Seurat+SingleR--从矩阵到细胞类型
单细胞数据标准分析Seurat流程： Seurat官网：https://satijalab.org/seurat/...
Harmony整合单细胞数
Seurat结合Harmony的整合流程参考：单细胞测序分析: Seurat V3联合harmony进行单细胞数...
跟着Cell学单细胞转录组分析(二):单细胞转录组测序文件的读入
分析单细胞转录组测序的软件和方法有很多，最流行的莫过于Seurat包，可以完成单细胞分析整个流程，我们整个教程也是...
Seurat使用教程（v3.0）
Seurat使用教程（v3.0） Seurat是一个分析单细胞转录组数据的R包，用于QC，分析和探索单细胞RNA-...
实用干货 | 运用Seurat包做单细胞数据的合并分析
之前的文章中为大家介绍了Seurat这款单细胞分析软件（实用干货|Seurat-单细胞转录组测序分析R包教学[ht...
Seurat包学习笔记（三）：Analysis of spati
本次教程中，我们将学习如何使用Seurat3处理空间转录组数据。整体的分析流程类似于Seurat的单细胞RNA-s...
单细胞转录组Seurat
引言：因为在尝试单细胞转录组分析流程，学习了Seurat的数据分析步骤，记录主要的分析过程，及一些心得。 Seur...
答读者问(6)：单细胞TPM矩阵如何分析？
问题一、有的文章只提供TPM的单细胞表达矩阵，可以用seurat分析吗？二、分析流程和用count矩阵有什么不...
sc-RAN-seq 数据分析||Seurat新版教程:Guid
Seurat - Guided Clustering Tutorial Seurat，一个单细胞数据分析工具箱十...
实用干货|Seurat-空间转录组数据基本分析教学
上一篇文章对用Seurat分析单细胞数据进行了比较详细的介绍（实用干货|Seurat-单细胞转录组测序分析R包教学...