https://github.com/morgankain/RRV_HostVectorCompetence
Raw File
Tip revision: be7e87c3c4c8af0420a8dd42cdcff5586fdbad90 authored by Morgan Kain on 25 May 2021, 16:23:11 UTC
Merge pull request #1 from morgankain/add-license-1
Tip revision: be7e87c
13_data.community.R
################################################################################
### Set up the mosquito and host communities for competence and R0 estimates ###
################################################################################

## Note: Because we do not have mosquito physiological responses and blood meals for all of the species in the community
 ## we simply scale up the mosquitoes that we do have these data for to assume that they are the whole mosquito community.
  ## For Brisbane the mosquitoes we do have data for are ~ 90% of the community, and for Cairns a bit over 70%

## Convert to proportions
mosq_prop <- mosq_prop %>% mutate(prop = count / sum(count))

## Pull out only the mosquito species that we have BOTH physiological data and mosquito biting data
mosq_prop_for_R0 <- mosq_prop %>% 
  filter(mos_species %in% as.character(unique(h_to_m_trans.pred$species))) %>%
  filter(mos_species %in% as.character(mosq_blood.prop$mos_species))

## Re-scale
mosq_prop_for_R0$prop <- mosq_prop_for_R0$prop / sum(mosq_prop_for_R0$prop)

## All the hosts that we have titer data for
physiology_species.h <- as.character(unique(tot_titer$host))
physiology_species.m <- as.character(unique(h_to_m_emp$species))

## All the hosts that show up in the community of interest
ecological_species.h <- as.character(unique(host_prop$host_species))
ecological_species.m <- as.character(unique(mosq_blood.prop$mos_species))

## Connections here between hosts with physiological responses and ecological data are the following:
 ## Note, some simplifications are needed here to align the various data sources that we have
  ## A potentially minor problem with the simple average taken here is that titer profiles aren't getting weighted
   ## appropriately by abundance (e.g. within marsupials); however we lack that level of detailed abundance data

 ## bushtail possum -> possum
 ## agile wallaby   -> within marsupial 
 ## grey kangaroo   -> within marsupial 
 ## bandicoot       -> dropped
 ## marsupial mouse -> dropped 
 ## horse           -> horse
 ## sheep           -> sheep
 ## pig             -> dropped because they are not allowed in urban settings
 ## cow             -> cattle
 ## grey h f fox    -> flying fox
 ## rat             -> rat
 ## rabbit          -> rabbit
 ## chicken         -> within bird
 ## black duck      -> within bird 
 ## little corella  -> within bird
 ## human           -> human

#####
## Steps to adjust h_to_m_trans, which already has titer incorporated
#####
 ## 1) Select only the mosquitoes in the community for which we have biting preferences
mosq.which        <- which(dimnames(h_to_m_trans_all_samps)[[3]] %in% ecological_species.m)
mosq.ordered      <- dimnames(h_to_m_trans_all_samps)[[3]][mosq.which]
mosq_prop_for_R0  <- mosq_prop_for_R0[match(mosq.ordered, as.character(mosq_prop_for_R0$mos_species)), ]

 ## 2) Average marsupials (actually macropods given the removal of bandicoot and marsupial mouse)
marsupials        <- c("agile_wallaby", "grey_kangaroo")
marsupials.data   <- h_to_m_trans_all_samps[, which(dimnames(h_to_m_trans_all_samps)[[2]] %in% marsupials), , ]
marsupials.data.s <- apply(marsupials.data, c(1, 3, 4), mean)
marsupials.data.s <- marsupials.data.s[, mosq.which, ]
## Do the same manipulation for AUC for the first panel of Figure 2: host. Not used in the rest of the model
marsupials.data.AUC   <- host_titer_AUC_all_samps[which(dimnames(host_titer_AUC_all_samps)[[1]] %in% marsupials), ]
marsupials.data.s.AUC <- colMeans(marsupials.data.AUC)

 ## 3) Average birds
birds             <- c("chicken", "black_duck", "little_corella")
birds.data        <- h_to_m_trans_all_samps[, which(dimnames(h_to_m_trans_all_samps)[[2]] %in% birds), , ]
birds.data.s      <- apply(birds.data, c(1, 3, 4), mean)
birds.data.s      <- birds.data.s[, mosq.which, ]
## Do the same manipulation for AUC for the first panel of Figure 2: host. Not used in the rest of the model
birds.data.AUC   <- host_titer_AUC_all_samps[which(dimnames(host_titer_AUC_all_samps)[[1]] %in% birds), ]
birds.data.s.AUC <- colMeans(birds.data.AUC)

h_to_m_trans_adj_to_com                <- array(dim = c(length(inf_days)
                                                        ## total length - summaries + cat and dog which don't have data because they never show titer
                                                          , length(ecological_species.h)
                                                          , length(mosq.which)))
h_to_m_trans_all_samps_adj_to_com      <- array(dim = c(length(inf_days)
                                                        ## total length - summaries + cat and dog which don't have data because they never show titer
                                                        , length(ecological_species.h)
                                                        , length(mosq.which), n_samps))
