Comparison of Harmony Algorithm

Comparison of Harmony Algorithm

Today I want to compare four types of harmony algorithm

Type1(harmony-pytorch)

from harmony import harmonize
import scanpy.external as sce
import scanpy as sc

adata=sc.read("/home/yxkang/test_methods/dataset/bct/bct_raw.h5ad")
print(adata)
# adata = adata[adata.obs["BATCH"].isin(["vis","spk"])] # vis spk and wal
# n_sample=1000
# adata = sc.pp.subsample(adata,n_obs=n_sample,copy=True)
sc.pp.filter_cells(adata,min_genes=200)
sc.pp.filter_genes(adata,min_cells=1)
sc.pp.normalize_total(adata, target_sum=1e4)
sc.pp.log1p(adata)
sc.pp.highly_variable_genes(adata, n_top_genes=2000, subset = True)#seurat_v3 neeed count based 
sc.pp.scale(adata,max_value=10.0)
sc.tl.pca(adata)

Z = harmonize(adata.obsm['X_pca'], adata.obs, batch_key ="BATCH")

adata.obsm['X_harmony'] = Z

sc.pp.neighbors(adata,use_rep="X_harmony")
#sc.tl.louvain(adata,resolution=3.0)
sc.tl.umap(adata)
sc.pl.umap(adata,color=['BATCH','celltype'])

the results are as follows:

Please pay attention to the cost time( 7 min 28s). Although the integrated result is normal, the running time is abnormal. I don’t know why.

Type2(harmony python)

import harmonypy as hm
import scanpy as sc 

adata=sc.read("/home/yxkang/test_methods/dataset/bct/bct_raw.h5ad")
print(adata)
# adata = adata[adata.obs["BATCH"].isin(["vis","spk"])] # vis spk and wal
# n_sample=1000
# adata = sc.pp.subsample(adata,n_obs=n_sample,copy=True)
sc.pp.filter_cells(adata,min_genes=200)
sc.pp.filter_genes(adata,min_cells=1)
sc.pp.normalize_total(adata, target_sum=1e4)
sc.pp.log1p(adata)
sc.pp.highly_variable_genes(adata, n_top_genes=2000, subset = True)#seurat_v3 neeed count based 
sc.pp.scale(adata,max_value=10.0)
sc.tl.pca(adata)

harmony_object = hm.run_harmony(adata.obsm['X_pca'], adata.obs, "BATCH")

# Store Harmony embeddings
adata.obsm['X_harmony'] = harmony_object.Z_corr.T

sc.pp.neighbors(adata,use_rep="X_harmony")
#sc.tl.louvain(adata,resolution=3.0)
sc.tl.umap(adata)
sc.pl.umap(adata,color=['BATCH','celltype'])

the result is as follows:

This type of use runs fast

Type3(scanpy external)

import scanpy.external as sce
import scanpy as sc

adata=sc.read("/home/yxkang/test_methods/dataset/bct/bct_raw.h5ad")
print(adata)
# adata = adata[adata.obs["BATCH"].isin(["vis","spk"])] # vis spk and wal
# n_sample=1000
# adata = sc.pp.subsample(adata,n_obs=n_sample,copy=True)
sc.pp.filter_cells(adata,min_genes=200)
sc.pp.filter_genes(adata,min_cells=1)
sc.pp.normalize_total(adata, target_sum=1e4)
sc.pp.log1p(adata)
sc.pp.highly_variable_genes(adata, n_top_genes=2000, subset = True)#seurat_v3 neeed count based 
sc.pp.scale(adata,max_value=10.0)
sc.tl.pca(adata)

sce.pp.harmony_integrate(adata, 'BATCH')
sc.pp.neighbors(adata,use_rep="X_pca_harmony")
sc.tl.umap(adata)
sc.pl.umap(adata,color=['BATCH','celltype'])

The result is as follows:

Comparison of Harmony Algorithm

Today I want to compare four types of harmony algorithm

Type1(harmony-pytorch)

from harmony import harmonize
import scanpy.external as sce
import scanpy as sc

