forked from bio-Pixel/panVC
-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathDataIntegration_step2_scRNA_ClusteredCelltype.R
55 lines (48 loc) · 1.97 KB
/
DataIntegration_step2_scRNA_ClusteredCelltype.R
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
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
#' @param dataset Input Seurat object for preclustered mono-cell type dataset
#' @param path Output path
ClusteredCelltype <- function(dataset, path) {
library(dplyr)
library(Seurat)
library(ggplot2)
library(doMC)
library(harmony)
library(RColorBrewer)
registerDoMC(50)
dir.create(paste0(path, gsub(".rds", "", x)))
a <- NormalizeData(dataset, normalization.method = "LogNormalize", scale.factor = 10000)
a <- FindVariableFeatures(a, selection.method = "vst", nfeatures = 2000)
batch <- c("percent.mt", "percent.rp", "S.Score", "G2M.Score")
if (length(unique([email protected]$Tech)) > 1) {
batch <- c(batch, "Tech")
}
if (length(unique([email protected]$SampleID)) > 1) {
batch <- c(batch, "SampleID")
}
a <- SCTransform(a, vars.to.regress = batch, verbose = T)
a <- RunPCA(a, features = VariableFeatures(object = a))
a <- RunUMAP(a, dims = 1:30, label = T)
source("scRNA_Harmony.R")
group.by.vars = c()
if (length(unique([email protected]$Tech)) > 1) {
group.by.vars <- c(group.by.vars, "Tech")
}
if (length(unique([email protected]$SampleID)) > 1) {
group.by.vars <- c(group.by.vars, "SampleID")
}
if(length(group.by.vars)>1){
a <- BatchRemove_Harmony(
seurat_object = a, cores = 20, group.by.vars = group.by.vars,
colors = c(brewer.pal(8, "Set1"), colorRampPalette(brewer.pal(8, "Accent"))(length(unique([email protected]$SampleID)))), max.dim = 30,
out_dir = paste0(path, gsub(".rds", "", x))
)
}
# Calculation of marker genes
source("Seurat_marker.R")
if (length(unique([email protected]$Tech)) == 1) {
logfc.threshold <- ifelse(unique([email protected]$Tech) == "smart-seq2", 1, 0.25)
} else {
logfc.threshold <- 0.25
}
a <- Seurat_marker(seurat_object = a, FindNeighbor_reduction = "harmony", resolution = 2, cores = 10, top_n = 50, out_dir = paste0(path, gsub(".rds", "", x)))
return(a)
})