3.P639-1_Data_Cleaning.qmd

---
title: "P639-1 : Cheetah Samples | Data Cleaning Nasal Lavage"
output: 
  html_document:
    toc: yes
    toc_float:
      collapsed: no
  pdf_document:
    toc: yes
editor_options:
  chunk_output_type: console
header-includes:
  \usepackage{float}
---

```{=html}
<style type="text/css">
body{ /* Normal  */
      font-size: 14px;
  }
h1 { /* Header 1 */
  font-size: 28px;
}
h2 { /* Header 2 */
    font-size: 24px;
}
h3 { /* Header 3 */
  font-size: 20px;
}
h4 { /* Header 4 */
  font-size: 16px;
}
</style>
```

```{r}
#| include = FALSE
knitr::opts_chunk$set(fig.height=15, fig.width = 15, echo=FALSE,message=FALSE, warning=FALSE)
```

# Project Summary

How does the interferon response change post baseline and baseline post rhinovirus infection, and in particular the kinetics, how does it relate to viral clearance, patient history.

The primary questions are: How do type 1 and 2 interferon responses on the level of the interferon relate back to host characteristics.

```{r}
#| label = "setup",
#| echo = FALSE,
#| message = FALSE,
#| warning = FALSE,
#| cache = FALSE
library(knitr)
library(tidyverse)

library(ggthemes)
#library(edgeR)
library(RColorBrewer)
library(kableExtra)
library(viridis)
library(ggbeeswarm)
library(gplots)
library(upstartr)
library(janitor)
library(cowplot)
library(edgeR)
library(ggpubr)

library(ComplexHeatmap)
library(limma)
library(data.table)
library(multcomp)

#library(bRi)
library(miscHelpers)
library(RNAseQC)
library(countSubsetNorm)
library(apird)

library(ggrepel)
library(ggbeeswarm)
library(randomcoloR)
library(dendextend)
library(corrplot)
library(patchwork)
#library(ggsankey)
# Change package priority to dplyr, not a great way to do it but works for(if it aint broke dont fix it)
select <- dplyr::select
rename <- dplyr::rename
filter <- dplyr::filter
distinct <- dplyr::distinct
mutate <- dplyr::mutate


opts_chunk$set(
  fig.width=6, fig.height=4.25, cache = TRUE,
  echo=FALSE, warning=FALSE, message=FALSE)
options(stringsAsFactors = FALSE)

theme_set(
  theme_bw(20) +
    theme(panel.grid.major = element_blank(),
          panel.grid.minor = element_blank(),
          panel.border = element_rect(colour="black", fill = NA, size = 1),
          axis.text = element_text(colour="black"),
          axis.ticks = element_line(colour="black"),
          axis.text.x = element_text(angle=0)))
update_geom_defaults("point", list(shape=16))
grDevices::pdf.options(useDingbats = FALSE)
options(ggrepel.max.overlaps = Inf)
#Create 'not in' operator
`%notin%` <- Negate(`%in%`)
# test_plot <- ggplot2::ggplot(mapping = ggplot2::aes(x = 1:10, y = 1:10)) + ggplot2::geom_point()
# ggplot(mapping = aes(x = 1:10, y = 1:10)) + geom_point()

col.vec <- c(
  "none" = "#969696",
  "percent aligned" = "#CC79A7",
  "CV coverage, percent aligned" = "steelblue3",
  
  "total seqs, percent aligned" = "#feb24c",
  
  "CV coverage, total seqs, percent aligned" = "#fc4e2a",
  "total seqs" = "darkolivegreen",
  "CV coverage, total seqs" = "aquamarine"
  )
 rasterResolutionDpi <- 300
 
 
 clean_column_names <- function(col_names) {
  col_names %>%
    str_replace_all("\\.\\.\\.+", "_") %>%  # Replace multiple dots with a single underscore
    str_replace_all("\\.", "_") %>%         # Replace remaining dots with underscores
    str_replace_all("_$", "") %>%           # Remove trailing underscores
    str_replace_all("__+", "_") %>%         # Remove double underscores
    str_replace_all("\\s+", "_") %>%        # Replace spaces with underscores
    str_replace_all("[^A-Za-z0-9_]", "")    # Remove any remaining special characters
}
```

```{r}
#| label = "setFilenamesAndGlobalVariables"


projectNumber_lavage <- "P639-1"
dirBoxBase <- 
  file.path("~", "Library", "CloudStorage", "Box-Box")
dirRoot <-
  file.path(
    dirBoxBase, "Altman_Lab", # Box version
    # "~", "Documents", "Projects_local", # local version
    projectNumber_lavage)
dirPlots_lavage <- file.path("Figs",projectNumber_lavage)
dirresults_lavage <- file.path("Results",projectNumber_lavage)
dirDataclean <- file.path("Data_clean")
dirDataraw <- file.path("Data_raw")


```