adata=sc.read("/home/yxkang/test_methods/dataset/bct/bct_raw.h5ad")
print(adata)
# adata = adata[adata.obs["BATCH"].isin(["vis","spk"])] # vis spk and wal
# n_sample=1000
# adata = sc.pp.subsample(adata,n_obs=n_sample,copy=True)
sc.pp.filter_cells(adata,min_genes=200)
sc.pp.filter_genes(adata,min_cells=1)
sc.pp.normalize_total(adata, target_sum=1e4)
sc.pp.log1p(adata)
sc.pp.highly_variable_genes(adata, n_top_genes=2000, subset = True)#seurat_v3 neeed count based 
sc.pp.scale(adata,max_value=10.0)
sc.tl.pca(adata)

Z = harmonize(adata.obsm['X_pca'], adata.obs, batch_key ="BATCH")

adata.obsm['X_harmony'] = Z

sc.pp.neighbors(adata,use_rep="X_harmony")
#sc.tl.louvain(adata,resolution=3.0)
sc.tl.umap(adata)
sc.pl.umap(adata,color=['BATCH','celltype'])

the results are as follows:

Please pay attention to the cost time( 7 min 28s). Although the integrated result is normal, the running time is abnormal. I don’t know why.

Type2(harmony python)

import harmonypy as hm
import scanpy as sc 

adata=sc.read("/home/yxkang/test_methods/dataset/bct/bct_raw.h5ad")
print(adata)
# adata = adata[adata.obs["BATCH"].isin(["vis","spk"])] # vis spk and wal
# n_sample=1000
# adata = sc.pp.subsample(adata,n_obs=n_sample,copy=True)
sc.pp.filter_cells(adata,min_genes=200)
sc.pp.filter_genes(adata,min_cells=1)
sc.pp.normalize_total(adata, target_sum=1e4)
sc.pp.log1p(adata)
sc.pp.highly_variable_genes(adata, n_top_genes=2000, subset = True)#seurat_v3 neeed count based 
sc.pp.scale(adata,max_value=10.0)
sc.tl.pca(adata)

harmony_object = hm.run_harmony(adata.obsm['X_pca'], adata.obs, "BATCH")

# Store Harmony embeddings
adata.obsm['X_harmony'] = harmony_object.Z_corr.T

sc.pp.neighbors(adata,use_rep="X_harmony")
#sc.tl.louvain(adata,resolution=3.0)
sc.tl.umap(adata)
sc.pl.umap(adata,color=['BATCH','celltype'])

the result is as follows:

This type of use runs fast

Type 4(R harmony)

rm(list=ls())
method="harmony"
suppressPackageStartupMessages({
  library(Seurat)
  library(harmony)
  library(SingleCellExperiment)
  library(ggplot2)
})

start_time <- Sys.time()
####################### parameter setting ###############################
setwd("/Users/yxkang/Desktop/Medium/Harmony/")
data <- readRDS("/Users/yxkang/Desktop/dataset/bct/bct_raw.rds")
####################### parameter setting ###############################

print(table(colData(data)$BATCH , colData(data)$celltype))
print("=====================================")
print(data)
data_seurat=CreateSeuratObject(counts = counts(data),meta.data = as.data.frame(colData(data)))
data_seurat <- NormalizeData(data_seurat, verbose = FALSE)
data_seurat <- FindVariableFeatures(data_seurat, selection.method = "vst", nfeatures = 2000, verbose = FALSE)

# Run the standard workflow for visualization and clustering
data_seurat <- ScaleData(data_seurat, verbose = FALSE)
data_seurat <- RunPCA(data_seurat, npcs = 30, verbose = F)
#data_seurat <- RunUMAP(data_seurat, reduction = "pca", dims = 1:30, verbose = F)
#DimPlot(data_seurat,reduction = "umap",group.by = "BATCH") + plot_annotation(title = "data before integration")

data_seurat <- data_seurat %>% RunHarmony("BATCH", plot_convergence = F,max.iter.harmony=50)
#data_seurat <- RunTSNE(data_seurat, reduction = "harmony", dims = 1:30, verbose = F)
data_seurat <- RunUMAP(data_seurat, reduction = "harmony", dims = 1:30, verbose = F)
# data_seurat <- FindNeighbors(data_seurat, reduction = "harmony", dims = 1:30,verbose=FALSE)
# data_seurat <- FindClusters(data_seurat,verbose=FALSE,resolution = 0.4)

# p1=DimPlot(data_seurat, reduction = "tsne", group.by = "BATCH", label.size = 10)+ggtitle("Integrated Batch")
# p2=DimPlot(data_seurat, reduction = "tsne", group.by = "celltype",label.size = 10)+ggtitle("Integrated Celltype")
# p= p1 + p2 
# print(p)
# print("harmony tsne done")
# ggsave("harmony_pancreas_tsne.png",p)

