大家好，今天周六，就简单些，跑个harmony去除批次效应的全流程代码，以后再介绍其他整合单细胞的方法。
harmony官方教程：https://htmlpreview.github.io/？https://github.com/immunogenomics/harmony/blob/master/doc/quickstart.html
1 加载r包，我加载的是 自己数据
.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"))library(Seurat)load("~/silicosis/silicosis_cluster_merge.rds")head(All.merge@meta.data)
可以看到，组别信息放在了stim列
2 然后，先对seurat对象走一遍标准流程
subset_data= All.mergeoptions(warn = -1)subset_data = subset_data %>%Seurat::NormalizeData(verbose = FALSE) %>%FindVariableFeatures(selection.method = "vst", nfeatures = 2000) %>%ScaleData(verbose = FALSE,vars.to.regress = 'percent.mt') %>%RunPCA(npcs = 30, verbose = FALSE)ElbowPlot(subset_data, ndims = 30)VlnPlot(subset_data, features = c("nFeature_RNA", "nCount_RNA", "percent.mt"), ncol = 3)
3 接着，使用harmony去除批次效应
#2 去批次效应---------library('harmony')#subset_data$stim=subset_data$orig.identsubset_data <- subset_data %>% RunHarmony("stim", plot_convergence = TRUE)harmony_embeddings <- Embeddings(subset_data, 'harmony')#######################clusterdims = 1:30subset_data <- subset_data %>%RunUMAP(reduction = "harmony", dims = dims) %>%RunTSNE(reduction = "harmony", dims = dims) %>%FindNeighbors(reduction = "harmony", dims = dims)subset_data=FindClusters(subset_data)DimPlot(subset_data)DimPlot(subset_data,split.by = 'stim')head(subset_data@meta.data)
4.最后就是对每个亚群进行命名啦
如果你感兴趣，你可以看一下我之前的亚群命名直播
看完记得顺手点个“在看”哦！
