---
title: "COVID19: Differential abundance based on symptom onset"
output: html_notebook
---

Performing differential abundance testing using a negative binomial GLM. Cell counts per donor sample are used as input, and the total number of cells 
captured (after QC) are used to normalize the model counts. This notebook is concerned with using the initial clustering results and their differential 
abundance according to symptom onset in hospitalised patients.


```{r, warning=FALSE, message=FALSE}
library(ggplot2)
library(ggsci)
library(cowplot)
library(RColorBrewer)
library(reshape2)
library(colorspace)
library(ggthemes)
library(edgeR)
library(scales)
library(ggrepel)
library(dplyr)
library(kableExtra)
```


```{r}
covid.meta <- read.csv("~/Dropbox/COVID19/Data/Metadata FINAL 10122020.csv",
                        header=TRUE, stringsAsFactors=FALSE)

covid.meta$Days_from_onset[covid.meta$Days_from_onset %in% c("Not_known", "Healthy")] <- NA
covid.meta$Days_from_onset <- as.numeric(covid.meta$Days_from_onset)
```

Plot days since symptom onset by disease severity and Site.

```{r, fig.height=2.95, fig.width=3.95,}
covid.meta$OrderedSeverity <- ordered(covid.meta$D0_status_summary,
                                      levels=c("Healthy", "Asymptomatic", "Mild", "Moderate", "Severe", "Critical"))

ggplot(covid.meta[covid.meta$Collection_Day %in% c("D0") &
                      !covid.meta$D0_status_summary %in% c("LPS", "Non_covid", "Asymptomatic", "Healthy") &
                      !is.na(covid.meta$Days_from_onset), ],
       aes(x=OrderedSeverity, y=Days_from_onset)) +
    geom_boxplot() +
    labs(x="Disease severity", y="Days from\nsymptom onset") +
    theme_cowplot()  +
    theme(aspect=1,
          axis.text.x=element_text(angle=90, vjust=0.5, hjust=1)) +
    ggsave("~/Dropbox/COVID19/plot.dir/SymptomOnset_severity-boxplot.pdf",
           height=2.95, width=3.95, useDingbats=FALSE) +
    NULL
```


```{r, fig.height=2.95, fig.width=3.95,}
covid.meta$OrderedSeverity <- ordered(covid.meta$D0_status_summary,
                                      levels=c("Healthy", "Asymptomatic", "Mild", "Moderate", "Severe", "Critical"))

ggplot(covid.meta[covid.meta$Collection_Day %in% c("D0") &
                      !covid.meta$D0_status_summary %in% c("LPS", "Non_covid", "Asymptomatic", "Healthy") &
                      !is.na(covid.meta$Days_from_onset), ],
       aes(x=Site, y=Days_from_onset)) +
    geom_boxplot() +
    labs(x="Disease severity", y="Days from\nsymptom onset") +
    theme_cowplot()  +
    theme(aspect=1,
          axis.text.x=element_text(angle=90, vjust=0.5, hjust=1)) +
    ggsave("~/Dropbox/COVID19/plot.dir/SymptomOnset_Site-boxplot.pdf",
           height=2.95, width=3.95, useDingbats=FALSE) +
    NULL
```



```{r, warning=FALSE, message=FALSE}
all.meta <- read.table("~/Dropbox/COVID19/Data/COVID19_scMeta-data.tsv",
                       sep="\t", header=TRUE, stringsAsFactors=FALSE)
rownames(all.meta) <- all.meta$CellID

doublet.remove <- read.table("~/Dropbox/COVID19/Data/Annotations/All_doublets.tsv",
                             sep="\t", header=FALSE, stringsAsFactors=FALSE)

# remove BGCV01_CV0209 and CV0198
all.meta <- all.meta[!all.meta$sample_id %in% c("BGCV01_CV0902"), ]
all.meta <- all.meta[!all.meta$patient_id %in% c("CV0198"), ]
all.meta$Days_from_onset[all.meta$Days_from_onset %in% c("Not_known", "Healthy")] <- NA
all.meta$Days_from_onset <- as.numeric(all.meta$Days_from_onset)

initial_clustering <- read.table("~/Dropbox/COVID19/Data/combined_dec_MetaData.txt",
                                 sep="\t", header=TRUE, stringsAsFactors=FALSE, row.names=1)
rownames(initial_clustering) <- gsub(rownames(initial_clustering), pattern="^b", replacement="")
initial_clustering$initial_clustering <- gsub(initial_clustering$initial_clustering, pattern="^b", replacement="")
initial_clustering$CellID <- rownames(initial_clustering)

all.meta <- merge(all.meta, initial_clustering[, c("CellID", "initial_clustering")], by='CellID')

n.cell.vecc <- table(all.meta$sample_id)

all.meta <- all.meta[!all.meta$CellID %in% doublet.remove$V1, ]
```


