-
Notifications
You must be signed in to change notification settings - Fork 2
/
Copy pathprocess_TNK_final.R
executable file
·328 lines (240 loc) · 18.6 KB
/
process_TNK_final.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
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
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
## Author: Arja Ray
library(assertthat) # For sanity checks
library(ggplot2) # For pretty plots
library(cowplot) # For pretty plots
library(dittoSeq) # For pretty, colorblind friendly plots
library(dplyr) # For inline modification of matrices
library(grid) # For plotting multiple plots in one frame
library(gridExtra) # For plotting multiple plots in one frame
library(reshape2) # For "melting" dataframesb
library(scales) # To access break formatting functions
library(reshape)
library(ggrepel)
library(RColorBrewer)
library(pheatmap)
library(Seurat)
# Post loading the original data with the subsetted "TNK" object
# Looking at data at the original resolution
DimPlot(TNK, reduction = 'umap',group.by = "subset_louvain_res0.4", label = T)
tnk_Markers_0.4 <- read.table(file = "merged_colossal_paper_final_res0.6_TNK_subset_louvain_res0.4_markers.tsv", sep = '\t', header = TRUE,
stringsAsFactors = FALSE)
tnk_Markers_0.4_padj0.1 <- tnk_Markers_0.4[which(tnk_Markers_0.4$p_val_adj<0.1),]
tnk_Markers_0.4_padj0.1 <- tnk_Markers_0.4_padj0.1[order(tnk_Markers_0.4_padj0.1$avg_logFC,decreasing = TRUE),]
tnk_Markers_0.4_padj0.1 <- tnk_Markers_0.4_padj0.1[order(tnk_Markers_0.4_padj0.1$cluster,decreasing = FALSE),]
tnk_Markers_0.4_padj0.1_Top5 <- tnk_Markers_0.4_padj0.1 %>% group_by(cluster) %>% top_n(n = 5, wt = avg_logFC)
tnk_Markers_0.4_padj0.1_Top5_unique <- unique(tnk_Markers_0.4_padj0.1_Top5$gene)
pdf("DotPlot_TNK_res0.4.pdf", width = 10, height = 12)
DotPlot(TNK, features= tnk_Markers_0.4_padj0.1_Top5_unique, cols= 'RdYlBu', group.by = "subset_louvain_res0.4",
assay = "RNA", dot.scale = 5) + theme(axis.text.x=element_text(angle=45, hjust = 1, size = 10),
axis.text.y=element_text(size = 10), text = element_text(size = 14)) + coord_flip()
dev.off()
# cleaning up Neut and RBC contaminants
tnk_clean <- subset(TNK, subset = subset_louvain_res0.4 %in% c(0:8, 10))
tnk_clean <- RunUMAP(tnk_clean, dims = 1:30, reduction='harmony',
n.neighbors = 30, min.dist = 0.3, spread = 1, verbose = FALSE)
tnk_clean <- FindNeighbors(tnk_clean, dims = 1:30, k.param = 20,
verbose = FALSE, reduction = 'harmony')
tnk_clean_0.8 <- FindClusters(tnk_clean, verbose = TRUE,
algorithm = 1, resolution = 0.8, random.seed = 21212)
pdf("UMAP_tnk_clean_res0.8.pdf", width = 10, height = 10)
DimPlot(tnk_clean_0.8, reduction = 'umap', label = T)
dev.off()
tnk_clean_Markers_0.8 <- FindAllMarkers(tnk_clean_0.8, only.pos = TRUE,
min.pct = 0.25, logfc.threshold = 0.40, test.use="poisson",
assay = "RNA", latent.vars = "LIBRARY")
tnk_clean_Markers_0.8_padj0.1 <- tnk_clean_Markers_0.8[which(tnk_clean_Markers_0.8$p_val_adj<0.1),]
tnk_clean_Markers_0.8_padj0.1 <- tnk_clean_Markers_0.8_padj0.1[order(tnk_clean_Markers_0.8_padj0.1$avg_logFC,decreasing = TRUE),]
tnk_clean_Markers_0.8_padj0.1 <- tnk_clean_Markers_0.8_padj0.1[order(tnk_clean_Markers_0.8_padj0.1$cluster,decreasing = FALSE),]
tnk_clean_Markers_0.8_padj0.1_Top5 <- tnk_clean_Markers_0.8_padj0.1 %>% group_by(cluster) %>% top_n(n = 5, wt = avg_logFC)
tnk_clean_Markers_0.8_padj0.1_Top5_unique <- unique(tnk_clean_Markers_0.8_padj0.1_Top5$gene)
pdf("DotPlot_tnk_clean_res0.6.pdf", width = 10, height = 12)
DotPlot(tnk_clean_0.6, features= tnk_clean_Markers_0.6_padj0.1_Top5_unique, cols= 'RdYlBu',
assay = "RNA", dot.scale = 5) + theme(axis.text.x=element_text(angle=45, hjust = 1, size = 10),
axis.text.y=element_text(size = 10), text = element_text(size = 14)) + coord_flip()
dev.off()
# Adding the original identities at this resolution to the metadata
[email protected]$clean_louvain_res0.8 <- Idents(tnk_clean_0.8)
[email protected]$clean_louvain_res0.6 <- Idents(tnk_clean_0.6)
tnk_clean2 <- subset(tnk_clean, subset = clean_louvain_res0.8 %in% c(0:16))
tnk_clean2 <- RunUMAP(tnk_clean2, dims = 1:30, reduction='harmony',
n.neighbors = 30, min.dist = 0.3, spread = 1, verbose = FALSE)
tnk_clean2 <- FindNeighbors(tnk_clean2, dims = 1:30, k.param = 20,
verbose = FALSE, reduction = 'harmony')
tnk_clean2_0.8 <- FindClusters(tnk_clean2, verbose = TRUE,
algorithm = 1, resolution = 0.8, random.seed = 21212)
pdf("UMAP_tnk_clean2_res0.8.pdf", width = 10, height = 10)
DimPlot(tnk_clean2_0.8, reduction = 'umap', label = T)
dev.off()
tnk_clean2_Markers_0.8 <- FindAllMarkers(tnk_clean2_0.8, only.pos = TRUE,
min.pct = 0.25, logfc.threshold = 0.40, test.use="poisson",
assay = "RNA", latent.vars = "LIBRARY")
tnk_clean2_Markers_0.8_padj0.1 <- tnk_clean2_Markers_0.8[which(tnk_clean2_Markers_0.8$p_val_adj<0.1),]
tnk_clean2_Markers_0.8_padj0.1 <- tnk_clean2_Markers_0.8_padj0.1[order(tnk_clean2_Markers_0.8_padj0.1$avg_logFC,decreasing = TRUE),]
tnk_clean2_Markers_0.8_padj0.1 <- tnk_clean2_Markers_0.8_padj0.1[order(tnk_clean2_Markers_0.8_padj0.1$cluster,decreasing = FALSE),]
tnk_clean2_Markers_0.8_padj0.1_Top5 <- tnk_clean2_Markers_0.8_padj0.1 %>% group_by(cluster) %>% top_n(n = 5, wt = avg_logFC)
tnk_clean2_Markers_0.8_padj0.1_Top5_unique <- unique(tnk_clean2_Markers_0.8_padj0.1_Top5$gene)
pdf("DotPlot_tnk_clean2_res0.8.pdf", width = 10, height = 12)
DotPlot(tnk_clean2_0.8, features= tnk_clean_Markers_0.6_padj0.1_Top5_unique, cols= 'RdYlBu',
assay = "RNA", dot.scale = 5) + theme(axis.text.x=element_text(angle=45, hjust = 1, size = 10),
axis.text.y=element_text(size = 10), text = element_text(size = 14)) + coord_flip()
dev.off()
# annotating object for removal and re-Harmony
[email protected]$clean_louvain_res0.8 <- Idents(tnk_clean2_0.8)
tnk_clean3 <- subset(tnk_clean2, subset = clean_louvain_res0.8 %in% c(0:4, 6:16))
write.table([email protected], file = "TNK_clean_metadata.tsv", sep = '\t', qoute = FALSE)
# Post-Harmony reclean object
load("merged_colossal_paper_final_res0.6_TNK.RData")
tnk_reclean <- RunUMAP(TNK, dims = 1:30, reduction='harmony',
n.neighbors = 30, min.dist = 0.3, spread = 1, verbose = FALSE)
tnk_reclean <- FindNeighbors(tnk_reclean, dims = 1:30, k.param = 20,
verbose = FALSE, reduction = 'harmony')
tnk_reclean_0.8 <- FindClusters(tnk_reclean, verbose = TRUE,
algorithm = 1, resolution = 0.8, random.seed = 21212)
pdf("UMAP_tnk_reclean_res0.5.pdf", width = 10, height = 10)
DimPlot(tnk_reclean_0.8, reduction = 'umap', label = T)
dev.off()
tnk_reclean_Markers_0.8 <- FindAllMarkers(tnk_reclean_0.8, only.pos = TRUE,
min.pct = 0.25, logfc.threshold = 0.40, test.use="poisson",
assay = "RNA", latent.vars = "LIBRARY")
tnk_reclean_Markers_0.8_padj0.1 <- tnk_reclean_Markers_0.8[which(tnk_reclean_Markers_0.8$p_val_adj<0.1),]
tnk_reclean_Markers_0.8_padj0.1 <- tnk_reclean_Markers_0.8_padj0.1[order(tnk_reclean_Markers_0.8_padj0.1$avg_logFC,decreasing = TRUE),]
tnk_reclean_Markers_0.8_padj0.1 <- tnk_reclean_Markers_0.8_padj0.1[order(tnk_reclean_Markers_0.8_padj0.1$cluster,decreasing = FALSE),]
tnk_reclean_Markers_0.8_padj0.1_Top5 <- tnk_reclean_Markers_0.8_padj0.1 %>% group_by(cluster) %>% top_n(n = 5, wt = avg_logFC)
tnk_reclean_Markers_0.8_padj0.1_Top5_unique <- unique(tnk_reclean_Markers_0.8_padj0.1_Top5$gene)
pdf("DotPlot_tnk_reclean_res0.8.pdf", width = 10, height = 12)
DotPlot(tnk_reclean_0.8, features= tnk_reclean_Markers_0.8_padj0.1_Top5_unique, cols= 'RdYlBu',
assay = "RNA", dot.scale = 5) + theme(axis.text.x=element_text(angle=45, hjust = 1, size = 10),
axis.text.y=element_text(size = 10), text = element_text(size = 14)) + coord_flip()
dev.off()
# Final Cleanup-removing one platelet cluster
[email protected]$clean_louvain_res0.8 <- Idents(tnk_reclean_0.8)
tnk_reclean2 <- subset(tnk_reclean_0.8, subset = clean_louvain_res0.8 %in% c(0:12, 14:16)) # removing one platelet cluster
pdf("UMAP_tnk_reclean2_res0.8.pdf", width = 10, height = 10)
DimPlot(tnk_reclean2, reduction = 'umap', label = T)
dev.off()
tnk_reclean2 <- RunUMAP(tnk_reclean2, dims = 1:30, reduction='harmony',
n.neighbors = 30, min.dist = 0.5, spread = 1, verbose = FALSE)
[email protected]$reclean_louvain_0.8 <- Idents(tnk_reclean2)
annotations_all <- c("Naive CD4", "CD56hi NK", "Cytotoxic CD8", "Memory CD4", "Memory CD8", "Memory CD4", "CD56lo NK", "Naive CD8", "Treg",
"MAIT", "GammaDelta", "Cycling CD8", "Cycling CD4", "Cytotoxic CD8", "ISG+", "Naive CD4")
names(annotations_all) <- levels([email protected]$reclean_louvain_0.8)
tnk_reclean2 <- RenameIdents(tnk_reclean2, annotations_all)
[email protected]$annotated_clusters <- Idents(tnk_reclean2)
DimPlot(tnk_reclean2, reduction = 'umap', group.by = "annotated_clusters", label = T)
pdf("ViolinPlot_tnk_reclean2_Ident_Markers_Original.pdf", width = 10, height = 6)
VlnPlot(tnk_reclean2, features = c("IL7R", "LTB", "CD3E","CD3D", "CD8A","CD8B" ,"PRF1", "GNLY"), group.by = "reclean_louvain_0.8", ncol = 4, pt.size = 0)
dev.off()
write.table([email protected], "TNK_Reclean2_metadata.tsv", sep = '\t', quote = FALSE)
save (tnk_reclean2, tnk_reclean_Markers_0.8, file = "tnk_cleaned.RData")
mylevels_reclean2 <- c("Naive CD4", "Memory CD4", "Treg", "Cycling CD4", "ISG+", "Naive CD8",
"Memory CD8", "Cytotoxic CD8", "Cycling CD8", "CD56lo NK", "CD56hi NK", "MAIT", "GammaDelta")
mygenelist_reclean2 <- c("LTB", "IL7R", "VIM", "CCR7", "LEF1", "SELL", "TCF7", "JUNB", "FYB",
"AQP3", "PLP2", "S100A4", "CD52", "NOSIP",
"FOXP3", "RTKN2", "CD27", "IL32", "IL2RA",
"ACTG1", "COTL1", "HIST1H4C", "STMN1", "HMGB2", "TUBB", "MKI67",
"ISG15", "MX1", "IFI6", "IFIT3", "MT2A",
"CD8A", "CD8B",
"CCL5", "GZMK","DUSP2",
"GZMH","FGFBP2","GNLY",
"NKG7","TYROBP","SPON2",
"NCAM1", "PTGDS", "FCER1G",
"KLRB1", "TRAV1-2","DUSP1",
"TRDV2", "TRGV9", "TRDC")
[email protected]$annotated_clusters <- factor([email protected]$annotated_clusters, levels = mylevels_reclean2)
Idents(tnk_reclean2) <- factor (Idents(tnk_reclean2), levels = mylevels_reclean2)
pdf("DotPlot_tnk_reclean2_reorder_mygenelist.pdf", width = 9, height = 12)
DotPlot(tnk_reclean2, features= mygenelist_reclean2, cols= 'RdYlBu', group.by = "annotated_clusters",
assay = "RNA") + theme(axis.text.x=element_text(angle=45, hjust = 1, size = 10),
axis.text.y=element_text(size = 10), text = element_text(size = 14)) + coord_flip()
dev.off()
COMET_10X_CLINICAL_SCORES = read.csv("COMET_10X_CLINICAL_SCORES_PAPER.csv")
[email protected]$Age <- COMET_10X_CLINICAL_SCORES$Age[match([email protected]$SAMPLE.by.SNPs,COMET_10X_CLINICAL_SCORES$SAMPLE.by.SNPs)]
[email protected]$Day_after_onset <- COMET_10X_CLINICAL_SCORES$Day_after_onset[match([email protected]$SAMPLE.by.SNPs,COMET_10X_CLINICAL_SCORES$SAMPLE.by.SNPs)]
[email protected]$ICU_vs_FLOOR <- COMET_10X_CLINICAL_SCORES$ICU_vs_FLOOR[match([email protected]$SAMPLE.by.SNPs,COMET_10X_CLINICAL_SCORES$SAMPLE.by.SNPs)]
[email protected]$Other_infection_type <- COMET_10X_CLINICAL_SCORES$Other_infection_type[match([email protected]$SAMPLE.by.SNPs,COMET_10X_CLINICAL_SCORES$SAMPLE.by.SNPs)]
[email protected]$Death <- COMET_10X_CLINICAL_SCORES$Death[match([email protected]$SAMPLE.by.SNPs,COMET_10X_CLINICAL_SCORES$SAMPLE.by.SNPs)]
[email protected]$NIH_score <- COMET_10X_CLINICAL_SCORES$NIH_score[match([email protected]$SAMPLE.by.SNPs,COMET_10X_CLINICAL_SCORES$SAMPLE.by.SNPs)]
[email protected]$Sampling_score <- COMET_10X_CLINICAL_SCORES$Sampling_score[match([email protected]$SAMPLE.by.SNPs,COMET_10X_CLINICAL_SCORES$SAMPLE.by.SNPs)]
[email protected]$Overall_score <- COMET_10X_CLINICAL_SCORES$Overall_score[match([email protected]$SAMPLE.by.SNPs,COMET_10X_CLINICAL_SCORES$SAMPLE.by.SNPs)]
[email protected]$Qualitative_score <- COMET_10X_CLINICAL_SCORES$Qualitative_score[match([email protected]$SAMPLE.by.SNPs,COMET_10X_CLINICAL_SCORES$SAMPLE.by.SNPs)]
# Adding IF-responsive gene signature scores and plotting
ISG_genes = c('MT2A', 'ISG15', 'LY6E', 'IFIT1', 'IFIT2', 'IFIT3', 'IFITM1', 'IFITM3', 'IFI44L', "IFI6", "MX1", "IFI27") #final ISG signature
# Get the cell values from the Seurat object
gene_df <- get_genes_s3(ISG_genes, tnk_reclean2)
# Define the score vector. In this case, we call the score for a cell the mean value of all genes in the vector
score <- colMeans(gene_df)
[email protected]$ISG_score <- score[match(colnames(tnk_reclean2), names(score))]
pdf("ViolinPlot_ISG_score_Covid_All.pdf", width = 8, height = 10)
VlnPlot(tnk_reclean2, features = "ISG_score", group.by = "covid_status", pt.size = 0, cols = mycolorscovid)
dev.off()
pdf("ViolinPlot_ISG_score_Severity_All.pdf", width = 8, height = 10)
VlnPlot(tnk_reclean2, features = "ISG_score", group.by = "Qualitative_score", pt.size = 0, cols = mycolorsseverity)
dev.off()
pdf("ViolinPlot_ISG_score_Severity_splitBycovid.pdf", width = 8, height = 10)
VlnPlot(tnk_reclean2, features = "ISG_score", group.by = "covid_status", split.by = "Qualitative_score", pt.size = 0, cols = mycolorsseverity)
dev.off()
pdf("UMAP_tnk_reclean2_Severity.pdf", width = 10, height = 10)
DimPlot(tnk_reclean2, reduction = 'umap', group.by = "Qualitative_score", label = F, cols = mycolorsseverity)
#dittoDimPlot(tnk_reclean, var = 'reclean_louvain_previous', reduction.use = 'umap', do.label = F)
dev.off()
pdf("UMAP_tnk_reclean2_Severity.pdf", width = 10, height = 10)
dittoDimPlot(tnk_reclean2, var = 'Qualitative_score', reduction.use = 'umap', do.label = F, colors = mycolorsseverity)
dev.off()
# annotating for Correlation/phEMD analysis
annotations_final <- c("TNK_Naive_CD4", "TNK_Memory_CD4", "TNK_Treg", "TNK_Cycling_CD4",
"TNK_ISG+", "TNK_Naive_CD8", "TNK_Memory_CD8", "TNK_Cytotoxic_CD8", "TNK_Cycling_CD8",
"TNK_CD56lo_NK", "TNK_CD56hi_NK", "TNK_MAIT", "TNK_GammaDelta")
names(annotations_final) <- levels([email protected]$annotated_clusters)
[email protected]$harmony_final_clusters <- annotations_final[match([email protected]$annotated_clusters, names(annotations_final))]
write.table([email protected], "TNK_Final_metadata.tsv", sep = '\t', quote = FALSE)
save (tnk_reclean2, file = "TNK_Final.RData")
## Calculating frequencies across subtypes with Covid status and Severity
# By Covid Status
y <- table([email protected][,c('SAMPLE.by.SNPs', 'annotated_clusters')])
#melt table into one column by id='SAMPLE.by.SNPs; call 'x'
ClustersByCovidStatus <- melt(y, id='SAMPLE.by.SNPs')
#make data.frame from tnk that contains SAMPLE.by.SNPs event #s; call 'temp'
temp <- table([email protected]$SAMPLE.by.SNPs)
temp <- data.frame(temp)
#add a column to x called 'total' that contains total # of cells per SAMPLE.by.SNPs
ClustersByCovidStatus$total <- temp$Freq[match(ClustersByCovidStatus$SAMPLE.by.SNPs, temp$Var1)]
#add a column to x called 'freq' that gives the fraction of cells of each seurat_clusters per total # of cells in sample
ClustersByCovidStatus$freq <- ClustersByCovidStatus$value/ClustersByCovidStatus$total
#generate temp2 dataframe that has covid_status matched to SAMPLE.by.SNPs
temp2 <- [email protected][,c('SAMPLE.by.SNPs', 'covid_status')]
temp2 <- data.frame(temp2)
#add column to x called 'covid_status' that adds pos/neg/ctrl to x matched to SAMPLE.by.SNPs
ClustersByCovidStatus$covid_status <- temp2$covid_status[match(ClustersByCovidStatus$SAMPLE.by.SNPs, temp2$SAMPLE.by.SNPs)]
#make 'seurat_clusters' numbers factors, not values; or else ggplot gets confused
ClustersByCovidStatus$annotated_clusters <- factor(ClustersByCovidStatus$annotated_clusters)
write.csv(ClustersByCovidStatus, file = "Freq_Covid_Status.csv")
# By Severity
z <- table([email protected][,c('SAMPLE.by.SNPs', 'annotated_clusters')])
#melt table into one column by id='SAMPLE.by.SNPs; call 'x'
ClustersBySeverity <- melt(z, id='SAMPLE.by.SNPs')
#make data.frame from tnk that contains SAMPLE.by.SNPs event #s; call 'temp'
temp3 <- table([email protected]$SAMPLE.by.SNPs)
temp3 <- data.frame(temp3)
#add a column to x called 'total' that contains total # of cells per SAMPLE.by.SNPs
ClustersBySeverity$total <- temp3$Freq[match(ClustersBySeverity$SAMPLE.by.SNPs, temp3$Var1)]
#add a column to x called 'freq' that gives the fraction of cells of each seurat_clusters per total # of cells in sample
ClustersBySeverity$freq <- ClustersBySeverity$value/ClustersBySeverity$total
#ggplot(ClustersByICU, aes(x=broad_clusters, y=freq))+geom_point(stat = 'identity')
temp4 <- [email protected][,c('SAMPLE.by.SNPs', 'Qualitative_score')]
temp4 <- data.frame(temp4)
ClustersBySeverity$Qualitative_score <- temp4$Qualitative_score[match(ClustersBySeverity$SAMPLE.by.SNPs, temp4$SAMPLE.by.SNPs)]
ClustersBySeverity$annotated_clusters <- factor(ClustersBySeverity$annotated_clusters)
write.csv(ClustersBySeverity, file = "Freq_Severity.csv")
## Calculating ISG scores for different groups
ISG_score_ctrl <- [email protected]$ISG_score[[email protected]$covid_status == 'ctrl']
ISG_score_neg <- [email protected]$ISG_score[[email protected]$covid_status == 'neg']
ISG_score_pos <- [email protected]$ISG_score[[email protected]$covid_status == 'pos']
ISG_score_Covid <- cbind(ISG_score_ctrl, ISG_score_neg, ISG_score_pos)
write.csv(ISG_score_Covid, file = "ISG_score_Covid.csv")
ISG_score_neg_mild <- [email protected]$ISG_score[[email protected]$covid_status == 'neg' & [email protected]$Qualitative_score == 'MILD']
ISG_score_neg_severe <- [email protected]$ISG_score[[email protected]$covid_status == 'neg' & [email protected]$Qualitative_score == 'SEVERE']
ISG_score_pos_mild <- [email protected]$ISG_score[[email protected]$covid_status == 'pos' & [email protected]$Qualitative_score == 'MILD']
ISG_score_pos_severe <- [email protected]$ISG_score[[email protected]$covid_status == 'pos' & [email protected]$Qualitative_score == 'SEVERE']
ISG_score_Severity <- cbind(ISG_score_neg_mild, ISG_score_neg_severe, ISG_score_pos_mild, ISG_score_pos_severe)
write.csv(ISG_score_Severity, file = "ISG_score_Severity.csv")
##--END