自從seurat V5更新之后呢，很多小伙伴观堂，初學(xué)者居多吧让网，都有點(diǎn)不適應(yīng)，再加上網(wǎng)上有些人的“煽風(fēng)點(diǎn)火”师痕，導(dǎo)致大家望而卻步溃睹，好像這次更新非常可怕一樣胰坟。其實(shí)不然因篇，seurat的更新在我看來(lái)并沒(méi)有多大的變化，不必望而生畏笔横。此外竞滓，他的更新也是非常好的，首先第一點(diǎn)如他官網(wǎng)上所述吹缔，數(shù)據(jù)結(jié)構(gòu)發(fā)生了很大的改變商佑，這樣在運(yùn)行的時(shí)候不會(huì)耗費(fèi)太多的內(nèi)容。其次我認(rèn)為最好的地方就是數(shù)據(jù)整合這里厢塘，將目前比較優(yōu)秀的方法通過(guò)一句代碼實(shí)現(xiàn)茶没，非常方便肌幽。這次我們演示一下它的基本分析，其實(shí)很簡(jiǎn)單礁叔，也不會(huì)有太大的問(wèn)題牍颈。

首先我們下載安裝相關(guān)的軟件，讀入數(shù)據(jù)琅关，數(shù)據(jù)讀取沒(méi)有什么變化煮岁！

#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)

數(shù)據(jù)質(zhì)控什么的和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質(zhì)控
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

接下來(lái)就是數(shù)據(jù)整合了，也是他更新的地方涣易。提供了很多方法画机，例如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")

然后就是細(xì)胞注釋了：我的建議還是手動(dòng)荚醒！

#================================================================================
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一樣：

#所有細(xì)胞類型中兩組的差異
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的基本分析了，就這么簡(jiǎn)單隆嗅，沒(méi)有什么難的地方界阁。其他詳細(xì)內(nèi)容請(qǐng)?jiān)诠倬W(wǎng)觀看，給出的步驟很詳細(xì)了胖喳！詳細(xì)請(qǐng)參考：https://mp.weixin.qq.com/s/s0FlOruxzPEYcXfwfJd33Q