```{r}
cell.xtab <- as.data.frame(xtabs( ~ sample_id +  initial_clustering, data=all.meta))
cell.cast <- dcast(cell.xtab, sample_id ~ initial_clustering, value.var='Freq')
rownames(cell.cast) <- cell.cast$sample_id
cell.cast <- cell.cast[, -1]
cell.freq <- as.data.frame(t(sapply(rownames(cell.cast), FUN=function(X) as.numeric(cell.cast[X, ]/n.cell.vecc[X]),
                                    simplify=TRUE)))

rownames(cell.freq) <- rownames(cell.cast)
colnames(cell.freq) <- colnames(cell.cast)
cell.freq$sample_id <- rownames(cell.cast)
cell.melt <- melt(cell.freq, id.vars=c('sample_id'))
colnames(cell.melt) <- c("sample_id", "CellType", "Freq")
cell.melt$CellType <- as.character(cell.melt$CellType)

cell.freq.merge <- merge(cell.melt, covid.meta, by='sample_id')
```


## Cluster-based DA analysis: days since symptom onset

```{r}
# set up testing model
rownames(covid.meta) <- covid.meta$sample_id
init.meta <- covid.meta[!covid.meta$D0_status_summary %in% c("Non_covid", "LPS", "Healthy"), ]
init.meta$OrderedSeverity <- ordered(init.meta$D0_status_summary,
                                      levels=c("Asymptomatic", "Mild", "Moderate", "Severe", "Critical"))
init.meta <- init.meta[init.meta$D0_status_summary %in% c("Mild", "Moderate", "Severe", "Critical"), ]

init.meta$Days_from_onset[init.meta$Days_from_onset %in% c("Not_known")] <- NA
init.meta$Days_from_onset <- as.numeric(init.meta$Days_from_onset)
init.meta <- init.meta[!is.na(init.meta$Days_from_onset), ]

init.model <- model.matrix(~ Sex + Age + Site + Days_from_onset, 
                            data=init.meta[init.meta$Collection_Day %in% c("D0"), ])

# count cells
cell.freq.tab <- t(table(all.meta$sample_id[all.meta$Collection_Day %in% c("D0") &
                                                !is.na(all.meta$Days_from_onset) &
                                                !all.meta$D0_status_summary %in% c("LPS", "Non_covid", "Healthy", "Asymptomatic")],
                         all.meta$initial_clustering[all.meta$Collection_Day %in% c("D0") &
                                                         !is.na(all.meta$Days_from_onset) &
                                                         !all.meta$D0_status_summary %in% c("LPS", "Non_covid", "Healthy", "Asymptomatic")]))

cell.freq.tab <- cell.freq.tab[!rownames(cell.freq.tab) %in% c("Doublet", "Doublets", "Doublets:Bcell", "Doublets:Platelet"), ]
test.samps <- intersect(colnames(cell.freq.tab), names(n.cell.vecc))
cell.freq.tab <- cell.freq.tab[, test.samps]
init.model <- init.model[colnames(cell.freq.tab), ]

init.dge <- DGEList(cell.freq.tab, lib.size=log(n.cell.vecc[test.samps]))

#estimate dispersions
init.dge <- estimateDisp(init.dge, design=init.model)
init.linear.fit <- glmQLFit(init.dge, init.model, robust=TRUE)
init.res <- as.data.frame(topTags(glmQLFTest(init.linear.fit, coef=4), sort.by='none', n=Inf))
init.res$CellType <- rownames(init.res)

init.res$Sig <- as.numeric(init.res$FDR < 0.1)
init.res$Diff <- sign(init.res$logFC)
init.res$Diff[init.res$FDR > 0.1] <- 0
table(init.res$Diff)
```



