---

title: "COVID19: T cell differential abundance"

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.

```{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)

```

```{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)

cell.freq.merge <- read.table("~/Dropbox/COVID19/Data/Tcell_proportions.tsv",

                              sep="\t", header=TRUE, stringsAsFactors=FALSE)

cell.freq.merge <- cell.freq.merge[!cell.freq.merge$D0_status_summary %in% c("Non_covid"), ]

cell.freq.merge <- cell.freq.merge[!cell.freq.merge$patient_id %in% c("CV0198"), ]

# remove doublets, etc

cell.freq.merge <- cell.freq.merge[!cell.freq.merge$CellType %in% c("Doublets:Bcell", "Doublets:Platelet", "Doublets", "NK", "ILCs"), ]

```

```{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

# 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)

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

```

```{r}

tcell.annots <- read.table("~/Dropbox/COVID19/Data/Tcell_annotations_ext.tsv",

                           sep="\t", header=TRUE, stringsAsFactors=FALSE)

tcell.df.merge <- merge(tcell.annots, all.meta, by='CellID')

```

## Cluster-baed DA analysis: Infected vs. controls

```{r}

# set up testing model

rownames(covid.meta) <- covid.meta$sample_id

tcell.meta <- covid.meta[!covid.meta$D0_status_summary %in% c("Non_covid", "LPS"), ]

tcell.meta$OrderedSeverity <- ordered(tcell.meta$D0_status_summary,

                                      levels=c("Healthy", "Asymptomatic", "Mild", "Moderate", "Severe", "Critical"))

tcell.meta$Infected <- ordered(tcell.meta$Status,

                               levels=c("Healthy", "Covid"))

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

# count cells

cell.freq.tab <- t(table(tcell.df.merge$sample_id[tcell.df.merge$Collection_Day %in% c("D0") &

                                                        !tcell.df.merge$D0_status_summary %in% c("LPS", "Non_covid")],

                         tcell.df.merge$Sub.Annotation[tcell.df.merge$Collection_Day %in% c("D0") &

                                                        !tcell.df.merge$D0_status_summary %in% c("LPS", "Non_covid")]))

cell.freq.tab <- cell.freq.tab[!rownames(cell.freq.tab) %in% c("NK", "ILCs", "Doublets", "Doublets:Bcell", "Doublets:Platelet"), ]

test.samps <- intersect(colnames(cell.freq.tab), names(n.cell.vecc))

cell.freq.tab <- cell.freq.tab[, test.samps]

tcell.model <- tcell.model[colnames(cell.freq.tab), ]

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

#estimate dispersions

tcell.dge <- estimateDisp(tcell.dge, design=tcell.model)

tcell.linear.fit <- glmQLFit(tcell.dge, tcell.model, robust=TRUE)

tcell.res <- as.data.frame(topTags(glmQLFTest(tcell.linear.fit, coef=4), sort.by='none', n=Inf))

tcell.res$CellType <- rownames(tcell.res)

tcell.res$Sig <- as.numeric(tcell.res$FDR < 0.1)

tcell.res$Diff <- sign(tcell.res$logFC)

tcell.res$Diff[tcell.res$FDR > 0.1] <- 0

table(tcell.res$Diff)

```

```{r, warning=FALSE, message=FALSE, fig.height=2.95, fig.width=4.15}

mx.lfc <- max(abs(tcell.res$logFC))

eps <- mx.lfc * 0.05

infect.resplot.labels <- tcell.res$CellType

infect.resplot.labels[tcell.res$Sig == 0] <- ""

ggplot(tcell.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=infect.resplot.labels),

                     fill='white') +

    labs(x="log fold change", y=expression(paste("-log"[10], " FDR"))) +

    ggsave("~/Dropbox/COVID19/plot.dir/Tcell_edgeR_volcano_caseVscontrol-linear.pdf",

           height=2.95, width=4.15, useDingbats=FALSE) +

    NULL

```

## Cluster-baed DA analysis: COVID19 severity

