Seurat结合Harmony的整合流程

参考：
单细胞测序分析: Seurat V3联合harmony进行单细胞数据整合分析
2021-07-21 harmony整合不同平台的单细胞数据

方法1

使用Harmony之前，需要构建一个Seurat对象， 进行一个标准的Seurat分析（包括PCA）。
Seurat中使用Harmony的流程与常规流程的区别是：可以所有细胞创建一个Seurat 对象，不需要为每个数据集创建一个Seurat 对象（Seurat list）。

R语言中%>%的含义是什么呢，管道函数啦，就是把左件的值发送给右件的表达式，并作为右件表达式函数的第一个参数。

library(devtools)
install_github("immunogenomics/harmony")
library(Seurat)
library(cowplot)
library(harmony)

#################    创建Seurat对象等一系列计算     ########################

pbmc <- CreateSeuratObject(counts = cbind(stim.sparse, ctrl.sparse), project = "PBMC", min.cells = 5) %>%
    Seurat::NormalizeData(verbose = FALSE) %>%
    FindVariableFeatures(selection.method = "vst", nfeatures = 2000) %>% 
    ScaleData(verbose = FALSE) %>% 
    RunPCA(pc.genes = pbmc@var.genes, npcs = 20, verbose = FALSE)

##1.  %>%管道函数的作用：将前一步的结果直接传参给下一步的函数，从而省略了中间的赋值步骤，可以大量减少内存中的对象，节省内存
##   例如：anscombe_tidy <- anscombe %>%mutate(observation = seq_len(n()))
##  等价于 anscombe_tidy=mutate(anscombe,observation = seq_len(n()))
##2.  上述函数进行了 创建Seurat对象 >-数据标准化>- 计算高变基因>- 数据缩放 >-PCA降维

##############  更改维度信息   #########################
#在object的metadata中定义细胞ID信息，变量名为stim.
pbmc@meta.data$stim <- c(rep("STIM", ncol(stim.sparse)), rep("CTRL", ncol(ctrl.sparse)))
##1.函数形式：rep(x, time = , length = , each = ,)
##    x：代表的是你要进行复制的对象，可以是一个向量或者是一个因子。
##    times：代表的是复制的次数，只能为正数。负数以及NA值都会为错误值。复制是指的是对整个向量进行复制。
##    each：代表的是对向量中的每个元素进行复制的次数。
##    length.out：代表的是最终输出向量的长度。 
##2.ncol(stim.sparse) 计算stim.sparse有多少列
##3.我们只需将两个数据集分开即可，所有将stim定义为两个细胞系的名字




#现在还没有使用Harmony矫正主成分分析结果的数据，数据集之间还是有很大的差异。
options(repr.plot.height = 5, repr.plot.width = 12)
##options为环境设置函数，修改绘图区域大小
p1 <- DimPlot(object = pbmc, reduction = "pca", pt.size = .1, group.by = "stim", do.return = TRUE)
##DimPlot绘制降维图。object为创建的   Seurat对象；reduction为降维方法，如果没有指定，首先搜索umap，然后是tsne，然后是pca；
##pt.size为调节绘图点的大小；group.by为分组依据。
p2 <- VlnPlot(object = pbmc, features = "PC_1", group.by = "stim", do.return = TRUE, pt.size = .1)
##绘制小提琴图。object为创建的 Seurat对象；features用于绘图的特征，可以是基因表达数，可以为PC得分。
plot_grid(p1,p2)
##将数个图拼接

###############   经运行Harmony校正后可视化结果   #################
#运行Harmony进行数据整合(矫正批次效应)
输入：使用Harmony，需要一个Seurat 对象和指定metadata信息中需要整合的变量名。
输出：返回一个Seurat对象，以及矫正之后的Harmony 信息
plot_convergenc参数设置为TRUE，保证Harmony 在运行中每一次迭代都在使矫正越累越好。

options(repr.plot.height = 2.5, repr.plot.width = 6)
pbmc <- pbmc %>% 
    RunHarmony("stim", plot_convergence = TRUE)
##函数公式：RunHarmony(object, group.by.vars, ...)
##object为管道对象，必须计算过PCA；

##获取Harmony 矫正之后的信息，使用Embeddings()函数
harmony_embeddings <- Embeddings(pbmc, 'harmony')
harmony_embeddings[1:5, 1:5]

##查看数据Harmony整合之后的前两个维度上数据是不是很好的整合，最好是很好的整合结果。
options(repr.plot.height = 5, repr.plot.width = 12)
p1 <- DimPlot(object = pbmc, reduction = "harmony", pt.size = .1, group.by = "stim", do.return = TRUE)
p2 <- VlnPlot(object = pbmc, features = "harmony_1", group.by = "stim", do.return = TRUE, pt.size = .1)
plot_grid(p1,p2)