dimnames(h_to_m_trans_adj_to_com)[[2]]           <- ecological_species.h
dimnames(h_to_m_trans_all_samps_adj_to_com)[[2]] <- ecological_species.h

## Create a smaller matrix for AUC as well
host_titer_AUC_all_samps_adj_to_com                <- matrix(nrow = length(ecological_species.h), ncol = n_samps, data = 0)
dimnames(host_titer_AUC_all_samps_adj_to_com)[[1]] <- ecological_species.h

### And then subset all of the mosquito data to these species
dimnames(h_to_m_trans_all_samps_adj_to_com)[[3]]   <- mosq.ordered

## The summaries and cat and dog
h_to_m_trans_all_samps_adj_to_com[, dimnames(h_to_m_trans_all_samps_adj_to_com)[[2]] == "bird", , ]      <- birds.data.s
h_to_m_trans_all_samps_adj_to_com[, dimnames(h_to_m_trans_all_samps_adj_to_com)[[2]] == "macropod", , ]  <- marsupials.data.s 
h_to_m_trans_all_samps_adj_to_com[, dimnames(h_to_m_trans_all_samps_adj_to_com)[[2]] == "cat", , ]       <- 0
h_to_m_trans_all_samps_adj_to_com[, dimnames(h_to_m_trans_all_samps_adj_to_com)[[2]] == "dog", , ]       <- 0

## The rest of the species (select the estimates for all the mosquitoes we are retaining and make the name change if necessary)
h_to_m_trans_all_samps_adj_to_com[, dimnames(h_to_m_trans_all_samps_adj_to_com)[[2]] == "human", , ]       <- 
  h_to_m_trans_all_samps[, dimnames(h_to_m_trans_all_samps)[[2]] == "human", mosq.which, ]
h_to_m_trans_all_samps_adj_to_com[, dimnames(h_to_m_trans_all_samps_adj_to_com)[[2]] == "possum", , ]      <- 
  h_to_m_trans_all_samps[, dimnames(h_to_m_trans_all_samps)[[2]] == "brushtail_possum", mosq.which, ]
h_to_m_trans_all_samps_adj_to_com[, dimnames(h_to_m_trans_all_samps_adj_to_com)[[2]] == "flying_fox", , ]  <- 
  h_to_m_trans_all_samps[, dimnames(h_to_m_trans_all_samps)[[2]] == "grey_h_f_fox", mosq.which, ]
h_to_m_trans_all_samps_adj_to_com[, dimnames(h_to_m_trans_all_samps_adj_to_com)[[2]] == "cattle", , ]      <- 
  h_to_m_trans_all_samps[, dimnames(h_to_m_trans_all_samps)[[2]] == "cow", mosq.which, ]
h_to_m_trans_all_samps_adj_to_com[, dimnames(h_to_m_trans_all_samps_adj_to_com)[[2]] == "horse", , ]       <- 
  h_to_m_trans_all_samps[, dimnames(h_to_m_trans_all_samps)[[2]] == "horse", mosq.which, ]
h_to_m_trans_all_samps_adj_to_com[, dimnames(h_to_m_trans_all_samps_adj_to_com)[[2]] == "sheep", , ]       <- 
  h_to_m_trans_all_samps[, dimnames(h_to_m_trans_all_samps)[[2]] == "sheep", mosq.which, ]
h_to_m_trans_all_samps_adj_to_com[, dimnames(h_to_m_trans_all_samps_adj_to_com)[[2]] == "rabbit", , ]      <- 
  h_to_m_trans_all_samps[, dimnames(h_to_m_trans_all_samps)[[2]] == "rabbit", mosq.which, ]
h_to_m_trans_all_samps_adj_to_com[, dimnames(h_to_m_trans_all_samps_adj_to_com)[[2]] == "rat", , ]         <- 
  h_to_m_trans_all_samps[, dimnames(h_to_m_trans_all_samps)[[2]] == "rat", mosq.which, ]

## fill in for AUC 
host_titer_AUC_all_samps_adj_to_com[dimnames(host_titer_AUC_all_samps_adj_to_com)[[1]] == "bird", ]       <- birds.data.s.AUC
host_titer_AUC_all_samps_adj_to_com[dimnames(host_titer_AUC_all_samps_adj_to_com)[[1]] == "macropod", ]   <- marsupials.data.s.AUC 
host_titer_AUC_all_samps_adj_to_com[dimnames(host_titer_AUC_all_samps_adj_to_com)[[1]] == "cat", ]        <- 0
host_titer_AUC_all_samps_adj_to_com[dimnames(host_titer_AUC_all_samps_adj_to_com)[[1]] == "dog", ]        <- 0
host_titer_AUC_all_samps_adj_to_com[dimnames(host_titer_AUC_all_samps_adj_to_com)[[1]] == "human", ]      <- 
  host_titer_AUC_all_samps[dimnames(host_titer_AUC_all_samps)[[1]] == "human", ]
