上面那个帖子介绍了FindIntegrationAnchors和IntegrateData，说也适合整合两个datasets只有部分重叠。那么，从这个角度出发，也适合整合两个跨物种的数据集，利用重叠的homolog基因来实现重叠的部分，所以我们拿arabidopsis和rice的单细胞数据来试试看（随便从paper里面选了一组拟南芥根和水稻根的数据作为测试而已）。
注：我是用cellranger获得了比对矩阵文件。然后根据arabidopsis和rice的ortholog基因，把rice ortholog的基因ID换成了arabidopsis，方便后面构造arabidopsis和rice的可比文件。
library(Seurat)
library(SeuratData)
library(patchwork)
library(harmony)
#library(rliger)
library(RColorBrewer)
library(tidyverse)
library(reshape2)
library(ggsci)
library(ggstatsplot)
#常规的构建seurat对象
arabidopsis.data <- Read10X(data.dir = "arabidopsis/filtered_feature_bc_matrix/")
rice.data <- Read10X(data.dir = "rice/filtered_feature_bc_matrix/")
arabidopsis <- CreateSeuratObject(counts = arabidopsis.data, project = "arabidopsis", min.cells = 3, min.features = 200)
rice <- CreateSeuratObject(counts = rice.data, project = "rice", min.cells = 3, min.features = 200)
#做了一些简单的过滤
arabidopsis[["percent.mg"]] <- PercentageFeatureSet(arabidopsis, pattern ="^ATMG")
arabidopsis[["percent.cg"]] <- PercentageFeatureSet(arabidopsis, pattern ="^ATCG")
arabidopsis <- subset(arabidopsis, subset = nFeature_RNA >500& nFeature_RNA <5000& percent.mg <1 & percent.cg <3)
rice<- subset(rice, subset = nFeature_RNA >500& nFeature_RNA <7000)
#获得拟南芥和水稻共有基因的
common_genes <- rownames(arabidopsis)[rownames(arabidopsis) %in% rownames(rice)]
> length(common_genes)
[1] 9042
#提取共有的部分进行对比
arabidopsis <- arabidopsis[rownames(arabidopsis) %in% common_genes,]
rice <- rice[rownames(rice) %in% common_genes,]
#下面的和上面一个帖子一样了
wb_list <- list()
wb_list[["arabidopsis"]] <- arabidopsis
wb_list[["rice"]] <- rice
wb_list <- lapply(X = wb_list, FUN = function(x) {
x <- NormalizeData(x)
x <- FindVariableFeatures(x, selection.method = "vst", nfeatures = 2000)
})
wb_anchors <- FindIntegrationAnchors(object.list = wb_list, dims = 1:30)
wb_seurat <- IntegrateData(anchorset = wb_anchors, dims = 1:30)
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)
p1<-DimPlot(wb_seurat,reduction = "umap") + scale_color_npg()+
plot_annotation(title = "arabidopsis and rice, before integration")
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)
p2<-DimPlot(wb_seurat, reduction = "umap") +
scale_color_npg()+
plot_annotation(title = "arabidopsis and rice, after integration")
p1|p2
#我们可以明显看出合并之前和合并之后的差异。但是我个人觉得貌似合并的效果也没有那么好，不确定是真实的细胞类型的差异，还是算法带来的误差差异。
#下面就是常规的聚类和可视化
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()
#我们也可以查看每个cluster下面，不同物种之间的差异。
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()
后面，我们也可以根据marker基因来注释每个cluster，同时查看不同细胞类型在物种间的差异，例如上面一个帖子讲的那样。
网友评论