Harmony整合单细胞数
2022-04-12 本文已影响0人
oceanandshore
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
后续如何分析，暂时不知道，看到其他教程再更新。