```{r, warning=FALSE, message=FALSE, fig.height=2.95, fig.width=4.15}
mx.lfc <- max(abs(init.res$logFC))
eps <- mx.lfc * 0.05

init.resplot.labels <- init.res$CellType
init.resplot.labels[init.res$Sig == 0] <- ""

ggplot(init.res, aes(x=logFC, y=-log10(FDR), fill=as.character(Sig))) +
    geom_point(shape=21, size=4) +
    theme_cowplot() +
    scale_fill_manual(values=c("black", "red")) +
    scale_x_continuous(limits=c(-mx.lfc - eps, mx.lfc + eps), oob=squish) +
    guides(fill=guide_legend(title="FDR < 0.1")) +
    geom_text_repel(aes(label=init.resplot.labels),
                     fill='white') +
    labs(x="log fold change", y=expression(paste("-log"[10], " FDR"))) +
    ggsave("~/Dropbox/COVID19/plot.dir/All_edgeR_volcano_daysOnset-linear.pdf",
           height=2.95, width=4.15, useDingbats=FALSE) +
    NULL
```


Show a plot of just the DA clusters.

```{r, warning=FALSE, message=FALSE, fig.height=5.15, fig.width=9.95}
da.clusters <- unique(c(init.resplot.labels))

group.cols <- c("#2ca02c", "#1f77b4", "#9467bd", "#fed976", "#fd8d3c", "#e31a1c", "#800026", "#252525")
names(group.cols) <- c("Healthy", "LPS", "Non_covid", "Asymptomatic", "Mild", "Moderate", "Severe", "Critical")

# censor outliers to the 95th quantile for each cell type.
plotting.df <- cell.freq.merge[cell.freq.merge$Collection_Day %in% c("D0") &
                                   !cell.freq.merge$D0_status_summary %in% c("LPS", "Non_covid") &
                                   !is.na(cell.freq.merge$Days_from_onset) &
                                   cell.freq.merge$CellType %in% da.clusters, ]

q95.ct <- sapply(unique(plotting.df$CellType),
                 FUN=function(x) quantile(plotting.df[plotting.df$CellType %in% x, ]$Freq, c(0.05, 0.95), na.rm=TRUE)[2])
q95.df <- data.frame("CellType"=gsub(names(q95.ct), pattern="\\.95%", replacement=""), "Q95"=q95.ct)

plotting.df <- merge(plotting.df, q95.df, by='CellType')
plotting.df[plotting.df$Freq >= plotting.df$Q95, ]$Freq <- plotting.df[plotting.df$Freq >= plotting.df$Q95, ]$Q95

plotting.df$D0_status_summary <- ordered(plotting.df$D0_status_summary,
                                         levels=c("Healthy", "Asymptomatic", "Mild", "Moderate", "Severe", "Critical"))

ggplot(plotting.df[!plotting.df$D0_status_summary %in% c("Healthy", "Asymptomatic"), ],
       aes(x=Days_from_onset, y=Freq)) +
    #geom_boxplot(outlier.size=0.5, coef=1.5) +
    geom_point(aes(colour=D0_status_summary), size=1) +
    stat_smooth(method=MASS::rlm) +
    theme_cowplot() +
    facet_wrap(~CellType, scales="free_y", nrow=3) +
    scale_colour_manual(values=group.cols) +
    expand_limits(y=c(0)) +
    theme(aspect=1/1.3,
          legend.key.size=unit(0.4, "cm"),
          strip.background=element_rect(fill='white', colour='white'),
          strip.text=element_text(size=12)) +
    labs(x="Severity", y="Proportion") +
    guides(colour=guide_legend(title="Severity")) +
    ggsave("~/Dropbox/COVID19/plot.dir/All_proportions-DA_daysOnset-points.pdf",
           height=5.15, width=9.95, useDingbats=FALSE) +
    NULL
```