#下游分析
下游分析都是基于Harmony矫正之后的PCA降维数据，不是基因表达数据和直接的PCA降维数据。设置reduction = 'harmony'，后续分析就会调用Harmony矫正之后的PCA降维数据。
要使用校正后的Harmony embeddings而不是PC进行后续分析，请设置reduction ='harmony'。

使用Harmony 矫正之后的数据，UMAP 和 Nearest Neighbor分析。
pbmc <- pbmc %>% 
    RunUMAP(reduction = "harmony", dims = 1:20) %>% 
    FindNeighbors(reduction = "harmony", dims = 1:20) %>% 
    FindClusters(resolution = 0.5) %>% 
    identity()

##UMAP 结果
在UMAP embedding中，我们可以看到更复杂的结构。由于我们使用harmony embeddings，因此UMAP embeddings混合得很好。
options(repr.plot.height = 4, repr.plot.width = 10)
DimPlot(pbmc, reduction = "umap", group.by = "stim", pt.size = .1, split.by = 'stim')

##聚类分析
options(repr.plot.height = 4, repr.plot.width = 6)
DimPlot(pbmc, reduction = "umap", label = TRUE, pt.size = .1)

方法2

参考：
2021-05-26 单细胞分析之harmony与Seurat
解读-Harmony原理介绍和官网教程

