## Background

Pertussis (whooping cough) is a common lung infection cause by the bacteria B. Pertussis. It can affect anyone but it is most deadly fro infants (under 1 year of age)


## CDC trackind data 
The CDC tracks the number of Pertussis cases:

```{r, echo=FALSE}
cdc <- data.frame(year = c(1922L,1923L,1924L,1925L,
                                          1926L,1927L,1928L,1929L,1930L,1931L,
                                          1932L,1933L,1934L,1935L,1936L,
                                          1937L,1938L,1939L,1940L,1941L,1942L,
                                          1943L,1944L,1945L,1946L,1947L,
                                          1948L,1949L,1950L,1951L,1952L,
                                          1953L,1954L,1955L,1956L,1957L,1958L,
                                          1959L,1960L,1961L,1962L,1963L,
                                          1964L,1965L,1966L,1967L,1968L,1969L,
                                          1970L,1971L,1972L,1973L,1974L,
                                          1975L,1976L,1977L,1978L,1979L,1980L,
                                          1981L,1982L,1983L,1984L,1985L,
                                          1986L,1987L,1988L,1989L,1990L,
                                          1991L,1992L,1993L,1994L,1995L,1996L,
                                          1997L,1998L,1999L,2000L,2001L,
                                          2002L,2003L,2004L,2005L,2006L,2007L,
                                          2008L,2009L,2010L,2011L,2012L,
                                          2013L,2014L,2015L,2016L,2017L,2018L,
                                          2019L,2020L,2021L,2022L,2023L, 2024L, 2025L),
         cases = c(107473,164191,165418,152003,
                                          202210,181411,161799,197371,
                                          166914,172559,215343,179135,265269,
                                          180518,147237,214652,227319,103188,
                                          183866,222202,191383,191890,109873,
                                          133792,109860,156517,74715,69479,
                                          120718,68687,45030,37129,60886,
                                          62786,31732,28295,32148,40005,
                                          14809,11468,17749,17135,13005,6799,
                                          7717,9718,4810,3285,4249,3036,
                                          3287,1759,2402,1738,1010,2177,2063,
                                          1623,1730,1248,1895,2463,2276,
                                          3589,4195,2823,3450,4157,4570,
                                          2719,4083,6586,4617,5137,7796,6564,
                                          7405,7298,7867,7580,9771,11647,
                                          25827,25616,15632,10454,13278,
                                          16858,27550,18719,48277,28639,32971,
                                          20762,17972,18975,15609,18617,
                                          6124,2116,3044,7063, 22538, 21996)
       )

```

> Q. I want a plot of year vs cases 

```{r}
library(ggplot2)
ggplot(cdc)+ aes(year,cases) + geom_point() + geom_line()
```

> Q. Add annotation lines for the major miletstone of wP vaccination roll-out (1946) and the switch to the aP vaccine (1996).

```{r}
ggplot(cdc) +
  aes(year, cases) +
  geom_point() +
  geom_line() + 
  geom_vline(xintercept = 1946, col= "blue", lty=2) +
  geom_vline(xintercept = 1996, col= "red", lty= 2) + 
  geom_vline(xintercept =2020, col="gray", lty=2)
```

Cases decreased dramatically after 1946, largely due to the introduction and widespread use of vaccines. However, cases began to rise again in the early 2000s. One possible reason is that the newer vaccine may not provide protection that lasts as long or is as strong as the older version. Over time, the bacteria may also have evolved, making it better able to survive despite vaccination. In addition, improvements in diagnostic testing and disease surveillance have made it easier to detect and report cases that may have previously gone unnoticed. Another contributing factor could be declining vaccination rates in some populations. The data also shows a correlation between peaks in cases and countries that use this particular vaccine, although the newer vaccine tends to show smaller peaks compared to earlier outbreaks.


## Exploring CMI-PD Poject

The CMI-PD project < https://www.cmi-pb.org > mission is to provide the scientific community with a comprehensive, high-quality and freely accessible resource of Pertussis booster vaccination.

Basically, make available a large dataset on the immune response to Pertussis. They used a "booster" vaccination as a proxy for Pertussis infection

