NNI_paper/03_integration.R at main · JinmiaoChenLab/NNI_paper · GitHub

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
55
56
57
setwd("/scbio7/mengwei/NNI_brain/revision/01_integration/")
library(Seurat)
library(pbmcapply)
library(stringr)
library(FastIntegration)
library(SeuratData)
library(SeuratDisk)
options(bitmapType = "cairo")
options(future.globals.maxSize = 500000 * 1024^2)
library(future)
library(harmony)
plan("multisession", workers = 30)


rna.list = readRDS("../../data/rna_list.rds")

rna.list = pbmclapply(
  1:length(rna.list), function(i) {
    return(SCTransform(rna.list[[i]], verbose = T))
  }, mc.cores = 20
)

saveRDS(rna.list, "rna_list.rds", compress = F)
rna.list = readRDS("rna_list.rds")
features = SelectIntegrationFeatures(rna.list)
rna.list = PrepSCTIntegration(object.list = rna.list, anchor.features = features)
rna.list = lapply(X = rna.list, FUN = RunPCA, features = features)

anchors = FindIntegrationAnchors(object.list = rna.list, normalization.method = "SCT",
                                 anchor.features = features, dims = 1:30, reduction = "rpca", k.anchor = 20)
saveRDS(anchors, "anchors.rds", compress = F)

anchors = readRDS("anchors.rds")
sct = IntegrateData(anchorset = anchors, normalization.method = "SCT", dims = 1:30)
saveRDS(sct, "rpca.rds", compress = F)

rna.list = merge(rna.list[[1]], rna.list[2:20])

rna.list = RunPCA(rna.list, features = features)
rna.list = RunUMAP(rna.list, dims = 1:30)

rna.list = RunHarmony(rna.list, group.by.vars = "sample", reduction = "pca", verbose = T)
rna.list = RunUMAP(rna.list, dims = 1:50, reduction = "harmony")

saveRDS(rna.list, "harmony.rds", compress = F)

rna.data = readRDS("rpca.rds")
rna.data = RunPCA(rna.data)
rna.data = RunUMAP(rna.data, reduction = "pca", dims = 1:30)
saveRDS(rna.data, "rpca.rds", compress = F)

rna = readRDS("/scbio7/mengwei/NNI_brain/data/rna_annotated.rds")

rna.data@meta.data[Cells(rna), "cell.type"] = rna$cell.type


DimPlot(rna.data, group.by = "cell.type")