Seurat v3的标准整合流程会返回一个新的assay存储了数据调整的expr counts <- scRNA[["integrated"]]@data；相比之下，Harmony整合方法更像是一种新的聚类降维方式，基于原矩阵的PCA数据，得到一个新的harmony降维矩阵，[ 2 dimensional reductions calculated: pca, harmony]，接下来的流程则基于用harmony 降维数据（代替默认的pca降维矩阵）来进行聚类分群。
options(warn = -1)
library(Seurat)
library(dplyr)
library(ggplot2)
library(BuenColors)
library(ggtree)
######################################################################################
#Seurat v3的标准整合流程，可选是否对样本间去除批次效应
stdIntg <- function(project,group = "orig.ident"){
project.list <- SplitObject(project, split.by = group)
for (i in 1:length(project.list)) {
project.list[[i]] <- NormalizeData(project.list[[i]], normalization.method = "LogNormalize", scale.factor = 10000)
project.list[[i]] <- FindVariableFeatures(project.list[[i]], selection.method = "vst",nfeatures = 3000)
}
features <- SelectIntegrationFeatures(object.list = project.list,nfeatures = 3000)
anchors <- FindIntegrationAnchors(object.list = project.list,anchor.features = features)
integrated <- IntegrateData(anchorset = anchors,dims = 1:20)
DefaultAssay(integrated) <- "integrated"
return(integrated)
}
######################################################################################
set.seed(123)
setwd("~")
scRNA <- readRDS("*.rds")
sub_Tmp <- subset(scRNA,subset = Majortype == "Excitatory neuron")
sub_Tmp <- subset(sub_Tmp,downsample=100) #减小数据集大小对‘seurat_clusters’默认分类变量抽样100个细胞进行测试
rm(scRNA); gc()
##对数据进行整合
if(T){
sub_Tmp <- stdIntg(sub_Tmp)
sub_Tmp <- ScaleData(sub_Tmp, features = rownames(sub_Tmp))
sub_Tmp <- RunPCA(sub_Tmp, features = VariableFeatures(sub_Tmp), npcs = 40, verbose = F)
sub_Tmp <- FindNeighbors(sub_Tmp, dims = 1:20)
sub_Tmp <- FindClusters(sub_Tmp, resolution = 1)
sub_Tmp <- RunUMAP(sub_Tmp, dims = 1:20)
}
#层次聚类
Idents(sub_Tmp) <- sub_Tmp$seurat_clusters #设置聚类层级
top50 <- FindAllMarkers(object = sub_Tmp, only.pos = T, min.pct = 0.25, logfc.threshold = 0.25) %>% group_by(cluster) %>% top_n(n = 50, wt = avg_log2FC )
sub_Tmp<-BuildClusterTree(sub_Tmp,features = top50$gene) #对top50的marker基因进行层次聚类
Tool(object = sub_Tmp, slot = 'BuildClusterTree')
tree <- sub_Tmp@tools$BuildClusterTree
tree$tip.label <- paste0("ex_", tree$tip.label)
View(fortify(tree))
##画出层次聚类图
ggtree::ggtree(tree, aes(x, y),size=0.8) +
ggtree::geom_tree() +
ggtree::theme_tree() +
ggtree::geom_tiplab(vjust= 2,hjust = 0.5,align = T) +
ggtree::geom_tippoint(color = jdb_color_maps[1:length(tree$tip.label)], shape = 15, size = 10) +
coord_flip() +scale_x_reverse()
