写在前面
之前我们介绍了Seurat、Harmony，rliger三个包，用于3'和5'数据合并的方法。🤒
但有时候我们会遇到两个datasets只有部分重叠，这和之前介绍的方法就有一点不同了。🤨
用到的包
rm(list = ls())
library(Seurat)
library(SeuratDisk)
library(SeuratWrappers)
library(patchwork)
library(harmony)
library(rliger)
library(RColorBrewer)
library(tidyverse)
library(reshape2)
library(ggsci)
library(ggstatsplot)
示例数据
这里我们提供1个3’ PBMC dataset和1个whole blood dataset。😉
umi_gz <- gzfile("./GSE149938_umi_matrix.csv.gz",'rt')
umi <- read.csv(umi_gz,check.names = F,quote = "")
matrix_3p <- Read10X_h5("./3p_pbmc10k_filt.h5",use.names = T)
创建Seurat对象。🧐
srat_wb <- CreateSeuratObject(t(umi),project = "whole_blood")
srat_3p <- CreateSeuratObject(matrix_3p,project = "pbmc10k_3p")
rm(umi_gz)
rm(umi)
rm(matrix_3p)
srat_wb
srat_3p
修改metadata
为了方便后续分析，这里我们对metadata进行一下注释修改。 😁
colnames(srat_wb@meta.data)[1] <- "cell_type"
srat_wb@meta.data$orig.ident <- "whole_blood"
srat_wb@meta.data$orig.ident <- as.factor(srat_wb@meta.data$orig.ident)
head(srat_wb[[]])
基础质控
做一下标准操作，计算线粒体基因和核糖体基因。🥳
srat_wb <- SetIdent(srat_wb,value = "orig.ident")
srat_wb[["percent.mt"]] <- PercentageFeatureSet(srat_wb, pattern = "^MT-")
srat_wb[["percent.rbp"]] <- PercentageFeatureSet(srat_wb, pattern = "^RP[SL]")
srat_3p[["percent.mt"]] <- PercentageFeatureSet(srat_3p, pattern = "^MT-")
srat_3p[["percent.rbp"]] <- PercentageFeatureSet(srat_3p, pattern = "^RP[SL]")
p1 <- VlnPlot(srat_wb, ncol = 4,
features = c("nFeature_RNA","nCount_RNA","percent.mt","percent.rbp"))
p2 <- VlnPlot(srat_3p, ncol = 4,
features = c("nFeature_RNA","nCount_RNA","percent.mt","percent.rbp"))
p1/p2
交集基因
whole blood dataset使用的是Cell Ranger GRCh38-2020A进行注释，与3’ PBMC dataset差的比较多，所以我们先看一下有多少共同基因吧。🤩
# table(rownames(srat_3p) %in% rownames(srat_wb))
common_genes <- rownames(srat_3p)[rownames(srat_3p) %in% rownames(srat_wb)]
length(common_genes)
过滤基因
我们设置一下过滤条件，把一些表达过低或过高的细胞去掉，以及一些线粒体基因过高的细胞（细胞状态不佳）。✌️
srat_3p <- subset(srat_3p, subset = nFeature_RNA > 500 & nFeature_RNA < 5000 & percent.mt < 15)
srat_wb <- subset(srat_wb, subset = nFeature_RNA > 1000 & nFeature_RNA < 6000)
srat_3p <- srat_3p[rownames(srat_3p) %in% common_genes,]
srat_wb <- srat_wb[rownames(srat_wb) %in% common_genes,]
数据整合
8.1 合并为list
wb_list <- list()
wb_list[["pbmc10k_3p"]] <- srat_3p
wb_list[["whole_blood"]] <- srat_wb
8.2 Normalization与特征基因
for (i in 1:length(wb_list)) {
wb_list[[i]] <- NormalizeData(wb_list[[i]], verbose = F)
wb_list[[i]] <- FindVariableFeatures(wb_list[[i]], selection.method = "vst", nfeatures = 2000, verbose = F)
}
8.3 寻找Anchors并整合数据
wb_anchors <- FindIntegrationAnchors(object.list = wb_list, dims = 1:30)
wb_seurat <- IntegrateData(anchorset = wb_anchors, dims = 1:30)
rm(wb_list)
rm(wb_anchors)
整合效果可视化
9.1 整合前
DefaultAssay(wb_seurat) <- "RNA"
wb_seurat <- NormalizeData(wb_seurat, verbose = F)
wb_seurat <- FindVariableFeatures(wb_seurat, selection.method = "vst", nfeatures = 2000, verbose = F)
wb_seurat <- ScaleData(wb_seurat, verbose = F)
wb_seurat <- RunPCA(wb_seurat, npcs = 30, verbose = F)
wb_seurat <- RunUMAP(wb_seurat, reduction = "pca", dims = 1:30, verbose = F)
DimPlot(wb_seurat,reduction = "umap") +
scale_color_npg()+
plot_annotation(title = "10k 3' PBMC and whole blood, before integration")
9.2 整合后
DefaultAssay(wb_seurat) <- "integrated"
wb_seurat <- ScaleData(wb_seurat, verbose = F)
wb_seurat <- RunPCA(wb_seurat, npcs = 30, verbose = F)
wb_seurat <- RunUMAP(wb_seurat, reduction = "pca", dims = 1:30, verbose = F)
DimPlot(wb_seurat, reduction = "umap") +
scale_color_npg()+
plot_annotation(title = "10k 3' PBMC and white blood cells, after integration")
降维与聚类
10.1 聚类可视化
wb_seurat <- FindNeighbors(wb_seurat, dims = 1:30, k.param = 10, verbose = F)
wb_seurat <- FindClusters(wb_seurat, verbose = F)
ncluster <- length(unique(wb_seurat[[]]$seurat_clusters))
mycol <- colorRampPalette(brewer.pal(8, "Set2"))(ncluster)
DimPlot(wb_seurat,label = T, reduction = "umap",
cols = mycol, repel = T) +
NoLegend()
10.2 具体查看及可视化
count_table <- table(wb_seurat@meta.data$seurat_clusters,
wb_seurat@meta.data$orig.ident)
count_table
#### 可视化
count_table %>%
as.data.frame() %>%
ggbarstats(x = Var2,
y = Var1,
counts = Freq)+
scale_fill_npg()
<img src="https://upload-images.jianshu.io/upload_images/24475539-a8a83e37395b6820.png" alt="面包" style="zoom:25%;" />
<center>最后祝大家早日不卷!~</center>
需要示例数据的小伙伴，在公众号回复
Merge2获取吧！点个在看吧各位~ ✐.ɴɪᴄᴇ ᴅᴀʏ 〰
<center> <b>📍 往期精彩 <b> </center>
📍 <font size=1>🤨 Google | 谷歌翻译崩了我们怎么办！？（附完美解决方案）</font>
📍 <font size=1>🤩 scRNA-seq | 吐血整理的单细胞入门教程</font>
📍 <font size=1>🤔 Reticulate | 如何在Rstudio中优雅地调用Python！？</font>
📍 <font size=1>🤣 NetworkD3 | 让我们一起画个动态的桑基图吧~</font>
📍 <font size=1>🤩 RColorBrewer | 再多的配色也能轻松搞定！~</font>
📍 <font size=1>🧐 rms | 批量完成你的线性回归</font>
📍 <font size=1>🤒 CMplot | 连Nature上的曼哈顿图都卷起来啦</font>
📍 <font size=1>🤩 CMplot | 完美复刻Nature上的曼哈顿图</font>
📍 <font size=1>🤠 Network | 高颜值动态网络可视化工具</font>
📍 <font size=1>🤗 boxjitter | 完美复刻Nature上的高颜值统计图</font>
📍 <font size=1>🤫 linkET | 完美解决ggcor安装失败方案（附教程）</font>
📍 <font size=1>......</font>
本文由mdnice多平台发布