```{r}
#| label = "set_directory",
#| cache = FALSE
if (!dir.exists(dirPlots_lavage)) dir.create(dirPlots_lavage)
if (!dir.exists(dirresults_lavage)) dir.create(dirresults_lavage)
#if (!dir.exists(dirDataclean)) dir.create(dirDataclean)
#if (!dir.exists(dirDataraw)) dir.create(dirDataraw)
```

# 1. Load in data

```{r}
#setwd("~/Library/CloudStorage/Box-Box/BRI_Projects/P506-1/")

design <- read_csv("Data_raw/P639-1_design_raw.csv") %>% arrange(libid) %>% 
  mutate(visit = ifelse(is.na(visit),"Unknown",visit)) %>% 
  mutate(site = ifelse(is.na(site),"Unknown",site))


counts_raw <- read_csv("Data_raw/P639-1_counts_raw.csv") %>% relocate(c(geneName,design$libid))


```

```{r}
design %>% 
  count(visit,site) %>% arrange(site) %>% 

kable(align="l",
      caption="Total libraries") %>% 
  kable_styling(bootstrap_options = "striped", full_width = FALSE)


design %>% 
  count(visit) %>% 

kable(align="l",
      caption="Total libraries") %>% 
  kable_styling(bootstrap_options = "striped", full_width = FALSE)
```

## 2.1 Load key.

## 2.2 Filter gene key to protein coding (pc) genes that occur in the count data set and have valid hgnc symbols.

```{r}
#| label = "key",
#| message = FALSE

if(file.exists(file.path(dirDataraw,"key_pc.csv"))){
  key.pc <- read_csv("Data_raw/key_pc.csv")
} else {
  library(biomaRt)
  ensembl <- useEnsembl(biomart = "genes", dataset = "hsapiens_gene_ensembl")
  full_key <- biomaRt::getBM(attributes = c("ensembl_gene_id","hgnc_symbol",	
                                "transcript_biotype"), mart = ensembl) %>% 
              #rename variable to match count data
              dplyr::rename(geneName = ensembl_gene_id,gene_biotype = transcript_biotype) 
  
  key.pc <- full_key %>% 
    # Keep only valid hgnc symbols
    filter(!is.na(hgnc_symbol)) %>% 
    # Keep protein coding genes only
    filter(gene_biotype == "protein_coding") %>% 
    # Remove duplicate entries
    distinct(geneName, .keep_all=TRUE) %>% 
    # Keep only genes found in dataset
    filter(geneName %in% counts_raw$geneName) %>% 
    arrange(geneName)
  
  write_csv(key.pc,"Data_raw/key_pc.csv")
}
```

## 2.3 Filter counts to protien coding genes

```{r}
counts.pc <- counts_raw %>% 
  filter(geneName %in% key.pc$geneName) %>% 
  arrange(geneName)
```

# 3. RNAseq Quality Control

## 3.1 Quality control by sequencing and alignment metrics

```{r}
#| label = "setQcCuts"
# set QC thresholds
# Usually start with lenient cutoffs and then increase or decrese as needed 
pct_aligned_cutoff = 75
fastq_total_reads_cutoff = 1E6
median_cv_coverage_cutoff = 1
pct_aligned_cutoff_2 = 50
```

```{r}
pal.subjectId <-
  big_colorblind_pal(
    n_distinct(design$subjectId),
    shuffle_colors = TRUE,
    drop_yellow = TRUE, drop_black = TRUE) %>%
  setNames(
    design %>%
      dplyr::arrange(subjectId) %>%
      dplyr::pull(subjectId) %>%
      unique())

pal.visit <-
  big_colorblind_pal(
    n_distinct(design$visit),
    shuffle_colors = TRUE,
    drop_yellow = TRUE, drop_black = TRUE) %>%
  setNames(
    design %>%
      dplyr::arrange(visit) %>%
      dplyr::pull(visit) %>%
      unique())

pal.site <-
  big_colorblind_pal(
    n_distinct(design$site),
    shuffle_colors = TRUE,
    drop_yellow = TRUE, drop_black = TRUE) %>%
  setNames(
    design %>%
      dplyr::arrange(site) %>%
      dplyr::pull(site) %>%
      unique())

pal.site <-
  paletteer::paletteer_d("MoMAColors::Rattner") %>% as.character() %>%
  setNames(
    design %>%
      dplyr::arrange(site) %>%
      dplyr::pull(site) %>%
      unique())

```

# 3.2 Look at Median CV coverage

