NNI_paper/08_deg_MAST.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
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
setwd("/scbio7/mengwei/NNI_brain/revision/03_Disease_DEG/")
library(Seurat)
library(pbmcapply)
library(stringr)
library(FastIntegration)
library(SeuratData)
library(SeuratDisk)
library(ggplot2)
library(dplyr)
library(ggpubr)
`%!in%` <- Negate('%in%')
options(bitmapType = "cairo")

rna.data = readRDS("../01_integration/rna_annotated.rds")

rna.sub = subset(rna.data, group %in% c("CTRL", "PDD"))
rna.sub = SetIdent(rna.sub, value = "cell.type")
a = data.frame(table(rna.sub$cell.type, rna.sub$group))
a$Var2 = as.character(a$Var2)
a$Var1 = as.character(a$Var1)

a = a[which(a$Freq > 30),]
a = a[which(a$Var1 %in% names(which(table(a$Var1) > 1))),]
a = as.character(unique(a$Var1))

DefaultAssay(rna.sub) = "RNA"
rna.sub = NormalizeData(rna.sub)

mm = pbmclapply(
  a, function(ct) {
    m = FindMarkers(rna.sub, ident.1 = "PDD", ident.2 = "CTRL", group.by = "group", logfc.threshold = 0.5, subset.ident = ct, test.use = "MAST", latent.vars = c("sex", "sample"))
    m$gene = rownames(m)
    m$disease = "PDD"
    m$cell_type = ct
    return(m)
  }, mc.cores = 20
)
mm = do.call(rbind, mm)
mm = mm[which(mm$p_val_adj < 0.05),]
write.table(mm, "DEG_PDD_MAST.txt", sep = "\t", quote = F, row.names = F, col.names = T)


##########################################################################################

rna.sub = subset(rna.data, group %in% c("CTRL", "DLB"))
rna.sub = SetIdent(rna.sub, value = "cell.type")
a = data.frame(table(rna.sub$cell.type, rna.sub$group))
a$Var2 = as.character(a$Var2)
a$Var1 = as.character(a$Var1)

a = a[which(a$Freq > 30),]
a = a[which(a$Var1 %in% names(which(table(a$Var1) > 1))),]
a = as.character(unique(a$Var1))

DefaultAssay(rna.sub) = "RNA"
rna.sub = NormalizeData(rna.sub)
mm = pbmclapply(
  a, function(ct) {
    m = FindMarkers(rna.sub, ident.1 = "DLB", ident.2 = "CTRL", group.by = "group", logfc.threshold = 0.5, subset.ident = ct, test.use = "MAST", latent.vars = c("sex", "sample"))
    m$gene = rownames(m)
    m$disease = "DLB"
    m$cell_type = ct
    return(m)
  }, mc.cores = 20
)
mm = do.call(rbind, mm)
mm = mm[which(mm$p_val_adj < 0.05),]
write.table(mm, "DEG_DLB_MAST.txt", sep = "\t", quote = F, row.names = F, col.names = T)

##########################################################################################

rna.sub = subset(rna.data, group %in% c("CTRL", "ADD"))
rna.sub = SetIdent(rna.sub, value = "cell.type")
a = data.frame(table(rna.sub$cell.type, rna.sub$group))
a$Var2 = as.character(a$Var2)
a$Var1 = as.character(a$Var1)

a = a[which(a$Freq > 30),]
a = a[which(a$Var1 %in% names(which(table(a$Var1) > 1))),]
a = as.character(unique(a$Var1))

DefaultAssay(rna.sub) = "RNA"
rna.sub = NormalizeData(rna.sub)
mm = pbmclapply(
  a, function(ct) {
    m = FindMarkers(rna.sub, ident.1 = "ADD", ident.2 = "CTRL", group.by = "group", logfc.threshold = 0.5, subset.ident = ct, test.use = "MAST", latent.vars = c("sex", "sample"))
    m$gene = rownames(m)
    m$disease = "ADD"
    m$cell_type = ct
    return(m)
  }, mc.cores = 20
)
mm = do.call(rbind, mm)
mm = mm[which(mm$p_val_adj < 0.05),]
table(mm$cell_type)
write.table(mm, "DEG_ADD_MAST.txt", sep = "\t", quote = F, row.names = F, col.names = T)