ww_report_template.Rmd

---
knit: (function(input, ...) {
    rmarkdown::render(
      input,
      output_file = paste0("ww_report", '_', Sys.Date(), '.pdf'
      ),
      envir = globalenv()
    )
  })

output: pdf_document

header-includes:
- \usepackage{pdflscape}
- \newcommand{\blandscape}{\begin{landscape}}
- \newcommand{\elandscape}{\end{landscape}}
---

```{r setup, include=FALSE}
#SARS-CoV-2 Wastewater Epidemiology Lab Regional Report Generator
knitr::opts_chunk$set(echo = FALSE, warning = FALSE, message = FALSE)
library(readxl)
library(lubridate)
library(plyr)
library(tidyverse)
library(knitr)
library(kableExtra)
library(googledrive)
```


```{r sampleID}
###### INSERT INFORMATION HERE:
region <- "CSC"
sites <- c("BMI","MH")
site_labels <- c("Bath and Millhaven Combined", "Millhaven")
names(site_labels) <- sites
date <- Sys.Date() #Sys.Date() will give you today's date or use a date in quotes ex. "2021-02-01"
comments <- NULL
authorizedby <- "Robert Delatolla" #your name here inside "quotes"
analysis <- "Solids"
######
```

![](img/Logo_address.png) ![](img/cheologo.png)

```{r Addresses}
addresses.df <- read_excel("data/Client List.xlsx")
address <- filter(addresses.df, Region == region)
address_counter <- 1:nrow(address)

for(i in address_counter){
  if(is.na(address$Name[i]) == TRUE) {
  clientName <- NULL
  clientTitle <- NULL
  clientOffice <- NULL
  clientAddress1 <- NULL
  clientAddress2 <- NULL
  clientAddress3 <- NULL
  } 
  
  else {
  clientName <- address$Name[i]
  clientTitle <- address$Title[i]
  clientOffice <- address$Office_Institution[i]
  clientAddress1 <- address$Address_line_1[i]
  clientAddress2 <- address$Address_line_2[i]
  clientAddress3 <- address$Address_line_3[i]
  }
}
```
`r clientName`  
`r clientTitle`  
`r clientOffice`  
`r clientAddress1`
`r clientAddress2`
`r clientAddress3`


# Wastewater Trend Analysis Report: SARS-CoV-2

```{r variables}
##read file and define variables for text portion of report
drive_download(as_id("1bJzYQ_V63of4Ym15RsiSlT9uA7R9OSth"), path = "data/report_data.csv", overwrite = TRUE)

data.df <- read.csv(file.path(getwd(), "data/report_data.csv"))

datetoday <- Sys.Date()

variant.df <- data.df %>%
  mutate(sampleDate = as.Date(sampleDate)) %>%
  filter(grepl(paste(sites, collapse = "|"), siteID),
         testB117 == TRUE)

detection <- ifelse(last(variant.df$detectB117 == TRUE), "Detected", "Not detected")

```  
Report Date: `r datetoday` &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; Authorized by: `r authorizedby`  
Comments:  
  
\noindent\rule{16cm}{0.4pt}
### Longitudinal data for 5 weeks ending `r date`
### Method Type: `r analysis` RT-qPCR &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; Region: `r region`
### B.1.1.7 variant: `r detection` on `r last(variant.df$sampleDate)`.


## SARS-CoV-2 Detection in Wastewater by CDC N1 & N2 RT-qPCR