```{r}
#| label = "cv.dups",
#| echo = FALSE
if (!dir.exists(file.path(dirPlots_lavage,"cleaning"))) dir.create(file.path(dirPlots_lavage,"cleaning"),recursive = TRUE)
median_cv_plot <- design %>% 
ggplot(aes(median_cv_coverage, pct_aligned)) +
           geom_point(aes(color=site,shape = visit),
                      size=3) +
  #Label points outside cutoffs
  geom_text_repel(data=filter(design,
                            median_cv_coverage > median_cv_coverage_cutoff | 
                            pct_aligned < pct_aligned_cutoff),
            aes(label=libid), show.legend = FALSE) +
  # Add cutoff lines
  geom_vline(xintercept=median_cv_coverage_cutoff, linetype=2) +
  geom_hline(yintercept=pct_aligned_cutoff, linetype=2) +
  #Beautify
  theme_classic() +
  #theme(legend.position = "none") +
  labs(x="Median CV coverage", y="Percent aligned") + scale_color_manual(values = pal.site)
  # lims(x=c(0,1.2), y=c(0,1))

median_cv_plot
#ggsave("Figs/cleaning/median_cv_coverage.pdf",plot = median_cv_plot,width = 15)
ggsave(file.path(dirPlots_lavage,"cleaning/median_cv_coverage.png"),plot = median_cv_plot,height = 10)



median_cv_plot_site_only <- design %>% 
ggplot(aes(median_cv_coverage, pct_aligned)) +
           geom_point(aes(color=site),
                      size=3) +
  #Label points outside cutoffs
  geom_text_repel(data=filter(design,
                            median_cv_coverage > median_cv_coverage_cutoff | 
                            pct_aligned < pct_aligned_cutoff),
            aes(label=libid), show.legend = FALSE) +
  # Add cutoff lines
  geom_vline(xintercept=median_cv_coverage_cutoff, linetype=2) +
  geom_hline(yintercept=pct_aligned_cutoff, linetype=2) +
  #Beautify
  theme_classic_custom +
  #theme(legend.position = "none") +
  labs(x="Median CV coverage", y="Percent aligned") + scale_color_manual(values = pal.site) + ggtitle("Median CV by Percent Aligned (Site)")
  # lims(x=c(0,1.2), y=c(0,1))

median_cv_plot_site_only
#ggsave("Figs/cleaning/median_cv_coverage.pdf",plot = median_cv_plot,width = 15)
ggsave(file.path(dirPlots_lavage,"cleaning/median_cv_coverage_site_only.png"),plot = median_cv_plot_site_only,height = 10)



```

### 3.3 Plot read counts

```{r}
library(scales)
P639_design_w_QC <- design %>%  mutate(col.group = ifelse(fastq_total_reads >= 1E6, "1,000,000+",
                            ifelse(fastq_total_reads >= 5E5, "500,000+",
                                   ifelse(fastq_total_reads >= 2.5E5, "250,000+",
                                          ifelse(fastq_total_reads >= 1E5, "100,000+",
                                                 "Less than 100,000")))),
          col.group = factor(col.group, levels=c("1,000,000+","500,000+","250,000+",
                                                 "100,000+","Less than 100,000")))

total_reads_QC_plot <- P639_design_w_QC %>% ggplot(aes(x=reorder_within(libid, by=fastq_total_reads,within = col.group), 
           y=fastq_total_reads, fill=site,shape=visit))  +
    geom_col() +

    # Add cutoff line
    geom_hline(yintercept = fastq_total_reads_cutoff) +

    # Beautify
    theme_classic()  +
    theme(axis.text.x = element_text(size=rel(0.75),
                                     angle = 90, hjust = 1)) +
    labs(x="libID", y="Total reads", fill="") +
    scale_x_reordered() +
    scale_y_continuous(trans = 'log10',
    breaks = trans_breaks("log10", function(x) 10^x),
    labels = trans_format("log10", math_format(10^.x))) + scale_fill_manual(values = pal.site)

QC_combined_plot <- median_cv_plot+total_reads_QC_plot+plot_layout(nrow=1,guides = "collect")

ggsave(file.path(dirPlots_lavage,paste("cleaning/total_reading_medianCV_comb.pdf",sep = "")),plot = QC_combined_plot,height = 10,width = 18)

total_reads_QC_plot_failed <- P639_design_w_QC %>% subset(col.group!="1,000,000+") %>% ggplot(aes(x=reorder_within(libid, by=fastq_total_reads,within = col.group), 
           y=fastq_total_reads, fill=site,shape=visit))  +
    geom_col() +

    # Add cutoff line
    geom_hline(yintercept = fastq_total_reads_cutoff) +

    # Beautify
    theme_classic()  +
    theme(axis.text.x = element_text(size=rel(0.75),
                                     angle = 90, hjust = 1)) +
    labs(x="libID", y="Total reads", fill="") +
    scale_x_reordered() +
    scale_y_continuous(trans = 'log10',
    breaks = trans_breaks("log10", function(x) 10^x),
    labels = trans_format("log10", math_format(10^.x))) + scale_fill_manual(values = pal.site)
QC_combined_plot_site_only <- median_cv_plot_site_only+total_reads_QC_plot_failed+plot_layout(nrow=1,guides = "collect")
ggsave(file.path(dirPlots_lavage,paste("cleaning/total_reading_medianCV_comb_failed_by_site.pdf",sep = "")),plot = QC_combined_plot_site_only,height = 7,width = 15)
QC_combined_plot_site_only_pct_50 <- median_cv_plot_site_only_pct_alig_50+total_reads_QC_plot_failed+plot_layout(nrow=1,guides = "collect")
ggsave(file.path(dirPlots_lavage,paste("cleaning/total_reading_medianCV_comb_failed_by_site_pct_aligned_50.pdf",sep = "")),plot = QC_combined_plot_site_only_pct_50,height = 10,width = 10)
```