Many of these changes seem to be driven by critically ill patients.


```{r, warning=FALSE, message=FALSE}
write.table(init.res,
            file="~/Dropbox/COVID19/Data/SymptomOnset_resTable.txt",
            sep="\t", quote=FALSE, row.names=FALSE)

# print P-values to html table
clon.lm.pvalues <- data.frame("CellType"=init.res$CellType, "logFC"=init.res$logFC,
                              "Pvalue"=init.res$PValue,
                              "FDR"=init.res$FDR)

kbl(clon.lm.pvalues) %>% kable_paper(full_width=FALSE) %>% 
    save_kable("~/Dropbox/COVID19/plot.dir/DA_symptomOnset_trend_pvals.html",
               self_contained=TRUE)

kable(clon.lm.pvalues)
```


## Cluster-based DA analysis: days since symptom onset - excluding critical patients

```{r}
# set up testing model
rownames(covid.meta) <- covid.meta$sample_id
sanscrit.meta <- covid.meta[!covid.meta$D0_status_summary %in% c("Non_covid", "LPS", "Healthy", "Critical"), ]
sanscrit.meta$Days_from_onset[sanscrit.meta$Days_from_onset %in% c("Not_known")] <- NA
sanscrit.meta$Days_from_onset <- as.numeric(sanscrit.meta$Days_from_onset)
sanscrit.meta <- sanscrit.meta[!is.na(sanscrit.meta$Days_from_onset), ]

sanscrit.model <- model.matrix(~ Sex + Age + Site + Days_from_onset, 
                            data=sanscrit.meta[sanscrit.meta$Collection_Day %in% c("D0"), ])

# count cells
cell.freq.tab <- t(table(all.meta$sample_id[all.meta$Collection_Day %in% c("D0") &
                                                !is.na(all.meta$Days_from_onset) &
                                                !all.meta$D0_status_summary %in% c("LPS", "Non_covid", "Healthy", "Asymptomatic", "Critical")],
                         all.meta$initial_clustering[all.meta$Collection_Day %in% c("D0") &
                                                         !is.na(all.meta$Days_from_onset) &
                                                         !all.meta$D0_status_summary %in% c("LPS", "Non_covid", "Healthy", "Asymptomatic", "Critical")]))

cell.freq.tab <- cell.freq.tab[!rownames(cell.freq.tab) %in% c("Doublet", "Doublets", "Doublets:Bcell", "Doublets:Platelet"), ]
test.samps <- intersect(colnames(cell.freq.tab), names(n.cell.vecc))
cell.freq.tab <- cell.freq.tab[, test.samps]
sanscrit.model <- sanscrit.model[colnames(cell.freq.tab), ]

sanscrit.dge <- DGEList(cell.freq.tab, lib.size=log(n.cell.vecc[test.samps]))

#estimate dispersions
sanscrit.dge <- estimateDisp(sanscrit.dge, design=sanscrit.model)
sanscrit.linear.fit <- glmQLFit(sanscrit.dge, sanscrit.model, robust=TRUE)
sanscrit.res <- as.data.frame(topTags(glmQLFTest(sanscrit.linear.fit, coef=4), sort.by='none', n=Inf))
sanscrit.res$CellType <- rownames(sanscrit.res)

sanscrit.res$Sig <- as.numeric(sanscrit.res$FDR < 0.1)
sanscrit.res$Diff <- sign(sanscrit.res$logFC)
sanscrit.res$Diff[sanscrit.res$FDR > 0.1] <- 0
table(sanscrit.res$Diff)
```



