-
Notifications
You must be signed in to change notification settings - Fork 1
/
figure3.Rmd
232 lines (179 loc) · 9.3 KB
/
figure3.Rmd
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
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
```{r}
library(Seurat)
library(ggplot2)
library(patchwork)
library(ggpubr)
library(dplyr)
library(Seurat)
library(patchwork)
library(sctransform)
library(ggplot2)
library(pheatmap)
library(Seurat)
library(ggplot2)
library(Matrix)
library(RColorBrewer)
library(dplyr)
library(scales)
library(data.table)
library(stats)
library("Nebulosa")
```
# Explore the mapping + ref_query results
```{r}
sample <- readRDS('/Users/jamrute/Documents/Graduate_School/Thesis_Lab/Lavine_Projects/Primary/COVID19_PBMC/azimuth/ref_query/cleaned_withPlatelets_query.rds')
```
# Plot heatmaps of the top Diff exp genes
```{r}
Idents(sample) <- "predicted.celltype.l2"
currCell <- subset(sample, idents = c("Plasmablast"))
Idents(currCell) <- "condition"
currCell_subset <- subset(currCell, idents = c("Control", "Alive_D7", "Dead_D7"))
currCell_subset$condition <- factor(currCell_subset$condition, levels = c("Control", "Alive_D7", "Dead_D7"), ordered = TRUE)
DGE <- read.csv2("/Users/jamrute/Documents/Graduate_School/Thesis_Lab/Lavine_Projects/Primary/COVID19_PBMC/azimuth/ref_query/single_cell_DGE/AliveD7_DeadD7/DE_Plasmablast.csv", header=TRUE, sep=',', row.names = 1)
DGE <- filter(DGE, as.numeric(DGE$p_val_adj) < 0.05 & as.numeric(DGE$avg_log2FC) > 0.5)
genes <- grep("RPS", DGE$gene, invert=TRUE, value = TRUE)
genes <- grep("RPL", genes, invert=TRUE, value = TRUE)
DefaultAssay(currCell_subset) <- "SCT"
currCell_subset.averageexpression <- AverageExpression(currCell_subset, features = genes, group.by = "condition", assays = "SCT")
currCell_subset.averageexpression <- as.matrix(currCell_subset.averageexpression)
```
```{r}
pheatmap(currCell_subset.averageexpression[[1]], scale="row", features = genes,col=colorspace::diverge_hsv(240), cexCol=0.5, cellwidth=8, cluster_rows=TRUE, fontsize_row=6, fontsize_col=6, cluster_cols = FALSE, legend = FALSE, cellheight = 8)
```
# B-cell UMAPs
```{r}
Idents(sample) <- "predicted.celltype.l2"
currCell <- subset(sample, idents = c("B naive", "B intermediate", "B memory", "Plasmablast"))
Idents(currCell) <- "condition"
currCell_subset <- subset(currCell, idents = c("Control", "Alive_D7", "Dead_D7"))
currCell_subset$condition <- factor(currCell_subset$condition, levels = c("Control", "Alive_D7", "Dead_D7"), ordered = TRUE)
```
```{r}
Idents(currCell_subset) <- "condition"
cd14_mono_control <- subset(currCell_subset, idents = c("Control"))
cd14_mono_day0 <- subset(currCell_subset, idents = c("Alive_D7"))
cd14_mono_day7 <- subset(currCell_subset, idents = c("Dead_D7"))
```
```{r}
DimPlot(cd14_mono_control, group.by = "predicted.celltype.l2") + xlim(-6, 6) + ylim(8, 15)
DimPlot(cd14_mono_day0, group.by = "predicted.celltype.l2") + xlim(-6, 6) + ylim(8, 15)
DimPlot(cd14_mono_day7, group.by = "predicted.celltype.l2") + xlim(-6, 6) + ylim(8, 15)
```
```{r}
DefaultAssay(currCell_subset) <- "predicted_ADT"
FeaturePlot(currCell_subset, "IgD") + scale_color_viridis(option = "C") + xlim(-6, 6) + ylim(7, 15)
FeaturePlot(currCell_subset, "CD73") + scale_color_viridis(option = "C") + xlim(-6, 6) + ylim(7, 15)
FeaturePlot(currCell_subset, "CD45RA") + scale_color_viridis(option = "C") + xlim(-6, 6) + ylim(7, 15)
FeaturePlot(currCell_subset, "CD71") + scale_color_viridis(option = "C") + xlim(-6, 6) + ylim(7, 15)
```
# Z-scores for key gene signatures
```{r}
Idents(sample) <- "predicted.celltype.l2"
currCell_subset <- subset(sample, idents = c("Plasmablast"))
Idents(currCell_subset) <- "condition"
currCell_subset <- subset(currCell_subset, idents = c("Control", "Alive_D7", "Dead_D7"))
DefaultAssay(currCell_subset) <- "SCT"
expdata <- GetAssayData(currCell_subset)
Pop1<-c("IFIT3","ISG15","IFI6","MX1","IFI44L","ISG20","IFITM1","IFI27")
pops<-list(Pop1)
#Z-Scores
z_scores<-NULL
for (i in 1:length(pops)) {
genes <- pops[[i]]
zz <- which(tolower(rownames(expdata)) %in% tolower(genes))
av <- numeric(ncol(expdata))
geneExp <- as.matrix(expdata[zz, ])
geneExp <- t(scale(t(geneExp)))
geneExp[is.nan(geneExp)] <- 0
z_scores <- rbind(z_scores,(av + colSums(geneExp) / length(zz)))
}
#write.table(z_scores, "./z_scores.tsv", sep="\t", quote=F, col.names=NA)
[email protected]$pop1_z<-z_scores[1,]
cd14_mono_AvgZ<-NULL
cd14_mono_AvgZ<-rbind(cd14_mono_AvgZ,(aggregate(currCell_subset$pop1_z, by=list([email protected]), mean)[,2]))
Idents(currCell_subset) <- "condition"
cd14_mono_control <- subset(currCell_subset, idents = c("Control"))
cd14_mono_day0 <- subset(currCell_subset, idents = c("Alive_D7"))
cd14_mono_day7 <- subset(currCell_subset, idents = c("Dead_D7"))
```
```{r}
FeaturePlot(object=cd14_mono_control, features = "pop1_z",pt.size=.5) + scale_color_gradientn(colors=c("blue","turquoise2","yellow","red","red4"), oob=scales::squish, limits=c(0,1)) + xlim(3, 6) + ylim(9, 12)
FeaturePlot(object=cd14_mono_day0, features = "pop1_z",pt.size=.5) + scale_color_gradientn(colors=c("blue","turquoise2","yellow","red","red4"), oob=scales::squish, limits=c(0,1)) + xlim(3, 6) + ylim(9, 12)
FeaturePlot(object=cd14_mono_day7, features = "pop1_z",pt.size=.5) + scale_color_gradientn(colors=c("blue","turquoise2","yellow","red","red4"), oob=scales::squish, limits=c(0,1)) + xlim(3, 6) + ylim(9, 12)
```
```{r}
Idents(currCell_subset) <- "predicted.celltype.l2"
currCell_subset2 <- subset(currCell_subset, idents = c("Plasmablast"))
Idents(currCell_subset2) <- "condition"
cd14_mono_control <- subset(currCell_subset2, idents = c("Control"))
cd14_mono_day0 <- subset(currCell_subset2, idents = c("Alive_D7"))
cd14_mono_day7 <- subset(currCell_subset2, idents = c("Dead_D7"))
write.csv(cd14_mono_control$pop1_z, file ="/Users/jamrute/Desktop/ISG_plasmablast_control.csv", quote = FALSE)
write.csv(cd14_mono_day0$pop1_z, file ="/Users/jamrute/Desktop/ISG_plasmablast_aliveD7.csv", quote = FALSE)
write.csv(cd14_mono_day7$pop1_z, file ="/Users/jamrute/Desktop/ISG_plasmablast_deceasedD7.csv", quote = FALSE)
```
# Get the percent Plasmablasts
```{r}
Idents(sample) <- "predicted.celltype.l2"
currCell <- subset(sample, idents = c("B naive", "B intermediate", "B memory", "Plasmablast"))
Idents(currCell) <- "condition"
currCell_subset <- subset(currCell, idents = c("Control", "Alive_D7", "Dead_D7"))
currCell_subset$condition <- factor(currCell_subset$condition, levels = c("Control", "Alive_D7", "Dead_D7"), ordered = TRUE)
```
```{r}
Idents(currCell_subset) <- "predicted.celltype.l2"
prop.table(table(Idents(currCell_subset), currCell_subset$diseaseStatus), margin = 2)*100
```
# IL6 Signaling
```{r}
genes <- c("CEBPB","CSNK2A1","ELK1","FOS","GRB2","HRAS","IL6","IL6R","IL6ST","JAK1","JAK2","JAK3","JUN","MAP2K1","MAPK3","PTPN11","RAF1","SHC1","SOS1","SRF","STAT3")
sample$condition <- factor(sample$condition, levels = c("Control", "Alive_D0", "Dead_D0", "Alive_D7", "Dead_D7"), ordered = TRUE)
```
```{r}
DefaultAssay(sample) <- "SCT"
sample.averageexpression <- AverageExpression(sample, features = genes, group.by = "condition", assays = "SCT")
sample.averageexpression <- as.matrix(sample.averageexpression)
pheatmap(sample.averageexpression[[1]], scale="row", features = genes,col=colorspace::diverge_hsv(240), cexCol=0.5, cellwidth=8, cluster_rows=TRUE, fontsize_row=6, fontsize_col=6, cluster_cols = FALSE, legend = FALSE, cellheight = 8)
```
```{r}
DefaultAssay(sample) <- "SCT"
expdata <- GetAssayData(sample)
Pop1<-genes
pops<-list(Pop1)
#Z-Scores
z_scores<-NULL
for (i in 1:length(pops)) {
genes <- pops[[i]]
zz <- which(tolower(rownames(expdata)) %in% tolower(genes))
av <- numeric(ncol(expdata))
geneExp <- as.matrix(expdata[zz, ])
geneExp <- t(scale(t(geneExp)))
geneExp[is.nan(geneExp)] <- 0
z_scores <- rbind(z_scores,(av + colSums(geneExp) / length(zz)))
}
#write.table(z_scores, "./z_scores.tsv", sep="\t", quote=F, col.names=NA)
[email protected]$pop1_z<-z_scores[1,]
cd14_mono_AvgZ<-NULL
cd14_mono_AvgZ<-rbind(cd14_mono_AvgZ,(aggregate(sample$pop1_z, by=list([email protected]), mean)[,2]))
Idents(sample) <- "condition"
control <- subset(sample, idents = c("Control"))
alive_day0 <- subset(sample, idents = c("Alive_D0"))
dead_day0 <- subset(sample, idents = c("Dead_D0"))
alive_day7 <- subset(sample, idents = c("Alive_D7"))
dead_day7 <- subset(sample, idents = c("Dead_D7"))
```
```{r}
FeaturePlot(object=control, features = "pop1_z") + scale_color_gradientn(colors=c("blue","turquoise2","yellow","red","red4"), oob=scales::squish, limits=c(0,0.5))
FeaturePlot(object=alive_day0, features = "pop1_z") + scale_color_gradientn(colors=c("blue","turquoise2","yellow","red","red4"), oob=scales::squish, limits=c(0,0.5))
FeaturePlot(object=dead_day0, features = "pop1_z") + scale_color_gradientn(colors=c("blue","turquoise2","yellow","red","red4"), oob=scales::squish, limits=c(0,0.5))
FeaturePlot(object=alive_day7, features = "pop1_z") + scale_color_gradientn(colors=c("blue","turquoise2","yellow","red","red4"), oob=scales::squish, limits=c(0,0.5))
FeaturePlot(object=dead_day7, features = "pop1_z") + scale_color_gradientn(colors=c("blue","turquoise2","yellow","red","red4"), oob=scales::squish, limits=c(0,0.5))
```
```{r}
write.csv(control$pop1_z, file ="/Users/jamrute/Desktop/control.csv", quote = FALSE)
write.csv(alive_day0$pop1_z, file ="/Users/jamrute/Desktop/aliveD0.csv", quote = FALSE)
write.csv(dead_day0$pop1_z, file ="/Users/jamrute/Desktop/deadD0.csv", quote = FALSE)
write.csv(alive_day7$pop1_z, file ="/Users/jamrute/Desktop/aliveD7.csv", quote = FALSE)
write.csv(dead_day7$pop1_z, file ="/Users/jamrute/Desktop/deadD7.csv", quote = FALSE)
```