#saveRDS(data_seurat,file="harmony.rds")

p3=DimPlot(data_seurat, reduction = "umap", group.by = "BATCH", label.size = 10)+ggtitle("Integrated Batch")
p4=DimPlot(data_seurat, reduction = "umap", group.by = "celltype",label.size = 10)+ggtitle("Integrated Celltype")
p= p3 + p4
print(p)
print("harmony umap done")
ggsave("harmony_bct_umap.png",p)
#saveRDS(data_seurat,file="harmony.rds")
end_time <- Sys.time()
# Calculate execution time
execution_time <- end_time - start_time
sprintf("runing harmony cost: %ds", round(execution_time))

the result is as follows:

it cost only 16s

Type5(R harmony standardalone)

rm(list=ls())
method="harmony"
suppressPackageStartupMessages({
  library(Seurat)
  library(harmony)
  library(SingleCellExperiment)
  library(ggplot2)
})

####################### parameter setting ###############################
#setwd("/Users/yxkang/Desktop/Medium/Harmony/")
data <- readRDS("/project/MultiSampleIstar/dataset/bct/bct_raw.rds")
####################### parameter setting ###############################

print(table(colData(data)$BATCH , colData(data)$celltype))
print("=====================================")
print(data)
data_seurat=CreateSeuratObject(counts = counts(data),meta.data = as.data.frame(colData(data)))
data_seurat <- NormalizeData(data_seurat, verbose = FALSE)
data_seurat <- FindVariableFeatures(data_seurat, selection.method = "vst", nfeatures = 2000, verbose = FALSE)

# Run the standard workflow for visualization and clustering
data_seurat <- ScaleData(data_seurat, verbose = FALSE)
data_seurat <- RunPCA(data_seurat, npcs = 30, verbose = F)
#data_seurat <- RunUMAP(data_seurat, reduction = "pca", dims = 1:30, verbose = F)
#DimPlot(data_seurat,reduction = "umap",group.by = "BATCH") + plot_annotation(title = "data before integration")

############################################################################################
############################################################################################
############################################################################################
############################################################################################
start_time <- Sys.time()

pca_mat = data_seurat@reductions$pca@cell.embeddings
meta_data = data_seurat@meta.data

saveRDS(pca_mat,file="./pca_mat.rds")
saveRDS(meta_data,file="./meta_data.rds")

harmony_mat = RunHarmony(pca_mat,meta_data, "BATCH",plot_convergence = F,max.iter.harmony=50)
data_seurat[["harmony"]] <- CreateDimReducObject(
  embeddings = harmony_mat,        # N_cells × N_dim  matrix
  key = "harmony_",                # prefix for each dimension
  assay = DefaultAssay(data_seurat) # \
)

end_time <- Sys.time()
# Calculate execution time
execution_time <- end_time - start_time
sprintf("runing harmony cost: %ds", round(execution_time))
############################################################################################
############################################################################################
############################################################################################
############################################################################################

#data_seurat <- RunTSNE(data_seurat, reduction = "harmony", dims = 1:30, verbose = F)
data_seurat <- RunUMAP(data_seurat, reduction = "harmony", dims = 1:30, verbose = F)
# data_seurat <- FindNeighbors(data_seurat, reduction = "harmony", dims = 1:30,verbose=FALSE)
# data_seurat <- FindClusters(data_seurat,verbose=FALSE,resolution = 0.4)

# p1=DimPlot(data_seurat, reduction = "tsne", group.by = "BATCH", label.size = 10)+ggtitle("Integrated Batch")
# p2=DimPlot(data_seurat, reduction = "tsne", group.by = "celltype",label.size = 10)+ggtitle("Integrated Celltype")
# p= p1 + p2 
# print(p)
# print("harmony tsne done")
# ggsave("harmony_pancreas_tsne.png",p)

#saveRDS(data_seurat,file="harmony.rds")

p3=DimPlot(data_seurat, reduction = "umap", group.by = "BATCH", label.size = 10)+ggtitle("Integrated Batch")
p4=DimPlot(data_seurat, reduction = "umap", group.by = "celltype",label.size = 10)+ggtitle("Integrated Celltype")
p= p3 + p4
print(p)
print("harmony umap done")
ggsave("harmony_bct_umap.png",p)
#saveRDS(data_seurat,file="harmony.rds")