```{r, warning=FALSE, message=FALSE, fig.height=2.95, fig.width=4.15}
mx.lfc <- max(abs(sanscrit.res$logFC))
eps <- mx.lfc * 0.05

sanscrit.resplot.labels <- sanscrit.res$CellType
sanscrit.resplot.labels[sanscrit.res$Sig == 0] <- ""

ggplot(sanscrit.res, aes(x=logFC, y=-log10(FDR), fill=as.character(Sig))) +
    geom_point(shape=21, size=4) +
    theme_cowplot() +
    scale_fill_manual(values=c("black", "red")) +
    scale_x_continuous(limits=c(-mx.lfc - eps, mx.lfc + eps), oob=squish) +
    guides(fill=guide_legend(title="FDR < 0.1")) +
    geom_text_repel(aes(label=sanscrit.resplot.labels),
                     fill='white') +
    labs(x="log fold change", y=expression(paste("-log"[10], " FDR"))) +
    ggsave("~/Dropbox/COVID19/plot.dir/All_edgeR_volcano_daysOnset-noCritical-linear.pdf",
           height=2.95, width=4.15, useDingbats=FALSE) +
    NULL
```


Show a plot of just the DA clusters.

```{r, warning=FALSE, message=FALSE, fig.height=5.15, fig.width=9.95}
da.clusters <- unique(c(sanscrit.resplot.labels))

group.cols <- c("#2ca02c", "#1f77b4", "#9467bd", "#fed976", "#fd8d3c", "#e31a1c", "#800026", "#252525")
names(group.cols) <- c("Healthy", "LPS", "Non_covid", "Asymptomatic", "Mild", "Moderate", "Severe", "Critical")

# censor outliers to the 95th quantile for each cell type.
plotting.df <- cell.freq.merge[cell.freq.merge$Collection_Day %in% c("D0") &
                                   !cell.freq.merge$D0_status_summary %in% c("LPS", "Non_covid") &
                                   !is.na(cell.freq.merge$Days_from_onset) &
                                   cell.freq.merge$CellType %in% da.clusters, ]

q95.ct <- sapply(unique(plotting.df$CellType),
                 FUN=function(x) quantile(plotting.df[plotting.df$CellType %in% x, ]$Freq, c(0.05, 0.95), na.rm=TRUE)[2])
q95.df <- data.frame("CellType"=gsub(names(q95.ct), pattern="\\.95%", replacement=""), "Q95"=q95.ct)

plotting.df <- merge(plotting.df, q95.df, by='CellType')
plotting.df[plotting.df$Freq >= plotting.df$Q95, ]$Freq <- plotting.df[plotting.df$Freq >= plotting.df$Q95, ]$Q95

plotting.df$D0_status_summary <- ordered(plotting.df$D0_status_summary,
                                         levels=c("Healthy", "Asymptomatic", "Mild", "Moderate", "Severe"))

ggplot(plotting.df[!plotting.df$D0_status_summary %in% c("Healthy", "Asymptomatic", "Critical"), ],
       aes(x=Days_from_onset, y=Freq)) +
    #geom_boxplot(outlier.size=0.5, coef=1.5) +
    geom_point(aes(colour=D0_status_summary), size=1) +
    stat_smooth(method=MASS::rlm) +
    theme_cowplot() +
    facet_wrap(~CellType, scales="free_y", nrow=3) +
    scale_colour_manual(values=group.cols) +
    expand_limits(y=c(0)) +
    theme(aspect=1/1.3,
          legend.key.size=unit(0.4, "cm"),
          strip.background=element_rect(fill='white', colour='white'),
          strip.text=element_text(size=12)) +
    labs(x="Severity", y="Proportion") +
    guides(colour=guide_legend(title="Symptom duration")) +
    ggsave("~/Dropbox/COVID19/plot.dir/All_proportions-DA_daysOnset-noCritical-points.pdf",
           height=5.15, width=9.95, useDingbats=FALSE) +
    NULL
```

Many of these changes seem to be driven by critically ill patients.


```{r, warning=FALSE, message=FALSE}
write.table(sanscrit.res,
            file="~/Dropbox/COVID19/Data/SymptomOnset-noCritical_resTable.txt",
            sep="\t", quote=FALSE, row.names=FALSE)

# print P-values to html table
clon.lm.pvalues <- data.frame("CellType"=sanscrit.res$CellType, "logFC"=sanscrit.res$logFC,
                              "Pvalue"=sanscrit.res$PValue,
                              "FDR"=sanscrit.res$FDR)

kbl(clon.lm.pvalues) %>% kable_paper(full_width=FALSE) %>% 
    save_kable("~/Dropbox/COVID19/plot.dir/DA_symptomOnset-noCritical_trend_pvals.html",
               self_contained=TRUE)

kable(clon.lm.pvalues)
```


