自从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
网友评论