They make their data available as JSON format API. We can read this into R the `read_json()` function from the **jsonlite** package: 

```{r}
library(jsonlite)
subject <- read_json("https://www.cmi-pb.org/api/v5_1/subject", simplifyVector=TRUE)
head(subject)
```

> Q. How many aP and wP individuals are there in this `subject` table?

```{r}
table(subject$infancy_vac)
```

> Q. How many male/female are there?

```{r}
table(subject$biological_sex)
```
 > What is the breakdown of `biological_sex` race subjects?
 
> Is this representtaive of the US population?

```{r}
table(subject$race, subject$biological_sex)
```

We can now read more tales form the CMI-PB database 

```{r}
specimen <- read_json("https://www.cmi-pb.org/api/v5_1/specimen", simplifyVector=TRUE)
ab_titer <- read_json("https://www.cmi-pb.org/api/v5_1/plasma_ab_titer", simplifyVector=TRUE)
```

```{r}
head(specimen)
```

```{r}
head(ab_titer)
```

To make sense of all this data we need to "join" (a.k.a "merge" or "link") all these tables together. Only then will you know that a give Ab measurement (from the `ab-titer` table) was collected on a certain data from a certain date (from the `specimen` table) wP or aP subject (from the `subject` table)

We can use **dplyr** and the `*join()` fmaily of functions to do this. 

```{r}
library(dplyr)
meta <- inner_join(subject, specimen)
head(meta)
```

Let's do one more `inner_join` to join the `ab_titer` with all our `meta` data.  

```{r}
abdata <- inner_join(ab_titer, meta)
head (abdata)
```

> Q. How many different Ab "isotype" values are in this dataset?

```{r}
table(abdata$isotype)
```

> Q. How many different "antigen" values are measured? 

```{r}
table(abdata$antigen)
```

Let's focus on IgG isotype

```{r}
igg <- abdata |> filter(isotype=="IgG")
head(igg)
```

Make a plot `MFI_normalised` values for all `antigen` values 

ˆ
```{r}
ggplot(igg)+ aes(MFI_normalised, antigen) + geom_boxplot()
```

The antigens "PT", "FIM2/3", and "FHA" appear to have the wide range of values.

> Q. Is there a difference for these repsonses between aP and wP individuals? 

```{r}
ggplot(igg)+ aes(MFI_normalised, antigen, col=infancy_vac) + geom_boxplot() 
```

```{r}
ggplot(igg)+ aes(MFI_normalised, antigen) + geom_boxplot() + facet_wrap(~infancy_vac)
```

> Q. Is there a difference with time (i.e. before booster shot vs after booster shot)? 

```{r}
ggplot(igg)+ aes(MFI_normalised, antigen, col=infancy_vac) + geom_boxplot() + facet_wrap(~visit)
```


```{r}
# Filter 2021 dataset 
abdata.21 <- abdata %>% filter(dataset == "2021_dataset")

abdata.21 %>% 
  filter(isotype == "IgG", antigen == "PT") %>%
  ggplot(aes(x = planned_day_relative_to_boost,
             y = MFI_normalised,
             col = infancy_vac,
             group = subject_id)) +
  
    # Individual subject trajectories (thin, faint)
    geom_line(alpha = 0.25, linewidth = 0.6) +
    geom_point(alpha = 0.4, size = 1.5) +
    
    # Smooth bold average trendlines for aP and wP
    stat_summary(fun = mean,
                 geom = "smooth",
                 se = FALSE,
                 linewidth = 2.8,
                 span = 0.5,
                 aes(group = infancy_vac, col = infancy_vac)) +
    
    # Reference lines
    geom_vline(xintercept = 0, linetype = "dashed") +
    geom_vline(xintercept = 14, linetype = "dashed") +
    
    labs(title = "CMI-PB 2021 dataset IgG PT",
         subtitle = "Dashed lines at day 0 (pre boost) and day 14 (post boost)",
         x = "Day relative to boost",
         y = "Normalised MFI") +
    
    theme_minimal(base_size = 14)


```
