Compare the log fold changes with and without the critical patients.

```{r}
colnames(init.res) <- c(paste0("Full.", colnames(init.res)[c(1:5)]), "CellType", "Full.Sig", "Full.Diff")
colnames(sanscrit.res) <- c(paste0("woCritical.", colnames(sanscrit.res)[c(1:5)]), "CellType", "woCritical.Sig", "woCritical.Diff")
```


```{r, warning=FALSE, message=FALSE, fig.height=2.95, fig.width=2.95}
compare.model <- merge(init.res, sanscrit.res, by='CellType')

compare.model$Sigs <- 0
compare.model$Sigs[compare.model$Full.Sig == 1 & compare.model$woCritical.Sig == 1] <- 1

cell.labels <- compare.model$CellType
cell.labels[compare.model$Sigs == 0] <- ""

ggplot(compare.model, aes(x=Full.logFC, y=woCritical.logFC)) +
    geom_hline(yintercept=0, lty=2) +
    geom_vline(xintercept=0, lty=2) +
    geom_point() +
    theme_cowplot() +
    labs(x="log Fold Change: Full model",
         y="log Fold Change: w/o Critical") +
    geom_text_repel(aes(label=cell.labels),
                    force=50) +
    ggsave("~/Dropbox/COVID19/plot.dir/SymptomOnset_compare-points.pdf",
           height=2.95, width=2.95, useDingbats=FALSE) +
    NULL
```


## Sensitivity analysis

Remove the outliers w.r.t. time since symptom onset to test the sensitivity.

## Cluster-based DA analysis: days since symptom onset - excluding > 24 days onset

```{r, fig.height=2.95, fig.width=3.95,}

ggplot(covid.meta[covid.meta$Collection_Day %in% c("D0") &
                      covid.meta$Days_from_onset <= 24 &
                      !covid.meta$D0_status_summary %in% c("LPS", "Non_covid", "Asymptomatic", "Healthy") &
                      !is.na(covid.meta$Days_from_onset), ],
       aes(x=OrderedSeverity, y=Days_from_onset)) +
    geom_boxplot() +
    labs(x="Disease severity", y="Days from\nsymptom onset") +
    theme_cowplot()  +
    theme(aspect=1,
          axis.text.x=element_text(angle=90, vjust=0.5, hjust=1)) +
    ggsave("~/Dropbox/COVID19/plot.dir/SymptomOnset_severity-sensitivityAnalysis-boxplot.pdf",
           height=2.95, width=3.95, useDingbats=FALSE) +
    NULL
```

This shows the distribution of symptom duration w.r.t. disease severity.

