-
Notifications
You must be signed in to change notification settings - Fork 1
/
4.5_Fig5_Global_spread.R
443 lines (396 loc) · 23.7 KB
/
4.5_Fig5_Global_spread.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
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
#======load packages======
library(readxl)
library(patchwork)
library(dplyr)
library(rgdal)
library(ggpattern)
library(ggplot2)
library(sf)
library(stringr)
library(ggpubr)
#==read global map==
worldmap <- st_read("World_map.shp")
st_crs(worldmap)
worldmap <- st_transform(worldmap, CRS("+proj=longlat +ellps=WGS84 +towgs84=0,0,0,0,0,0,0 +no_defs"))
worldmap$iso_a3[worldmap$name == "Taiwan"] <- "CHN"
nine <- st_read("nine.shp")
nine <- st_transform(nine, CRS("+proj=longlat +ellps=WGS84 +towgs84=0,0,0,0,0,0,0 +no_defs"))
WHO_country<-read_excel("Data/WHO_country.xlsx",sheet=1)
#==read sequence data==
mydata.variant4 <- readRDS("F:/Output/Total_seq_full_20211031_after_dedu2_2.rds")
mydata.variant4 <- mydata.variant4[,c(3,4,7)]
names(mydata.variant4)[3] <- "lineage"
mydata.variant4 <- mydata.variant4 %>% filter(!(lineage == "B.1.1.7" & date_collect < as.Date("2020-09-20")))
mydata.variant4 <- mydata.variant4 %>% filter(!(lineage == "B.1.351" & date_collect < as.Date("2020-05-11")))
mydata.variant4 <- mydata.variant4 %>% filter(!(lineage == "P.1" & date_collect < as.Date("2020-11-03")))
mydata.variant4 <- mydata.variant4 %>% filter(!(lineage == "B.1.617.2" & date_collect < as.Date("2020-10-23")))
mydata.variant4 <- mydata.variant4 %>% filter(!(lineage == "C.37" & date_collect < as.Date("2020-12-22")))
mydata.variant4 <- mydata.variant4 %>% filter(!(lineage == "B.1.621" & date_collect < as.Date("2021-01-11")))
mydata.variant1 <- mydata.variant4 %>% filter(!is.na(date_collect))
all(!is.na(mydata.variant1$country)) # T
mydata.variant1$period[mydata.variant1$date_collect <= as.Date("2020-12-31") &
mydata.variant1$date_collect >= as.Date("2020-01-01")] <- "2020"
mydata.variant1$period[mydata.variant1$date_collect <= as.Date("2021-03-31") &
mydata.variant1$date_collect >= as.Date("2021-01-01")] <- "Jan-Mar"
mydata.variant1$period[mydata.variant1$date_collect <= as.Date("2021-06-30") &
mydata.variant1$date_collect >= as.Date("2021-04-01")] <- "Apr-Jun"
mydata.variant1$period[mydata.variant1$date_collect <= as.Date("2021-10-31") &
mydata.variant1$date_collect >= as.Date("2021-07-01")] <- "Jul-Oct"
mydata.variant1 <- mydata.variant1 %>% filter(date_collect >= as.Date("2020-01-01")) %>% filter(date_collect <= as.Date("2021-10-31"))
tmp <- mydata.variant1 %>% group_by(country, period, lineage) %>% dplyr::summarise(case = n())
tmp1 <- mydata.variant1 %>% group_by(country, period) %>% dplyr::summarise(seq = n())
tmp <- left_join(tmp, tmp1)
tmp2 <- left_join(tmp, WHO_country[,c(2:4)], by = c("country" = "GBD_location"))
all(!is.na(tmp2$ISO3)) #T
#==stat==
res <- tmp2 %>% group_by(period,lineage) %>%
summarise(var = sum(case))
res1 <- res %>% group_by(period) %>% summarise(var1 = sum(var))
res2 <- full_join(res, res1)
res2$prop <- round(res2$var * 100/ res2$var1, digits = 1)
#==plot global map==
first <- "2020"
second <- "Jan-Mar"
third <- "Apr-Jun"
four <- "Jul-Oct"
tmp2 <- tmp2 %>% filter(seq > 10) # only those countries with more than 10 sequences in each period are calculated
worldmap1 <- left_join(worldmap, tmp2[(tmp2$period == first & tmp2$lineage == "B.1.1.7"),], by=c("iso_a3" = "ISO3"))
worldmap2 <- left_join(worldmap,tmp2[(tmp2$period == second & tmp2$lineage == "B.1.1.7"),], by=c("iso_a3" = "ISO3"))
worldmap3 <- left_join(worldmap,tmp2[(tmp2$period == third & tmp2$lineage == "B.1.1.7"),], by=c("iso_a3" = "ISO3"))
worldmap4 <- left_join(worldmap,tmp2[(tmp2$period == four & tmp2$lineage == "B.1.1.7"),], by=c("iso_a3" = "ISO3"))
worldmap5 <- left_join(worldmap,tmp2[(tmp2$period == first & tmp2$lineage == "B.1.351"),], by=c("iso_a3" = "ISO3"))
worldmap6 <- left_join(worldmap, tmp2[(tmp2$period == second & tmp2$lineage == "B.1.351"),], by=c("iso_a3" = "ISO3"))
worldmap7 <- left_join(worldmap,tmp2[(tmp2$period == third & tmp2$lineage == "B.1.351"),], by=c("iso_a3" = "ISO3"))
worldmap8 <- left_join(worldmap,tmp2[(tmp2$period == four & tmp2$lineage == "B.1.351"),], by=c("iso_a3" = "ISO3"))
worldmap9 <- left_join(worldmap,tmp2[(tmp2$period == first & tmp2$lineage == "P.1"),], by=c("iso_a3" = "ISO3"))
worldmap10 <- left_join(worldmap,tmp2[(tmp2$period == second & tmp2$lineage == "P.1"),], by=c("iso_a3" = "ISO3"))
worldmap11 <- left_join(worldmap,tmp2[(tmp2$period == third & tmp2$lineage == "P.1"),], by=c("iso_a3" = "ISO3"))
worldmap12 <- left_join(worldmap,tmp2[(tmp2$period == four & tmp2$lineage == "P.1"),], by=c("iso_a3" = "ISO3"))
worldmap13 <- left_join(worldmap, tmp2[(tmp2$period == first & tmp2$lineage == "B.1.617.2"),], by=c("iso_a3" = "ISO3"))
worldmap14 <- left_join(worldmap,tmp2[(tmp2$period == second & tmp2$lineage == "B.1.617.2"),], by=c("iso_a3" = "ISO3"))
worldmap15 <- left_join(worldmap, tmp2[(tmp2$period == third & tmp2$lineage == "B.1.617.2"),], by=c("iso_a3" = "ISO3"))
worldmap16 <- left_join(worldmap,tmp2[(tmp2$period == four & tmp2$lineage == "B.1.617.2"),], by=c("iso_a3" = "ISO3"))
worldmap17 <- left_join(worldmap, tmp2[(tmp2$period == first & tmp2$lineage == "C.37"),], by=c("iso_a3" = "ISO3"))
worldmap18 <- left_join(worldmap,tmp2[(tmp2$period == second & tmp2$lineage == "C.37"),], by=c("iso_a3" = "ISO3"))
worldmap19 <- left_join(worldmap,tmp2[(tmp2$period == third & tmp2$lineage == "C.37"),], by=c("iso_a3" = "ISO3"))
worldmap20 <- left_join(worldmap,tmp2[(tmp2$period == four & tmp2$lineage == "C.37"),], by=c("iso_a3" = "ISO3"))
worldmap21 <- left_join(worldmap, tmp2[(tmp2$period == first & tmp2$lineage == "Ref"),], by=c("iso_a3" = "ISO3"))
worldmap22 <- left_join(worldmap,tmp2[(tmp2$period == second & tmp2$lineage == "Ref"),], by=c("iso_a3" = "ISO3"))
worldmap23 <- left_join(worldmap,tmp2[(tmp2$period == third & tmp2$lineage == "Ref"),], by=c("iso_a3" = "ISO3"))
worldmap24 <- left_join(worldmap,tmp2[(tmp2$period == four & tmp2$lineage == "Ref"),], by=c("iso_a3" = "ISO3"))
worldmap25 <- left_join(worldmap, tmp2[(tmp2$period == first & tmp2$lineage == "B.1.621"),], by=c("iso_a3" = "ISO3"))
worldmap26 <- left_join(worldmap, tmp2[(tmp2$period == second & tmp2$lineage == "B.1.621"),], by=c("iso_a3" = "ISO3"))
worldmap27 <- left_join(worldmap, tmp2[(tmp2$period == third & tmp2$lineage == "B.1.621"),], by=c("iso_a3" = "ISO3"))
worldmap28 <- left_join(worldmap, tmp2[(tmp2$period == four & tmp2$lineage == "B.1.621"),], by=c("iso_a3" = "ISO3"))
#==Global map==
theme_map1 <- theme(axis.ticks = element_blank(),
axis.line = element_blank(),
axis.text = element_blank(),
legend.justification=c(0,0),
legend.position = c(-0.03,0.05) ,
legend.background = element_blank(),
legend.direction = "horizontal",
legend.key.size = unit(0.35,"cm"),
legend.key.width = unit(0.22,"cm"),
legend.key.height = unit(0.18,"cm"),
legend.title = element_text(size = 5),
legend.text = element_text(size = 5,lineheight=4),
legend.spacing = unit(4, "cm"), legend.spacing.x = NULL, # the spacing between legends (unit)
legend.spacing.y = NULL,#the spacing between legends (unit)
axis.title = element_blank(),
plot.margin = margin(0.35, 0.35, 0.35, 0.35, "cm"),
panel.background = element_rect(fill = "white"),
legend.box.background =element_blank(),
legend.box.margin= margin(0, 0, 0, 0, "cm"),
panel.spacing = unit(0,"cm"),
panel.border = element_rect(fill='transparent',colour="transparent"),
legend.key = element_blank(),
plot.title = element_text(hjust = 0.5, size = 6, vjust = 0))
theme_map <- theme(axis.ticks = element_blank(),
axis.line = element_blank(),
axis.text = element_blank(),
legend.justification=c(0,0),
legend.position = "none" ,
legend.background = element_blank(),
legend.key.size = unit(0.4,"cm"),
legend.key.width = unit(0.15,"cm"),
legend.title = element_text(size = 5),
legend.text = element_text(size = 5,lineheight=4),
legend.spacing = unit(4, "cm"), legend.spacing.x = NULL, # the spacing between legends (unit)
legend.spacing.y = NULL,#the spacing between legends (unit)
axis.title = element_blank(),
plot.margin = margin(0.35, 0.35, 0.35, 0.35, "cm"),
panel.background = element_rect(fill = "white"),
legend.box.background =element_blank(),
legend.box.margin= margin(0, 0, 0, 0, "cm"),
panel.spacing = unit(0,"cm"),
panel.border = element_rect(fill='transparent',colour="transparent"),
legend.key = element_blank(),
plot.title = element_text(hjust = 0.5, size = 6, vjust = 0))
ggplot() +
geom_sf(data = worldmap21, aes(fill = case/seq), size = 0.12, color = "black")+
geom_sf(data = nine, color="black",size = 0.06)+
coord_sf(xlim = c(-170, 170), ylim = c(-57, 90))+
theme_map+
scale_fill_gradient(low = "white",high = "#005757", na.value="grey85",
breaks = seq(0, 1, 0.25), labels = c("0","","50","","100"),limits = c(0, 1))+
ggtitle("Jan-Dec, 2020")-> p21
ggplot() +
geom_sf(data = worldmap22, aes(fill = case/seq), size = 0.12, color = "black")+
geom_sf(data = nine, color="black",size = 0.06)+
coord_sf(xlim = c(-170, 170), ylim = c(-57, 90)) +
theme_map+
scale_fill_gradient(low = "white",high = "#005757", na.value="grey85",
breaks = seq(0, 1, 0.25), labels = c("0","","50","","100"),limits = c(0, 1))+
ggtitle("Jan-Mar, 2021")-> p22
ggplot() +
geom_sf(data = worldmap23, aes(fill = case/seq), size = 0.12, color = "black")+
geom_sf(data = nine, color="black",size = 0.06)+
coord_sf(xlim = c(-170, 170), ylim = c(-57, 90)) +
theme_map+
scale_fill_gradient(low = "white",high = "#005757", na.value="grey85",
breaks = seq(0, 1, 0.25), labels = c("0","","50","","100"),limits = c(0, 1))+
ggtitle("Apr-Jun, 2021")-> p23
ggplot() +
geom_sf(data = worldmap24, aes(fill = case/seq), size = 0.12, color = "black")+
geom_sf(data = nine, color="black",size = 0.06)+
coord_sf(xlim = c(-170, 170), ylim = c(-57, 90)) +
theme_map1+
scale_fill_gradient("Non-variant",
low = "white",high = "#005757", na.value="grey85",
breaks = seq(0, 1, 0.25), labels = c("0","","50","","100"),limits = c(0, 1))+
guides(fill = guide_colorbar(frame.colour = "black",
title.position="top",
ticks.colour = "black"),
size = guide_legend(title.position="top"))+
ggtitle("Jul-Oct, 2021")-> p24
p24
ggplot() +
geom_sf(data = worldmap1, aes(fill = case/seq), size = 0.12, color = "black")+
geom_sf(data = nine, color="black",size = 0.06)+
coord_sf(xlim = c(-170, 170), ylim = c(-57, 90)) +
theme_map+
scale_fill_gradient(low = colorRampPalette(c(colorRampPalette(c("white", "#21610B"))(20)[2], "#21610B"))(20)[2],
high = "#21610B", na.value="grey85",
breaks = seq(0, 1, 0.25), labels = c("0","","50","","100"),limits = c(0, 1))-> p1
ggplot() +
geom_sf(data = worldmap2, aes(fill = case/seq), size = 0.12, color = "black")+
geom_sf(data = nine, color="black",size = 0.06)+
coord_sf(xlim = c(-170, 170), ylim = c(-57, 90)) +
theme_map+
scale_fill_gradient(low = colorRampPalette(c(colorRampPalette(c("white", "#21610B"))(20)[2], "#21610B"))(20)[2],
high = "#21610B", na.value="grey85",
breaks = seq(0, 1, 0.25), labels = c("0","","50","","100"),limits = c(0, 1))-> p2
ggplot() +
geom_sf(data = worldmap3, aes(fill = case/seq), size = 0.12, color = "black")+
geom_sf(data = nine, color="black",size = 0.06)+
coord_sf(xlim = c(-170, 170), ylim = c(-57, 90)) +
theme_map+
scale_fill_gradient(low = colorRampPalette(c("white", "#21610B"))(20)[2],high = "#21610B", na.value="grey85",
breaks = seq(0, 1, 0.25), labels = c("0","","50","","100"),limits = c(0, 1))-> p3
ggplot() +
geom_sf(data = worldmap4, aes(fill = case/seq), size = 0.12, color = "black")+
geom_sf(data = nine, color="black",size = 0.06)+
coord_sf(xlim = c(-170, 170), ylim = c(-57, 90)) +
theme_map1+
guides(fill = guide_colorbar(frame.colour = "black",
title.position="top",
ticks.colour = "black"),
size = guide_legend(title.position="top"))+
scale_fill_gradient("Alpha",
low = colorRampPalette(c("white", "#21610B"))(20)[2],high = "#21610B", na.value="grey85",
breaks = seq(0, 1, 0.25), labels = c("0","","50","","100"),limits = c(0, 1))-> p4
ggplot() +
geom_sf(data = worldmap5, aes(fill = case/seq), size = 0.12, color = "black")+
geom_sf(data = nine, color="black",size = 0.06)+
coord_sf(xlim = c(-170, 170), ylim = c(-57, 90)) +
theme_map+
scale_fill_gradient(low = colorRampPalette(c("white", "#868A08"))(20)[2],high = "#868A08", na.value="grey85",
breaks = seq(0, 1, 0.25), labels = c("0","","50","","100"),limits = c(0, 1)) -> p5
ggplot() +
geom_sf(data = worldmap6, aes(fill = case/seq), size = 0.12, color = "black")+
geom_sf(data = nine, color="black",size = 0.06)+
coord_sf(xlim = c(-170, 170), ylim = c(-57, 90)) +
theme_map+
scale_fill_gradient(low = colorRampPalette(c("white", "#868A08"))(20)[2],high = "#868A08", na.value="grey85",
breaks = seq(0, 1, 0.25), labels = c("0","","50","","100"),limits = c(0, 1)) -> p6
ggplot() +
geom_sf(data = worldmap7, aes(fill = case/seq), size = 0.12, color = "black")+
geom_sf(data = nine, color="black",size = 0.06)+
coord_sf(xlim = c(-170, 170), ylim = c(-57, 90)) +
theme_map+
scale_fill_gradient(low = colorRampPalette(c("white", "#868A08"))(20)[2],high = "#868A08", na.value="grey85",
breaks = seq(0, 1, 0.25), labels = c("0","","50","","100"),limits = c(0, 1)) -> p7
ggplot() +
geom_sf(data = worldmap8, aes(fill = case/seq), size = 0.12, color = "black")+
geom_sf(data = nine, color="black",size = 0.06)+
coord_sf(xlim = c(-170, 170), ylim = c(-57, 90)) +
theme_map1+
guides(fill = guide_colorbar(frame.colour = "black",
title.position="top",
ticks.colour = "black"),
size = guide_legend(title.position="top"))+
scale_fill_gradient("Beta",
low = colorRampPalette(c("white", "#868A08"))(20)[2],high = "#868A08", na.value="grey85",
breaks = seq(0, 1, 0.25), labels = c("0","","50","","100"),limits = c(0, 1)) -> p8
ggplot() +
geom_sf(data = worldmap9, aes(fill = case/seq), size = 0.12, color = "black")+
geom_sf(data = nine, color="black",size = 0.06)+
coord_sf(xlim = c(-170, 170), ylim = c(-57, 90)) +
theme_map+
scale_fill_gradient(low = colorRampPalette(c("white", "#380B61"))(20)[2],high = "#380B61", na.value="grey85",
breaks = seq(0, 1, 0.25), labels = c("0","","50","","100"),limits = c(0, 1)) -> p9
ggplot() +
geom_sf(data = worldmap10, aes(fill = case/seq), size = 0.12, color = "black")+
geom_sf(data = nine, color="black",size = 0.06)+
coord_sf(xlim = c(-170, 170), ylim = c(-57, 90)) +
theme_map+
scale_fill_gradient(low = colorRampPalette(c("white", "#380B61"))(20)[2],high = "#380B61", na.value="grey85",
breaks = seq(0, 1, 0.25), labels = c("0","","50","","100"),limits = c(0, 1)) -> p10
ggplot() +
geom_sf(data = worldmap11, aes(fill = case/seq), size = 0.12, color = "black")+
geom_sf(data = nine, color="black",size = 0.06)+
coord_sf(xlim = c(-170, 170), ylim = c(-57, 90)) +
theme_map+
scale_fill_gradient(low = colorRampPalette(c("white", "#380B61"))(20)[2],high = "#380B61", na.value="grey85",
breaks = seq(0, 1, 0.25), labels = c("0","","50","","100"),limits = c(0, 1)) -> p11
ggplot() +
geom_sf(data = worldmap12, aes(fill = case/seq), size = 0.12, color = "black")+
geom_sf(data = nine, color="black",size = 0.06)+
coord_sf(xlim = c(-170, 170), ylim = c(-57, 90)) +
theme_map1+
guides(fill = guide_colorbar(frame.colour = "black",
title.position="top",
ticks.colour = "black"),
size = guide_legend(title.position="top"))+
scale_fill_gradient("Gamma",
low = colorRampPalette(c("white", "#380B61"))(20)[2],high = "#380B61", na.value="grey85",
breaks = seq(0, 1, 0.25), labels = c("0","","50","","100"),limits = c(0, 1)) -> p12
ggplot() +
geom_sf(data = worldmap13, aes(fill = case/seq), size = 0.12, color = "black")+
geom_sf(data = nine, color="black",size = 0.06)+
coord_sf(xlim = c(-170, 170), ylim = c(-57, 90)) +
theme_map+
scale_fill_gradient(low = colorRampPalette(c("white", "#8A0808"))(20)[2],high = "#8A0808", na.value="grey85",
breaks = seq(0, 1, 0.25), labels = c("0","","50","","100"),limits = c(0, 1)) -> p13
ggplot() +
geom_sf(data = worldmap14, aes(fill = case/seq), size = 0.12, color = "black")+
geom_sf(data = nine, color="black",size = 0.06)+
coord_sf(xlim = c(-170, 170), ylim = c(-57, 90)) +
theme_map+
scale_fill_gradient(low = colorRampPalette(c("white", "#8A0808"))(20)[2],high = "#8A0808", na.value="grey85",
breaks = seq(0, 1, 0.25), labels = c("0","","50","","100"),limits = c(0, 1)) -> p14
ggplot() +
geom_sf(data = worldmap15, aes(fill = case/seq), size = 0.12, color = "black")+
geom_sf(data = nine, color="black",size = 0.06)+
coord_sf(xlim = c(-170, 170), ylim = c(-57, 90)) +
theme_map+
scale_fill_gradient(low = colorRampPalette(c("white", "#8A0808"))(20)[2],high = "#8A0808", na.value="grey85",
breaks = seq(0, 1, 0.25), labels = c("0","","50","","100"),limits = c(0, 1)) -> p15
ggplot() +
geom_sf(data = worldmap16, aes(fill = case/seq), size = 0.12, color = "black")+
geom_sf(data = nine, color="black",size = 0.06)+
coord_sf(xlim = c(-170, 170), ylim = c(-57, 90)) +
theme_map1+
guides(fill = guide_colorbar(frame.colour = "black",
title.position="top",
ticks.colour = "black"),
size = guide_legend(title.position="top"))+
scale_fill_gradient("Delta",
low = colorRampPalette(c("white", "#8A0808"))(20)[2],high = "#8A0808", na.value="grey85",
breaks = seq(0, 1, 0.25), labels = c("0","","50","","100"),limits = c(0, 1)) -> p16
#==output==
tiff(file = "Output0106/Fig5.tiff", width = 8.2, height = 5.8,
units = "in", compression = "lzw", res = 600)
ggarrange(p21,p22,p23,p24,p1,p2,p3,p4,p5,p6,p7,p8,p9,p10,p11,p12,
p13,p14,p15, p16, nrow =5, ncol = 4, heights = c(1.1,1,1,1,1))
dev.off()
pdf(file = "Output0106/Fig5.pdf", width = 8.2, height = 5.8)
ggarrange(p21,p22,p23,p24,p1,p2,p3,p4,p5,p6,p7,p8,p9,p10,p11,p12,
p13,p14,p15, p16, nrow =5, ncol = 4, heights = c(1.1,1,1,1,1))
dev.off()
#==for VOI==
ggplot() +
geom_sf(data = worldmap17[worldmap17$name != "Greenland",], aes(fill = case/seq), size = 0.2, color = "black")+
geom_sf(data = nine, color="black",size = 0.15)+
coord_sf(xlim = c(-170, -40), ylim = c(-57, 90)) +
theme_map+
labs(title = "Jan-Dec, 2020")+
theme(plot.title = element_text(hjust = 0.5, size = 6, vjust = 0))+
scale_fill_gradient(low = "white",high = "#004B97", na.value="grey85",
breaks = seq(0, 1, 0.25), labels = c("0","","50","","100"),limits = c(0, 1)) -> p17
ggplot() +
geom_sf(data = worldmap18[worldmap18$name != "Greenland",], aes(fill = case/seq), size = 0.2, color = "black")+
geom_sf(data = nine, color="black",size = 0.15)+
coord_sf(xlim = c(-170, -40), ylim = c(-57, 90)) +
theme_map+
labs(title = "Jan-Mar, 2021")+
theme(plot.title = element_text(hjust = 0.5, size = 6, vjust = 0))+
scale_fill_gradient(low = "white",high = "#004B97", na.value="grey85",
breaks = seq(0, 1, 0.25), labels = c("0","","50","","100"),limits = c(0, 1)) -> p18
ggplot() +
geom_sf(data = worldmap19[worldmap19$name != "Greenland",], aes(fill = case/seq), size = 0.2, color = "black")+
geom_sf(data = nine, color="black",size = 0.15)+
coord_sf(xlim = c(-170, -40), ylim = c(-57, 90)) +
theme_map+
labs(title = "Apr-Jun, 2021")+
theme(plot.title = element_text(hjust = 0.5, size = 6, vjust = 0))+
scale_fill_gradient(low = "white",high = "#004B97", na.value="grey85",
breaks = seq(0, 1, 0.25), labels = c("0","","50","","100"),limits = c(0, 1)) -> p19
ggplot() +
geom_sf(data = worldmap20[worldmap20$name != "Greenland",], aes(fill = case/seq), size = 0.2, color = "black")+
geom_sf(data = nine, color="black",size = 0.15)+
coord_sf(xlim = c(-170, -40), ylim = c(-57, 90)) +
theme_map1+
ggtitle("Jul-Oct, 2021")+
theme(plot.title = element_text(hjust = 0.5, size = 6, vjust = 0),
legend.key.width = unit(0.4,"cm"))+
guides(fill = guide_colorbar(frame.colour = "black",
title.position="top",
ticks.colour = "black"),
size = guide_legend(title.position="top"))+
scale_fill_gradient("Prevalence of Lambda (%)",
low = "white",high = "#004B97", na.value="grey85",
breaks = seq(0, 1, 0.25), labels = c("0","25","50","75","100"),limits = c(0, 1)) -> p20
p20
ggplot() +
geom_sf(data = worldmap25[worldmap25$name != "Greenland",], aes(fill = case/seq), size = 0.2, color = "black")+
geom_sf(data = nine, color="black",size = 0.15)+
coord_sf(xlim = c(-170, -40), ylim = c(-57, 90)) +
theme_map+
scale_fill_gradient(low = "white",high = "#DE5258", na.value="grey85",
breaks = seq(0, 1, 0.25), labels = c("0","","50","","100"),limits = c(0, 1)) -> p25
ggplot() +
geom_sf(data = worldmap26[worldmap26$name != "Greenland",], aes(fill = case/seq), size = 0.2, color = "black")+
geom_sf(data = nine, color="black",size = 0.15)+
coord_sf(xlim = c(-170, -40), ylim = c(-57, 90)) +
theme_map+
scale_fill_gradient(low = "white",high = "#DE5258", na.value="grey85",
breaks = seq(0, 1, 0.25), labels = c("0","","50","","100"),limits = c(0, 1)) -> p26
ggplot() +
geom_sf(data = worldmap27[worldmap27$name != "Greenland",], aes(fill = case/seq), size = 0.2, color = "black")+
geom_sf(data = nine, color="black",size = 0.15)+
coord_sf(xlim = c(-170, -40), ylim = c(-57, 90)) +
theme_map+
theme(plot.title = element_text(hjust = 0.5, size = 7, vjust = 0))+
scale_fill_gradient(low = "white",high = "#DE5258", na.value="grey85",
breaks = seq(0, 1, 0.25), labels = c("0","","50","","100"),limits = c(0, 1)) -> p27
ggplot() +
geom_sf(data = worldmap28[worldmap28$name != "Greenland",], aes(fill = case/seq),size = 0.2, color = "black")+
geom_sf(data = nine, color="black",size = 0.15)+
coord_sf(xlim = c(-170, -40), ylim = c(-57, 90)) +
theme_map1+
theme(legend.key.width = unit(0.4,"cm"))+
guides(fill = guide_colorbar(frame.colour = "black",
title.position="top",
ticks.colour = "black"),
size = guide_legend(title.position="top"))+
scale_fill_gradient("Prevalence of Mu (%)",
low = "white",high = "#DE5258", na.value="grey85",
breaks = seq(0, 1, 0.25), labels = c("0","25","50","75","100"),limits = c(0, 1)) -> p28
p28
#==output==
tiff(file = "Output0106/Extended/ED_Fig8.tif", width = 7, height = 4,
units = "in", compression = "lzw", res = 300)
ggarrange(p17,p18,p19, p20,p25,p26,p27,p28,nrow =2, ncol = 4, heights = c(1,0.9))
dev.off()