# 3.3 Plot total counts per library. Libraries outside cutoffs are labeled as "questionable". The minimum total sequences cutoff set above is indicated by a horizontal line.


# PCA plots for QC

```{r}
#| echo = FALSE,
#| warning = FALSE
#Calculate PCA for all data.
PCA.all <- counts.pc %>% 
  column_to_rownames("geneName") %>% 
  #Convert to log counts per million
  cpm(., log=TRUE) %>% 
  t() %>% 
  #Calc PCA
  prcomp()

PC1.label <- paste("PC1 (", summary(PCA.all)$importance[2,1]*100, "%)", sep="")
PC2.label <-paste("PC2 (", summary(PCA.all)$importance[2,2]*100, "%)", sep="")

# Extract PC values
PCA.all.dat <- as.data.frame(PCA.all$x) %>% 
  rownames_to_column("libid") %>% 
  # Select PCs for plotting
  dplyr::select(libid, PC1:PC3) %>% 
  # Merge with metadata
  left_join(P639_design_w_QC, by="libid")

plot_pca_unnorm_by_qc <- ggplot(PCA.all.dat, aes(PC1, PC2)) +
           geom_point(aes(color=site,shape=visit),alpha=0.8,
                      size=3) + 
  #Beautify
  theme_classic() +
  theme(legend.title = element_blank(),
        legend.position = "bottom") +
  labs(x=PC1.label, y=PC2.label, title="Un-normalized logCPM") + scale_color_manual(values = pal.site) +
  coord_fixed(ratio=1) 
plot_pca_unnorm_by_qc
ggsave(file.path(dirPlots_lavage,"cleaning/pca_unnorm_by_qc.pdf"),plot = plot_pca_unnorm_by_qc)

plot_pca_unnorm_by_qc_by_site <- ggplot(PCA.all.dat, aes(PC1, PC2)) +
           geom_point(aes(color=site),alpha=0.8,
                      size=3) + 
  #Beautify
  theme_classic() +
  theme(legend.title = element_blank(),
        legend.position = "bottom") +
  geom_text_repel(data=filter(PCA.all.dat,
                            median_cv_coverage > median_cv_coverage_cutoff | 
                            pct_aligned < pct_aligned_cutoff |
                             fastq_total_reads < fastq_total_reads_cutoff),
            aes(label=libid), show.legend = FALSE) + 
  labs(x=PC1.label, y=PC2.label, title="Un-normalized logCPM (Site)") + scale_color_manual(values = pal.site) +
  coord_fixed(ratio=1) 
plot_pca_unnorm_by_qc_by_site

ggsave(file.path(dirPlots_lavage,"cleaning/pca_unnorm_by_qc_by_site.pdf"),plot = plot_pca_unnorm_by_qc_by_site,height = 10,width = 10)
```

## Check sample validity by sex chromosome gene expression, kinship comparison, and stimulant PCA

```{r}
#| label = "inferSexByRna",
#| dependson = "applyQcCuts"
## infer sex based on RNA-seq reads
P639_design_w_QC$logXyRatio <-
  logXYratio(counts.pc %>% relocate(c(geneName,P639_design_w_QC$libid)) %>% column_to_rownames("geneName"), gene_ID = "ensembl",use_annotables = FALSE)
# plot histogram, and use it to determine break point

logXyRatioThreshold <- 6.0
pal.sex <-
  c("F"="orange", "M"="blue")
P639_design_w_QC$sexByRna <-
  ifelse(P639_design_w_QC$logXyRatio >= logXyRatioThreshold, "F", "M") %>%
  factor(levels=c("F", "M"))
```

```{r}
#| label = "plotLogXyRatio",
#| dependson = "inferSexByRna",
#| fig.width = 7.25,
#| fig.height = 4.25
plot.tmp <-
  ggplot(
    P639_design_w_QC,
    mapping = aes(x = logXyRatio, fill = sexByRna)) +
  geom_histogram(color = "black", position = "dodge") +
  geom_vline(xintercept = logXyRatioThreshold, linetype = "dashed") +
  labs(x = "log ratio of X to Y chromosome reads") +
  scale_fill_manual("sex by RNA", values = pal.sex, na.value="gray60")

print(plot.tmp)

```