```{r plot}
enddate <- as.Date(date)
startdate <- as.Date(enddate - days(35))
#startdate <- as.POSIXct(ymd(enddate) - days(35))

data.df <- filter(data.df, grepl(paste(sites, collapse = "|"), siteID))

plotdata.df <- filter(data.df,
                      sampleDate <= as.Date(enddate) &
                      sampleDate >= as.Date(startdate))

plotdata.df$N1_N2_avg <- (plotdata.df$covN1_nPMMoV_meanMnNr+
                            plotdata.df$covN2_nPMMoV_meanMnNr)/2

plotdata.df$N1_N2_avg_mL <- plotdata.df$N1_N2_avg_L/1000

if(enddate-startdate < 15)
{
  xbreaks <- "1 day"
}
if(enddate-startdate > 14 | enddate-startdate < 57)
{
  xbreaks <- "1 week"
}
if(enddate-startdate > 56)
{
  xbreaks <- "1 month"
}

plotdata.df$confidence <- ifelse(
  is.na(plotdata.df$viral_confidence), "0",
  ifelse(plotdata.df$viral_confidence == FALSE, "1", "2"))

plotdata.df$norm_confidence <- ifelse(
  (is.na(plotdata.df$N1_N2_avg)) & (plotdata.df$qualityFlag == FALSE), "0",
  ifelse(
      ((plotdata.df$N1_N2_avg == 0) & (plotdata.df$qualityFlag == FALSE)), "0",
  ifelse(
    (plotdata.df$qualityFlag == TRUE), "1",
    "2")
  )
  )

N1_N2_plotdata <- filter(plotdata.df, !(is.na(plotdata.df$qualityFlag)))

N1_N2_plotdata$N1_N2_avg_mL <- ifelse(is.na(N1_N2_plotdata$N1_N2_avg_mL), 0, 
                                   N1_N2_plotdata$N1_N2_avg_mL)

N1_N2_plotdata$N1_N2_avg <- ifelse(is.na(N1_N2_plotdata$N1_N2_avg), 0, 
                                   N1_N2_plotdata$N1_N2_avg)

N1_N2_plotdata$sampleDate <- as.POSIXct(N1_N2_plotdata$sampleDate)

N1_N2_plotdata <- arrange(N1_N2_plotdata, siteID, sampleDate)

variant_plotdata <- plotdata.df %>%
  filter(!(is.na(plotdata.df$fractionB117))) %>%
  mutate(pctB117 = fractionB117*100)
variant_plotdata$sampleDate <- as.POSIXct(variant_plotdata$sampleDate)
#variant_plotdata <- filter(variant_plotdata, !is.na(percent))

#variant_plotdata$allelecount <- paste(variant_plotdata$S_SN501Y_WT_alleles, "-", variant_plotdata$S_SN501Y_Var_alleles, sep = "")

N1ymax <- round_any(max(N1_N2_plotdata$N1_N2_avg_mL), 10000, f = ceiling)

viralload_plot <- ggplot(data = N1_N2_plotdata, aes(x = sampleDate, y = N1_N2_avg_mL))+
  theme_classic()+
  geom_point(aes(shape = confidence, colour = confidence), size = 3.0)+
  geom_line(aes(group = 1), stat = "identity", size = 1, color = "black")+
  geom_hline(yintercept = 1, color = "darkorange2")+
  scale_shape_manual(values = c("0" = 18, "1" = 16, "2" = 16),
                     labels = c("0" = "Non detection", "1" = "Low confidence", "2" = "High Confidence"),
                     name = "Confidence:")+
  scale_color_manual(values = c("0" = "black", "1" = "lightcoral", "2" = "seagreen3"),
                     labels = c("0" = "Non detection", "1" = "Low confidence", "2" = "High Confidence"),
                     name = "Confidence:")+
  scale_x_datetime(date_breaks = waiver(), date_minor_breaks = "1 day", 
                   date_labels = "%Y %b %d")+
  scale_y_continuous(trans = 'log10', limits = c(1, N1ymax))+
  labs(x = "Date Sampled", y = "Viral Load (cp/mL)")+
  facet_wrap(vars(siteID), labeller = labeller(siteID = site_labels))+
  theme(plot.title = element_text(size = 16, hjust = 0.5),
        legend.position = "top", 
        legend.direction = "horizontal",
        axis.text.x = element_text(angle = 90, hjust = 1),
        panel.border = element_rect(colour = "black", fill = NA),
        strip.background = element_blank())

nrm_copies_plot <- ggplot(data = N1_N2_plotdata, aes(x = sampleDate, y = N1_N2_avg))+
  theme_classic()+
  geom_point(aes(shape = norm_confidence, colour = norm_confidence), size = 3.0)+
  geom_line(aes(group = 1), stat = "identity", size = 1, color = "black")+
  scale_shape_manual(values = c("0" = 18, "1" = 16, "2" = 16),
                     labels = c("0" = "Non detection", "1" = "Low confidence", "2" = "High Confidence"),
                     name = "Confidence:")+
  scale_color_manual(values = c("0" = "black", "1" = "lightcoral", "2" = "seagreen3"),
                     labels = c("0" = "Non detection", "1" = "Low confidence", "2" = "High Confidence"),
                     name = "Confidence:") +
  scale_x_datetime(date_breaks = waiver(), date_minor_breaks = "1 day", 
                   date_labels = "%Y %b %d")+
  #scale_y_continuous(trans = 'log10', limits = c(1, N1ymax))+
  labs(x = "Date Sampled", y = "Normalized viral copies")+
  facet_wrap(vars(siteID), labeller = labeller(siteID = site_labels))+
  theme(plot.title = element_text(size = 16, hjust = 0.5),
        legend.position = "top", 
        legend.direction = "horizontal",
        axis.text.x = element_text(angle = 90, hjust = 1),
        panel.border = element_rect(colour = "black", fill = NA),
        strip.background = element_blank())

plot(viralload_plot)
```

