424 lines (331 with data), 18.5 kB
---
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
```
Datasets
Models