Remove samples that fail QC

```{r}
final_summary <- P639_design_w_QC %>%
  count(visit, site) %>%
  pivot_wider(
    names_from = visit,
    values_from = n
  ) %>%
  select(site, `Screening/Enrollment`, `Visit 1`,
         `Cold A Clinic`, Unknown) %>%
  # Add row totals
  rowwise() %>%
  mutate(Total = sum(c_across(`Screening/Enrollment`:Unknown), na.rm = TRUE)) %>%
  ungroup() %>%
  # Add column totals as a new row
  bind_rows(
    summarise(.,
              across(`Screening/Enrollment`:Total, ~ sum(.x, na.rm = TRUE)) %>%
                mutate(site = "TOTAL") %>%
                select(site, everything()))
  )

final_summary

P639_design_w_QC_remove_failed <- P639_design_w_QC %>% subset(pct_aligned>=pct_aligned_cutoff&fastq_total_reads>=fastq_total_reads_cutoff&median_cv_coverage<=median_cv_coverage_cutoff)
#P639_design_w_QC_remove_failed_pct_aligned_50 <- P639_design_w_QC %>% subset(pct_aligned>=pct_aligned_cutoff_2&fastq_total_reads>=fastq_total_reads_cutoff&median_cv_coverage<=median_cv_coverage_cutoff)
failed_samples <- P639_design_w_QC %>% subset(libid%notin%P639_design_w_QC_remove_failed$libid)
#failed_samples_pct_aligned_50 <- P639_design_w_QC %>% subset(libid%notin%P639_design_w_QC_remove_failed_pct_aligned_50$libid)

counts.pc_w_remove_failed <- counts.pc[,c("geneName",P639_design_w_QC_remove_failed$libid)]

total_samples_by_site <- P639_design_w_QC %>% count(site) %>% rename("Total"="n")
total_samples_by_site <- P639_design_w_QC %>% count(site,siteId,siteNumber) %>% rename("Total"="n")
```

Write the design w QC for pass fail by site
```{r}
P639_design_w_QC %>% select(libid,subjectId,site,siteId,siteNumber,barcode,visit,dateCollected,pct_aligned,fastq_total_reads,median_cv_coverage) %>% mutate(site = ifelse(siteId == "Boston","Boston Childrens",ifelse(siteId =="?","TBD",site))) %>% 
  mutate(QC = ifelse(pct_aligned>=pct_aligned_cutoff&fastq_total_reads>=fastq_total_reads_cutoff&median_cv_coverage<=median_cv_coverage_cutoff,"pass","fail")) %>% write_csv("Data_clean/P639_Nasal_Lavage_Design_by_QC_Site.csv")
```

Of the total number of libraries `r length(P639_design_w_QC$libid)` `r length(failed_samples$libid)` Failed

```{r}
P639_design_w_QC_remove_failed %>% 
  count(visit,site) %>% arrange(site) %>% 

kable(align="l",
      caption="Total libraries") %>% 
  kable_styling(bootstrap_options = "striped", full_width = FALSE)


P639_design_w_QC_remove_failed %>% 
  count(site) %>% rename("Samples"="n") %>% left_join(total_samples_by_site ) %>% arrange(site) %>% 

kable(align="l",
      caption="Total libraries") %>% 
  kable_styling(bootstrap_options = "striped", full_width = FALSE)
```

The `r length(failed_samples$libid)` that failed QC

```{r}
failed_samples %>% 
  count(visit,site,siteNumber,siteId) %>% arrange(site) %>% 

kable(align="l",
      caption="Total libraries") %>% 
  kable_styling(bootstrap_options = "striped", full_width = FALSE)

failed_samples %>% 
  count(site) %>% rename("Samples"="n") %>% left_join(total_samples_by_site ) %>% arrange(site) %>%
  mutate(Percent_failed = round((Samples/Total)*100,digits = 2)) %>% 

kable(align="l",
      caption="Total libraries") %>% 
  kable_styling(bootstrap_options = "striped", full_width = FALSE)





failed_summary <- failed_samples %>%
  count(visit, site) %>%
  pivot_wider(
    names_from = visit,
    values_from = n
  ) %>%
  select(site, `Screening/Enrollment`,
         `Cold A Clinic`, Unknown) %>%
  # Add row totals
  rowwise() %>%
  mutate(Total = sum(c_across(`Screening/Enrollment`:Unknown), na.rm = TRUE)) %>%
  ungroup() %>%
  # Add column totals as a new row
  bind_rows(
    summarise(.,
              across(`Screening/Enrollment`:Total, ~ sum(.x, na.rm = TRUE)) %>%
                mutate(site = "TOTAL") %>%
                select(site, everything()))
  )

failed_summary
```