host_titer_AUC_all_samps_adj_to_com[dimnames(host_titer_AUC_all_samps_adj_to_com)[[1]] == "possum", ]     <- 
  host_titer_AUC_all_samps[dimnames(host_titer_AUC_all_samps)[[1]] == "brushtail_possum", ]
host_titer_AUC_all_samps_adj_to_com[dimnames(host_titer_AUC_all_samps_adj_to_com)[[1]] == "flying_fox", ] <- 
  host_titer_AUC_all_samps[dimnames(host_titer_AUC_all_samps)[[1]] == "grey_h_f_fox", ]
host_titer_AUC_all_samps_adj_to_com[dimnames(host_titer_AUC_all_samps_adj_to_com)[[1]] == "cattle", ]     <- 
  host_titer_AUC_all_samps[dimnames(host_titer_AUC_all_samps)[[1]] == "cow", ]
host_titer_AUC_all_samps_adj_to_com[dimnames(host_titer_AUC_all_samps_adj_to_com)[[1]] == "horse", ]      <- 
  host_titer_AUC_all_samps[dimnames(host_titer_AUC_all_samps)[[1]] == "horse", ]
host_titer_AUC_all_samps_adj_to_com[dimnames(host_titer_AUC_all_samps_adj_to_com)[[1]] == "sheep", ]      <- 
  host_titer_AUC_all_samps[dimnames(host_titer_AUC_all_samps)[[1]] == "sheep", ]
host_titer_AUC_all_samps_adj_to_com[dimnames(host_titer_AUC_all_samps_adj_to_com)[[1]] == "rabbit", ]     <- 
  host_titer_AUC_all_samps[dimnames(host_titer_AUC_all_samps)[[1]] == "rabbit", ]
host_titer_AUC_all_samps_adj_to_com[dimnames(host_titer_AUC_all_samps_adj_to_com)[[1]] == "rat", ]        <- 
  host_titer_AUC_all_samps[dimnames(host_titer_AUC_all_samps)[[1]] == "rat", ]

#####
## And subset prop_inf for the community of interest
#####
prop_inf_for_R0 <- prop_inf
prop_inf_for_R0 <- prop_inf_for_R0 %>% mutate(
  host = mapvalues(host
                   , from = c(
  "brushtail_possum"  
, "agile_wallaby"     
, "grey_kangaroo"         
, "bandicoot"        
, "marsupial_mouse"      
, "horse"             
, "sheep"              
, "pig"              
, "cow"                 
, "grey_h_f_fox"      
, "rat"                
, "rabbit"            
, "chicken"             
, "black_duck"        
, "little_corella"        
, "human"              
, "cat"                    
, "dog"
) , to = c(
  "possum"  
, "macropod"     
, "macropod"         
, "marsupial"        
, "marsupial_mouse"      
, "horse"             
, "sheep"              
, "pig"              
, "cattle"                 
, "flying_fox"      
, "rat"                
, "rabbit"            
, "bird"             
, "bird"        
, "bird"        
, "human"              
, "cat"                    
, "dog"                     
)
))

## Last adjustment for proportion that become infected from the infection experiments
prop_inf_for_R0 <- prop_inf_for_R0 %>% 
  group_by(host) %>%
  summarize(tot_inf = sum(tot_inf), tot_vir = sum(tot_vir)) %>%
  mutate(num_inf = tot_vir / tot_inf)

prop_inf_for_R0[prop_inf_for_R0$host == "cat" | prop_inf_for_R0$host == "dog", 4] <- 0

prop_inf_for_R0 <- prop_inf_for_R0[match(host_prop$host_species, prop_inf_for_R0$host), ]

## Grab the info for host seroprevalence for Brisbane for use in both Brisbane and Cairns since we do not have this data for Cairns
host_sero <- host_sero[match(prop_inf_for_R0$host, host_sero$host), ]

#################################################################
### Supplemental Titer plot #2: averaged titer curves (to be made, maybe...??)
#####

####
## Plot the probabilities of infection over time from the titre curves and the mosquito infection probabilities scaled
## to consider the proportion of infected hosts that show a viraemic response
####

h_to_m_trans_all_samps_adj_to_com.gg        <- melt(h_to_m_trans_all_samps_adj_to_com)
names(h_to_m_trans_all_samps_adj_to_com.gg) <- c("day", "host", "mosquito", "sample", "probability")
h_to_m_trans_all_samps_adj_to_com.gg        <- h_to_m_trans_all_samps_adj_to_com.gg %>%
  group_by(host, day) %>% 
  summarise(
    lwr = quantile(probability, 0.20)
  , est = quantile(probability, 0.500)
  , upr = quantile(probability, 0.80)
    )

inf_prob_plots <- ggplot(h_to_m_trans_all_samps_adj_to_com.gg, aes(day, est)) + 
  geom_line() + 
  geom_ribbon(aes(x = day, ymin = lwr, ymax = upr), alpha = 0.25) +
  facet_wrap(~host)
back to top