# EMIT Main Quarantine Paper Data Curation - Producing authoritative, publishable dataset for the EMIT Quarantine study
# Program Objective: working with the 8 cleaned datasets to create a single datatable of data to use for analysis and publication:
# Note: The 11 cleaned dfs produced by the "EMIT_Quarantine_Main_Cleaning.R" script can be found in the directory:
# /Users/jbueno/Box Sync/EMIT/EMIT_Data_Analysis_Jake/UK vs. UMD data and analysis/UK Quarantine Study Data and Notes/Curated Data/Cleaned Data"
# Author: Jacob Bueno de Mesquita
# Date: September 21-23, 2018, December 2018; February 2019
# Summary: Merged all of the EMIT Quarantine Data into a single file called "QuarantineMergedData.csv" in the "/Users/jbueno/Box Sync/EMIT/EMIT_Data_Analysis_Jake/EMIT_Quarantine/Curated Data/Analytical Datasets" directory

#### Load required packages and set working directory ####
library(tidyverse)
library(RcppRoll)
library(readxl)
library(knitr)
library(data.table)
library(lubridate)

setwd("/Users/jbueno/Box Sync/EMIT/EMIT_Data_Analysis_Jake/EMIT_Quarantine/")

sessionInfo()

#### Reading in the cleaned datasets and exploring them ####

pcr <- read.csv("Curated Data/Cleaned Data/CDCpcr.csv")
head(pcr)
str(pcr)
summary(pcr)

microneut <- read.csv("Curated Data/Cleaned Data/CDCserology.csv")
head(microneut)
str(microneut)
summary(microneut)

HAI <- read.csv("Curated Data/Cleaned Data/GLASGOWserology.csv")
head(HAI)
str(HAI)
summary(HAI)

VolunteerList <- read.csv("Curated Data/Cleaned Data/DonorRecipientSubjectIDlist.csv")
head(VolunteerList)
str(VolunteerList)
summary(VolunteerList)

Demog <- read.csv("Curated Data/Cleaned Data/QuarantineDemographics.csv")
head(Demog)
str(Demog)
summary(Demog)

Symptoms <- read.csv("Curated Data/Cleaned Data/Symptoms for all volunteers Q1-Q3.csv")
head(Symptoms)
str(Symptoms)
summary(Symptoms)

PhysicalExam <- read.csv("Curated Data/Cleaned Data/PhysicalExam.csv")
head(PhysicalExam)
str(PhysicalExam)
summary(PhysicalExam)

Calendar <- read.csv("Curated Data/Cleaned Data/QuarantineCalendar.csv")
head(Calendar)
str(Calendar)
summary(Calendar)

# UKg2mean <- read.csv("Curated Data/Cleaned Data/UKg2mean.csv")
# head(UKg2mean)
# str(UKg2mean)
# summary(UKg2mean)
# Note that in the UKg2mean df, the pcr positive data for subject 109 on study day 4 for the coarse aerosol fraction has been eliminated (I believe this is because 109 was not considered to be an infected donor - likely because the NPS was negative on all study days or no more than 1 study day). It is not ideal to scrub data like this in the raw datafiles that are read in for analysis. But this raw data file was used directly in the analysis aiming to get at GM and GSD for the aerosol RNA copies in fine and coarse aerosols and email correspondence from Don Milton indicated that the aerosol pcr data for 109 should be ignored, so for the purpose of that analysis this file was edited. It should also be noted that in that analysis of GM and GSD, in the instances where there were PCR replicates for a particular sample that were qualitatively different, i.e., one was positive and the other was negative, the negative was treated as a non-detected sample, and its value was imputed as LOQ * 1/sqrt(2). The LOQ for flu A in the EMIT study was taken to be 2,000 RNA copies/sample. It should also be noted here that in models of aerosol shedding that we may later construct, we will most likely use tobit regression to account for censored data below the quantification limit and thus not rely on the rudimentary LOQ * 1/sqrt(2) method. 