# 4. Create DGEList object

```{r}
#| label = "dgelist"
file.tmp <- file.path(dirDataclean, paste(projectNumber_lavage,"_dat_clean_unnormalized.RDS",sep=""))
if(!file.exists(file.tmp)){
dat.pc <- DGEList(
  #count table. move gene names to column names
  counts=as.matrix(column_to_rownames(counts.pc_w_remove_failed,
                                      "geneName")),
  #metadata
  samples=P639_design_w_QC_remove_failed %>% arrange(libid),
  #keep genes in count table
  genes=key.pc)

saveRDS(dat.pc, 
     file=file.tmp)
}else{
  dat.pc <- readRDS(file.tmp)
}
```

# 5. Rare gene filtering

The raw gene sets contain highly variable, low abundance/rare genes.

## Threshold for filtering :

Present in minimum percent of samples = 10 minimum CPM = 0.8

```{r}
#| eval = FALSE
dat.pc.abund<-RNAetc::filter_rare(dat = dat.pc, min.pct = 10,min.CPM = 0.8,gene.var = "geneName",plot = TRUE)
```

# 6. Normalize the data

## 6.1 Calculate factors to scale library sizes.

```{r}
#| label = "norm",
#| eval = FALSE
dat.pc.abund.norm <- calcNormFactors(dat.pc.abund)
```

## 6.2 Normalize with voom

```{r}
#| label = "voom3"
file.tmp <- file.path(dirDataclean, paste(projectNumber_lavage,"_dat_voom_normalized.RDS",sep=""))
if(!file.exists(file.tmp)){
dat.voom <- voomWithQualityWeights(
                          dat.pc.abund.norm,
                          design=model.matrix(~visit,
                                 data=dat.pc.abund.norm$samples),
                          plot=TRUE)

saveRDS(dat.voom,file=file.tmp)
}else{
  dat.voom <- readRDS(file=file.tmp)
}
```

# Add clinical data
```{r}
file.tmp <- file.path(dirDataclean, paste(projectNumber_lavage,"_dat_voom_normalized_w_clinical.RDS",sep=""))
if(!file.exists(file.tmp)){
dat.voom_w_clinical <- dat.voom
clinical_data <- readxl::read_xlsx(file.path(dirDataraw,"EOS_FENO_AST.xlsx")) %>% rename("subjectId"="ID")
dat.voom_w_clinical$targets <- dat.voom_w_clinical$targets %>% left_join(clinical_data)
saveRDS(dat.voom_w_clinical,file=file.tmp)
}else{
  dat.voom <- readRDS(file=file.tmp)
}
```

```{r}
pal.subjectId <-
  big_colorblind_pal(
    n_distinct(dat.voom$targets$subjectId),
    shuffle_colors = TRUE,
    drop_yellow = TRUE, drop_black = TRUE) %>%
  setNames(
    dat.voom$targets %>%
      dplyr::arrange(subjectId) %>%
      dplyr::pull(subjectId) %>%
      unique())

pal.visit <-
  big_colorblind_pal(
    n_distinct(dat.voom$targets$visit),
    shuffle_colors = TRUE,
    drop_yellow = TRUE, drop_black = TRUE) %>%
  setNames(
    dat.voom$targets %>%
      dplyr::arrange(visit) %>%
      dplyr::pull(visit) %>%
      unique())

pal.site <-
  big_colorblind_pal(
    n_distinct(dat.voom$targets$site),
    shuffle_colors = TRUE,
    drop_yellow = TRUE, drop_black = TRUE) %>%
  setNames(
    dat.voom$targets %>%
      dplyr::arrange(site) %>%
      dplyr::pull(site) %>%
      unique())

```

# 8. PCA final data set

## 8.1 PCA by Treatment