这个方法读取和前面的不一样
dir = c('Data/GSE139324) sample_name <- c('HNC01PBMC') scRNAlist

## 1、安装
library(harmony)
install_github("immunogenomics/harmony")
library(devtools)

## 2、创建Seurat对象
library(Seurat)
library(tidyverse)
library(patchwork)
rm(list=ls())
dir = c('Data/GSE139324/GSE139324_SUB/GSM4138110', 
        'Data/GSE139324/GSE139324_SUB/GSM4138111', 
        'Data/GSE139324/GSE139324_SUB/GSM4138128',
        'Data/GSE139324/GSE139324_SUB/GSM4138129',
        'Data/GSE139324/GSE139324_SUB/GSM4138148',
        'Data/GSE139324/GSE139324_SUB/GSM4138149',
        'Data/GSE139324/GSE139324_SUB/GSM4138162',
        'Data/GSE139324/GSE139324_SUB/GSM4138163',
        'Data/GSE139324/GSE139324_SUB/GSM4138168',
        'Data/GSE139324/GSE139324_SUB/GSM4138169')       
#GSE139324是存放数据的目录

sample_name <- c('HNC01PBMC', 'HNC01TIL', 'HNC10PBMC', 'HNC10TIL', 'HNC20PBMC',
                 'HNC20TIL',  'PBMC1', 'PBMC2', 'Tonsil1', 'Tonsil2')
scRNAlist <- list()
for(i in 1:length(dir)){
  counts <- Read10X(data.dir = dir[i])

#这里怎么还有去线粒体的步骤呢？
 scRNAlist[[i]] <- CreateSeuratObject(counts, project=sample_name[i], min.cells=3, min.features = 200)
  scRNAlist[[i]][["percent.mt"]] <- PercentageFeatureSet(scRNAlist[[i]], pattern = "^MT-")
  scRNAlist[[i]] <- subset(scRNAlist[[i]], subset = percent.mt < 10) 
}   
saveRDS(scRNAlist, "scRNAlist.rds")

##==harmony整合多样本==##
library(harmony)
scRNAlist <- readRDS("scRNAlist.rds")
##PCA降维
scRNA_harmony <- merge(scRNAlist[[1]], y=c(scRNAlist[[2]], scRNAlist[[3]], scRNAlist[[4]], scRNAlist[[5]], 
                                           scRNAlist[[6]], scRNAlist[[7]], scRNAlist[[8]], scRNAlist[[9]], scRNAlist[[10]]))
 scRNA_harmony <- NormalizeData(scRNA_harmony) %>% FindVariableFeatures() %>% ScaleData() %>% RunPCA(verbose=FALSE)
##整合
system.time({scRNA_harmony <- RunHarmony(scRNA_harmony, group.by.vars = "orig.ident")})
#   用户   系统   流逝 
# 34.308  0.024 34.324

#user  system elapsed 
#62.90    4.06   72.08 

#降维聚类
scRNA_harmony <- RunUMAP(scRNA_harmony, reduction = "harmony", dims = 1:30)
scRNA_harmony <- FindNeighbors(scRNA_harmony, reduction = "harmony", dims = 1:30) %>% FindClusters()
##作图
#group_by_cluster
plot1 = DimPlot(scRNA_harmony, reduction = "umap", label=T) 
#group_by_sample
plot2 = DimPlot(scRNA_harmony, reduction = "umap", group.by='orig.ident') 
#combinate
plotc <- plot1+plot2
plotc

自己做

参考：
2021-05-26 单细胞分析之harmony与Seurat
单细胞转录组高级分析一：多样本合并与批次校正
看了几个帖子都说是用Harmony分析的话，先用Seurat跑到PCA结束，然后再用Harmony分析。但是我自己比对了下这几个帖子的分析流程，感觉不太一样，就跟着帖子分析吧。最难的地方还是读取多样本合并Seurat这一步，就是这步很多教程说的不是很清楚，循坏看不太懂。弄好合并对象之后跟着做就行。

整体思路：

1.读取数据；
2.创建Seurat对象（这里看教程上面的代码应该是先对每个样本创建Seurat对象，然后再 merge合并成一个）
3.整合成一个之后就开始用教程的代码。代码如下：

library(harmony)
install_github("immunogenomics/harmony")
library(devtools)
library(Seurat)
library(tidyverse)
library(patchwork)
library(tidyverse)
rm(list=ls())

#1.数据准备，下载好的数据都是压缩文件。读取文件方式现在只会标准的三个文件的读取方式。
#下载好数据之后，可以用下面的代码整理成标准三文件格式，
#但是教程上面的还有改名这一步，他们的代码暂时看不懂，先用下面的吧。
#改名以后再学，或者后面修图的时候改下名字。
#压缩文件下载好之后，解压成文件夹，改这个就行"GSE139324_RAW" 


fs=list.files('D:/临时/GSE139324_RAW/','^GSM')
fs
library(tidyverse)
samples=str_split(fs,'_',simplify = T)[,1]

lapply(unique(samples),function(x){
  y=fs[grepl(x,fs)]
  folder=paste0("GSE139324_RAW/", str_split(y[1],'_',simplify = T)[,1])
  dir.create(folder,recursive = T)
  #为每个样本创建子文件夹
  file.rename(paste0("GSE139324_RAW/",y[1]),file.path(folder,"barcodes.tsv.gz"))
  #重命名文件，并移动到相应的子文件夹里
  file.rename(paste0("GSE139324_RAW/",y[2]),file.path(folder,"features.tsv.gz"))
  file.rename(paste0("GSE139324_RAW/",y[3]),file.path(folder,"matrix.mtx.gz"))
})


#2.数据读取有两种方法：一种是直接全部读入，创建对象；另一种方法是先对每个样本创建对象，再将所有对象合并为最终的对象。

library(Seurat)
samples=list.files("GSE139324_RAW/")
samples
dir <- file.path('./GSE139324_RAW',samples)
names(dir) <- samples

#这里用合并方法2
scRNAlist <- list()
for(i in 1:length(dir)){
  print(i)
  counts <- Read10X(data.dir = dir[i])
  scRNAlist[[i]] <- CreateSeuratObject(counts, min.cells=3, min.features = 200) # 这里记得改参数
}
scRNA2 <- merge(scRNAlist[[1]], y=c(scRNAlist[[2]], scRNAlist[[3]], 
                                    scRNAlist[[4]], scRNAlist[[5]], scRNAlist[[6]], scRNAlist[[7]], 
                                    scRNAlist[[8]], scRNAlist[[9]], scRNAlist[[10]]))
dim(scRNA2)   #查看基因数和细胞总数
table(scRNA2@meta.data$orig.ident)  #查看每个样本的细胞数

# 在参考的第二个示例中，在合并对象之后有一个Warning message，其他教程没看到，暂时不知道怎么处理 
Warning message:
  In CheckDuplicateCellNames(object.list = objects) :
  Some cell names are duplicated across objects provided. Renaming to enforce unique cell names.

#3. 开始跑流程
scRNA2 <- NormalizeData(scRNA2) %>% FindVariableFeatures() %>% ScaleData() %>% RunPCA(verbose=FALSE)


#开始整合
system.time({scRNA2 <- RunHarmony(scRNA2, group.by.vars = "orig.ident")})

#降维聚类
scRNA2 <- FindNeighbors(scRNA2, reduction = "harmony", dims = 1:15) %>% FindClusters(resolution = 0.5)


#降维可视化
scRNA2 <- RunUMAP(scRNA2, reduction = "harmony", dims = 1:15)

#画图看看整合效果
#group_by_cluster
plot1 = DimPlot(scRNA2, reduction = "umap", label=T) 
#group_by_sample
plot2 = DimPlot(scRNA2, reduction = "umap", group.by='orig.ident') 
#combinate
plotc <- plot1+plot2
plotc

后续如何分析，暂时不知道，看到其他教程再更新。