```{r}
# set up testing model
rownames(covid.meta) <- covid.meta$sample_id
symoutlier.meta <- covid.meta[!covid.meta$D0_status_summary %in% c("Non_covid", "LPS", "Healthy"), ]
symoutlier.meta$Days_from_onset[symoutlier.meta$Days_from_onset %in% c("Not_known")] <- NA
symoutlier.meta$Days_from_onset <- as.numeric(symoutlier.meta$Days_from_onset)
symoutlier.meta <- symoutlier.meta[!is.na(symoutlier.meta$Days_from_onset), ]
symoutlier.meta <- symoutlier.meta[symoutlier.meta$Days_from_onset <= 24, ]

symoutlier.model <- model.matrix(~ Sex + Age + Site + Days_from_onset, 
                            data=symoutlier.meta[symoutlier.meta$Collection_Day %in% c("D0"), ])

# count cells
cell.freq.tab <- t(table(all.meta$sample_id[all.meta$Collection_Day %in% c("D0") &
                                                !is.na(all.meta$Days_from_onset) &
                                                all.meta$Days_from_onset <= 24 &
                                                !all.meta$D0_status_summary %in% c("LPS", "Non_covid", "Healthy")],
                         all.meta$initial_clustering[all.meta$Collection_Day %in% c("D0") &
                                                         !is.na(all.meta$Days_from_onset) &
                                                         all.meta$Days_from_onset <= 24 &
                                                         !all.meta$D0_status_summary %in% c("LPS", "Non_covid", "Healthy")]))

cell.freq.tab <- cell.freq.tab[!rownames(cell.freq.tab) %in% c("Doublet", "Doublets", "Doublets:Bcell", "Doublets:Platelet"), ]
test.samps <- intersect(colnames(cell.freq.tab), names(n.cell.vecc))
cell.freq.tab <- cell.freq.tab[, test.samps]
symoutlier.model <- symoutlier.model[colnames(cell.freq.tab), ]

symoutlier.dge <- DGEList(cell.freq.tab, lib.size=log(n.cell.vecc[test.samps]))

#estimate dispersions
symoutlier.dge <- estimateDisp(symoutlier.dge, design=symoutlier.model)
symoutlier.linear.fit <- glmQLFit(symoutlier.dge, symoutlier.model, robust=TRUE)
symoutlier.res <- as.data.frame(topTags(glmQLFTest(symoutlier.linear.fit, coef=4), sort.by='none', n=Inf))
symoutlier.res$CellType <- rownames(symoutlier.res)

symoutlier.res$Sig <- as.numeric(symoutlier.res$FDR < 0.1)
symoutlier.res$Diff <- sign(symoutlier.res$logFC)
symoutlier.res$Diff[symoutlier.res$FDR > 0.1] <- 0
table(symoutlier.res$Diff)
```



```{r, warning=FALSE, message=FALSE, fig.height=2.95, fig.width=4.15}
mx.lfc <- max(abs(symoutlier.res$logFC))
eps <- mx.lfc * 0.05

symoutlier.resplot.labels <- symoutlier.res$CellType
symoutlier.resplot.labels[symoutlier.res$Sig == 0] <- ""

ggplot(symoutlier.res, aes(x=logFC, y=-log10(FDR), fill=as.character(Sig))) +
    geom_point(shape=21, size=4) +
    theme_cowplot() +
    scale_fill_manual(values=c("black", "red")) +
    scale_x_continuous(limits=c(-mx.lfc - eps, mx.lfc + eps), oob=squish) +
    guides(fill=guide_legend(title="FDR < 0.1")) +
    geom_text_repel(aes(label=symoutlier.resplot.labels),
                     fill='white') +
    labs(x="log fold change", y=expression(paste("-log"[10], " FDR"))) +
    ggsave("~/Dropbox/COVID19/plot.dir/All_edgeR_volcano_daysOnset-SensitivityAnalysis-linear.pdf",
           height=2.95, width=4.15, useDingbats=FALSE) +
    NULL
```


Show a plot of just the DA clusters.