# However, in the UKg2_all df (read-in below), the pcr copies for 109 are kept in the copy.num.replicate variable and in the final.copies.replicate variable. The process is to multiply the copy.num.replicate variable by the calibration factor (cfactor) variable to get out the final.copies.replicate variable. 
# If subject 109 is, in fact, not to be considered as a positive donor, and thus excluded from analysis of aerosol viral shedding, there will be a selection step in the analysis script to handle this. Thus, we will keep the data from 109 as it is in this UKg2_all datafile. 

UKg2_all <- read.csv("Curated Data/Cleaned Data/UKg2_all.csv")
head(UKg2_all)
str(UKg2_all)
summary(UKg2_all)

UKcough <- read.csv("Curated Data/Cleaned Data/UKCoughData.csv")
head(UKcough)
str(UKcough)
summary(UKcough)

#### Merging CDCserology.csv (microneut) with DonorRecipientSubjectIDlist.csv (VolunteerList) ####
#microneut1 <- microneut %>%
  #filter(VisitType == "Q baseline" & Seroconversion == 1) %>%
  #rename(SeroconversionBeforeQ = Seroconversion) #%>%
  #select(-X, -Randomization, -DrawDate, -VisitType, -QuarantineNumber, -Microneutralization.Titer.to.A.Wisconsin.67.2005)
VolunteerList_MN <- VolunteerList %>%
  left_join(microneut, by = c("SubjectID" = "SubjectID", "Randomization" = "Randomization")) %>%
  #left_join(microneut1, by = c("SubjectID" = "SubjectID")) #%>%
  #filter(VisitType == "F/up") %>%
  select(-X.x, -X.y, -QuarantineNumber.y) %>%
  rename(QuarantineNumber = QuarantineNumber.x)

# The commented out lines above can be used if you want to generate a single row per subjectID with the serology seroconverted beforeQ and during Q results included
# Decided to leave this because of the raw data containing draw dates and MN titer measurements
# If we didn't care about the draw dates, we might present these findings in the wide format with one row per subjectID, but instead we will present the full data here and long format appears better suited for this df

#### Merging VolunteerListMN (generated in previous steps) with GLASGOWserology.csv (HAI data) ####

VolList_MN_HAI <- VolunteerList_MN %>%
  left_join(HAI, by = c("SubjectID" = "SubjectID")) %>%
  select(-QuarantineNumber.y, -X, -EMIT) %>%
  rename(QuarantineNumber = QuarantineNumber.x, RandomizationDIRCR = Randomization.x, RandomizationDR = Randomization.y)

#### Merging VolList_MN_HAI (generated in previous steps) with QuarantineDemographics.csv (Demog data) ####

VolList_MN_HAI_Demog <- VolList_MN_HAI %>%
  left_join(Demog, by = c("SubjectID" = "SubjectID")) %>%
  select(-Randomization, -X)

#### Merging VolList_MN_HAI_Demog (generated in previous steps) with PhysicalExam.csv (PhysicalExam data) + CDCpcr.csv (pcr data) ####
# First merge the pcr data with the PhysicalExam data 
# Can do this easily becasue we have StudyDay data for both of these dfs

# First prepare the data to merge VisitName (StudyDay) variable
PhysicalExam <- PhysicalExam %>%
  mutate(VisitName = str_remove_all(VisitName, "Quarantine Day "), 
         VisitName = str_remove_all(VisitName, " for Donors Only"), 
         VisitName = str_remove_all(VisitName, " for Recipients Only"), 
         VisitName = str_remove_all(VisitName, "Quarantine Discharge Day "),
         VisitName = strtoi(VisitName, base=0))

PhysExam_PCR <- pcr %>%
  full_join(PhysicalExam, by = c("SubjectID" = "SubjectID", "StudyDay" = "VisitName")) %>%
  select(-X.x, -X.y, -Randomization.x, -Randomization.y)

# Now merge this PhysExam_PCR df with the VolList_MN_HAI_Demog df

VolList_MN_HAI_Demog_PhysExam_PCR <- VolList_MN_HAI_Demog %>%
  full_join(PhysExam_PCR, by = c("SubjectID" = "SubjectID"))

