-
Notifications
You must be signed in to change notification settings - Fork 2
/
Copy pathFigure 2 and S3 code - script.R
190 lines (144 loc) · 7.94 KB
/
Figure 2 and S3 code - script.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
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
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
library(dplyr)
library(ggplot2)
library(ggpattern)
library(reshape2)
library(cowplot)
library(Seurat)
library(dplyr)
library(RColorBrewer)
load("F:/4-7-20m-Final-1030/Figure 2+S2/Add_ADvsPso_New_1102.Rda")
####################### Figure 2A and B
rash_lfc = as.data.frame(foldchange[,6])
colnames(rash_lfc) = "Rash vs N"
rash_p = as.matrix(pvals_prop[,6])
vis_coef <- rash_lfc
pvals_adj <- p.adjust(rash_p, method = "BH")
vis_pvals_adj <- abs(log10(pvals_adj+1e-10))
#knitr::kable( vis_coef )
vis_coef <- as.data.frame(vis_coef)
vis_coef$CellType = rownames(vis_coef)
vis_coef$CellType <- factor(vis_coef$CellType, levels = cell_identity$Identities)
#vis_coef$CellType <- factor(vis_coef$CellType, levels = 1:41)
vis_coef = melt(vis_coef)
colnames(vis_coef)[2:3] = c("Comparison", "coefs")
vis_coef$IsSig = abs(vis_pvals_adj) >= abs(log10(0.05+1e-10))
new_meta = meta %>% filter(dis %in% c("AD1","AE", "Pso", "LP", "BP", "NML")) %>%
group_by(new_level, identity) %>% summarise(n = n()) %>%
mutate(freq = n / sum(n))
new_meta$new_level2 <- ""
new_meta$new_level2 <- ifelse(new_meta$new_level=="NML",
"NML(7)", "Rash(24)")
vis_coef$p_val <- ""
vis_coef$p_val <- ifelse(vis_coef$IsSig=="TRUE",
"P-value < 0.05", "P-value> 0.05")
################################### save for tables ##################################
write.xlsx(vis_coef, "F:/4-7-20m-Final-1030/Figure 2+S2/Final/vis_coef_Rash_N.xlsx")
write.xlsx(rash_lfc, "F:/4-7-20m-Final-1030/Figure 2+S2/Final/rash_lfc_Rash_N.xlsx")
write.xlsx(as.matrix(new_meta), "F:/4-7-20m-Final-1030/Figure 2+S2/Final/new_meta_Rash_N.xlsx")
write.xlsx(as.matrix(pvals_adj), "F:/4-7-20m-Final-1030/Figure 2+S2/Final/p_value_adj_Rash_N.xlsx")
########################################################################################
p1 = ggplot(data = new_meta, aes(x = identity, y = freq,
group = new_level, fill = new_level2)) +
geom_bar(stat= "identity", position = "dodge")+
theme_classic()+ylab("proportion") +
theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)) +
labs(fill = "") + xlab("CD45+ Cell Types")+ylab("Proportion of CD45+ Cells")+
scale_fill_manual(values = colorRampPalette(brewer.pal(12, "Paired"))(2))
p2 = ggplot(vis_coef, aes(x = CellType, y = coefs, fill = p_val)) +
geom_bar(stat = "identity") +
scale_fill_manual(values = c("red4", "gray")) +
theme_classic() +
labs(fill = "") + xlab("CD45+ Cell Types")+ylab("Log2FC Rash vs NML")+
theme(axis.text.x=element_text(angle=90,hjust=1,vjust=0.5))
pp <- list(p1, p2)
plot_grid(plotlist=pp, ncol=1, align='v')
####################### Figure 2C
our <- readRDS("F:/4-7-20m-Final-1030/our_no_218_203.rds")
Idents(our) <- our$ID
our <- RenameIdents(our, `1` = "Tcm", `2` = "Trm1", `3` = "eTreg1", `4` = "Trm2", `5` = "CTLex", `6` = "CTLem",
`7` = "Tet", `8` = "Tmm1", `9` = "ILC/NK", `10` = "NK", `11` = "Trm3", `12` = "cmTreg",
`13` = "Tmm2", `14` = "eTreg2", `15` = "Tmm3", `16` = "ILC2", `17` = "Tn",
`18` = "moDC1", `19` = "LC1", `20` = "InfMono", `21` = "Mac1", `22` = "Mono",
`23` = "LC2", `24` = "moDC2", `25` = "moDC3", `26` = "DC1", `27` = "DC2",
`28` = "LC3", `29` = "Mac2", `30` = "Plasma", `31` = "B", `32` = "Mac3",
`33` = "DC3", `34` = "migDC", `35` = "Mac4", `36` = "Mast", `37` = "Trm-c",
`38` = "Treg-c", `39` = "Mast-c", `40` = "ILC/NK-c", `41` = "CTL-c")
our[["Ident2"]] <- Idents(object = our)
Idents(our) <- our$ID
plot_df <- as.data.frame.matrix([email protected])
plot_df$sample_name <- plot_df$donor
test_df <- plot_df %>% group_by(donor, Ident2, dis) %>%
summarise(n = n())
test_df <- test_df %>% group_by(donor, dis) %>% mutate(Freq = n / sum(n))
test_df$sample_name_condition <- paste0(test_df$donor,' ', test_df$dis)
library(RColorBrewer)
test_df$dis <- factor(test_df$dis, levels = c("NML", "AD1", "Pso", "AE", "LP", "BP"))
#test_df <- subset(test_df, dis!= "AE")
test_df$Ident2 <- factor(test_df$Ident2, levels = c("Tcm", "Trm1", "eTreg1", "Trm2", "CTLex", "CTLem", "Tet",
"Tmm1", "ILC/NK", "NK", "Trm3", "cmTreg", "Tmm2", "eTreg2",
"Tmm3", "ILC2", "Tn", "moDC1", "LC1", "InfMono", "Mac1",
"Mono", "LC2", "moDC2", "moDC3", "DC1", "DC2", "LC3",
"Mac2", "Plasma", "B", "Mac3", "DC3", "migDC", "Mac4",
"Mast", "Trm-c", "Treg-c", "Mast-c", "ILC/NK-c", "CTL-c"))
write.xlsx(as.matrix(test_df), "F:/4-7-20m-Final-1030/Figure 2+S2/Final/Figure2C.xlsx")
ggplot(data=test_df, aes(x=sample_name_condition, y = Freq, fill=Ident2)) +
geom_bar(stat="identity") +
facet_grid(cols = vars(dis),
labeller = label_wrap_gen(width = 16,multi_line = TRUE),
scales="free", space = "free") +
xlab("Individual Sample") + ylab("Proportion of CD45+ Cells for each sample") +
scale_fill_manual(values = colorRampPalette(brewer.pal(12, "Paired"))(41)) +
theme_bw()+
theme(legend.position="bottom", axis.text.x = element_blank(), axis.ticks.x = element_blank())+
guides(fill=guide_legend(nrow=5,byrow=TRUE))
######################################### S2AB ###############################################
dis_lfc = as.data.frame(foldchange[,c(1,3,7)])
colnames(dis_lfc) = colnames(foldchange)[c(1,3,7)]
dis_p = as.matrix(pvals_prop[,c(1,3, 7)])
vis_coef <- dis_lfc
pvals_adj <- p.adjust(dis_p, method = "BH")
vis_pvals_adj <- abs(log10(pvals_adj+1e-10))
#knitr::kable( vis_cp.adjustoef )
vis_coef <- as.data.frame(vis_coef)
vis_coef$CellType = rownames(vis_coef)
vis_coef$CellType <- factor(vis_coef$CellType, levels = cell_identity$Identities)
#vis_coef$CellType <- factor(vis_coef$CellType, levels = 1:41)
vis_coef = melt(vis_coef)
colnames(vis_coef)[2:3] = c("Comparison", "coefs")
vis_coef$IsSig = abs(vis_pvals_adj) >= abs(log10(0.05+1e-10))
new_meta = meta %>% filter(dis %in% c("AD1","Pso", "NML")) %>%
group_by(dis, identity) %>% summarise(n = n()) %>%
mutate(freq = n / sum(n))
new_meta$dis = factor(new_meta$dis, levels = c("AD1","Pso", "NML"))
p1 = ggplot(data = new_meta, aes(x = identity, y = freq, group = dis, fill = dis)) +
geom_bar(stat = "identity",position = "dodge")+
theme_classic()+ylab("Proportion of CD45+ Cells") + xlab(" ") +
theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))
p2 = ggplot(vis_coef, aes(x = CellType, y = coefs, fill = Comparison, col = IsSig)) +
geom_col_pattern(
pattern = ifelse(vis_coef$IsSig, "none", "stripe"),
pattern_density = 0.35,
pattern_colour = 'black',
position = "dodge"
) +
theme_classic() +
theme(axis.text.x=element_text(angle=90,hjust=1,vjust=0.5)) +
ylab("Log2FC Rash vs NML of log(P)") +
xlab("")+
guides(col = guide_legend(override.aes =
list(
pattern = c("none", "stripe")
)
)) +
guides(fill = guide_legend(override.aes =
list(
pattern = c("none", "none", "none")
)
))
pp <- list(p1, p2)
plot_grid(plotlist=pp, ncol=1, align='v')
write.xlsx(vis_coef, "F:/4-7-20m-Final-1030/Figure 2+S2/Final/vis_coef_AD_Pso_N.xlsx")
write.xlsx(dis_lfc, "F:/4-7-20m-Final-1030/Figure 2+S2/Final/dis_lfc_AD_Pso_N.xlsx")
write.xlsx(as.matrix(new_meta), "F:/4-7-20m-Final-1030/Figure 2+S2/Final/new_meta_AD_Pso_N.xlsx")
write.xlsx(dis_p, "F:/4-7-20m-Final-1030/Figure 2+S2/Final/dis_p_AD_Pso_N.xlsx")
write.xlsx(as.matrix(pvals_adj), "F:/4-7-20m-Final-1030/Figure 2+S2/Final/pvals_adj_AD_Pso_N.xlsx")