```{r}

# set up testing model

rownames(covid.meta) <- covid.meta$sample_id

tcell.meta <- covid.meta[!covid.meta$D0_status_summary %in% c("Non_covid", "LPS", "Healthy"), ]

tcell.meta$OrderedSeverity <- ordered(tcell.meta$D0_status_summary,

                                      levels=c("Asymptomatic", "Mild", "Moderate", "Severe", "Critical"))

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

# count cells

cell.freq.tab <- t(table(tcell.df.merge$sample_id[tcell.df.merge$Collection_Day %in% c("D0") &

                                                        !tcell.df.merge$D0_status_summary %in% c("LPS", "Non_covid", "Healthy")],

                         tcell.df.merge$Sub.Annotation[tcell.df.merge$Collection_Day %in% c("D0") &

                                                        !tcell.df.merge$D0_status_summary %in% c("LPS", "Non_covid", "Healthy")]))

cell.freq.tab <- cell.freq.tab[!rownames(cell.freq.tab) %in% c("NK", "ILCs", "Doublets", "Doublets:Bcell", "Doublets:Platelet"), ]

test.samps <- intersect(colnames(cell.freq.tab), names(n.cell.vecc))

cell.freq.tab <- cell.freq.tab[, test.samps]

tcell.model <- tcell.model[colnames(cell.freq.tab), ]

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

#estimate dispersions

tcell.dge <- estimateDisp(tcell.dge, design=tcell.model)

tcell.linear.fit <- glmQLFit(tcell.dge, tcell.model, robust=TRUE)

tcell.res <- as.data.frame(topTags(glmQLFTest(tcell.linear.fit, coef=4), sort.by='none', n=Inf))

tcell.res$CellType <- rownames(tcell.res)

tcell.res$Sig <- as.numeric(tcell.res$FDR < 0.1)

tcell.res$Diff <- sign(tcell.res$logFC)

tcell.res$Diff[tcell.res$FDR > 0.1] <- 0

table(tcell.res$Diff)

```

```{r, warning=FALSE, message=FALSE, fig.height=2.95, fig.width=4.15}

mx.lfc <- max(abs(tcell.res$logFC))

eps <- mx.lfc * 0.05

tcell.resplot.labels <- tcell.res$CellType

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

ggplot(tcell.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=tcell.resplot.labels),

                     fill='white') +

    labs(x="log fold change", y=expression(paste("-log"[10], " FDR"))) +

    ggsave("~/Dropbox/COVID19/plot.dir/Tcell_edgeR_volcano_caseOnly-linear.pdf",

           height=2.95, width=4.15, useDingbats=FALSE) +

    NULL

```

Quadratic changes.

```{r}

tcell.quad.res <- as.data.frame(topTags(glmQLFTest(tcell.linear.fit, coef=5), sort.by='none', n=Inf))

tcell.quad.res$CellType <- rownames(tcell.quad.res)

tcell.quad.res$Sig <- as.numeric(tcell.quad.res$FDR < 0.1)

tcell.quad.res$Diff <- sign(tcell.quad.res$logFC)

tcell.quad.res$Diff[tcell.quad.res$FDR > 0.1] <- 0

table(tcell.quad.res$Diff)

```

```{r, warning=FALSE, message=FALSE, fig.height=2.95, fig.width=4.15}

mx.lfc <- max(abs(tcell.quad.res$logFC))

eps <- mx.lfc * 0.05

tcell.quad.resplot.labels <- tcell.quad.res$CellType

tcell.quad.resplot.labels[tcell.quad.res$Sig == 0] <- ""

ggplot(tcell.quad.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=tcell.quad.resplot.labels),

                    force=10,

                    fill='white') +

    labs(x="log fold change", y=expression(paste("-log"[10], " FDR"))) +

    ggsave("~/Dropbox/COVID19/plot.dir/Tcell_edgeR_volcano_caseOnly-quadratic.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=4.15, fig.width=9.95}

da.clusters <- setdiff(unique(c(tcell.resplot.labels, tcell.quad.resplot.labels)), "CD4.Th17")

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") &

                           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"), ],

       aes(x=D0_status_summary, y=Freq, fill=D0_status_summary)) +

    geom_boxplot(outlier.size=0.5, coef=1.5) +

    theme_cowplot() +

    facet_wrap(~CellType, scales="free_y", nrow=3) +

    scale_fill_manual(values=group.cols) +

    expand_limits(y=c(0)) +

    theme(aspect=1/1.3,

          axis.text.x=element_blank(),

          axis.ticks.x=element_blank(),

          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(fill=guide_legend(title="Severity")) +

    ggsave("~/Dropbox/COVID19/plot.dir/Tcells_proportions-DA_CasesOnly-boxplot.pdf",

           height=4.15, width=9.95, useDingbats=FALSE) +

    NULL

```

Also show the equivalent plot of healthy vs. infected.

