使用Harmony算法移除批次效应

安装补充包

读取数据

1
2
3
4
5
6
7
8
9
10
11
12
13
14
root_path = "."
f_read10x <- function(dirN, ...){
tp_samples <- list.files(file.path(root_path, dirN))
tp_dir <- file.path(root_path, dirN, tp_samples)
names(tp_dir) <- tp_samples
scRNA <- list()
for (lc_ba in names(tp_dir)){
counts <- Read10X(data.dir = tp_dir[[lc_ba]])
scRNA[[lc_ba]] <- CreateSeuratObject(counts, project = lc_ba, ...)
scRNA[[lc_ba]][["percent.mt"]] <- PercentageFeatureSet(scRNA[[lc_ba]], pattern = "^MT-")
scRNA[[lc_ba]][["percent.ERCC"]] <- PercentageFeatureSet(scRNA[[lc_ba]], pattern = "^ERCC-")
}
scRNA
}

质量控制

1
2
3
4
5
6
7
8
9
10
options(repr.plot.width=12, repr.plot.height=6)
options(ggrepel.max.overlaps = Inf)
library(ggplot2)
for (scRNAo in scRNA){
print(VlnPlot(scRNAo, features = c("nFeature_RNA", "nCount_RNA", "percent.mt"), ncol = 3) + labs(title = scRNAo@project.name))
}

scRNA$APAP1 <- subset(scRNA$APAP1, subset = nFeature_RNA > 200 & nFeature_RNA < 3000)

scRNA <- Reduce(function(...) merge(...), scRNA)

预处理函数

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
f_ScaleData_RunPCA <- function(scRNA){
scRNA <- FindVariableFeatures(scRNA, selection.method = "vst", nfeatures = 2000)
lc_all.genes <- rownames(scRNA)
scRNA <- ScaleData(scRNA, features = lc_all.genes)

scRNA <- RunPCA(scRNA, features = VariableFeatures(object = scRNA))
print(ElbowPlot(scRNA, ndims = 40))
scRNA
}
require(harmony)
f_RunHarmony <- function(scRNA, dims=1:30, batchN="batch"){
scRNA = scRNA %>% RunHarmony(batchN, plot_convergence = TRUE, max.iter.harmony = 30)
scRNA <- scRNA %>% RunUMAP(reduction = "harmony", dims = dims)
scRNA
}

分群函数

1
2
3
4
5
f_FindNeighbors <- function(scRNA, resolution = 0.5){
scRNA <- scRNA %>% FindNeighbors(reduction = "harmony") %>% FindClusters(resolution = resolution)
scRNA[[paste0('h_resolution_', resolution)]] <- Idents(scRNA)
scRNA
}

查看分群质量

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
f_metaG2G <- function(metaG, matrixN=F){
res <- list()
alltype <- unique(metaG[[1]])
for(type in alltype){
res[[type]] <- rownames(metaG)[metaG[[1]] == type]
if (matrixN){
res[[type]] <- gsub('-','.',res[[type]])
}
}
res
}

f_br_cluster_f <- function(sObject, lc_groupN){
lc_filter <- unlist(unique(sObject[[lc_groupN]]))
lc_filter <- lc_filter[!is.na(lc_filter)]
lc_filter
}

f_br_cluster <- function(sObject, lc_groupN, lc_labelN, lc_prop = F){
lc_g <- f_metaG2G(sObject[[lc_groupN]])
lc_l <- sObject[[lc_labelN]]
lc_l[[1]] <- as.character(lc_l[[1]])
res <- data.frame(row.names = f_br_cluster_f(lc_l, lc_labelN))
if(lc_prop){
for(Nm in names(lc_g)){
tmp <- prop.table(table(lc_l[lc_g[[Nm]],]))
res[[Nm]] <- tmp[rownames(res)]
}
}else{
for(Nm in names(lc_g)){
tmp <- table(lc_l[lc_g[[Nm]],])
res[[Nm]] <- tmp[rownames(res)]
}
}
res[is.na(res)] = 0
res
}

标注函数

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
f_add_annotation <- function(sObject, lc_ids,lc_metaName){
Idents(sObject) <- sObject[[lc_metaName]]
names(lc_ids) <- levels(sObject)
sObject <- RenameIdents(sObject, lc_ids)
sObject[[lc_metaName]] <- Idents(sObject)
sObject
}

f_mergeMeta_pre <- function(scRNA, metaN, baseMetaN){
scRNA[[metaN]] <- as.character(scRNA[[baseMetaN]][[1]])
scRNA
}
f_mergeMeta <- function(scRNA, metaN, metaData){
scRNA[[metaN]][rownames(metaData),] <- as.character(metaData[[1]])
scRNA
}
f_mergeMeta_end <- function(scRNA, metaN){
scRNA[[metaN]][[1]] <- as.factor(scRNA[[metaN]][[1]])
scRNA
}

使用Harmony算法移除批次效应
https://occdn.limour.top/2232.html
Author
Limour
Posted on
August 13, 2022
Licensed under