#### Merge VolList_MN_HAI_Demog_PhysExam_PCR df (generated in previous steps) with QuarantineCalendar.csv (Calendar data) ####
## The purpose of this merge is to get add dates to the merged file So that it can be merged with the Symptom data (which only has dates, and a funky study day classification)

# First get QuarantineCalendar.csv (Calendar data) ready for this merge

Calendar1 <- Calendar %>%
  select(Qdays, Q1dates) %>%
  mutate(Qdays = recode(Qdays, `D-2` = "D-3", `D-1a` = "D-2", `D-1b` = "D-1")) %>%
  filter(!is.na(Q1dates))
  
Calendar2 <- Calendar %>%
  select(Qdays, Q2dates) %>%
  mutate(Qdays = recode(Qdays, `D-2` = "D-3", `D-1a` = "D-2", D0 = "D-1", D0repeat = "D0")) %>%
  filter(!is.na(Q2dates)) %>%
  arrange(Q2dates)

Calendar3 <- Calendar %>%
  select(Qdays, Q3dates) %>%
  mutate(Qdays = recode(Qdays, `D-1a` = "D-1")) %>%
  filter(!is.na(Q3dates))

Calendar <- Calendar1 %>%
  full_join(Calendar2, by = c("Qdays")) %>%
  full_join(Calendar3, by = c("Qdays"))

CalendarLong <- Calendar %>%
  gather(key = QuarantineNumber, value = Date, Q1dates, Q2dates, Q3dates, na.rm = TRUE, convert = TRUE) %>%
  arrange(Date) %>%
  mutate(QuarantineNumber = str_remove_all(QuarantineNumber, "dates"),
         QuarantineNumber = str_remove(QuarantineNumber, "Q")) %>%
  mutate(QuarantineNumber = as.numeric(QuarantineNumber))

CalendarLong <- CalendarLong %>%
  mutate(Qdays = str_remove(Qdays, "D")) %>%
  mutate(Qdays = as.numeric(Qdays))

# Need to add on Subject IDs for each of the StudyDays
CalendarLong <- VolunteerList %>%
  left_join(CalendarLong, by = c("QuarantineNumber" = "QuarantineNumber")) %>%
  select(SubjectID, QuarantineNumber, Qdays, Date)

#Now can merge the Calendar (date and study day) df with the growing big df
VolList_MN_HAI_Demog_PhysExam_PCR_Date <- VolList_MN_HAI_Demog_PhysExam_PCR %>%
  full_join(CalendarLong, by = c("SubjectID" = "SubjectID", "StudyDay" = "Qdays", "QuarantineNumber" = "QuarantineNumber"))

#### Merging VolList_MN_HAI_Demog_PhysExam_PCR_Date df with Symptoms df ####

# First get Symptom date variable from factor to character
Symptoms <- Symptoms %>%
  mutate(SDC_date = as.character(SDC_date))

VolList_MN_HAI_Demog_PhysExam_PCR_Date_Symptoms <- VolList_MN_HAI_Demog_PhysExam_PCR_Date %>%
  full_join(Symptoms, by = c("SubjectID" = "SubjectID", "Date" = "SDC_date")) %>%
  select(-X, -Randomization) %>%
  rename(Randomization_DorIRorCR = RandomizationDIRCR,
         Randomization_DorR = RandomizationDR,
         Microneut_DrawDate = DrawDate,
         Microneut_VisitType = VisitType,
         Microneut_Seroconvert = Seroconversion, 
         HAI_dayminus2 = HAI.2, 
         HAI_day28 = HAI28,
         HAI_dayminus2_recodeNDA = HAI.2recodeNDA,
         HAI_day28_recodeNDA = HAI28recodeNDA,
         HAI_dayminus2_recodeNDA_x4 = HAI.2recodeNDA.4, 
         HAI_Seroconversion = SeroconversionHAI,
         Sx_Date = Date)

#### Need to merge in the pcr copy number data that exists for the G-II fine and coarse samples ####
VolList_MN_HAI_Demog_PhysExam_PCR_Date_Symptoms_G2copies <- VolList_MN_HAI_Demog_PhysExam_PCR_Date_Symptoms %>%
  full_join(UKg2_all, by = c("SubjectID" = "subject.id", "StudyDay" = "sample.day")) %>%
  select(-X)