```{r, warning=FALSE, message=FALSE, fig.height=4.15, fig.width=9.95}

infect.da.clusters <- setdiff(unique(c(infect.resplot.labels)), "CD4.Th17")

group.cols <- c("#2ca02c", "#252525")

names(group.cols) <- c("Healthy", "COVID-19")

# 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") &

                           cell.freq.merge$CellType %in% infect.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$Status <- factor(plotting.df$Status,

                             levels=c("Healthy", "Covid"),

                             labels=c("Healthy", "COVID-19"))

ggplot(plotting.df,

       aes(x=Status, y=Freq, fill=Status)) +

    geom_boxplot(outlier.size=0.5, coef=1.5) +

    theme_cowplot() +

    facet_wrap(~CellType, scales="free_y", nrow=3) +

    scale_fill_manual(values=group.cols) +

    expand_limits(y=c(0)) +

    theme(aspect=1/1.3,

          axis.text.x=element_blank(),

          axis.ticks.x=element_blank(),

          legend.key.size=unit(0.4, "cm"),

          strip.background=element_rect(fill='white', colour='white'),

          strip.text=element_text(size=12)) +

    labs(x="", y="Proportion") +

    guides(fill=guide_legend(title="Status")) +

    ggsave("~/Dropbox/COVID19/plot.dir/Tcells_proportions-DA_CaseVsControls-boxplot.pdf",

           height=4.15, width=9.95, useDingbats=FALSE) +

    NULL

```

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

```{r}

# set up testing model

rownames(covid.meta) <- covid.meta$sample_id

tcell.meta <- covid.meta[!covid.meta$D0_status_summary %in% c("Non_covid", "LPS", "Healthy"), ]

tcell.meta$OrderedSeverity <- ordered(tcell.meta$D0_status_summary,

                                      levels=c("Asymptomatic", "Mild", "Moderate", "Severe", "Critical"))

tcell.meta <- tcell.meta[tcell.meta$D0_status_summary %in% c("Mild", "Moderate", "Severe", "Critical"), ]

tcell.meta$Days_from_onset[tcell.meta$Days_from_onset %in% c("Not_known")] <- NA

tcell.meta$Days_from_onset <- as.numeric(tcell.meta$Days_from_onset)

tcell.meta <- tcell.meta[!is.na(tcell.meta$Days_from_onset), ]

tcell.model <- model.matrix(~ Sex + Age + Days_from_onset, 

                            data=tcell.meta[tcell.meta$Collection_Day %in% c("D0"), ])

# count cells

cell.freq.tab <- t(table(tcell.df.merge$sample_id[tcell.df.merge$Collection_Day %in% c("D0") &

                                                      !is.na(tcell.df.merge$Days_from_onset) &

                                                        !tcell.df.merge$D0_status_summary %in% c("LPS", "Non_covid", "Healthy", "Asymptomatic")],

                         tcell.df.merge$Sub.Annotation[tcell.df.merge$Collection_Day %in% c("D0") &

                                                           !is.na(tcell.df.merge$Days_from_onset) &

                                                        !tcell.df.merge$D0_status_summary %in% c("LPS", "Non_covid", "Healthy", "Asymptomatic")]))

cell.freq.tab <- cell.freq.tab[!rownames(cell.freq.tab) %in% c("NK", "ILCs", "Doublets", "Doublets:Bcell", "Doublets:Platelet"), ]

test.samps <- intersect(colnames(cell.freq.tab), names(n.cell.vecc))

cell.freq.tab <- cell.freq.tab[, test.samps]

tcell.model <- tcell.model[colnames(cell.freq.tab), ]

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

#estimate dispersions

tcell.dge <- estimateDisp(tcell.dge, design=tcell.model)

tcell.linear.fit <- glmQLFit(tcell.dge, tcell.model, robust=TRUE)

tcell.res <- as.data.frame(topTags(glmQLFTest(tcell.linear.fit, coef=4), sort.by='none', n=Inf))

tcell.res$CellType <- rownames(tcell.res)

tcell.res$Sig <- as.numeric(tcell.res$FDR < 0.1)

tcell.res$Diff <- sign(tcell.res$logFC)

tcell.res$Diff[tcell.res$FDR > 0.1] <- 0

table(tcell.res$Diff)

```

```{r, warning=FALSE, message=FALSE, fig.height=2.95, fig.width=4.15}

mx.lfc <- max(abs(tcell.res$logFC))

eps <- mx.lfc * 0.05

tcell.resplot.labels <- tcell.res$CellType

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

ggplot(tcell.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=tcell.resplot.labels),

                     fill='white') +

    labs(x="log fold change", y=expression(paste("-log"[10], " FDR"))) +

    ggsave("~/Dropbox/COVID19/plot.dir/Tcell_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=4.15, fig.width=9.95}

da.clusters <- setdiff(unique(c(tcell.resplot.labels, tcell.quad.resplot.labels)), "CD4.Th17")

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") &

                                   !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(size=1) +

    stat_smooth() +

    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,

          axis.text.x=element_blank(),

          axis.ticks.x=element_blank(),

          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(fill=guide_legend(title="Severity")) +

    ggsave("~/Dropbox/COVID19/plot.dir/Tcells_proportions-DA_daysOnset-points.pdf",

           height=4.15, width=9.95, useDingbats=FALSE) +

    NULL

```