Low confidence: SARS-CoV-2 signal was detected below the method’s limit of quantification.

High confidence: SARS-CoV-2 signal was detected and meets QA/QC requirements. 

Orange line represents the limit of quantification of the assay

```{r}
plot(nrm_copies_plot)
```

```{r variant}
variantplot <- ggplot(data = variant_plotdata, aes(x = sampleDate, y = pctB117))+
  theme_classic()+
  geom_point(size = 3.0)+
  geom_line(stat = "identity", size = 1)+
  #scale_color_discrete(labels = c("S_SN501Y_percent" = "N501Y", "S_Sdel_percent" = "Sdel"), name = "Allele:")+
  scale_x_datetime(date_breaks = waiver(), date_minor_breaks = "1 day", 
                   date_labels = "%Y %b %d")+
  scale_y_continuous(limits = c(0, 100))+
  labs(x = "Date Sampled", y = "Percent Variant (%)")+
  facet_wrap(vars(siteID))+
  theme(plot.title = element_text(size = 16, hjust = 0.5),
        legend.position = "top", 
        legend.direction = "horizontal",
        axis.text.x = element_text(angle = 90, hjust = 1),
        panel.border = element_rect(colour = "black", fill = NA),
        strip.background = element_blank())

#plot(variantplot)

```

NOTE:  
- This report is the only valid reference in case of differences with other transmitted documents.  
- The analysis results refer only to what was provided for testing.  
- The information provided in this report is for research use only.  

See individual sample report for variances and limitations to sample analysis.  
  
&nbsp;  

Robert Delatolla, Kamya Bhatnagar, Patrick M. D’Aoust, Élisabeth Mercier,
Alex MacKenzie, and Tyson Graber

Department of Civil Engineering, University of Ottawa & Children’s Hospital
of Eastern Ontario – Research Institute

\newpage

```{r spreadsheet}
# Upload ww_reports
drive_download(as_id("157D3pWH0r47xx4-rNua_Djbc606r_9vM"), path = "data/ww_reports.xlsx", overwrite = TRUE)

ww_reports <- read_excel(path = "data/ww_reports.xlsx", sheet = "7 - WWMeasure") %>%
  filter(grepl(paste(sites, collapse = "|"), sampleID)) %>%
  filter(!is.na(uWwMeasureID), !is.na(value)) %>%
  select(sampleID, labID, analysisDate, reportDate, type, value) %>%
  filter(type == "covN1" | type == "covN2") %>%
  group_by(sampleID) %>%
  mutate(mean = mean(value)) %>%
  select(-type, -value) %>%
  distinct() %>%
  mutate(viral_cp_mL = mean/1000) %>%
  select(-mean)

ww_reports$submitted_id <- NA
ww_reports$arr_tmp <- NA
ww_reports$condition <- NA
ww_reports$pcr_date <- NA

report.df <- rename(ww_reports,
                     "Sample ID" = sampleID,
                     "Submitted ID" = submitted_id,
                     "Location" = labID,
                     "Date Sampled" = analysisDate,
                     "Date Received" = reportDate,
                     "Arrival Temp" = arr_tmp,
                     "Sample Condition" = condition,
                     "qPCR Date" = pcr_date,
                     "Average N1 & N2 Copies/mL" = viral_cp_mL)

col_order <- c("Sample ID", "Submitted ID", "Location", "Date Sampled",
               "Date Received", "Arrival Temp", "Sample Condition", "qPCR Date",
               "Average N1 & N2 Copies/mL")

report.df <- report.df[,col_order]

report.df <- filter(report.df, as.Date(`Date Sampled`) >= as.Date(startdate),
                    as.Date(`Date Sampled`) <= as.Date(enddate))   
report.df <- arrange(report.df, `Sample ID`, desc(`Date Sampled`)) 
#report.df$`Average N1 Copies/mL` <- round(report.df$`Average N1 Copies/mL`, 2)
                                              
kable(report.df, booktabs = T, align = "c") %>%
kable_styling(latex_options = "scale_down") %>%
#kableExtra::landscape() %>%
row_spec(0, angle = 45)

#16cm line for portrait, 23cm line for landscape
```

\noindent\rule{16cm}{0.4pt}  

End of report