```{r}
#| echo = FALSE,
#| warning = FALSE

#Calculate PCA for voom norm data.
file.tmp <- file.path(dirDataclean, paste(projectNumber_lavage,"_PCA.RDS",sep=""))
if(!file.exists(file.tmp)){
PCA.voom <- as.data.frame(dat.voom$E) %>% 
  t() %>% 
  #Calc PCA
  prcomp(scale. = TRUE)
saveRDS(PCA.voom,file=file.tmp)
}else{
  PCA.voom <- readRDS(file.tmp)
}


PC1.label <- paste("PC1 (", summary(PCA.voom)$importance[2,1]*100, "%)", sep="")
PC2.label <-paste("PC2 (", summary(PCA.voom)$importance[2,2]*100, "%)", sep="")
PC3.label <-paste("PC3 (", summary(PCA.voom)$importance[2,3]*100, "%)", sep="")

# Extract PC values
PCA.voom.dat<- as.data.frame(PCA.voom$x) %>% 
  rownames_to_column("libid") %>%
  # Select PCs for plotting
  dplyr::select(libid, PC1:PC3) %>% 
  # Merge with metadata
  left_join(as.data.frame(dat.voom$targets),
            by="libid")

plot_normalized_site<-ggplot(PCA.voom.dat ,
                  aes(PC1, PC2)) +
           geom_point(aes(color=site),
                      size=3) +
  #Beautify
  theme_classic() +
  theme(
        legend.position = "right",legend.text = element_text(size=12)) + 
  labs(x=PC1.label, y=PC2.label, title="Voom normalized log2 CPM") +scale_color_manual(values = pal.site) +
  coord_fixed(ratio=1)
plot_normalized_site
ggsave(file.path(dirPlots_lavage,"cleaning/PCA_normalized_site.pdf"),plot = plot_normalized_site)

plot_normalized_visit<-ggplot(PCA.voom.dat ,
                  aes(PC1, PC2)) +
           geom_point(aes(color=visit),
                      size=3) +
  #Beautify
  theme_classic() +
  theme(
        legend.position = "right",legend.text = element_text(size=12)) + 
  labs(x=PC1.label, y=PC2.label, title="Voom normalized log2 CPM") +scale_color_manual(values = pal.visit) +
  coord_fixed(ratio=1)
plot_normalized_visit
ggsave(file.path(dirPlots_lavage,"cleaning/PCA_normalized_visit.pdf"),plot = plot_normalized_visit)

# Clinical variables
plot_normalized_EOS<-ggplot(PCA.voom.dat ,
                  aes(PC1, PC2)) +
           geom_point(aes(color=`Eosinophil Count (V0)`),
                      size=3) +
  #Beautify
  theme_classic() +
  theme(
        legend.position = "right",legend.text = element_text(size=12)) + 
  labs(x=PC1.label, y=PC2.label, title="Voom normalized log2 CPM") +#scale_color_manual(values = pal.visit) +
  coord_fixed(ratio=1)
plot_normalized_EOS
ggsave(file.path(dirPlots_lavage,"cleaning/PCA_normalized_EOS.pdf"),plot = plot_normalized_EOS)


```

# BY QC metrics
```{r}
plot_normalized_fastq_reads<-ggplot(PCA.voom.dat ,
                  aes(PC1, PC2)) +
           geom_point(aes(color=fastq_total_reads),
                      size=3) +
  #Beautify
  theme_classic() +
  theme(
        legend.position = "right",legend.text = element_text(size=12)) + 
  labs(x=PC1.label, y=PC2.label, title="Voom normalized log2 CPM") +scale_color_viridis(direction = -1) +
  coord_fixed(ratio=1)
plot_normalized_fastq_reads
ggsave(file.path(dirPlots_lavage,"cleaning/PCA_normalized_total_reads.pdf"),plot = plot_normalized_fastq_reads)


plot_normalized_pct_aligned<-ggplot(PCA.voom.dat ,
                  aes(PC1, PC2)) +
           geom_point(aes(color=pct_aligned),
                      size=3) +
  #Beautify
  theme_classic() +
  theme(
        legend.position = "right",legend.text = element_text(size=12)) + 
  labs(x=PC1.label, y=PC2.label, title="Voom normalized log2 CPM") +scale_color_viridis(direction = -1) +
  coord_fixed(ratio=1)
plot_normalized_pct_aligned
ggsave(file.path(dirPlots_lavage,"cleaning/PCA_normalized_pct_aligned.pdf"),plot = plot_normalized_pct_aligned)


plot_normalized_median_cv<-ggplot(PCA.voom.dat ,
                  aes(PC1, PC2)) +
           geom_point(aes(color=median_cv_coverage),
                      size=3) +
  #Beautify
  theme_classic() +
  theme(
        legend.position = "right",legend.text = element_text(size=12)) + 
  labs(x=PC1.label, y=PC2.label, title="Voom normalized log2 CPM") +scale_color_viridis(direction = 1) +
  coord_fixed(ratio=1)
plot_normalized_median_cv
ggsave(file.path(dirPlots_lavage,"cleaning/PCA_normalized_median_cv.pdf"),plot = plot_normalized_median_cv)


plot_normalized_mapped_reads<-ggplot(PCA.voom.dat ,
                  aes(PC1, PC2)) +
           geom_point(aes(color=mapped_reads),
                      size=3) +
  #Beautify
  theme_classic() +
  theme(
        legend.position = "right",legend.text = element_text(size=12)) + 
  labs(x=PC1.label, y=PC2.label, title="Voom normalized log2 CPM") +scale_color_viridis(direction = -1) +
  coord_fixed(ratio=1)
plot_normalized_mapped_reads
ggsave(file.path(dirPlots_lavage,"cleaning/PCA_normalized_median_cv.pdf"),plot = plot_normalized_median_cv)


plot_normalized_well<-ggplot(PCA.voom.dat ,
                  aes(PC1, PC2)) +
           geom_point(aes(color=well),
                      size=3) +
  #Beautify
  theme_classic() +
  theme(
        legend.position = "right",legend.text = element_text(size=12)) + 
  labs(x=PC1.label, y=PC2.label, title="Voom normalized log2 CPM") +#scale_color_viridis(direction = 1) +
  coord_fixed(ratio=1)
plot_normalized_well
ggsave(file.path(dirPlots_lavage,"cleaning/PCA_normalized_well.pdf"),plot = plot_normalized_well)


plot_normalized_collectionTime<-ggplot(PCA.voom.dat %>% mutate(collectionTime = as.POSIXct(collectionTime)),
                  aes(PC1, PC2)) +
           geom_point(aes(color=collectionTime),
                      size=3) +
  #Beautify
  theme_classic() +
  theme(
        legend.position = "right",legend.text = element_text(size=12)) + 
  labs(x=PC1.label, y=PC2.label, title="Voom normalized log2 CPM") +
  coord_fixed(ratio=1)
plot_normalized_collectionTime
ggsave(file.path(dirPlots_lavage,"cleaning/PCA_normalized_collectionTime.pdf"),plot = plot_normalized_collectionTime)


plot_normalized_subjectId<-ggplot(PCA.voom.dat,
                  aes(PC1, PC2)) +
           geom_point(aes(color=subjectId),
                      size=3) +
  #Beautify
  theme_classic() +
  theme(
        legend.position = "right",legend.text = element_text(size=12)) + 
  labs(x=PC1.label, y=PC2.label, title="Voom normalized log2 CPM") +
  coord_fixed(ratio=1)
plot_normalized_collectionTime
ggsave(file.path(dirPlots_lavage,"cleaning/PCA_normalized_collectionTime.pdf"),plot = plot_normalized_collectionTime)
```

