Skip to content
Permalink
master
Switch branches/tags

Name already in use

A tag already exists with the provided branch name. Many Git commands accept both tag and branch names, so creating this branch may cause unexpected behavior. Are you sure you want to create this branch?

MalariaResistance/ResistanceAnalysis.Rmd

Go to file
 
 
Cannot retrieve contributors at this time
1375 lines (1210 sloc) 60.8 KB
Raw Blame
Open in GitHub Desktop
---
title: "Resistance"
author: "hqx"
date: "2022-12-02"
output: html_document
---
```{r setup, include=FALSE}
knitr::opts_knit$set(root.dir = '~/Dropbox/research/current/resistance/MalariaResistance/')
knitr::opts_chunk$set(echo = TRUE)
```
```{r include=FALSE}
library(tidyverse)
library(ggplot2)
library(patchwork)
library(cowplot)
library(readxl)
library(betareg)
library(modelr)
library(ggthemes)
library(scales)
library(ggpattern)
```
## Resistance running results
Loading resistance results with varying strain diversity and Transmission Potential
```{r}
straRef<-data.frame(nstrains=ceiling(10^(seq(0.7,2.75,0.15))),rawStrains =10^(seq(0.7,2.75,0.15)))
#oneLR<-read.csv("v13_results_filtered.csv")
#oneLR<-read.csv("v14_results_all.csv")
#oneLR<-read.csv("v13_results_all.csv")
#oneLR<-read.csv("v13_622_results.csv")
#oneLR<-read.csv("v15_622_results.csv")
#oneLR<-read.csv("v625_all_results.csv")
oneLR<-read.csv("fullModel/v628_results_all.csv")
oneLR<-oneLR%>%mutate(PW_NU=round(PW_NU),
PW_SU=round(PW_SU),
PW_AU=round(PW_AU),
PW_ND=round(PW_ND),
PW_SD=round(PW_SD),
PW_AD=round(PW_AD),
PR_NU=round(PR_NU),
PR_SU=round(PR_SU),
PR_AU=round(PR_AU),
PR_ND=round(PR_ND),
PR_SD=round(PR_SD),
PR_AD=round(PR_AD))
oneLR<-oneLR%>%mutate(
m_N = P_acc_N/(NU+ND+1e-6),m_S =P_acc_S/(SU+SD+1e-6), m_A = P_acc_A/(AU+AD+1e-6),
PW = PW_NU+PW_SU+PW_AU+PW_ND+PW_SD+PW_AD,
PR = PR_NU+PR_SU+PR_AU+PR_ND+PR_SD+PR_AD,
P=PW+PR,N = NU+ND+SU+SD+AU+AD,
resistP = PR/P,
adjB = bavgnull*(0.5819519*(seasonality>0) + (seasonality==0)),
rNU = (PW_NU+PR_NU)/(NU+1e-6),
rND = (PW_ND+PR_ND)/(ND+1e-6),
rSU = (PW_SU+PR_SU)/(SU+1e-6),
rSD = (PW_SD+PR_SD)/(SD+1e-6),
rAU = (PW_AU+PR_AU)/(AU+1e-6),
rAD = (PW_AD+PW_AD)/(AD+1e-6),
I_NU = 1-exp(-rNU),
I_ND = 1-exp(-rND),
I_SU = 1-exp(-rSU),
I_SD = 1-exp(-rSD),
I_AU = 1-exp(-rAU),
I_AD = 1-exp(-rAD),
NN = NU+ND,
S = SU+SD,
A = AU+AD,
U = NU+SU+AU,
D = SD+AD+ND,
adjPrev = (I_NU*NU+I_ND*ND+I_SU*SU+I_SD*SD+I_AU*AU+I_AD*AD)/N
)
oneLR<-oneLR %>%mutate(infN = (1-1/nstrains)^m_N,
infS = (1-1/nstrains)^m_S,
infA = (1-1/nstrains)^m_A)
oneLR<-oneLR%>%left_join(straRef,by="nstrains")
```
Pick Peak Prevalence Region
```{r}
subDat<- oneLR %>% filter(wildOnly==1,d1==0,mixCost==0.9,cloneCost==0.1,varyCost==1,ResistanceEfficacy==0.006666667) %>%
group_by(nstrains) %>%
arrange(adjPrev,desc=T) %>%
slice_max(adjPrev, prop=0.3) %>%
dplyr::select(nstrains,adjB) %>%
mutate(prev_order = seq_along(nstrains))
oneLR <- oneLR %>% left_join(subDat )
oneLR <- oneLR %>% mutate(scenario = paste(d1,wildOnly))
oneLR$wildOnlyChar <- factor(oneLR$wildOnly, levels = c(1,0),labels = c("Wild-type Only", "Resistance Invasion"))
oneLR$treatment<-factor(oneLR$d1,levels=c(0,0.05,0.2),labels=c("No Treatment","Treatment: 63.2%","Treatment: 98.2%"))
oneLR$resEffect<-factor(oneLR$ResistanceEfficacy,c(0.006666667,0.013071895,0.034657359),
labels = c("Resistance Efficacy: 87.5%","Resistance Efficacy: 77.0%","Resistance Efficacy: 50%"))
oneLR$scenario <- factor(oneLR$scenario, levels = c("0 1", "0.05 1","0.2 1","0 0", "0.05 0","0.2 0"),
labels = c("Wild-type Only; No Treatment",
"Wild-type Only; Treatment: 63.2%",
"Wild-type Only; Treatment: 98.2%",
"Resistant Only; No Treatment",
"Resistant invasion; Treatment: 63.2%",
"Resistant invasion; Treatment: 98.2%"))
#policy<-policy %>% left_join(subDat)
```
### calculate adjPrev and resistP as a function of biting and strain diversity
Prevalence without drug or resistant types
```{r}
subWild<-oneLR%>%filter(wildOnly==1,d1==0,
cloneCost==0.1,mixCost==0.9,varyCost==1,ResistanceEfficacy==0.006666667)
p1<-ggplot(subWild, aes(adjB,rawStrains))+geom_tile(aes(fill=adjPrev))+
geom_tile(data=subWild%>%filter(P<5),fill="grey")+
geom_tile(data = subWild%>%filter(d1==0,prev_order<=10,prev_order>1),color = "grey",fill=NA,linewidth=0.4)+
geom_tile(data = subWild%>%filter(d1==0,prev_order==1),color = "white",fill=NA,linewidth=0.5)+
#geom_raster(data = oneLR%>%filter(wildOnly==0,d1==0,prev_order<=10),fill = "white",alpha=0.3)+
xlab("Transmission Potential") + ylab("Number of Strains")+
scale_fill_viridis_c(option="magma",name = "Prevalence",
limits = c(0, 1),values=seq(0,1,0.25))+
scale_x_continuous(trans = "log10")+
scale_y_continuous(trans = "log10")+
facet_grid(.~scenario)+
theme_cowplot()+
theme(legend.position="right", legend.key.height = unit(1,"cm"))
```
### Resistant allele percentage with clone cost = 0.1 and mix cost = 0.9
```{r}
oneLR<-oneLR%>%dplyr::select(-RW)
subP<-oneLR%>%filter(wildOnly==1,ResistanceEfficacy==0.006666667)%>%dplyr::select(No,d1,adjPrev)%>%mutate(RW=adjPrev)%>%dplyr::select(No,d1,RW)
oneLR<-oneLR%>%left_join(subP)
```
```{r}
#addRateName <- function(string){
# paste("Daily Treatment Rate:",string)
#}
p2<-ggplot(oneLR%>%filter(wildOnly==0,d1>0,varyCost==1,cloneCost == 0.1, mixCost==0.9,ResistanceEfficacy==0.006666667,adjPrev>0.01),
aes(adjB,rawStrains))+
geom_tile(aes(fill=resistP))+
#geom_tile(data = oneLR%>%filter(wildOnly==0,d1>0,varyCost==1,cloneCost == 0.1, mixCost==0.9,prev_order<=10),color = "white",fill=NA,size=0.5)+
xlab("Transmission Potential") + ylab("Number of Strains")+
scale_x_continuous(trans = "log10")+
scale_y_continuous(trans = "log10")+
scale_fill_viridis_c(option="magma",name = "Resistance\nfrequency",
limits = c(0, 1),values=seq(0,1,0.25))+
facet_grid(resEffect~treatment)+
theme_cowplot()+
theme(legend.position="right", legend.key.height = unit(1,"cm"))
```
### varyCost == 1, vary over d1, clone cost ==0.1, mix cost == 0.9
```{r}
oneLR <- oneLR %>% mutate(case = paste(cloneCost,mixCost,sep=";"))
oneLR$case <- factor(oneLR$case)
p3<-ggplot(oneLR%>%filter(wildOnly==0, d1>0,varyCost==1,P>1),
aes(adjPrev,resistP))+
geom_point(aes(color=case))+
#scale_colour_tableau("Classic Blue-Red 12")+
scale_color_manual(values=c('#80cdc1','#018571','#fdcc8a','#fc8d59','#d7301f','#cccccc','#969696','#525252',
'#bdc9e1','#74a9cf','#0570b0'),name=expression(paste("s"[single],";s"[mixed])))+
#scale_fill_manual(values=c('#f7f4f9','#e7e1ef','#d4b9da','#c994c7','#df65b0','#e7298a','#ce1256','#980043','#67001f'),name=expression(paste("s"[single],";s"[mixed])))+
#scale_color_manual(values=c('#d7191c','#fdae61','#74add1'))+
#geom_point(aes(color=nstrains))+
#scale_color_viridis_c(option="turbo",trans="log10",name = "number of\nstrains")+
#geom_point(aes(color=!is.na(prev_order)))+
#scale_color_viridis_d(option="turbo")+
facet_grid(treatment~resEffect)+
xlab("Prevalence")+ylab("Resistance Frequency")+
theme_cowplot()
```
### combine p1, p2, p3 to generate the fig2
```{r}
layout <- "
ABB
CCC
CCC
"
fig2<-p1+p2+p3+plot_layout(design=layout)+plot_annotation(tag_levels = 'A')
ggsave("Figures/fig2.tiff",fig2, width=11,height=8)
#ggsave("fig2.eps",fig2, width=11,height=8,device="eps")
ggsave("Figures/fig2.png",fig2, width=11,height=8)
ggsave("Figures/fig2.pdf",fig2, width=11,height=8)
```
Prevalence in wild-type only scenarios
```{r}
p23<-ggplot(oneLR%>%filter(wildOnly==1,
cloneCost==0.1,mixCost==0.9,ResistanceEfficacy==0.006666667), aes(adjB,rawStrains))+
geom_tile(aes(fill=adjPrev))+
geom_tile(data = oneLR%>%filter(wildOnly==1,
cloneCost==0.1,mixCost==0.9,ResistanceEfficacy==0.006666667,P<10),
fill="grey")+
xlab("Transmission Potential") + ylab("Number of Strains")+
scale_fill_viridis_c(option="magma",name = "Prevalence",
limits = c(0, 1),values=seq(0,1,0.25))+
scale_x_continuous(trans = "log10")+
scale_y_continuous(trans = "log10")+
facet_grid(.~scenario)+
theme_cowplot()+
theme(legend.position="right", legend.key.height = unit(1,"cm"),
plot.background = element_rect(fill="white"))
#ggsave("figS2_PW_prev.eps",p23,width=10,height=4,device="eps")
ggsave("Figures/figS2_PW_prev.tiff",p23,width=10,height=4)
ggsave("Figures/figS2_PW_prev.png",p23,width=10,height=4)
ggsave("Figures/figS2_PW_prev.pdf",p23,width=10,height=4)
```
### infectivity comparisons
```{r}
newLR<-oneLR %>% filter(wildOnly==0,varyCost==1,d1>0.05,ResistanceEfficacy==0.006666667,cloneCost==0.1,mixCost==0.9,P>1) %>%
dplyr::select(adjB,rawStrains,infN,infS,infA) %>%
pivot_longer(3:5,names_to="GI",values_to="infectivity")
newLR$GI_classes<-factor(newLR$GI,levels=c("infN","infS","infA"),
labels=c("G0","G1","G2"))
p16<-ggplot(newLR,
aes(adjB,rawStrains))+
geom_raster(aes(fill=infectivity))+
xlab("Transmission Potential") + ylab("Number of Strains")+
scale_x_continuous(trans = "log10")+
scale_y_continuous(trans = "log10")+
scale_fill_viridis_c(option="magma",name = "Infectivity",
limits = c(0, 1),values=seq(0,1,0.25))+
facet_grid(.~GI_classes)+
theme_cowplot()+
theme(legend.position="right", legend.key.height = unit(1,"cm"),
plot.background = element_rect(fill="white"))
ggsave("Figures/figS3_infectivity.tiff",p16,width=10,height=4)
ggsave("Figures/figS3_infectivity.png",p16,width=10,height=4)
ggsave("Figures/figS3_infectivity.pdf",p16,width=10,height=4)
```
```{r eval=FALSE, include=FALSE}
ggplot(oneLR%>%filter(wildOnly==0,d1>0,varyCost==0,cloneCost==0.6),
aes(adjB,rawStrains))+
geom_raster(aes(fill=resistP))+
xlab("Transmission Potential") + ylab("Number of Strains")+
scale_x_continuous(trans = "log10")+
scale_y_continuous(trans = "log10")+
scale_fill_viridis_c(option="turbo",name = "Resistance\nfrequency",
limits = c(0, 1),values=seq(0,1,0.25))+
facet_grid(~d1, labeller = labeller(d1 =addRateName))+
theme_cowplot()+
theme(legend.position="right", legend.key.height = unit(1,"cm"))
```
Resistance vs. Prevalence, cloneCost by Mix Cost
varyCost == 0
```{r eval=FALSE, include=FALSE}
addCostName <- function(string){
paste("Fixed Cost:",string)
}
ggplot(oneLR%>%filter(wildOnly==0, d1>0,varyCost==0),
aes(adjPrev,resistP))+
geom_point(aes(color=nstrains))+
scale_color_viridis_c(option="turbo",trans="log10",name = "number of\nstrains")+facet_grid(d1~mixCost,labeller=labeller(mixCost = addCostName,d1=addRateName))+
xlab("Prevalence")+ylab("Resistance Frequency")+
theme_cowplot()
```
varyCost == 1
```{r}
p141<-ggplot(oneLR%>%filter(wildOnly==0, d1==0.05,varyCost==1,P>1),
aes(adjPrev,resistP))+
geom_point(aes(color=adjB))+
scale_color_viridis_c(option="magma",trans="log10",name = "Transmission\n Potential")+
facet_grid(case~treatment+resEffect)+
xlab("Prevalence")+ylab("Resistance Frequency")+
theme_cowplot()
p142<-ggplot(oneLR%>%filter(wildOnly==0, d1==0.2,varyCost==1,P>1),
aes(adjPrev,resistP))+
geom_point(aes(color=adjB))+
scale_color_viridis_c(option="magma",trans="log10",name = "Transmission\n Potential")+
facet_grid(case~treatment+resEffect)+
xlab("Prevalence")+ylab("Resistance Frequency")+
theme_cowplot()
ggsave("Figures/figS4_prevResist.tiff",p141,width=10,height=10)
ggsave("Figures/figS4_prevResist.png",p141,width=10,height=10)
ggsave("Figures/figS4_prevResist.pdf",p141,width=10,height=10)
ggsave("Figures/figS5_prevResist.tiff",p142,width=10,height=10)
ggsave("Figures/figS5_prevResist.png",p142,width=10,height=10)
ggsave("Figures/figS5_prevResist.pdf",p142,width=10,height=10)
```
```{r eval=FALSE, include=FALSE}
ggplot(oneLR%>%filter(wildOnly==0, d1==0.05,varyCost==1,P>1,nstrains==80),
aes(log10(adjPrev*adjB/365),resistP))+
#geom_jitter(alpha=0.5,aes(color=as.factor(mixCost)))+
#scale_color_manual(values=c('#d7191c','#fdae61','#74add1'))+
geom_point(aes(color=nstrains))+
scale_color_viridis_c(option="turbo",trans="log10",name = "number of\nstrains")+
#geom_point(aes(color=!is.na(prev_order)))+
#scale_color_viridis_d(option="turbo")+
facet_grid(cloneCost~mixCost,
labeller="label_both")+
xlab("Prevalence")+ylab("Resistance Frequency")+
theme_cowplot()
```
## Empirical Prevalence vs. Resistance
Read in WWARN_ACT_Partner_Drug_Mol_Surveyor_Data for PfCRT.
```{r}
a<-read_excel("Empirical/WWARN_ACT_Partner_Drug_Mol_Surveyor_Data_2022.xls")
#a<-read_excel("WWARN_ACT_Partner_Drug_Mol_Surveyor_Data.xls")
#write.csv(unique(a%>%dplyr::select(lat,lon)),file="uniqueGeo.csv",row.names=F,col.names=c("latitude","longitude"))
a<-type_convert(a)
cont<-read.csv("Empirical/continents.csv")
a<-a%>%left_join(cont,by="country")
a$unique_id<-1:nrow(a)
a$lon<-round(a$lon,digits=5)
a$lat<-round(a$lat,digits=5)
geoPoint<-read.csv("Empirical/analysis_uniqueGeo.csv",na.strings = "null")
geoPoint$long<-round(geoPoint$long,digits = 5)
geoPoint$lat<-round(geoPoint$lat,digits = 5)
a_orig<-a
a<-a_orig
a<-a%>%left_join(geoPoint,by=c("lon"="long","lat"="lat"))
a<-a%>%mutate(yearDiff = abs((`study start`+`study end`)/2-year))
a<-a%>%group_by(unique_id)%>%mutate(minYearDiff = min(yearDiff))%>%
filter(yearDiff==minYearDiff) %>% slice_head(n=1)
#a%>%filter(country=="Cambodia")
a<-a%>%filter(!is.na(value))
a <- a%>%rename(study.start = `study start`)
a <- a%>%rename(study.end = `study end`,Prevalence = value)
```
Calculate the percentage of resistance:
"if y also assumes the extremes 0 and 1, a useful transformation in practice is (y · (n − 1) + 0.5)/n where n is the sample size"
Smithson M, Verkuilen J (2006). “A Better Lemon Squeezer? Maximum-Likelihood Regression with Beta-Distributed Dependent Variables.” Psychological Methods, 11(1), 54–71.
```{r}
a<-a%>%mutate(resisP1 = (present+(`Mixed present`-present)*0.5)/tested,
resisP2 = (present+(`Mixed present`-present))/(tested+(`Mixed present`-present)),
# only consider non-mixed inf
resisP3 = present/tested,
# beta regression
resisP4 = (resisP3*(tested-1)+0.5)/tested,
resisP5 = (resisP2*(tested-1)+0.5)/tested)
```
Perform beta regression of the data, Senegal was excluded because the country had different first line treatment from 1990-1994
```{r}
b<-a%>%filter(study.end<=2000,study.start>=1990,tested>20, country!="Senegal")
tesb<-betareg(resisP5~Prevalence,weight=tested,data= b)
dgrid<-b %>% data_grid(Prevalence)
dgrid<-dgrid %>% add_predictions(tesb)
```
Plot the empirical together with regression line
```{r}
p4<-ggplot(b,aes(Prevalence,resisP5))+
geom_point(aes(size = tested,shape=Continent,color=Continent),alpha=0.8) +
scale_size(name = "Sample Size")+
geom_line(aes(Prevalence,pred),data=dgrid,linetype=1,color="grey20")+
labs(
#x=expression(paste(italic("Pf"),"Prevalence"[2-10])),
x = expression(paste(italic("P. falciparum"), " Prevalence in Children Aged 2-10")),
y=expression(paste("Frequency of Resistant Genotypes (",italic("pfcrt"), " 76T)",sep=""))) +
scale_color_manual(values=c('#d7191c','#fdae61','#74add1','#525252'))+
#scale_color_colorblind()+
theme_cowplot()+
theme(axis.title = ggtext::element_markdown(),plot.background = element_rect(fill="white"))
ggsave2("Figures/fig3_empprev.tiff",p4,width=8,height=6)
ggsave2("Figures/fig3_empprev.pdf",p4,width=8,height=6)
ggsave2("Figures/fig3_empprev.png",p4,width=8,height=6)
```
### Empirical VC and nStrain
transmissibility: 0.5 churcher et al. 2017
VC
1. PNG: VC, 2.2-29.2; Infectivity: 0.1 Anopheles farauti, Timinao et al. 2021
2. Asia: VC, 0.06-2.6, Rhosenberg et al. 1990
2.29 4.73, RATTANARITHIKUl et al. 1996, southern thailand
2.37-6.5, 2000, Lao, sekong., Vythilingam etal. 2003
The Transmission Potential of An. dirus was 0.009-0.428, while that of An. minimus was 0.048-0.186. Lao. Toma et al. 2000
0.06-0.63 Turkey
China 0.5-0.7
Infectivity: 0.1 Coleman et al. 2004, Anopheles dirus, Thailand
3. south america: VC, An. Darlingi, 0.84-3.6, Brazil, 2022, zimmerman
infectivity: 0.35, KLEIN et al. 1991
4 Africa: VC 0.02 Bali, nigeria,
infectivity 10% Bousema et al . 2011
VC: 0.001-20
Diversity
1. PNG 290-1094, tessama et al. 2015
2. asia: 1100-1700 Calculated from Tonkin-Hill et al. 2021
3. South America: 113-351, Rougeron et al. 2017
4. Africa: 4000-20000 Otto, Tan, Day, chen
```{r}
vcd<-read_excel("Empirical/VC_D.xlsx",sheet="Combined")
#vcd <- vcd %>%rowwise%>% mutate(meanTP = mean(TP_min, TP_max), meanD = mean(D_min,D_max))
newcd <- rbind(setNames(vcd[,c("Continent","TP_min","D_min")],c("Continent","TP","D")),
setNames(vcd[,c("Continent","TP_min","D_max")],c("Continent","TP","D")),
setNames(vcd[,c("Continent","TP_max","D_max")],c("Continent","TP","D")),
setNames(vcd[,c("Continent","TP_max","D_min")],c("Continent","TP","D")))
newcd <- newcd%>%ungroup()%>%arrange(Continent)
newcd$TP[newcd$TP<min(subWild$adjB)]<-min(subWild$adjB)
newcd$D[newcd$D<min(subWild$rawStrains)]<-min(subWild$rawStrains)
```
```{r}
#addRateName <- function(string){
# paste("Daily Treatment Rate:",string)
#}
p35<-ggplot(oneLR%>%filter(wildOnly==0,d1>0,varyCost==1,cloneCost == 0.1,mixCost ==0.9, ResistanceEfficacy==0.006666667,adjPrev>0.01, treatment=="Treatment: 98.2%"),
aes(adjB,rawStrains))+
geom_tile(aes(fill=resistP))+
geom_polygon(data = newcd, aes(TP,D, color=Continent),linewidth=1.5,fill=NA)+
#geom_tile(data = oneLR%>%filter(wildOnly==0,d1>0,varyCost==1,cloneCost == 0.1, mixCost==0.9,prev_order<=10),color = "white",fill=NA,size=0.5)+
xlab("Transmission Potential") + ylab("Number of Strains")+
scale_x_continuous(trans = "log10")+
scale_y_continuous(trans = "log10")+
#scale_fill_viridis_c(option="magma",name = "Resistance\nfrequency",
# limits = c(0, 1),values=seq(0,1,0.25))+
scale_fill_gradient(high="grey90",low="grey20",
limits=c(0,1),name = "Simulated\nResistance\nFrequency")+
scale_color_manual(values=c('#d7191c','#fdae61','#74add1','#525252'))+
ggtitle("Empirical Parameter Ranges of Continents")+
#facet_grid(resEffect~treatment)+
theme_cowplot()+
theme(legend.position="right", legend.key.height = unit(1,"cm"))
```
```{r}
#fig3<-p4/p35+
# plot_annotation(tag_levels = 'A')
ggsave("Figures/fig3_empRange.pdf",p35, width=8,height=6)
ggsave("Figures/fig3_empRange.png",p35, width=8,height=6)
ggsave("Figures/fig3_empRange.tiff",p35, width=8,height=6)
```
## detailed analyses of percentage of compartments of P, N, and prev
```{r}
# individual class
library(ggpattern)
NProp<-oneLR%>%filter( prev_order==1,varyCost==1) %>%
select(No,NU,ND,SU,SD,AU,AD,rawStrains,adjB,mixCost,cloneCost,d1,wildOnly) %>%
pivot_longer(c(NU,ND,SU,SD,AU,AD), names_to = "NCompart", values_to = "NCVal")
NProp<-NProp %>% group_by(No,d1,wildOnly) %>% mutate(PrevProp = NCVal/sum(NCVal))
NProp <- NProp %>% mutate(HostImm = str_sub(NCompart,1,1),
DrugTreat = str_sub(NCompart,2,2))
NProp$DrugTreat <- factor(NProp$DrugTreat, levels = c("U", "D"),
labels = c("Untreated", "Treated"))
NProp$HostImm <- factor(NProp$HostImm, levels = c("N","S","A"),
labels = c("G0", "G1", "G2"))
NProp <- NProp %>% mutate(scenario = paste(d1,wildOnly))
NProp$wildOnlyChar <- factor(NProp$wildOnly, levels = c(1,0),labels = c("Wild-type Only", "Resistance Invasion"))
NProp$scenario <- factor(NProp$scenario, levels = c("0 1", "0.05 1","0.2 1","0 0", "0.05 0","0.2 0"),
labels = c("Wild-type Only; No Treatment",
"Wild-type Only; Daily Treatment Rate: 0.05",
"Wild-type Only; Daily Treatment Rate: 0.2",
"Resistant Only; No Treatment",
"Resistant invasion; Daily Treatment Rate: 0.05",
"Resistant invasion; Daily Treatment Rate: 0.2"))
NProp$NCompart <- factor(NProp$NCompart, levels = c("NU","ND","SU","SD","AU","AD"))
#ggplot(NProp%>%filter(cloneCost==0.1,mixCost==0.9)) +
# geom_col(aes(log10(nstrains),PrevProp, fill=NCompart),position = "stack")+
# scale_fill_brewer(palette = "Paired")+
# geom_line(data = oneLR%>%filter(prev_order==1,cloneCost==0.1,mixCost==0.9,varyCost==1),
# aes(log10(nstrains), resistP))+
# geom_line(data = oneLR%>%filter(prev_order==1,cloneCost==0.1,mixCost==0.9,varyCost==1),
# aes(log10(nstrains), adjPrev),color="blue")+
# facet_grid(wildOnly~d1,scales="free_y",labeller = "label_both")
p5<-ggplot(NProp%>%filter(cloneCost==0.1,mixCost==0.9)) +
geom_col_pattern(aes(x = rawStrains,y = PrevProp, fill = HostImm, pattern = DrugTreat),
color ="grey30",
size = 0.2,
pattern_fill = "black",
pattern_angle = 45,
pattern_density = 0.1,
pattern_spacing = 0.025,
pattern_key_scale_factor = 0.6) +
geom_line(data = oneLR%>%filter(prev_order==1,cloneCost==0.1,mixCost==0.9,varyCost==1),
aes(rawStrains, resistP),linetype = 2, color='#d7191c')+
scale_x_continuous(trans = "log10", name = "Number of Strains")+
scale_y_continuous(
# Features of the first axis
name = "Fraction of Hosts",
# Add a second axis and specify its features
sec.axis = sec_axis(~., name="Resistance Frequency")
)+
scale_fill_manual(values = c('#74add1','#ffffb2','#fdae61'), name = "Generalized\nImmunity") +
#scale_color_manual(values = c('#74add1','#d7191c','#fdae61')) +
scale_pattern_manual(values = c(Treated = "stripe", Untreated = "none"), name = "Drug")+
facet_grid(wildOnlyChar~d1,scales="free_y",labeller = labeller(d1=addRateName))+
guides(pattern = guide_legend(override.aes = list(fill = "white")),
fill = guide_legend(override.aes = list(pattern = "none")))+
theme_cowplot()+
theme(axis.title.y.right = element_text(color='#d7191c'),
axis.text.y.right = element_text(color='#d7191c'),
plot.background = element_rect(fill="white"))
#ggsave("fig4_propHost.tiff", p5, width=10,height=6)
ggsave2("fig4_propHost.pdf",p5, width=10,height=6)
ggsave2("fig4_propHost.png",p5, width=10,height=6)
ggsave2("fig4_propHost.tiff",p5, width=10,height=6)
```
# fix at one strain diversity
inspect the effect of increasing Transmission Potential
```{r eval=FALSE, include=FALSE}
# individual class
#library(ggpattern)
NProp<-oneLR%>%filter( nstrains==113,varyCost==1) %>%
select(No,NU,ND,SU,SD,AU,AD,rawStrains,adjB,mixCost,cloneCost,d1,wildOnly) %>%
pivot_longer(c(NU,ND,SU,SD,AU,AD), names_to = "NCompart", values_to = "NCVal")
NProp<-NProp %>% group_by(No,d1,wildOnly) %>% mutate(PrevProp = NCVal/sum(NCVal))
NProp <- NProp %>% mutate(HostImm = str_sub(NCompart,1,1),
DrugTreat = str_sub(NCompart,2,2))
NProp$DrugTreat <- factor(NProp$DrugTreat, levels = c("U", "D"),
labels = c("Untreated", "Treated"))
NProp$HostImm <- factor(NProp$HostImm, levels = c("N","S","A"),
labels = c("G0", "G1", "G2"))
NProp <- NProp %>% mutate(scenario = paste(d1,wildOnly))
NProp$wildOnlyChar <- factor(NProp$wildOnly, levels = c(1,0),labels = c("Wild-type Only", "Resistance Invasion"))
NProp$scenario <- factor(NProp$scenario, levels = c("0 1", "0.05 1","0.2 1","0 0", "0.05 0","0.2 0"),
labels = c("Wild-type Only; No Treatment",
"Wild-type Only; Daily Treatment Rate: 0.05",
"Wild-type Only; Daily Treatment Rate: 0.2",
"Resistant Only; No Treatment",
"Resistant invasion; Daily Treatment Rate: 0.05",
"Resistant invasion; Daily Treatment Rate: 0.2"))
NProp$NCompart <- factor(NProp$NCompart, levels = c("NU","ND","SU","SD","AU","AD"))
#ggplot(NProp%>%filter(cloneCost==0.1,mixCost==0.9)) +
# geom_col(aes(log10(nstrains),PrevProp, fill=NCompart),position = "stack")+
# scale_fill_brewer(palette = "Paired")+
# geom_line(data = oneLR%>%filter(prev_order==1,cloneCost==0.1,mixCost==0.9,varyCost==1),
# aes(log10(nstrains), resistP))+
# geom_line(data = oneLR%>%filter(prev_order==1,cloneCost==0.1,mixCost==0.9,varyCost==1),
# aes(log10(nstrains), adjPrev),color="blue")+
# facet_grid(wildOnly~d1,scales="free_y",labeller = "label_both")
ggplot(NProp%>%filter(cloneCost==0.1,mixCost==0.9)) +
geom_col_pattern(aes(x = adjB,y = PrevProp, fill = HostImm, pattern = DrugTreat),
color ="grey30",
size = 0.2,
pattern_fill = "black",
pattern_angle = 45,
pattern_density = 0.1,
pattern_spacing = 0.025,
pattern_key_scale_factor = 0.6) +
geom_line(data = oneLR%>%filter(nstrains==113,cloneCost==0.1,mixCost==0.9,varyCost==1),
aes(adjB, resistP),linetype = 2, color='#d7191c')+
scale_x_continuous(trans = "log10", name = "Transmission Potential")+
scale_y_continuous(
# Features of the first axis
name = "Fraction of Hosts",
# Add a second axis and specify its features
sec.axis = sec_axis(~., name="Resistance Frequency")
)+
scale_fill_manual(values = c('#74add1','#ffffb2','#fdae61'), name = "Generalized\nImmunity") +
#scale_color_manual(values = c('#74add1','#d7191c','#fdae61')) +
scale_pattern_manual(values = c(Treated = "stripe", Untreated = "none"), name = "Drug")+
facet_grid(wildOnlyChar~d1,scales="free_y",labeller = labeller(d1=addRateName))+
guides(pattern = guide_legend(override.aes = list(fill = "white")),
fill = guide_legend(override.aes = list(pattern = "none")))+
theme_cowplot()+
theme(axis.title.y.right = element_text(color='#d7191c'),
axis.text.y.right = element_text(color='#d7191c'))
```
```{r eval=FALSE, include=FALSE}
# individual class
#library(ggpattern)
NProp<-oneLR%>%filter( nstrains==113,varyCost==1) %>%
select(No,PW_NU,PW_ND,PW_SU,PW_SD,PW_AU,PW_AD, PR_NU,PR_ND,PR_SU,PR_SD,PR_AU,PR_AD,rawStrains,adjB,mixCost,cloneCost,d1,wildOnly,adjPrev) %>%
pivot_longer(c(PW_NU,PW_ND,PW_SU,PW_SD,PW_AU,PW_AD, PR_NU,PR_ND,PR_SU,PR_SD,PR_AU,PR_AD), names_to = "NCompart", values_to = "NCVal")
NProp<-NProp %>%mutate(Genotype = str_sub(NCompart,1,2),
HostImm = str_sub(NCompart,4,4),
DrugTreat = str_sub(NCompart,5,5))
NProp<-NProp %>% group_by(No,d1,wildOnly) %>% mutate(PrevProp = NCVal/sum(NCVal))
NProp$DrugTreat <- factor(NProp$DrugTreat, levels = c("U", "D"),
labels = c("Untreated", "Treated"))
NProp$HostImm <- factor(NProp$HostImm, levels = c("N","S","A"),
labels = c("G0", "G1", "G2"))
NProp <- NProp %>% mutate(scenario = paste(d1,wildOnly))
NProp$wildOnlyChar <- factor(NProp$wildOnly, levels = c(1,0),labels = c("Wild-type Only", "Resistance Invasion"))
NProp$scenario <- factor(NProp$scenario, levels = c("0 1", "0.05 1","0.2 1","0 0", "0.05 0","0.2 0"),
labels = c("Wild-type Only; No Treatment",
"Wild-type Only; Daily Treatment Rate: 0.05",
"Wild-type Only; Daily Treatment Rate: 0.2",
"Resistant Only; No Treatment",
"Resistant invasion; Daily Treatment Rate: 0.05",
"Resistant invasion; Daily Treatment Rate: 0.2"))
#NProp$NCompart <- factor(NProp$NCompart, levels = c("NU","ND","SU","SD","AU","AD"))
#ggplot(NProp%>%filter(cloneCost==0.1,mixCost==0.9)) +
# geom_col(aes(log10(nstrains),PrevProp, fill=NCompart),position = "stack")+
# scale_fill_brewer(palette = "Paired")+
# geom_line(data = oneLR%>%filter(prev_order==1,cloneCost==0.1,mixCost==0.9,varyCost==1),
# aes(log10(nstrains), resistP))+
# geom_line(data = oneLR%>%filter(prev_order==1,cloneCost==0.1,mixCost==0.9,varyCost==1),
# aes(log10(nstrains), adjPrev),color="blue")+
# facet_grid(wildOnly~d1,scales="free_y",labeller = "label_both")
ggplot(NProp%>%filter(cloneCost==0.1,mixCost==0.9,wildOnly==0,d1>0)) +
geom_col_pattern(aes(x = adjB,y = PrevProp, fill = HostImm, pattern = DrugTreat),
color ="grey30",
size = 0.2,
pattern_fill = "black",
pattern_angle = 45,
pattern_density = 0.05,
pattern_spacing = 0.035,
pattern_key_scale_factor = 0.6) +
geom_line(
aes(adjB, adjPrev),linetype = 2, color="blue")+
scale_x_continuous(trans = "log10", name = "Transmission Potential")+
scale_y_continuous(
# Features of the first axis
name = "Fraction of Parasites\n",
# Add a second axis and specify its features
sec.axis = sec_axis(~., name="Prevalence")
)+
scale_fill_manual(values = c('#74add1','#ffffb2','#fdae61'), name = "Generalized\nImmunity") +
#scale_color_manual(values = c('#74add1','#d7191c','#fdae61')) +
scale_pattern_manual(values = c(Treated = "stripe", Untreated = "none"), name = "Drug")+
facet_grid(Genotype~d1,scales="free_y",labeller = labeller(d1=addRateName,Genotype=c(PR="Resistant Parasites", PW="Wild-type Parasites")))+
guides(pattern = guide_legend(override.aes = list(fill = "white")),
fill = guide_legend(override.aes = list(pattern = "none")))+
theme_cowplot()+
theme(axis.title.y.right = element_text(color="blue"))
```
## Example of progression in the course of infection for the self-limiting property of resistance
```{r eval=FALSE, include=FALSE}
head(NProp%>%filter(cloneCost==0.1,mixCost==0.9,d1==0.05,rawStrains>112,rawStrains<120))
head(NProp%>%filter(cloneCost==0.1,mixCost==0.9,d1==0.05,rawStrains>19,rawStrains<30))
unique(oneLR$nstrains)
```
### Combine time series of high and low diversity
```{r}
tmSer1 <- read_csv("fullModel/v13_param_1506_ResisInvTimeSeries.csv")
tmSer1$run <- "lowDiv"
tmSer2 <- read_csv("fullModel/v13_param_3354_ResisInvTimeSeries.csv")
tmSer2$run <- "highDiv"
tmSer <- rbind(tmSer1,tmSer2)
tmSer$run <- factor(tmSer$run,levels=c("lowDiv","highDiv"))
hostDynm <- tmSer %>% filter(variable %in% c("NU","ND","SU","SD","AU","AD")) %>%
mutate(year = floor(time/365))
hostDynm <- hostDynm %>% group_by(run,year,variable) %>% summarize(value = mean(value)) %>%
mutate(Genotype = "Host") %>%
ungroup() %>% group_by(run,year) %>% mutate(H = sum(value),
prevProp = value/H,
HostImm = str_sub(variable,1,1),
DrugTreat = str_sub(variable,2,2))
hostDynm$DrugTreat <- factor(hostDynm$DrugTreat, levels = c("U", "D"),
labels = c("Untreated", "Treated"))
hostDynm$HostImm <- factor(hostDynm$HostImm, levels = c("N","S","A"),
labels = c("G0", "G1", "G2"))
# parasite
parasiteDynm <- tmSer %>% filter(variable %in% c("PW_NU","PW_ND","PW_SU","PW_SD","PW_AU","PW_AD",
"PR_NU","PR_ND","PR_SU","PR_SD","PR_AU","PR_AD")) %>%
mutate(year = floor(time/365))
parasiteDynm <- parasiteDynm %>% group_by(run,year,variable) %>%
summarize(value = mean(value)) %>% mutate(Genotype = str_sub(variable,1,2)) %>%
ungroup() %>% group_by(run,year,Genotype) %>% mutate(H = sum(value),
prevProp = value/H,
HostImm = str_sub(variable,4,4),
DrugTreat = str_sub(variable,5,5))
parasiteDynm$DrugTreat <- factor(parasiteDynm$DrugTreat, levels = c("U", "D"),
labels = c("Untreated", "Treated"))
parasiteDynm$HostImm <- factor(parasiteDynm$HostImm, levels = c("N","S","A"),
labels = c("G0", "G1", "G2"))
```
### Fig 5
```{r}
p6<-ggplot(hostDynm %>% filter(year<30)) +
geom_col_pattern(aes(x = year,y = value, fill = HostImm, pattern = DrugTreat),
color ="grey20",
width = 1,
linewidth=0.1,
pattern_fill = "black",
#pattern_color = NA,
pattern_angle = 45,
pattern_density = 0.1,
pattern_spacing = 0.02,
pattern_key_scale_factor = 0.6
)+
ylab("Host Population Size")+
scale_fill_manual(values = c('#74add1','#ffffb2','#fdae61'), name = "Generalized\nImmunity")+
#scale_color_manual(values = c('#74add1','#d7191c','#fdae61')) +
scale_pattern_manual(values = c(Treated = "stripe", Untreated = "none"),
name = "Drug")+
facet_grid(.~run,scales="free_y",
labeller = labeller(
run = c(highDiv = "High Diversity: 113 strains",
lowDiv = "Low Diversity: 20 strains")))+
guides(pattern = guide_legend(override.aes = list(fill = "white",color="black")),
fill = guide_legend(override.aes = list(pattern = "none")))+
theme_cowplot()+
theme(panel.grid.major.y = element_line(color="grey",linewidth = 0.5,linetype=2))
#ggsave("test.tiff",p6, width=6,height=12)
p8<-ggplot(parasiteDynm %>% filter(year<30)) +
geom_col_pattern(aes(x = year,y = value, fill = HostImm, pattern = DrugTreat),
color ="grey20",
width = 1,
linewidth=0.1,
pattern_fill = "black",
#pattern_color = NA,
pattern_angle = 45,
pattern_density = 0.1,
pattern_spacing = 0.02,
pattern_key_scale_factor = 0.6
)+
ylab("Parasite Population Size")+
scale_fill_manual(values = c('#74add1','#ffffb2','#fdae61'), name = "Generalized\nImmunity")+
#scale_color_manual(values = c('#74add1','#d7191c','#fdae61')) +
scale_pattern_manual(values = c(Treated = "stripe", Untreated = "none"),
name = "Drug")+
facet_grid(Genotype~run,switch="y",
labeller = labeller(Genotype=c(Host = "Host",PR="Resistant",
PW="Wild-type"),
run = c(highDiv = "High Diversity: 113 strains",
lowDiv = "Low Diversity: 20 strains")))+
guides(pattern = guide_legend(override.aes = list(fill = "white",color="black")),
fill = guide_legend(override.aes = list(pattern = "none")))+
theme_cowplot()+
theme(panel.grid.major.y = element_line(color="grey",linewidth = 0.5,linetype=2))
#ggsave("test.tiff",p6, width=6,height=12)
p7<-ggplot(parasiteDynm %>% filter(year<30,Genotype=="PW",run=="lowDiv")) +
geom_col_pattern(aes(x = year,y = value, fill = HostImm, pattern = DrugTreat),
color ="grey20",
width = 1,
linewidth=0.1,
pattern_fill = "black",
#pattern_color = NA,
pattern_angle = 45,
pattern_density = 0.1,
pattern_spacing = 0.02,
pattern_key_scale_factor = 0.6
)+
ylab("Wild-type")+
scale_fill_manual(values = c('#74add1','#ffffb2','#fdae61'), name = "Generalized\nImmunity")+
#scale_color_manual(values = c('#74add1','#d7191c','#fdae61')) +
scale_pattern_manual(values = c(Treated = "stripe", Untreated = "none"),
name = "Drug")+
theme_cowplot()+
theme(legend.position="none",
plot.background = element_rect(fill="white",color="black"))
#ggsave("test.tiff",p6, width=6,height=12)
layout<-"
AA
BB
BB
"
#fig4<-p6+inset_element(p7,left=0.1,bottom=0.12,right=0.5,top=0.333)
fig4<-p6+(p8+inset_element(p7,left=0.1,bottom=0.2,right=0.5,top=0.5))+
plot_annotation(tag_levels = 'A')+
plot_layout(design=layout,guides = 'collect')
ggsave("Figures/fig5_combinedPanels.pdf",fig4, width=10,height=12)
ggsave("Figures/fig5_combinedPanels.png",fig4, width=10,height=12)
ggsave("Figures/fig5_combinedPanels.tiff",fig4, width=10,height=12)
```
## Policy
```{r}
policy<-read.csv("policy/policy_all_results.csv")
head(policy)
straRef<-data.frame(nstrains=ceiling(10^(seq(0.7,2.75,0.15))),rawStrains = 10^(seq(0.7,2.75,0.15)))
policy<-policy%>%left_join(straRef,by="nstrains")
policy$treatment<-factor(policy$d1,levels=c(0,0.05,0.2),labels=c("No Treatment","Treatment: 63.2%","Treatment: 98.2%"))
policy <- policy %>% mutate(case = paste(cloneCost,mixCost,sep=";"))
policy<-policy%>%mutate(PR = round(PR_NU+PR_ND+PR_SU+PR_SD+PR_AU+PR_AD),
PW = round(PW_NU+PW_ND+PW_SU+PW_SD+PW_AU+PW_AD),
P = round(PR+PW),
N = NU+ND+SU+SD+AU+AD,
resistP = PR/P,
adjB = bavgnull*(0.5819519*(seasonality>0) + (seasonality==0)),
rNU = (PW_NU+PR_NU)/(NU+1),
rND = (PW_ND+PR_ND)/(ND+1),
rSU = (PW_SU+PR_SU)/(SU+1),
rSD = (PW_SD+PR_SD)/(SD+1),
rAU = (PW_AU+PR_AU)/(AU+1),
rAD = (PW_AD+PW_AD)/(AD+1),
I_NU = 1-exp(-rNU),
I_ND = 1-exp(-rND),
I_SU = 1-exp(-rSU),
I_SD = 1-exp(-rSD),
I_AU = 1-exp(-rAU),
I_AD = 1-exp(-rAD),
adjPrev = (I_NU*NU+I_ND*ND+I_SU*SU+I_SD*SD+I_AU*AU+I_AD*AD)/N
)
```
```{r}
subDat<- oneLR %>% filter(wildOnly==0,d1>0,mixCost==0.9,cloneCost==0.1,ResistanceEfficacy==0.006666667) %>%
group_by(d1,nstrains) %>%
arrange(adjPrev,desc=T) %>%
slice_max(adjPrev, prop=0.3) %>%
dplyr::select(d1,nstrains,adjB) %>%
mutate(prev_order = seq_along(nstrains))
policy <- policy %>% left_join(subDat)
```
```{r}
addRateName <- function(string){
paste("Daily Treatment Rate:",string)
}
policy<-policy %>% group_by(No,mixCost,cloneCost,varyCost,d1) %>%
mutate(changeRate = (resistP-lead(resistP))/(resistP))
p26<-ggplot(policy %>%filter(d1>0,time==0,varyCost==1),
aes(adjB,rawStrains))+
geom_tile(aes(fill=changeRate))+
xlab("Transmission Potential") + ylab("Number of Strains")+
scale_x_continuous(trans = "log10")+
scale_y_continuous(trans = "log10")+
scale_fill_viridis_c(option="magma",name = "Percentage of\nreduction\nin resistance")+
facet_grid(case~treatment)+
theme_cowplot()+
theme(legend.position="right", legend.key.height = unit(1,"cm"))
ggsave("Figures/figS5_policyRaster.tiff",p26, width=8,height=12)
ggsave("Figures/figS5_policyRaster.png",p26, width=8,height=12)
ggsave("Figures/figS5_policyRaster.pdf",p26, width=8,height=12)
```
```{r}
addTP <- function(string){
paste("TP: ",round(as.numeric(string),2),sep="")
}
uniqAdjb<- unique(policy$adjB)
ggplot(policy%>%filter(mixCost==0.9,cloneCost==0.1,varyCost==1,adjB %in% uniqAdjb[seq(3,30,3)],P>10,resistP<=1))+
aes(time/365, resistP, group=nstrains, color=nstrains)+
geom_point(size=2)+scale_color_viridis_c(option="magma",trans="log10",name = "Numberof\nStrains")+
geom_path(aes(group=nstrains))+
xlab("Year")+ylab("Frequency of Resistance")+
scale_x_continuous( breaks=breaks_pretty(n = 5))+
facet_grid(adjB~treatment,labeller=labeller(adjB=addTP))+
theme_cowplot()+
theme(
legend.title = element_text(size=10),
legend.text = element_text(size=10))
```
```{r}
meanRes <- policy%>% filter(d1>0,prev_order<10,
resistP<=1)%>%
group_by(time,mixCost,cloneCost,varyCost,treatment,nstrains)%>%
summarize(ss = n(),
meanRes = mean(resistP,na.rm=T),
sdRes = sd(resistP,na.rm=T))
p27<-ggplot(meanRes%>%filter(mixCost==0.9,cloneCost==0.1,varyCost==1))+
aes(time/365, meanRes, group=nstrains, color=nstrains)+
geom_point(size=2)+scale_color_viridis_c(option="magma",trans="log10",name = "Numberof\nStrains")+
geom_path(aes(group=nstrains))+
xlab("Year")+ylab("Frequency of Resistance")+
scale_x_continuous( breaks=breaks_pretty(n = 5))+
facet_grid(.~treatment)+
theme_cowplot()+
theme(
legend.title = element_text(size=10),
legend.text = element_text(size=10))
ggsave("Figures/fig5_policy.tiff",p27, width=8,height=4)
ggsave("Figures/fig5_policy.png",p27, width=8,height=4)
ggsave("Figures/fig5_policy.pdf",p27, width=8,height=4)
```
## Trajectory of resistance change after policy change
Policy in Kenya: CQ (-1998), SP(1998-2004), ACT(2004-) Hemming-Schroeder et al. (2018)
Policy in Ghana: CQ (-2004), ACT (2004-) Flegg et al.
Policy in Cambodia: mefloquine 1985, ACT(2000-) Delacollette
Policy in Thailand: SP 1973, quinine-tetra 1980,mefloquine 1985, ATC (1995-) Delacollette,Rasmussen
Brazil: CQ, before 1990, Quinine+tetracycline 1990-2006, 2007 ACT, Gama
PNG: Nsanzabana
```{r}
drugUsage<-read_excel("DrugUsage.xlsx",sheet="DrugPolicy")
colorScheme<-read_excel("DrugUsage.xlsx",sheet="colorScheme")
drugUsage$FirstLineDrug<-factor(drugUsage$FirstLineDrug,levels=colorScheme$FirstLineDrug)
drugUsage <- drugUsage %>% left_join(colorScheme)
anno<-read_excel("DrugUsage.xlsx",sheet="Annotate")
```
```{r}
a<-a %>% mutate(siteStudyUnique = str_c(`study Id`,lon,lat,sep="_"))
ccc<- a %>% filter(country %in% c("Kenya","Ghana","Cambodia","Thailand","Papua New Guinea","Brazil"),study.start>1980,tested>20) %>%
arrange(siteStudyUnique, study.start)
p13<-ggplot()+
geom_rect(data=drugUsage, aes(ymin=ymin,ymax=ymax,xmin=xmin,xmax=xmax,fill=pfcrt),alpha=0.7)+
geom_rect(data=anno, aes(ymin=annolow,ymax=annohigh,xmin=xmin,xmax=xmax),
color="black",fill=NA)+
geom_text(data = anno, aes(wordsx,wordsy, label=FirstLineDrug) )+
geom_path(data=ccc, aes((study.start+study.end)/2, resisP5, group=siteStudyUnique),color="grey20",size=0.5,linetype=2)+
#geom_jitter(data=ccc, aes((study.start+study.end)/2, resisP5,fill=Prevalence),shape=21, size=3,alpha=0.5)+
geom_point(data=ccc, aes((study.start+study.end)/2, resisP5),color="black",fill="black",shape=21, lwd=3,alpha=0.5)+
scale_fill_manual(values=c('#91bfdb','#ffffbf','#fc8d59'),name =expression(paste("selects for ", italic("pfcrt")," 76T")))+
scale_x_continuous(breaks=breaks_pretty(n=10))+
scale_y_continuous(breaks=c(0,0.25,0.5,0.75,1))+
labs(x="Year", y = "Frequency of Resistant Genotypes")+
facet_wrap(Continent~country,ncol=2,labeller = label_wrap_gen(multi_line=FALSE,width=40))+
#theme_clean(base_size=16)+
theme_cowplot()+
theme(legend.position = "top",axis.title.y = ggtext::element_markdown(),
legend.key=element_rect(colour="black"),
#legend.text = element_text(size=14)
panel.grid.major.y = element_line(linetype = "dotted",color="grey")
)
ggsave("fig6_empDrug_Traject_policy.pdf",p13,width=10,height=8)
ggsave("fig6_empDrug_Traject_policy.tiff",p13,width=10,height=8)
ggsave("fig6_empDrug_Traject_policy.png",p13,width=10,height=8)
```
## Comparison to a GI model
```{r}
GI<-read_csv("../GIonly/v13_GImean_results.csv")
GI<-GI%>%mutate(
m_N = P_acc_N/(NU+ND+1e-6),m_S =P_acc_S/(SU+SD+1e-6), m_A = P_acc_A/(AU+AD+1e-6),
PW = PW_NU+PW_SU+PW_AU+PW_ND+PW_SD+PW_AD,
PR = PR_NU+PR_SU+PR_AU+PR_ND+PR_SD+PR_AD,
P=PW+PR,N = NU+ND+SU+SD+AU+AD,
resistP = PR/P,
adjB = bavgnull*(0.5819519*(seasonality>0) + (seasonality==0)),
rNU = (PW_NU+PR_NU)/(NU+1e-6),
rND = (PW_ND+PR_ND)/(ND+1e-6),
rSU = (PW_SU+PR_SU)/(SU+1e-6),
rSD = (PW_SD+PR_SD)/(SD+1e-6),
rAU = (PW_AU+PR_AU)/(AU+1e-6),
rAD = (PW_AD+PW_AD)/(AD+1e-6),
I_NU = 1-exp(-rNU),
I_ND = 1-exp(-rND),
I_SU = 1-exp(-rSU),
I_SD = 1-exp(-rSD),
I_AU = 1-exp(-rAU),
I_AD = 1-exp(-rAD),
NN = NU+ND,
S = SU+SD,
A = AU+AD,
U = NU+SU+AU,
D = SD+AD+ND,
adjPrev = (I_NU*NU+I_ND*ND+I_SU*SU+I_SD*SD+I_AU*AU+I_AD*AD)/N
)
GI<-GI%>%left_join(straRef,by="nstrains")
```
### Prevalence
```{r}
GI <- GI %>% mutate(case = paste(cloneCost,mixCost,sep=";"))
p10<-ggplot(GI%>%filter(wildOnly==0,d1>0, varyCost==1,P>1) %>%
arrange(case,adjB),
aes(adjPrev,resistP))+
geom_point(aes(color=case))+
geom_path(alpha=0.6,aes(color=case))+
#scale_colour_tableau("Classic Blue-Red 12")+
#scale_color_manual(values=c('#fd8d3c','#e6550d','#a63603','#969696','#636363','#252525','#c6dbef','#6baed6','#08519c'),name=expression(paste("s"[single],";s"[mixed])))+
scale_color_manual(values=c('#fdcc8a','#fc8d59','#d7301f','#cccccc','#969696','#525252',
'#bdc9e1','#74a9cf','#0570b0'),name=expression(paste("s"[single],";s"[mixed])))+
#scale_color_manual(values=c('#d7191c','#fdae61','#74add1'))+
#geom_point(aes(color=nstrains))+
#scale_color_viridis_c(option="turbo",trans="log10",name = "number of\nstrains")+
#geom_point(aes(color=!is.na(prev_order)))+
#scale_color_viridis_d(option="turbo")+
xlim(0,1.05)+
facet_grid(~d1,
labeller=labeller(d1 =addRateName))+
xlab("Prevalence")+ylab("Resistance Frequency")+
theme_cowplot()+theme(plot.background = element_rect(fill="white"))
#ggsave("figS3.tiff",p10,width=11,height=6)
ggsave2("Figures/figS6_GIonlyModel.pdf",plot=p10,width=10,height=5)
ggsave2("Figures/figS6_GIonlyModel.tiff",plot=p10,width=10,height=5)
ggsave2("Figures/figS6_GIonlyModel.png",plot=p10,width=10,height=5)
```
### Detailed patterns
```{r}
#library(ggpattern)
# Host
NProp<-GI%>%filter(varyCost==1) %>%
select(No,NU,ND,SU,SD,AU,AD,rawStrains,adjB,mixCost,cloneCost,d1,wildOnly) %>%
pivot_longer(c(NU,ND,SU,SD,AU,AD), names_to = "NCompart", values_to = "NCVal")
NProp<-NProp %>% group_by(No,d1,wildOnly) %>% mutate(PrevProp = NCVal/sum(NCVal))
NProp <- NProp %>% mutate(HostImm = str_sub(NCompart,1,1),
DrugTreat = str_sub(NCompart,2,2))
NProp$DrugTreat <- factor(NProp$DrugTreat, levels = c("U", "D"),
labels = c("Untreated", "Treated"))
NProp$HostImm <- factor(NProp$HostImm, levels = c("N","S","A"),
labels = c("G0", "G1", "G2"))
NProp <- NProp %>% mutate(scenario = paste(d1,wildOnly))
NProp$wildOnly <- factor(NProp$wildOnly, levels = c(1,0),labels = c("Wild-type Only", "Resistance Invasion"))
NProp$scenario <- factor(NProp$scenario, levels = c("0 1", "0.05 1","0.2 1","0 0", "0.05 0","0.2 0"),
labels = c("Wild-type Only; No Treatment",
"Wild-type Only; Daily Treatment Rate: 0.05",
"Wild-type Only; Daily Treatment Rate: 0.2",
"Resistant Only; No Treatment",
"Resistant invasion; Daily Treatment Rate: 0.05",
"Resistant invasion; Daily Treatment Rate: 0.2"))
NProp$NCompart <- factor(NProp$NCompart, levels = c("NU","ND","SU","SD","AU","AD"))
#ggplot(NProp%>%filter(cloneCost==0.1,mixCost==0.9)) +
# geom_col(aes(log10(nstrains),PrevProp, fill=NCompart),position = "stack")+
# scale_fill_brewer(palette = "Paired")+
# geom_line(data = oneLR%>%filter(prev_order==1,cloneCost==0.1,mixCost==0.9,varyCost==1),
# aes(log10(nstrains), resistP))+
# geom_line(data = oneLR%>%filter(prev_order==1,cloneCost==0.1,mixCost==0.9,varyCost==1),
# aes(log10(nstrains), adjPrev),color="blue")+
# facet_grid(wildOnly~d1,scales="free_y",labeller = "label_both")
p26<-ggplot(NProp%>%filter(cloneCost==0.1,mixCost==0.9,wildOnly== "Resistance Invasion")) +
geom_col_pattern(aes(x = adjB,y = PrevProp, fill = HostImm, pattern = DrugTreat),
color ="grey30",
size = 0.2,
pattern_fill = "black",
pattern_angle = 45,
pattern_density = 0.1,
pattern_spacing = 0.025,
pattern_key_scale_factor = 0.6) +
geom_line(data = GI%>%filter(cloneCost==0.1,mixCost==0.9,varyCost==1,wildOnly==0),
aes(adjB, adjPrev),linetype = 2, color="blue")+
scale_x_continuous(trans = "log10", name = "Transmission Potential")+
scale_y_continuous(
# Features of the first axis
name = "Fraction of Hosts",
# Add a second axis and specify its features
sec.axis = sec_axis(~., name="Prevalence")
)+
scale_fill_manual(values = c('#74add1','#ffffb2','#fdae61'), name = "Generalized\nImmunity") +
#scale_color_manual(values = c('#74add1','#d7191c','#fdae61')) +
scale_pattern_manual(values = c(Treated = "stripe", Untreated = "none"), name = "Drug")+
facet_grid(.~d1,scales="free_y",labeller = labeller(d1=addRateName))+
guides(pattern = guide_legend(override.aes = list(fill = "white")),
fill = guide_legend(override.aes = list(pattern = "none")))+
theme_cowplot()+
theme(axis.title.y.right = element_text(color="blue"),
axis.text.y.right = element_text(color='blue'))
# theme(axis.title.y.right = element_text(color='#d7191c'))
```
```{r}
#library(ggpattern)
# Parasites
ParaProp<-GI%>%filter(varyCost==1) %>%
select(No,PW_NU,PW_ND,PW_SU,PW_SD,PW_AU,PW_AD, PR_NU,PR_ND,PR_SU,PR_SD,PR_AU,PR_AD,rawStrains,adjB,mixCost,cloneCost,d1,wildOnly,adjPrev) %>%
pivot_longer(c(PW_NU,PW_ND,PW_SU,PW_SD,PW_AU,PW_AD, PR_NU,PR_ND,PR_SU,PR_SD,PR_AU,PR_AD), names_to = "NCompart", values_to = "NCVal")
ParaProp<-ParaProp %>%mutate(Genotype = str_sub(NCompart,1,2),
HostImm = str_sub(NCompart,4,4),
DrugTreat = str_sub(NCompart,5,5))
ParaProp<-ParaProp %>% group_by(No,d1,wildOnly) %>% mutate(PrevProp = NCVal/sum(NCVal))
ParaProp$DrugTreat <- factor(ParaProp$DrugTreat, levels = c("U", "D"),
labels = c("Untreated", "Treated"))
ParaProp$HostImm <- factor(ParaProp$HostImm, levels = c("N","S","A"),
labels = c("G0", "G1", "G2"))
ParaProp <- ParaProp %>% mutate(scenario = paste(d1,wildOnly))
ParaProp$wildOnlyChar <- factor(ParaProp$wildOnly, levels = c(1,0),labels = c("Wild-type Only", "Resistance Invasion"))
ParaProp$scenario <- factor(ParaProp$scenario, levels = c("0 1", "0.05 1","0.2 1","0 0", "0.05 0","0.2 0"),
labels = c("Wild-type Only; No Treatment",
"Wild-type Only; Daily Treatment Rate: 0.05",
"Wild-type Only; Daily Treatment Rate: 0.2",
"Resistant Only; No Treatment",
"Resistant invasion; Daily Treatment Rate: 0.05",
"Resistant invasion; Daily Treatment Rate: 0.2"))
#ggplot(ParaProp%>%filter(cloneCost==0.1,mixCost==0.9)) +
# geom_col(aes(log10(nstrains),PrevProp, fill=NCompart),position = "stack")+
# scale_fill_brewer(palette = "Paired")+
# geom_line(data = oneLR%>%filter(prev_order==1,cloneCost==0.1,mixCost==0.9,varyCost==1),
# aes(log10(nstrains), resistP))+
# geom_line(data = oneLR%>%filter(prev_order==1,cloneCost==0.1,mixCost==0.9,varyCost==1),
# aes(log10(nstrains), adjPrev),color="blue")+
# facet_grid(wildOnly~d1,scales="free_y",labeller = "label_both")
p11<-ggplot(ParaProp%>%filter(cloneCost==0.1,mixCost==0.9,wildOnlyChar== "Resistance Invasion")) +
geom_col_pattern(aes(x = adjB,y = PrevProp, fill = HostImm, pattern = DrugTreat),
color ="grey30",
size = 0.2,
pattern_fill = "black",
pattern_angle = 45,
pattern_density = 0.1,
pattern_spacing = 0.025,
pattern_key_scale_factor = 0.6) +
scale_x_continuous(trans = "log10", name = "Transmission Potential")+
ylab("Fraction of Parasites")+
#scale_y_continuous(
# Features of the first axis
# name = "Fraction of Hosts",
# Add a second axis and specify its features
# sec.axis = sec_axis(~., name="Resistance Frequency")
#)+
scale_fill_manual(values = c('#74add1','#ffffb2','#fdae61'), name = "Generalized\nImmunity") +
#scale_color_manual(values = c('#74add1','#d7191c','#fdae61')) +
scale_pattern_manual(values = c(Treated = "stripe", Untreated = "none"), name = "Drug")+
facet_grid(Genotype~d1,scales="free_y",labeller = labeller(d1=addRateName,
Genotype = c(PR="Resistant Parasites", PW="Wild-type Parasites")))+
#guides(pattern = guide_legend(override.aes = list(fill = "white")),
# fill = guide_legend(override.aes = list(pattern = "none")))+
theme_cowplot()+
theme(legend.position = "none")
figS7 <- p26/p11+ plot_annotation(tag_levels = 'A')+plot_layout(heights=c(1,2))
ggsave("figS7_GIComp.tiff", figS7, height = 9, width=12)
ggsave("figS7_GIComp.png", figS7, height = 9, width=12)
ggsave("figS7_GIComp.pdf", figS7, height = 9, width=12)
```
### No Cost
```{r}
straRef<-data.frame(nstrains=ceiling(10^(seq(0.7,2.75,0.15))),rawStrains =10^(seq(0.7,2.75,0.15)))
oneLR<-read.csv("NoCost/noCost_allResults.csv")
oneLR<-oneLR%>%mutate(PW_NU=round(PW_NU),
PW_SU=round(PW_SU),
PW_AU=round(PW_AU),
PW_ND=round(PW_ND),
PW_SD=round(PW_SD),
PW_AD=round(PW_AD),
PR_NU=round(PR_NU),
PR_SU=round(PR_SU),
PR_AU=round(PR_AU),
PR_ND=round(PR_ND),
PR_SD=round(PR_SD),
PR_AD=round(PR_AD))
oneLR<-oneLR%>%mutate(
m_N = P_acc_N/(NU+ND+1e-6),m_S =P_acc_S/(SU+SD+1e-6), m_A = P_acc_A/(AU+AD+1e-6),
PW = PW_NU+PW_SU+PW_AU+PW_ND+PW_SD+PW_AD,
PR = PR_NU+PR_SU+PR_AU+PR_ND+PR_SD+PR_AD,
P=PW+PR,N = NU+ND+SU+SD+AU+AD,
resistP = PR/P,
adjB = bavgnull*(0.5819519*(seasonality>0) + (seasonality==0)),
rNU = (PW_NU+PR_NU)/(NU+1e-6),
rND = (PW_ND+PR_ND)/(ND+1e-6),
rSU = (PW_SU+PR_SU)/(SU+1e-6),
rSD = (PW_SD+PR_SD)/(SD+1e-6),
rAU = (PW_AU+PR_AU)/(AU+1e-6),
rAD = (PW_AD+PW_AD)/(AD+1e-6),
I_NU = 1-exp(-rNU),
I_ND = 1-exp(-rND),
I_SU = 1-exp(-rSU),
I_SD = 1-exp(-rSD),
I_AU = 1-exp(-rAU),
I_AD = 1-exp(-rAD),
NN = NU+ND,
S = SU+SD,
A = AU+AD,
U = NU+SU+AU,
D = SD+AD+ND,
adjPrev = (I_NU*NU+I_ND*ND+I_SU*SU+I_SD*SD+I_AU*AU+I_AD*AD)/N
)
oneLR<-oneLR %>%mutate(infN = (1-1/nstrains)^m_N,
infS = (1-1/nstrains)^m_S,
infA = (1-1/nstrains)^m_A)
oneLR<-oneLR%>%left_join(straRef,by="nstrains")
```
```{r}
oneLR <- oneLR %>% mutate(scenario = paste(d1,wildOnly))
oneLR$wildOnlyChar <- factor(oneLR$wildOnly, levels = c(1,0),labels = c("Wild-type Only", "Resistance Invasion"))
oneLR$treatment<-factor(oneLR$d1,levels=c(0,0.05,0.2),labels=c("No Treatment","Treatment: 63.2%","Treatment: 98.2%"))
oneLR$resEffect<-factor(oneLR$ResistanceEfficacy,c(0.006666667,0.013071895,0.034657359),
labels = c("Resistance Efficacy: 87.5%","Resistance Efficacy: 77.0%","Resistance Efficacy: 50%"))
oneLR$scenario <- factor(oneLR$scenario, levels = c("0 1", "0.05 1","0.2 1","0 0", "0.05 0","0.2 0"),
labels = c("Wild-type Only; No Treatment",
"Wild-type Only; Treatment: 63.2%",
"Wild-type Only; Treatment: 98.2%",
"Resistant Only; No Treatment",
"Resistant invasion; Treatment: 63.2%",
"Resistant invasion; Treatment: 98.2%"))
subWild<-oneLR%>%filter(wildOnly==1,d1==0,
cloneCost==0,mixCost==0,varyCost==1,ResistanceEfficacy==0.006666667)
p1<-ggplot(subWild, aes(adjB,rawStrains))+geom_tile(aes(fill=adjPrev))+
geom_tile(data=subWild%>%filter(P<5),fill="grey")+
#geom_tile(data = subWild%>%filter(d1==0,prev_order<=10,prev_order>1),color = "grey",fill=NA,linewidth=0.4)+
#geom_tile(data = subWild%>%filter(d1==0,prev_order==1),color = "white",fill=NA,linewidth=0.5)+
#geom_raster(data = oneLR%>%filter(wildOnly==0,d1==0,prev_order<=10),fill = "white",alpha=0.3)+
xlab("Transmission Potential") + ylab("Number of Strains")+
scale_fill_viridis_c(option="magma",name = "Prevalence",
limits = c(0, 1),values=seq(0,1,0.25))+
scale_x_continuous(trans = "log10")+
scale_y_continuous(trans = "log10")+
facet_grid(.~scenario)+
theme_cowplot()+
theme(legend.position="right", legend.key.height = unit(1,"cm"))
```
```{r}
oneLR <- oneLR %>% mutate(case = paste("single|mixed costs = ", cloneCost,";",mixCost,sep=""))
pGI2<-ggplot(oneLR%>%filter(wildOnly==0,d1>0,varyCost==1,ResistanceEfficacy==0.006666667,adjPrev>0.01),
aes(adjB,rawStrains))+
geom_tile(aes(fill=resistP))+
#geom_tile(data = oneLR%>%filter(wildOnly==0,d1>0,varyCost==1,cloneCost == 0.1, mixCost==0.9,prev_order<=10),color = "white",fill=NA,size=0.5)+
xlab("Transmission Potential") + ylab("Number of Strains")+
scale_x_continuous(trans = "log10")+
scale_y_continuous(trans = "log10")+
scale_fill_viridis_c(option="magma",name = "Resistance\nfrequency",
limits = c(0, 1),values=seq(0,1,0.25))+
facet_grid(case~treatment)+
theme_cowplot()+
theme(legend.position="right", legend.key.height = unit(1,"cm"))
pGI3<-ggplot(oneLR%>%filter(wildOnly==0, d1>0,varyCost==1,P>1),
aes(adjPrev,resistP))+
geom_point(aes(color=case))+
#scale_colour_tableau("Classic Blue-Red 12")+
scale_color_manual(values=c('#fdcc8a','#fc8d59','#d7301f','#cccccc','#969696','#525252',
'#bdc9e1','#74a9cf','#0570b0'),name=expression(paste("s"[single],";s"[mixed])))+
#scale_fill_manual(values=c('#f7f4f9','#e7e1ef','#d4b9da','#c994c7','#df65b0','#e7298a','#ce1256','#980043','#67001f'),name=expression(paste("s"[single],";s"[mixed])))+
#scale_color_manual(values=c('#d7191c','#fdae61','#74add1'))+
#geom_point(aes(color=nstrains))+
#scale_color_viridis_c(option="turbo",trans="log10",name = "number of\nstrains")+
#geom_point(aes(color=!is.na(prev_order)))+
#scale_color_viridis_d(option="turbo")+
facet_grid(treatment~resEffect)+
xlab("Prevalence")+ylab("Resistance Frequency")+
theme_cowplot()
```
### GI policy
```{r}
GIpolicy<-read.csv("GI_policy/GI_policy_results.csv")
GIpolicy<-GIpolicy%>%mutate(PW_NU=round(PW_NU),
PW_SU=round(PW_SU),
PW_AU=round(PW_AU),
PW_ND=round(PW_ND),
PW_SD=round(PW_SD),
PW_AD=round(PW_AD),
PR_NU=round(PR_NU),
PR_SU=round(PR_SU),
PR_AU=round(PR_AU),
PR_ND=round(PR_ND),
PR_SD=round(PR_SD),
PR_AD=round(PR_AD))
head(GIpolicy)
#straRef<-data.frame(nstrains=ceiling(10^(seq(0.7,2.75,0.15))),rawStrains = 10^(seq(0.7,2.75,0.15)))
#GIpolicy<-GIpolicy%>%left_join(straRef,by="nstrains")
GIpolicy$treatment<-factor(GIpolicy$d1,levels=c(0,0.05,0.2),labels=c("No Treatment","Treatment: 63.2%","Treatment: 98.2%"))
GIpolicy <- GIpolicy %>% mutate(case = paste(cloneCost,mixCost,sep=";"))
GIpolicy<-GIpolicy%>%mutate(PR = round(PR_NU+PR_ND+PR_SU+PR_SD+PR_AU+PR_AD),
PW = round(PW_NU+PW_ND+PW_SU+PW_SD+PW_AU+PW_AD),
P = round(PR+PW),
N = NU+ND+SU+SD+AU+AD,
resistP = PR/P,
adjB = bavgnull*(0.5819519*(seasonality>0) + (seasonality==0)),
rNU = (PW_NU+PR_NU)/(NU+1),
rND = (PW_ND+PR_ND)/(ND+1),
rSU = (PW_SU+PR_SU)/(SU+1),
rSD = (PW_SD+PR_SD)/(SD+1),
rAU = (PW_AU+PR_AU)/(AU+1),
rAD = (PW_AD+PW_AD)/(AD+1),
I_NU = 1-exp(-rNU),
I_ND = 1-exp(-rND),
I_SU = 1-exp(-rSU),
I_SD = 1-exp(-rSD),
I_AU = 1-exp(-rAU),
I_AD = 1-exp(-rAD),
adjPrev = (I_NU*NU+I_ND*ND+I_SU*SU+I_SD*SD+I_AU*AU+I_AD*AD)/N
)
GIpolicy$resistP[is.nan(GIpolicy$resistP)]<-0
```
```{r}
addRateName <- function(string){
paste("Daily Treatment Rate:",string)
}
p28<-ggplot(GIpolicy%>% filter(d1>0,varyCost==1,mixCost==0.9,cloneCost==0.1,adjPrev>0) )+
aes(time/365, resistP, group=No, color=adjB)+
geom_point(size=2)+scale_color_viridis_c(option="magma",trans="log10",name = "Transmission\nPotential")+
geom_path(aes(group=No))+
xlab("Year")+ylab("Frequency of Resistance")+
scale_x_continuous( breaks=breaks_pretty(n = 5))+
facet_grid(.~treatment)+
theme_cowplot()+
theme(
legend.title = element_text(size=10),
legend.text = element_text(size=10))
```
```{r}
ggsave("Figures/figS8_GIpolicy.tiff",p28, width=8,height=4)
ggsave("Figures/figS8_GIpolicy.png",p28, width=8,height=4)
ggsave("Figures/figS8_GIpolicy.pdf",p28, width=8,height=4)
```