```{r, warning=FALSE, message=FALSE, fig.height=5.15, fig.width=9.95}
da.clusters <- unique(c(symoutlier.resplot.labels))

group.cols <- c("#2ca02c", "#1f77b4", "#9467bd", "#fed976", "#fd8d3c", "#e31a1c", "#800026", "#252525")
names(group.cols) <- c("Healthy", "LPS", "Non_covid", "Asymptomatic", "Mild", "Moderate", "Severe", "Critical")

# censor outliers to the 95th quantile for each cell type.
plotting.df <- cell.freq.merge[cell.freq.merge$Collection_Day %in% c("D0") &
                                   !cell.freq.merge$D0_status_summary %in% c("LPS", "Non_covid") &
                                   !is.na(cell.freq.merge$Days_from_onset) &
                                   cell.freq.merge$CellType %in% da.clusters, ]

q95.ct <- sapply(unique(plotting.df$CellType),
                 FUN=function(x) quantile(plotting.df[plotting.df$CellType %in% x, ]$Freq, c(0.05, 0.95), na.rm=TRUE)[2])
q95.df <- data.frame("CellType"=gsub(names(q95.ct), pattern="\\.95%", replacement=""), "Q95"=q95.ct)

plotting.df <- merge(plotting.df, q95.df, by='CellType')
plotting.df[plotting.df$Freq >= plotting.df$Q95, ]$Freq <- plotting.df[plotting.df$Freq >= plotting.df$Q95, ]$Q95

plotting.df$D0_status_summary <- ordered(plotting.df$D0_status_summary,
                                         levels=c("Healthy", "Asymptomatic", "Mild", "Moderate", "Severe"))

ggplot(plotting.df[!plotting.df$D0_status_summary %in% c("Healthy") &
                       plotting.df$Days_from_onset <= 24, ],
       aes(x=Days_from_onset, y=Freq)) +
    #geom_boxplot(outlier.size=0.5, coef=1.5) +
    geom_point(aes(colour=D0_status_summary), size=1) +
    stat_smooth(method=MASS::rlm) +
    theme_cowplot() +
    facet_wrap(~CellType, scales="free_y", nrow=3) +
    scale_colour_manual(values=group.cols) +
    expand_limits(y=c(0)) +
    theme(aspect=1/1.3,
          legend.key.size=unit(0.4, "cm"),
          strip.background=element_rect(fill='white', colour='white'),
          strip.text=element_text(size=12)) +
    labs(x="Symptom duration", y="Proportion") +
    guides(colour=guide_legend(title="Symptom duration")) +
    ggsave("~/Dropbox/COVID19/plot.dir/All_proportions-DA_daysOnset-sensitivityAnalysis-points.pdf",
           height=5.15, width=9.95, useDingbats=FALSE) +
    NULL
```

This shows the relationship between symptom duration and cell type abundance, excluding the outlier patients.


```{r, warning=FALSE, message=FALSE}
write.table(symoutlier.res,
            file="~/Dropbox/COVID19/Data/SymptomOnset-sensitivityAnalysis_resTable.txt",
            sep="\t", quote=FALSE, row.names=FALSE)

# print P-values to html table
clon.lm.pvalues <- data.frame("CellType"=symoutlier.res$CellType, "logFC"=symoutlier.res$logFC,
                              "Pvalue"=symoutlier.res$PValue,
                              "FDR"=symoutlier.res$FDR)

kbl(clon.lm.pvalues) %>% kable_paper(full_width=FALSE) %>% 
    save_kable("~/Dropbox/COVID19/plot.dir/DA_symptomOnset-sensitivityAnalysis_trend_pvals.html",
               self_contained=TRUE)

kable(clon.lm.pvalues)
```


Compare the log fold changes with and without the critical patients.

```{r}
# colnames(init.res) <- c(paste0("Full.", colnames(init.res)[c(1:5)]), "CellType", "Full.Sig", "Full.Diff")
colnames(symoutlier.res) <- c(paste0("Ltday24.", colnames(symoutlier.res)[c(1:5)]), "CellType", "woCritical.Sig", "woCritical.Diff")
```


```{r, warning=FALSE, message=FALSE, fig.height=2.95, fig.width=3.35}
compare.model <- merge(init.res, symoutlier.res, by='CellType')

compare.model$Sigs <- 0
compare.model$Sigs[compare.model$Full.Sig == 1 & compare.model$woCritical.Sig == 1] <- 1

cell.labels <- compare.model$CellType
cell.labels[compare.model$Sigs == 0] <- ""

ggplot(compare.model, aes(x=Full.logFC, y=Ltday24.logFC)) +
    geom_hline(yintercept=0, lty=2) +
    geom_vline(xintercept=0, lty=2) +
    geom_point() +
    theme_cowplot() +
    labs(x="log Fold Change: Full model",
         y="log Fold Change: Symptoms\n<24 days") +
    geom_text_repel(aes(label=cell.labels),
                    force=50) +
    ggsave("~/Dropbox/COVID19/plot.dir/SymptomOnset_compare-sensitivityAnalysis-points.pdf",
           height=2.95, width=3.35, useDingbats=FALSE) +
    NULL
```