# 13. Co-variate correlation with PC



```{r}
#| fig.width = 8.5,
#| fig.height = 8.5,
#| message = FALSE
#PCA percents
pca <- stats::prcomp(t(dat.voom$E), scale.=TRUE)
pct <- summary(pca)$importance %>% t() %>% 
  as.data.frame() %>% 
  rownames_to_column("name") %>% 
  mutate(newname = paste(name, " (", signif(`Proportion of Variance`*100, 2), 
                      "%)", sep="")) %>% 
  select(name, newname)
#save(pct, file = "data_clean/PC_labs_voom_normalized.RData")
#PC values
pca.dat <- as.data.frame(pca$x) %>% 
  rownames_to_column() %>% 
  pivot_longer(-rowname) %>% 
  left_join(pct, by = "name") %>% 
  select(-name) %>% 
  pivot_wider(names_from = newname) %>% 
  arrange(rowname) %>% 
  column_to_rownames() %>%  select(1:13)

#Co-variates as numeric
meta.dat.cov <- dat.voom$targets %>% 
  arrange(libid) %>%
  select(subjectId, median_cv_coverage,fastq_total_reads,pct_aligned,visit,site,mapped_reads
         ) %>% 
  mutate_if(is.character, ~as.numeric(as.factor(.))) %>% 
  mutate_if(is.factor, ~as.numeric(.)) %>% 
  as.matrix()
rownames(meta.dat.cov) <- sort(dat.voom$targets$libid)

#Correlation
corr <- psych::corr.test(pca.dat, meta.dat.cov, method = "kendall", adjust = "none")

#heatmap
pdf(file.path(dirPlots_lavage,"cleaning/PCA_correlation.pdf"))
corrplot(t(corr$r), p.mat=t(corr$p), sig.level = c(0.01, 0.05), pch.cex = 0.9,
         insig = 'label_sig', method = "color", tl.col = 'black',
         col = rev(COL2('RdBu', 10)), cl.ratio = 0.2, 
         title = "Kendall correlation (r)\n **P < 0.01 *P < 0.05",
         mar=c(0,0,2,0))
invisible(dev.off())
```


```{r}
passed_summary <- dat.voom_lavage_wo_unknown$targets %>%
  count(visit, site) %>%
  pivot_wider(
    names_from = visit,
    values_from = n
  ) %>%
  select(site, Screening,
         `Cold A Clinic`) %>%
  # Add row totals
  rowwise() %>%
  mutate(Total = sum(c_across(Screening:`Cold A Clinic`), na.rm = TRUE)) %>%
  ungroup() %>%
  # Add column totals as a new row
  bind_rows(
    summarise(.,
              across(Screening:Total, ~ sum(.x, na.rm = TRUE)) %>%
                mutate(site = "TOTAL") %>%
                select(site, everything()))
  )

passed_summary
```



# R session

```{r}
sessionInfo()
```