# Previously, the UKg2mean df was merged into the larger data frame here, however the UKg2mean data takes the average of these replicates and uses 2000*1/sqrt(2) for the missing values, which were then averaged so that there is a single, averaged, RNA copy number value for each of the samples (subject.id-sample.day instance). This is fine when we are interested in computing the GM and GSD for the positive samples (there aren't too many: n<30). However, when we are interested in using models of aerosol shedding using tobit regression or other methods that account for observations below the limit of detection, we want the raw pcr data on each replicate. That is, we don't want to impute the values that we have for samples where one replicate was negative and another was positive. Rather, we want to have the full replicate data as it is. 
# Thus, we add the UKg2_all df, which has this raw, replicate data (however when I say raw, perhaps that is the wrong term, because these aren't really raw data for RNA copy number, but rather they have already been multiplied by a dilution factor and by the qPCR calibration factor). The qPCR calibration factor is included in this dataset and the final.copies.replicate variable is simply the product of the calibration factor times the copy.num.replicate variable. The copy.number.replicate variable looks like it is the variable that comes directly from the output of the MxPro pcr assay analysis program (converting ct to virus particle -- the assay was standardized on an EM quantified PR/8 Influenza A virus stock). However, I have not seen code from Jing that was used to process the raw MxPro pcr exported data and manipulate it to get out the copy.number.replicate variable. In looking at the data, though, it is clear that to get from copies.in to copy.num, a factor of 2,000 was used. This factor of 2,000 is likely intended by Jing to be a combination of the dilution factor for aerosols in this setup, which was 50, multiplied by the RNA copies to virus particle ratio, which was 250 for this setup. Jing was using 80 for this RNA copies to virus particle ratio. 50*80 = 2,000 and this is what was used. 
# Jake has since edited the EMIT_Quarantin_Main_Cleaning.R script to update the RNA copies to virus particle ratio from 80 to 250 (given data from Michael Grantham, the raw data is in the directories for this project in box.com). Thus, instead of using the factor of 2,000 to get from copies.in to copy.num, we will use 12,500. 
# However, new evidence shows that the RNA copy to virus particle ratios of 80 and 411 are the finalized, revised numbers and these should be used in the analysis. Thus, we have changed the RNA copies to virus particle ratios back to 80 and 411 for flu A and B, respectively. 

#### Merge in the cough count data that was collected during the half-hour G-II sample collection instances ####
VolList_MN_HAI_Demog_PhysExam_PCR_Date_Symptoms_G2copies_cough <- VolList_MN_HAI_Demog_PhysExam_PCR_Date_Symptoms_G2copies %>%
  full_join(UKcough, by = c("SubjectID" = "subject_id", "Sx_Date" = "visit_date")) %>%
  select(-X)

#### Sort so the data is in order of SubjectID and then StudyDay ####
VolList_MN_HAI_Demog_PhysExam_PCR_Date_Symptoms_G2copies_cough <- VolList_MN_HAI_Demog_PhysExam_PCR_Date_Symptoms_G2copies_cough %>%
  arrange(SubjectID, StudyDay)

#### Write the merged authoritative dataframe #####

write.csv(VolList_MN_HAI_Demog_PhysExam_PCR_Date_Symptoms_G2copies_cough, "Curated Data/Analytical Datasets/QuarantineMergedData.csv")

# To facilitate this analysis for the EMIT Natural Versus Artificial manuscript analyses, we will write out another copy of this dataset to the appropriate directory
# write.csv(VolList_MN_HAI_Demog_PhysExam_PCR_Date_Symptoms_G2copies_cough, "/Users/jbueno/Box Sync/EMIT/EMIT_Data_Analysis_Jake/Natural_vs_Artificial_Infection/Analytical Datasets/EMIT_Quarantine_data_full.csv")

#### Comment about this merging process ####
# This merge works, however it may appear that there is replication of certain variable values
# Perhaps there are other ways to make the data appear more elegantly merged
# But the good thing about this merged df is that it has all the data and works for analysis