ixxmu
commented
12 months ago seurat个性化细胞注释并把细分亚群放回总群 by 生信小博士
大家好，今天我们分享的是
1. 个性化细胞命名
2. 细分亚群
3. 把细分的亚群放回原来的总群
公众号内回复冒号后面的关键字获取直播代码：第四次直播
直播视频已经上传到b站up主： 生信小博士https://space.bilibili.com/1706877377
01
1. 个性化细胞命名
首先读入pbmc数据
.libPaths(c( "/home/data/t040413/R/x86_64-pc-linux-gnu-library/4.2","/home/data/t040413/R/yll/usr/local/lib/R/site-library","/home/data/refdir/Rlib/", "/usr/local/lib/R/library"))sessionInfo()getwd()dir.create('~/gzh/20231202_4_live_个性化注释')setwd('~/gzh/20231202_4_live_个性化注释')getwd()library(Seurat)pbmc=readRDS("~/gzh/featureplot_dotplot_vlnplot/pbmc3k_final.rds")DimPlot(pbmc,label=TRUE)pbmc$cell_type=pbmc$cell.type
未更改细胞id之前的umap图
然后我们进行单细胞个性化细胞命名
1# 1 选定自己想要命名的细胞----All.merge=pbmcplot <- DimPlot(All.merge, reduction = "umap")library(dplyr)plot$data %>%head()select.cells <- CellSelector(plot = plot) #运行到这一步，可以交互式的选择自己想要的细胞
head(select.cells)# subset_data=subset(All.merge,cells=select.cells)# DimPlot(subset_data,label = TRUE)## head(All.merge@meta.data)# head(subset_data@meta.data)All.merge=SetIdent(object = All.merge,cells = select.cells ,value = 'B')
改完名称之后，结果如下
这样就完成了单细胞个性化的细胞注释
02
2. 细分亚群
紧接着我们演示如何细分亚群
这里假如我们想对b细胞进行细分亚群
2 #2 细分亚群----bcell=subset(pbmc,idents = 'B')DimPlot(bcell,label = TRUE)|DimPlot(pbmc,label = TRUE)
细分亚群的代码如下
pbmc=bcell{2#2 将 QC 指标可视化为小提琴图----# [[ 运算符可以向对象 metadata 添加列。这是存储 QC 统计数据的好地方pbmc[["percent.mt"]] <- PercentageFeatureSet(pbmc, pattern = "^MT-")pbmc[["ribosomal.mt"]] <- PercentageFeatureSet(pbmc, pattern = "^RP[SL]")head(pbmc@meta.data)VlnPlot(pbmc, features = c("nFeature_RNA", "nCount_RNA", "percent.mt"), ncol = 3)# FeatureScatter 通常用于可视化特征与特征之间的关系，但也可以用于由对象计算的任何内容，即对象 metadata 中的列、PC 分数等。plot1 <- FeatureScatter(pbmc, feature1 = "nCount_RNA", feature2 = "percent.mt")plot2 <- FeatureScatter(pbmc, feature1 = "nCount_RNA", feature2 = "nFeature_RNA")plot1 + plot2# 显示前 5 个细胞的 QC 指标head(pbmc@meta.data, 5)#nFeature_RNA is the number of genes detected in each cell.# nCount_RNA is the total number of molecules detected within a cell.2.1 #2.1 qc---------pbmc <- subset(pbmc, subset = nFeature_RNA > 200 & nFeature_RNA < 2500 & percent.mt < 5)3 #3 数据归一化----pbmc <- NormalizeData(pbmc, normalization.method = "LogNormalize", scale.factor = 10000)4# 4 识别高变特征（特征选择）----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 + plot25 # 5 标准化数据------all.genes <- rownames(pbmc)pbmc <- ScaleData(pbmc, features = all.genes)5.1pbmc <- ScaleData(pbmc)5.2head(pbmc@meta.data)pbmc <- ScaleData(pbmc, vars.to.regress = c('ribosomal.mt', 'percent.mt'))#SCTransform()6#6 线性降维----pbmc <- RunPCA(pbmc, features = VariableFeatures(object = pbmc))# 通过几种不同的方式检查和可视化 PCA 结果print(pbmc[["pca"]], dims = 1:5, nfeatures = 5)VizDimLoadings(pbmc, dims = 1:2, reduction = "pca")7# 7 确定数据集的维度------# 注意：对于大数据集，此过程可能需要很长时间，为了方便起见，请注释掉。更多的近似技术，例如在 ElbowPlot() 中实现的技术，可用于减少计算时间pbmc <- JackStraw(pbmc, num.replicate = 100)pbmc <- ScoreJackStraw(pbmc, dims = 1:20)ElbowPlot(pbmc)8 #8 细胞聚类\分群-------pbmc <- FindNeighbors(pbmc, dims = 1:15)pbmc <- FindClusters(pbmc, resolution = 1)# 查看前 5 个细胞的 cluster IDshead(Idents(pbmc), 5)9#9 运行非线性降维 (UMAP/tSNE)----# 如果你还没有安装 UMAP，你可以通过 reticulate::py_install(packages = 'umap-learn')pbmc <- RunUMAP(pbmc, dims = 1:15)pbmc <- RunTSNE(pbmc, dims = 1:15)# 请注意，您可以设置 `label = TRUE` 或使用 LabelClusters 函数来帮助标记单个类群DimPlot(pbmc, reduction = "umap")bcell=pbmc}
最后我们就获得了分好群的b细胞。 我们发现分辨率为1时，b细胞可以分为四个亚群，但是否每个亚群都有生物学意义，自行判断哦！
03
3. 把细分的b细胞亚群放回原来的总群 （这一部分由于时间紧迫，没在视频里演示，这里我给出代码）
——现在的问题其实就是如何把b细胞的4个亚群放回总群里
#3 把细分的亚群放回原来的总群----pbmc=readRDS("~/gzh/featureplot_dotplot_vlnplot/pbmc3k_final.rds")dim(bcell)dim(pbmc)Idents(pbmc, cells = colnames(bcell)) <- Idents(bcell)DimPlot(pbmc,label = TRUE)
是不是特别简单~
如果你看不懂今天的教程，那你可以看下前三场视频热热身，再回来看哦~
文末友情宣传
看完记得顺手点个“赞”哦！
https://mp.weixin.qq.com/s/z2mPfyiETtQ_BqZA4E0YWQ