1145 lines (906 with data), 57.5 kB
---
title: "COVID19: T cell annotation - updated"
output: html_notebook
---
Annotating T cells subtypes.
```{r, warning=FALSE, message=FALSE}
library(ggplot2)
library(scran)
library(ggsci)
library(cowplot)
library(scattermore)
library(ComplexHeatmap)
library(RColorBrewer)
library(pheatmap)
library(scater)
library(reshape2)
library(colorspace)
library(dplyr)
library(ggthemes)
library(edgeR)
library(scales)
library(ggrepel)
library(viridis)
library(MASS)
```
I've taken all of the full clusters and selected the ones that I think contain T cells.
```{r, warning=FALSE, message=FALSE}
# tcell.umap <- read.table("~/Dropbox/COVID19/Data/Updated/Cambridge_Sanger_Newcastle_Tcell_UMAPs.tsv",
# sep="\t", header=TRUE, stringsAsFactors=FALSE)
# fix the Sanger barcodes
# tcell.umap$CellID[tcell.umap$DataSet %in% c("Sanger")] <- gsub(tcell.umap$CellID[tcell.umap$DataSet %in% c("Sanger")],
# pattern="^(CV[0-9]+_[A-Z0-9]+-CV[0-9]+_[A-Z0-9]+)_", replacement="")
tcell.clusters <- read.table("~/Dropbox/COVID19/Data/Updated/Cambridge_Sanger_Newcastle_Tcell_clusters.tsv",
sep="\t", header=TRUE, stringsAsFactors=FALSE)
# tcell.meta <- merge(tcell.umap, tcell.clusters, by='CellID')
```
```{r, warning=FALSE, message=FALSE}
all.meta <- read.table("~/Dropbox/COVID19/Data/Updated/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"), ]
```
To normalise the cell type proportions, I need to know how many cells were captured for each patient-sample combo.
```{r, warning=FALSE, message=FALSE}
n.cell.vecc <- table(all.meta$sample_id)
```
```{r, warning=FALSE, message=FALSE}
# tcell.merge <- merge(tcell.clusters, tcell.umap, by='CellID')
tcell.merge <- tcell.clusters
# tcell.merge <- merge(tcell.clusters, all.meta, by='CellID')
# tcell.merge <- merge(tcell.merge, tcell.umap, by='CellID')
```
```{r}
gene.clust <- read.table("~/Dropbox/COVID19/Data/Updated/Cambridge_Sanger_Newcastle-Tcells-cluster_gene-Ave.tsv",
sep="\t", header=TRUE, stringsAsFactors=FALSE)
protein.clust <- read.table("~/Dropbox/COVID19/Data/Updated/Cambridge_Sanger_Newcastle-Tcells-cluster_protein-Ave.tsv",
sep="\t", header=TRUE, stringsAsFactors=FALSE)
```
```{r, warning=FALSE, message=FALSE}
plot.proteins <- unique(c("CD19", "CD86", "CD11b", "CD56", "CD8", "CD4", "CD3", "CD16", "CD14", "CD335", "CD158", "SIGLEC7",
"FCER2", "TCR", "CTLA4", "CXCR5", "ICOSL", "CD40LG",
"CX3CR1", "CD79b", "CD20", "CD1c", "CD25", "CD45RO", "HLA-DR", "CD21", "CD5", "CD64", "CD127",
"CD34", "CD40", "CD123", "CD27", "CD1c", "TCR_Va7.2", "TCR_Vg2", "CD235ab", "CD10", "CD2", "CD71",
"TCR_Vg9", "TCR_Va24-Ja18", "CCR7", "CD45RA", "CD62L", "CD138", "CD38", "CD274", "CD279", "TIGIT"))
plot.genes <- unique(c("CD3G", "CD3E", "FCGR3A", "CD8A", "CD8B", "CD4", "CD80", "CD86", "CXCL10", "HLA.DRA", "HLA.DRB1", "VEGFA",
"CCR7", "CCR2", "CCR5", "CXCR3", "CD36", "FOXP3", "MS4A1", "LYZ", "GZMB", "CD68", "GATA3", "PECAM1", "RORC",
"FCGR2A", "FCGR2B", "FCGR3B", "TPO", "HBB", "NCR1", "NCR2", "F8", "HIF1A", "CTLA4", "CD28", "CXCR1",
"PDCD1", "TIM3", "LAG3", "TYMS", "PITX1", "TIGIT", "CXCL13", "BCL6", "TRAV1.2",
#"IGHD", "IGHM", "IGHA1", "IGHE", "IGHG1", "IGHG2", "IGHG3", "IGHG4",
"ITGA2B", "ITGAX", "TOX", "FYN", "NCAM1", "LRRC32", "SELL", "CXCR5", "CSF1", "IL7R", "CD40LG", "MYC", "MKI67",
"IFNG", "IL6", "THRB", "CD19", "CD38", "CD27", "F2", "F5", "THRB", "CD44", "IKZF2", "IKZF4"))
```
```{r, warning=FALSE, message=FALSE, fig.height=10.95, fig.width=9.95}
gene.plot.hm <- t(apply(t(gene.clust[, intersect(colnames(gene.clust), plot.genes)]), 1,
FUN=function(XP) XP/max(XP)))
gene.plot.hm[is.na(gene.plot.hm)] <- 0
colnames(gene.plot.hm) <- gene.clust$GeneID
pdf("~/Dropbox/COVID19/Updated_plots/Gene_heatmap-Tcells.pdf",
height=10.95, width=9.95, useDingbats=FALSE)
Heatmap(gene.plot.hm, name="Z-score")
dev.off()
Heatmap(gene.plot.hm, name="Z-score")
```
I will take these clusters and compare them to the profiles of the original annotations.
```{r, warning=FALSE, message=FALSE}
table(tcell.merge$Centre, tcell.merge$Sub.Cluster)
```
```{r, warning=FALSE, message=FALSE}
table(tcell.merge$initial_clustering, tcell.merge$Sub.Cluster)
```
```{r, warning=FALSE, message=FALSE, fig.height=10.95, fig.width=9.95}
rownames(protein.clust) <- protein.clust$Protein
protein.plot.hm <- protein.clust[intersect(protein.clust$Protein, plot.proteins),
setdiff(colnames(protein.clust), "Protein")]
#protein.plot.hm[protein.plot.hm < 0] <- 0
# gene.plot.hm <- t(apply(t(gene.clust[, intersect(colnames(gene.clust), plot.genes)]), 1,
# FUN=function(XP) XP/max(XP)))
# protein.plot.hm <- t(apply(t(protein.plot.hm), 1, FUN=function(XP) (XP - mean(XP))/sd(XP)))
protein.plot.hm[protein.plot.hm <= 0] <- 0
protein.plot.hm <- t(apply(protein.plot.hm, 1, FUN=function(XP) (XP - mean(XP))/sd(XP)))
protein.plot.hm[is.na(protein.plot.hm)] <- 0
rownames(protein.plot.hm) <- intersect(protein.clust$Protein, plot.proteins)
# pdf("~/Dropbox/COVID19/Updated_plots/protein_heatmap-celltypes.pdf",
# height=10.95, width=9.95, useDingbats=FALSE)
# Heatmap(protein.plot.hm, name="Z-score")
# dev.off()
Heatmap(protein.plot.hm, name="Z-score")
```
Compare the previous annotations.
```{r, warning=FALSE, message=FALSE}
old.annots <- read.csv("~/Dropbox/COVID19/Data/Annotations/annotation_200112.csv",
header=TRUE, stringsAsFactors=FALSE)
old.clusters <- read.table("~/Dropbox/COVID19/Data/Annotations/Cambridge_Sanger_Newcastle_Tcell_clusters.tsv",
sep="\t", header=TRUE, stringsAsFactors = TRUE)
colnames(old.clusters) <- c("CellID", "old.Sub.Cluster", "old.Louvain.Cluster")
tcell.annot.merge <- merge(tcell.merge, old.annots, by.x='CellID', by.y='X')
tcell.annot.merge <- merge(tcell.annot.merge, old.clusters, by='CellID')
```
```{r, warning=FALSE, message=FALSE}
table(tcell.annot.merge$full_clustering, tcell.annot.merge$Sub.Cluster)
```
```{r, warning=FALSE, message=FALSE}
table(tcell.annot.merge$old.Sub.Cluster, tcell.annot.merge$Sub.Cluster)
```
```{r, fig.height=4.95, fig.width=5.15}
tcell.clusts <- c("CD4.Prolif", "CD8.Prolif", "NK_prolif", "ILC1_3", "ILC2", "CD4.Naive", "CD8.Naive",
"CD4.Th17", "CD4.EM", "CD4.IL22", "CD4.Tfh", "CD4.CM", "CD4.Th1", "CD4.Th2", "Treg",
"gdT", "MAIT", "NKT", "CD8.EM", "CD8.TE", "NK_56hi", "NK_16hi")
comp.tab <- table(tcell.annot.merge$Sub.Cluster, tcell.annot.merge$full_clustering)
comp.tab <- apply(comp.tab[, tcell.clusts], 1, FUN=function(X) X/sum(X))
pheatmap(comp.tab, cluster_rows = FALSE)
```
This shows that the clustering only really identifies the major subsets.
```{r}
tcell.merge$Annotation <- "Unknown"
tcell.merge$Annotation[tcell.merge$Sub.Cluster %in% c()] <- "PoorQC"
tcell.merge$Annotation[tcell.merge$Sub.Cluster %in% c(15)] <- "NK"
tcell.merge$Annotation[tcell.merge$Sub.Cluster %in% c()] <- "Doublets"
tcell.merge$Annotation[tcell.merge$Sub.Cluster %in% c()] <- "ILCs"
tcell.merge$Annotation[tcell.merge$Sub.Cluster %in% c(9)] <- "gdT"
tcell.merge$Annotation[tcell.merge$Sub.Cluster %in% c(8, 10, 19)] <- "CD8.TE"
tcell.merge$Annotation[tcell.merge$Sub.Cluster %in% c(17)] <- "CD8.EM"
tcell.merge$Annotation[tcell.merge$Sub.Cluster %in% c(1, 29)] <- "Tfh-like"
tcell.merge$Annotation[tcell.merge$Sub.Cluster %in% c(2)] <- "Treg"
tcell.merge$Annotation[tcell.merge$Sub.Cluster %in% c(21)] <- "MAIT"
tcell.merge$Annotation[tcell.merge$Sub.Cluster %in% c(4, 5, 16, 18, 20, 22, 26)] <- "CD4.Naive"
tcell.merge$Annotation[tcell.merge$Sub.Cluster %in% c(6, 24)] <- "CD8.Naive"
tcell.merge$Annotation[tcell.merge$Sub.Cluster %in% c(3, 7, 12, 13, 14, 23, 27, 28)] <- "CD4.CM"
tcell.merge$Annotation[tcell.merge$Sub.Cluster %in% c(11)] <- "CD4.IL22"
tcell.merge$Annotation[tcell.merge$Sub.Cluster %in% c()] <- "Exhausted"
tcell.merge$Annotation[tcell.merge$Sub.Cluster %in% c(25)] <- "Proliferating"
table(tcell.merge$Annotation, tcell.merge$Sub.Cluster)
```
```{r, warning=FALSE, message=FALSE}
table(tcell.merge$Annotation, tcell.merge$Site)
```
Plot the markers genes for these groups: CD4, CD8, TCR V$\alpha$7.2, TCR V$\gamma$9, TCR V$\gamma$2, CD45RA, CCR7, CD45RO, CD27, CD28,
CTLA4, HLA-DR, FOXP3, CXCR5, GZMB, IFNG, TOX, PDCD1.
I think these annotations look OK.
```{r, warning=FALSE, message=FALSE}
# read in the gene expression
goi.df <- read.table("~/Dropbox/COVID19/Data/Updated/Cambridge_Sanger_Newcastle-Tcells_GeneGOI.tsv",
sep="\t", header=TRUE, stringsAsFactors=FALSE)
tcell.goi.merge <- merge(tcell.merge, goi.df, by='CellID')
```
```{r, warning=FALSE, message=FALSE}
# re-classify CD8.TE with high NCAM1 * FCGR3A expression
tcell.goi.merge$Annotation[tcell.goi.merge$NCAM1 > 0 & tcell.goi.merge$Annotation %in% c("CD8.TE", "NK") &
tcell.goi.merge$TRAV1.2 > 0] <- "NKT"
tcell.goi.merge$Annotation[tcell.goi.merge$FCGR3A > 0 & tcell.goi.merge$Annotation %in% c("CD8.TE", "NK") &
tcell.goi.merge$TRAV1.2 > 0] <- "NKT"
tcell.goi.merge$Annotation[tcell.goi.merge$FCGR3A > 0 & tcell.goi.merge$Annotation %in% c("CD8.TE") &
tcell.goi.merge$NCAM1 > 0] <- "NKT"
```
Remove doublets based on expression of markers for different cell lineages.
```{r}
# re-classify CD8.TE with high NCR1 * FCGR3A expression
bcell.markers <- tcell.goi.merge$IGHM > 0 | tcell.goi.merge$IGHG3 > 0 | tcell.goi.merge$IGHD > 0 | tcell.goi.merge$IGHG4 > 0 |
tcell.goi.merge$IGHG2 > 0 | tcell.goi.merge$IGHG1 > 0 | tcell.goi.merge$IGHA1 > 0
tcell.goi.merge$Annotation[bcell.markers] <- "Doublets"
table(tcell.goi.merge$Annotation, tcell.goi.merge$Site)
```
```{r, fig.height=9.75, fig.width=10.95}
n.cells <- length(unique(tcell.goi.merge$Sub.Cluster))
cell.cols <- colorRampPalette(pal_futurama()(12))(n.cells)
names(cell.cols) <- unique(tcell.goi.merge$Sub.Cluster)
tcell.goi.merge$Sub.Cluster <- factor(tcell.goi.merge$Sub.Cluster,
levels=c(1:29))
ggplot(tcell.goi.merge,
aes(x=UMAP1, y=UMAP2)) +
geom_scattermore(data=tcell.goi.merge[, c("UMAP1", "UMAP2")],
colour='grey80', alpha=0.5, size=0.1) +
geom_scattermore(aes(colour=Sub.Cluster)) +
theme_cowplot() +
facet_wrap(~Annotation) +
scale_colour_manual(values=cell.cols) +
guides(colour=guide_legend(override.aes = list(size=3)))
```
```{r, warning=FALSE, message=FALSE}
clust.model.matrix <- model.matrix(~ 0 + Annotation, data=tcell.goi.merge)
rownames(clust.model.matrix) <- tcell.goi.merge$CellID
gene.goi.by.cluster <- t(clust.model.matrix) %*% as.matrix(tcell.goi.merge[, setdiff(colnames(goi.df), colnames(tcell.merge))])
gene.ave.by.label <- as.data.frame(apply(gene.goi.by.cluster, 2, function(X) X/colSums(clust.model.matrix)))
rownames(gene.ave.by.label) <- gsub(rownames(gene.ave.by.label), pattern="Annotation", replacement="")
gene.goi.by.cluster.bin <- t(clust.model.matrix) %*% as.matrix(tcell.goi.merge[, setdiff(colnames(goi.df), colnames(tcell.merge))] > 0)
gene.bin.by.label <- as.data.frame(apply(gene.goi.by.cluster.bin, 2, function(X) X/sum(X)))
gene.bin.by.label[is.na(gene.bin.by.label)] <- 0
rownames(gene.bin.by.label) <- gsub(rownames(gene.bin.by.label), pattern="Annotation", replacement="")
```
```{r, message=FALSE, warning=FALSE}
gene.ave.label.df <- as.data.frame(apply(gene.ave.by.label, 2, FUN=function(XP) XP/max(XP)))
gene.ave.label.df$Annotation <- rownames(gene.ave.label.df)
gene.ave.label.melt <- melt(gene.ave.label.df, id.vars='Annotation')
colnames(gene.ave.label.melt) <- c("Annotation", "Gene", "AveExprs")
gene.bin.label.df <- gene.bin.by.label
gene.bin.label.df$Annotation <- rownames(gene.bin.label.df)
gene.bin.label.melt <- melt(gene.bin.label.df, id.vars='Annotation')
colnames(gene.bin.label.melt) <- c("Annotation", "Gene", "PropExprs")
dot.plot.df <- merge(gene.ave.label.melt, gene.bin.label.melt, by=c('Annotation', 'Gene'))
dot.plot.df$Gene <- as.character(dot.plot.df$Gene)
```
```{r, message=FALSE, warning=FALSE, fig.height=5.15, fig.width=8.25}
dot.plot.genes <- c("CD3E", "CD8A", "CD8B", "CD4", "CCR7", "CD27", "CD28", "TOX", "CTLA4", "IFNG", "PDCD1", "GZMB", "IFNG",
"TRAV1.2", "TRDV2", "TRGV9", "CD40LG", "RORC", "CD38", "NCAM1",
"CXCR5", "HLA.DRB1", "LAG3", "CD44", "FOXP3", "MKI67", "NCR1", "GARP", "GZMH", "CXCR4", "FCGR3A", "TCF7")
ggplot(dot.plot.df[dot.plot.df$Gene %in% dot.plot.genes, ],
aes(x=Gene, y=Annotation, colour=AveExprs, size=PropExprs*100)) +
geom_point() +
theme_cowplot() +
scale_colour_gradient(low='grey70', high='blue') +
theme(axis.text.x=element_text(angle=90, vjust=0.5, hjust=1, colour='black'),
axis.text.y=element_text(colour='black'),
panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
panel.background = element_blank(), axis.line = element_line(colour = "black")) +
guides(colour=guide_colorbar(title='Mean Expression'),
size=guide_legend(title='% Express')) +
# ggsave("~/Dropbox/COVID19/Updated_plots/Tcell_dotplot-genes.pdf",
# height=5.15, width=8.25, useDingbats=FALSE) +
NULL
```
Do the same for the proteins.
```{r, warning=FALSE, message=FALSE}
# read in the gene expression
prot.df <- read.table("~/Dropbox/COVID19/Data/Updated/Cambridge_Sanger_Newcastle-Tcells_ProteinGOI.tsv",
sep="\t", header=TRUE, stringsAsFactors=FALSE)
colnames(prot.df) <- c(paste0("Prot.", colnames(prot.df)[!colnames(prot.df) %in% c("CellID")]), "CellID")
tcell.prot.merge <- merge(tcell.goi.merge, prot.df, by=c('CellID'))
```
```{r, warning=FALSE, message=FALSE}
clust.model.matrix <- model.matrix(~ 0 + Annotation, data=tcell.prot.merge)
rownames(clust.model.matrix) <- tcell.prot.merge$CellID
goi.by.cluster <- t(clust.model.matrix) %*% as.matrix(tcell.prot.merge[, gsub(setdiff(colnames(prot.df), colnames(tcell.merge)),
pattern="\\.x", replacement="")])
ave.by.label <- as.data.frame(apply(goi.by.cluster, 2, function(X) X/colSums(clust.model.matrix)))
rownames(ave.by.label) <- gsub(rownames(ave.by.label), pattern="Annotation", replacement="")
colnames(ave.by.label) <- gsub(colnames(ave.by.label), pattern="Prot\\.", replacement="")
goi.by.cluster.bin <- t(clust.model.matrix) %*% as.matrix(tcell.prot.merge[, gsub(setdiff(colnames(prot.df), colnames(tcell.merge)),
pattern="\\.x", replacement="")] > 0)
bin.by.label <- as.data.frame(apply(goi.by.cluster.bin, 2, function(X) X/sum(X)))
bin.by.label[is.na(bin.by.label)] <- 0
rownames(bin.by.label) <- gsub(rownames(bin.by.label), pattern="Annotation", replacement="")
colnames(bin.by.label) <- gsub(colnames(bin.by.label), pattern="Prot\\.", replacement="")
```
```{r, message=FALSE, warning=FALSE}
ave.label.df <- as.data.frame(apply(ave.by.label, 2, FUN=function(XP) (XP-min(XP))/max(XP-min(XP))))
ave.label.df$Annotation <- rownames(ave.label.df)
ave.label.melt <- melt(ave.label.df, id.vars='Annotation')
colnames(ave.label.melt) <- c("Annotation", "Protein", "AveExprs")
bin.label.df <- bin.by.label
bin.label.df$Annotation <- rownames(bin.label.df)
bin.label.melt <- melt(bin.label.df, id.vars='Annotation')
colnames(bin.label.melt) <- c("Annotation", "Protein", "PropExprs")
protein.dot.plot.df <- merge(ave.label.melt, bin.label.melt, by=c('Annotation', 'Protein'))
protein.dot.plot.df$Protein <- as.character(protein.dot.plot.df$Protein)
```
```{r, message=FALSE, warning=FALSE, fig.height=5.15, fig.width=8.25}
dot.plot.prots <- c("CD3", "CD8", "CD4", "CCR7", "CD45RA", "CD45RO", "CD28", "CD27", "CTLA4", "CD56", "CD16", "ICOSL",
"CD274", "CXCR3", "CD40LG", "CD38", "ICOS", "PD1", "HAVCR2",
"CXCR5", "TCR_Va7.2", "TCR_Vg2", "TCR_Vg9", "TCR_Va24-Ja18", "CD62L", "CD44", "TCR_Va24.Ja18")
ggplot(protein.dot.plot.df[protein.dot.plot.df$Protein %in% dot.plot.prots, ],
aes(x=Protein, y=Annotation, colour=AveExprs, size=PropExprs*100)) +
geom_point() +
theme_cowplot() +
scale_colour_gradient(low='grey70', high=darken('red')) +
theme(axis.text.x=element_text(angle=90, vjust=0.5, hjust=1, colour='black'),
axis.text.y=element_text(colour='black'),
panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
panel.background = element_blank(), axis.line = element_line(colour = "black")) +
guides(colour=guide_colorbar(title='Mean Expression'),
size=guide_legend(title='% Express')) +
# ggsave("~/Dropbox/COVID19/Updated_plots/Tcell_dotplot-protein.pdf",
# height=5.15, width=8.25, useDingbats=FALSE) +
NULL
```
Write out the exhausted, proliferating and Tfh-like clusters for sub-clustering.
```{r}
write.table(tcell.prot.merge$CellID[tcell.prot.merge$Annotation %in% c("Proliferating")],
file="~/Dropbox/COVID19/Data/Updated/Proliferating_barcodes.tsv",
quote=FALSE, col.names=FALSE, row.names=FALSE, sep="\t")
write.table(tcell.prot.merge$CellID[tcell.prot.merge$Annotation %in% c("Tfh-like")],
file="~/Dropbox/COVID19/Data/Updated/Tfh_barcodes.tsv",
quote=FALSE, col.names=FALSE, row.names=FALSE, sep="\t")
write.table(tcell.prot.merge$CellID[tcell.prot.merge$Annotation %in% c("Exhausted")],
file="~/Dropbox/COVID19/Data/Updated/Exhausted_barcodes.tsv",
quote=FALSE, col.names=FALSE, row.names=FALSE, sep="\t")
write.table(tcell.prot.merge$CellID[tcell.prot.merge$Annotation %in% c("Doublets")],
file="~/Dropbox/COVID19/Data/Updated/Doublet_barcodes.tsv",
quote=FALSE, col.names=FALSE, row.names=FALSE, sep="\t")
```
Fill in the sub-annotations.
```{r, warning=FALSE, message=FALSE}
rownames(tcell.prot.merge) <- tcell.prot.merge$CellID
tcell.prot.merge$Sub.Annotation <- tcell.prot.merge$Annotation
# first, enforce the original annotations on all previously annotated cells
old.annots <- read.csv("~/Dropbox/COVID19/Data/Annotations/annotation_200112.csv",
header=TRUE, stringsAsFactors=FALSE)
rownames(old.annots) <- old.annots$X
annot.merge <- merge(tcell.prot.merge, old.annots, by.x='CellID', by.y='X', all.x=TRUE)
annot.merge$Sub.Annotation[!is.na(annot.merge$full_clustering)] <- annot.merge$full_clustering[!is.na(annot.merge$full_clustering)]
annot.merge <- annot.merge[!annot.merge$full_clustering %in% c("CD14_mono", "B_naive", "DC_prolif", "DC3", "HSC", "Mono_prolif",
"NK_16hi", "NK_56hi", "NK_prolif", "pDC", "Plasmablast",
"Platelets", "RBC"), ]
# re-classify CD4s within CD8 clusters
annot.merge$Sub.Annotation[(annot.merge$CD8A > 0 | annot.merge$CD8B > 0) &
annot.merge$CD4 <= 0 &
annot.merge$CXCR5 <= 0 &
annot.merge$Sub.Annotation %in% c("Tfh-like")] <- "CD8.TE"
# Tfh classification
annot.merge$Sub.Annotation[(annot.merge$CD8A <= 0 & annot.merge$CD8B <= 0) &
annot.merge$CD4 > 0 &
annot.merge$CXCR5 > 0 &
annot.merge$Sub.Annotation %in% c("CD8.TE")] <- "CD4.Tfh"
# fix the Treg annotations
annot.merge$Sub.Annotation[annot.merge$Sub.Annotation %in% c("CD4.CM") &
(annot.merge$CXCR5 > 0) & (annot.merge$CD40LG > 0)] <- "CD4.Tfh"
# final doublet classification
annot.merge$Sub.Annotation[annot.merge$IGHA1 > 0 | annot.merge$IGHD > 0 | annot.merge$IGHE > 0 | annot.merge$IGHG1 > 0 |
annot.merge$IGHG2 > 0 | annot.merge$IGHG3 > 0 | annot.merge$IGHG4 > 0] <- "Doublets"
# final NKT classification
annot.merge$Sub.Annotation[annot.merge$FCGR3A > 0 & annot.merge$NCAM1 > 0 &
annot.merge$Sub.Annotation %in% c("CD8.EM", "CD8.TE", "CD8.Naive")] <- "NKT"
annot.merge$Sub.Annotation[annot.merge$FCGR3A > 0 & annot.merge$NCR1 > 0 &
annot.merge$Sub.Annotation %in% c("CD8.EM", "CD8.TE", "CD8.Naive")] <- "NKT"
# fix the Treg annotations
annot.merge$Sub.Annotation[annot.merge$Sub.Annotation %in% c("CD4.CM", "CD4.Naive", "CD4.Tfh", "Treg") &
(annot.merge$IL22 > 0)] <- "CD4.IL22"
annot.merge$Sub.Annotation[annot.merge$Sub.Annotation %in% c("Treg") &
!(annot.merge$FOXP3 > 0 & annot.merge$IKZF2 > 0)] <- "CD4.IL22"
# Add in a Th17 phenotype
annot.merge$Sub.Annotation[annot.merge$Sub.Annotation %in% c("Treg", "CD4.CM", "CD4.Tfh", "CD4.EM", "CD4.Naive", "CD4.IL22") &
(annot.merge$CD8A <= 0 & annot.merge$CD8B <= 0) &
annot.merge$RORC > 0 &
(annot.merge$IL17A > 0 |
annot.merge$IL17F > 0 | annot.merge$IL21 > 0)] <- "CD4.Th17"
# Add in a Th1 phenotype
annot.merge$Sub.Annotation[annot.merge$Sub.Annotation %in% c("Treg", "CD4.CM", "CD4.EM", "CD4.Naive", "CD4.IL22") &
(annot.merge$CD8A <= 0 & annot.merge$CD8B <= 0) &
(annot.merge$IL2 > 0 | annot.merge$TNF > 0 |
annot.merge$LTA) &
annot.merge$TBX21 > 0] <- "CD4.Th1"
# Add in a Th2 phenotype
annot.merge$Sub.Annotation[annot.merge$Sub.Annotation %in% c("Treg", "CD4.CM", "CD4.EM", "CD4.Naive", "CD4.IL22") &
(annot.merge$CD8A <= 0 & annot.merge$CD8B <= 0) &
annot.merge$GATA3 > 0 &
(annot.merge$IL4 > 0 | annot.merge$IL5 > 0)] <- "CD4.Th2"
# merge the Tfh annotation
annot.merge$Sub.Annotation[annot.merge$Sub.Annotation %in% c("PD1-.Tfh-like", "PD1+.Tfh-like", "Tfh-like")] <- "CD4.Tfh"
# Add in the proliferating/activated/exhausted annotation
active.tcells <- read.csv("~/Dropbox/COVID19/Data/prolif_CD4CD8.csv",
header=FALSE, stringsAsFactors=FALSE)
annot.merge$Sub.Annotation[annot.merge$CellID %in% active.tcells$V1[active.tcells$V2 %in% c("prolif_CD4")]] <- "CD4.Prolif"
annot.merge$Sub.Annotation[annot.merge$CellID %in% active.tcells$V1[active.tcells$V2 %in% c("prolif_CD8")]] <- "CD8.Prolif"
# CD4 prolif
annot.merge$Sub.Annotation[annot.merge$Sub.Annotation %in% c("Proliferating") &
(annot.merge$CD8A <= 0 & annot.merge$CD8B <= 0)] <- "CD4.Prolif"
annot.merge$Sub.Annotation[annot.merge$Sub.Annotation %in% c("CD4.Prolif") &
!(annot.merge$TYMS > 0 | annot.merge$MKI67 > 0)] <- "CD4.CM"
# CD8A
annot.merge$Sub.Annotation[annot.merge$Sub.Annotation %in% c("Proliferating") &
(annot.merge$CD8A > 0 | annot.merge$CD8B > 0)] <- "CD8.Prolif"
annot.merge$Sub.Annotation[annot.merge$Sub.Annotation %in% c("CD8.Prolif") &
!(annot.merge$TYMS > 0 | annot.merge$MKI67 > 0)] <- "CD8.TE"
# table(annot.merge$Sub.Annotation, annot.merge$full_clustering)
table(annot.merge$Sub.Annotation, annot.merge$Site)
```
Compare to MNN-label transfers.
```{r, warning=FALSE, message=FALSE}
# only compare the MNN and new labels.
mnn.labels <- read.table("~/Dropbox/COVID19/Data/Updated/Tcell_MNN_labels.tsv",
sep="\t", header=TRUE, stringsAsFactors=FALSE)
mnn.compare <- merge(annot.merge, mnn.labels, by='CellID')
table(mnn.compare$Sub.Annotation.x, mnn.compare$Sub.Annotation.y)
```
```{r}
# write out the annotations, without the doublets, NK cells and ILCs
write.table(annot.merge[!annot.merge$Sub.Annotation %in% c("NK", "Prolif.NK", "ILCs", "B_naive", "CD14_mono",
"DC_prolif", "DC3", "HSC", "ILC1_3", "ILC2", "Mono_prolif",
"NK_16hi", "NK_56hi", "NK_prolif", "pDC", "Plasmablast",
"Platelets", "RBC",
"Doublets:Bcell", "Doublets:Platelet", "Doublets"),
c("CellID", "Sub.Annotation")],
file="~/Dropbox/COVID19/Data/Updated/Tcell_annotations_ext.tsv",
sep="\t", quote=FALSE, row.names=FALSE)
write.table(annot.merge[!annot.merge$Sub.Annotation %in% c("NK", "Prolif.NK", "ILCs", "B_naive", "CD14_mono",
"DC_prolif", "DC3", "HSC", "ILC1_3", "ILC2", "Mono_prolif",
"NK_16hi", "NK_56hi", "NK_prolif", "pDC", "Plasmablast",
"Platelets", "RBC",
"Doublets:Bcell", "Doublets:Platelet", "Doublets"),
c("CellID")],
file="~/Dropbox/COVID19/Data/Updated/Tcell_barcodes.tsv",
sep="\t", quote=FALSE, row.names=FALSE, col.names=FALSE)
```
Plot the proportions of T cells, normalised properly.
```{r}
cell.freq.tab <- table(annot.merge$sample_id,
annot.merge$Sub.Annotation)
cell.freq.df <- as.data.frame(sapply(rownames(cell.freq.tab), FUN=function(X) cell.freq.tab[X, ]/n.cell.vecc[X]))
cell.freq.df$CellType <- rownames(cell.freq.df)
cell.freq.df <- melt(cell.freq.df, id.vars=c("CellType"))
cell.freq.df$CellType <- as.character(cell.freq.df$CellType)
colnames(cell.freq.df) <- c("CellType", "sample_id", "Freq")
# cell.freq.df$Sex <- "Female"
# cell.freq.df$Sex[cell.freq.df$sample_id %in% unique(annot.merge$sample_id[annot.merge$Sex %in% c("M", "Male")])] <- "Male"
cell.freq.df$sample_id <- as.character(cell.freq.df$sample_id)
covid.meta <- read.csv("~/Dropbox/COVID19/Data/basic_COVID19_PBMC_metadata_160212.csv",
header=TRUE, stringsAsFactors=FALSE)
old.meta <- read.csv("~/Dropbox/COVID19/Data/Metadata FINAL 10122020.csv",
header=TRUE, stringsAsFactors=FALSE)
old.meta$sample_id[old.meta$sample_id %in% c("BGCV13_CV0201")] <- "BGCV06_CV0201"
old.meta$sample_id[old.meta$sample_id %in% c("BGCV06_CV0326")] <- "BGCV13_CV0326"
covid.meta$AgeRange <- covid.meta$Age
covid.meta$Status_on_day_collection_summary[covid.meta$Status_on_day_collection_summary %in% c("LPS_10hours")] <- "LPS"
covid.meta$Status_on_day_collection_summary[covid.meta$Status_on_day_collection_summary %in% c("LPS_90mins")] <- "LPS"
covid.meta <- merge(covid.meta[, !colnames(covid.meta) %in% c("Age")], old.meta[, c("sample_id", "patient_id", "Age")],
by='sample_id')
cell.freq.merge <- merge(cell.freq.df, covid.meta, by=c('sample_id'))
cell.freq.merge <- cell.freq.merge[!cell.freq.merge$Status_on_day_collection_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, fig.height=2.95, fig.width=9.15}
group.cols <- c("#2ca02c", "#1f77b4", "#9467bd", "#fed976", "#fd8d3c", "#e31a1c", "#800026", "#252525")
names(group.cols) <- c("Healthy", "LPS", "Non_covid", "Asymptomatic", "Mild", "Moderate", "Severe", "Critical")
cell.freq.merge$Status_on_day_collection_summary <- factor(cell.freq.merge$Status_on_day_collection_summary,
levels=c("Healthy", "Asymptomatic", "Mild",
"Moderate", "Severe", "Critical", "LPS", "Non_covid"))
```
```{r}
write.table(cell.freq.merge,
file="~/Dropbox/COVID19/Data/Updated/Tcell_proportions.tsv",
sep="\t", quote=FALSE, row.names=FALSE)
```
How do these compare to the previous frequencies?
```{r, fig.height=9.95, fig.width=10.95}
old.freqs <- read.table("~/Dropbox/COVID19/Data/Tcell_proportions.tsv",
sep="\t", header=TRUE, stringsAsFactors = FALSE)
old.freqs$sample_id[old.freqs$sample_id %in% c("BGCV13_CV0201")] <- "BGCV06_CV0201"
old.freqs$sample_id[old.freqs$sample_id %in% c("BGCV06_CV0326")] <- "BGCV13_CV0326"
old.freqs <- old.freqs[, c("sample_id", "CellType", "Freq")]
colnames(old.freqs) <- c("sample_id", "CellType", "old.Freq")
comp.freqs <- merge(cell.freq.merge, old.freqs, by=c('sample_id', 'CellType'))
comp.freqs$Switch <- comp.freqs$sample_id %in% c("BGCV06_CV0201", "BGCV13_CV0326")
ggplot(comp.freqs, aes(x=Freq, y=old.Freq, colour=Switch)) +
geom_point() +
theme_cowplot() +
scale_colour_manual(values=c("black", "red")) +
facet_wrap(~CellType, scales="free") +
NULL
```
```{r, warning=FALSE, message=FALSE, fig.height=5.95, fig.width=10.15}
ggplot(cell.freq.merge[cell.freq.merge$Collection_Day %in% c("D0") &
!cell.freq.merge$CellType %in% c("NK", "Prolif.NK", "ILCs", "B_naive", "CD14_mono",
"DC_prolif", "DC3", "HSC", "ILC1_3", "ILC2", "Mono_prolif",
"NK_16hi", "NK_56hi", "NK_prolif", "pDC", "Plasmablast",
"Platelets", "RBC",
"Doublets:Bcell", "Doublets:Platelet", "Doublets"), ],
aes(x=Status_on_day_collection_summary, y=Freq, fill=Status_on_day_collection_summary)) +
geom_boxplot() +
theme_cowplot() +
facet_wrap(~CellType, scales="free_y", ncol=4) +
scale_fill_manual(values=group.cols) +
expand_limits(y=c(0)) +
theme(axis.text.x=element_blank(),
strip.background=element_rect(fill='white', colour='white'),
strip.text=element_text(size=14)) +
labs(x="Severity Group", y="Proportion") +
guides(fill=guide_legend(title="Severity")) +
ggsave("~/Dropbox/COVID19/Updated_plots/Tcells_proportions-boxplot.pdf",
height=5.95, width=10.15, useDingbats=FALSE) +
NULL
```
For the figures these need to be a filled barchart (?!). Focus in on the Tfh-like cells.
```{r, warning=FALSE, message=FALSE, fig.height=2.15, fig.width=3.55}
ggplot(cell.freq.merge[cell.freq.merge$Collection_Day %in% c("D0") &
cell.freq.merge$CellType %in% c("CD4.Tfh"), ],
aes(x=Status_on_day_collection_summary, y=Freq, fill=Status_on_day_collection_summary)) +
geom_boxplot() +
theme_cowplot() +
facet_wrap(~CellType, scales="free_y") +
scale_fill_manual(values=group.cols) +
expand_limits(y=c(0)) +
theme(axis.text.x=element_blank(),
strip.background=element_rect(fill='white', colour='white'),
strip.text=element_text(size=14)) +
labs(x="Severity Group", y="Proportion") +
guides(fill=guide_legend(title="Severity")) +
ggsave("~/Dropbox/COVID19/Updated_plots/Tcells_proportions_Tfh-boxplot.pdf",
height=2.15, width=3.55, useDingbats=FALSE) +
NULL
```
Run a LM comparing categories.
```{r, warning=FALSE, message=FALSE, fig.height=2.15, fig.width=3.55}
ggplot(cell.freq.merge[cell.freq.merge$Collection_Day %in% c("D0") &
cell.freq.merge$CellType %in% c("CD4.Tfh"), ],
aes(x=Status_on_day_collection_summary, y=Freq, fill=Sex)) +
geom_boxplot() +
theme_cowplot() +
facet_wrap(~CellType, scales="free_y") +
scale_fill_colorblind() +
expand_limits(y=c(0)) +
theme(axis.text.x=element_blank(),
strip.background=element_rect(fill='white', colour='white'),
strip.text=element_text(size=14)) +
labs(x="Severity Group", y="Proportion") +
guides(fill=guide_legend(title="Severity")) +
ggsave("~/Dropbox/COVID19/Updated_plots/Tcells_proportions_Tfh-byGender_boxplot.pdf",
height=2.15, width=3.55, useDingbats=FALSE) +
NULL
```
Compare healthy vs. asymptomatic, and asymptomatic > critical.
```{r, warning=FALSE, message=FALSE}
healthy.tfh <- cell.freq.merge$Freq[cell.freq.merge$CellType %in% c("CD4.Tfh") &
cell.freq.merge$Status_on_day_collection_summary %in% c("Healthy") &
cell.freq.merge$Collection_Day %in% c("D0")]
asymp.tfh <- cell.freq.merge$Freq[cell.freq.merge$CellType %in% c("CD4.Tfh") &
cell.freq.merge$Status_on_day_collection_summary %in% c("Asymptomatic") &
cell.freq.merge$Collection_Day %in% c("D0")]
tfh.utest <- wilcox.test(x=healthy.tfh, y=asymp.tfh, alternative="two.sided")
tfh.utest
```
I should model the counts, normalized by the total numbers of captured cells.
```{r}
tfh.lm.list <- list()
covid.tfh.merge <- cell.freq.merge[!cell.freq.merge$Status_on_day_collection_summary %in% c("Healthy"), ]
covid.tfh.merge$OrderedStatus <- ordered(covid.tfh.merge$Status_on_day_collection_summary,
levels=c("Asymptomatic", "Mild", "Moderate", "Severe", "Critical"))
covid.tfh.merge$Days_from_onset <- as.numeric(covid.tfh.merge$Days_from_onset)
c.types <- c("CD4.Tfh")
for(i in seq_along(c.types)){
i.c <- c.types[i]
i.lm.df <- data.frame("sample_id"=rownames(cell.freq.tab), "Count"=cell.freq.tab[, i.c])
i.lm.df <- merge(i.lm.df, covid.tfh.merge[covid.tfh.merge$CellType %in% i.c, ], by='sample_id')
# i.tfh.lm <- glm(Count ~ Sex + Age + Site + OrderedStatus,
# data=i.lm.df, family=poisson(link="log"))
i.tfh.lm <- glm.nb(Count ~ Sex + Age + Site + OrderedStatus, data=i.lm.df, init.theta=1)
i.lm.res <- summary(i.tfh.lm)
tfh.lm.list[[i.c]] <- data.frame("CellType"=rep(i.c, nrow(i.lm.res$coefficients)),
"Beta"=i.lm.res$coefficients[, 1],
"SE"=i.lm.res$coefficients[, 2],
"STAT"=i.lm.res$coefficients[, 3],
"Pvalue"=i.lm.res$coefficients[, 4],
"Predictor"=rownames(i.lm.res$coefficients))
}
tfh.lm.df <- do.call(rbind.data.frame, tfh.lm.list)
tfh.lm.df$Predictor <- gsub(tfh.lm.df$Predictor, pattern="OrderedStatus", replacement="Severity")
tfh.lm.df$Predictor <- gsub(tfh.lm.df$Predictor, pattern="Male", replacement="")
tfh.lm.df$Predictor <- gsub(tfh.lm.df$Predictor, pattern="Site", replacement="")
tfh.lm.df <- tfh.lm.df[!tfh.lm.df$Predictor %in% c("(Intercept)"), ]
tfh.lm.df$P.Adjust <- p.adjust(tfh.lm.df$Pvalue)
```
```{r, warning=FALSE, message=FALSE, fig.height=3.95, fig.width=4.95}
ggplot(tfh.lm.df, aes(x=Predictor, y=Beta, fill=Pvalue < 0.05)) +
geom_hline(yintercept=0, lty=2, col='grey70') +
geom_errorbar(aes(ymin=Beta-SE, ymax=Beta+SE)) +
geom_point(shape=22, size=3) +
theme_cowplot() +
theme(strip.background=element_rect(fill='white', colour='white'),
strip.text=element_text(size=14)) +
coord_flip() +
scale_fill_manual(values=c("blue", "orange")) +
facet_wrap(~CellType, scales="free_x") +
ggsave("~/Dropbox/COVID19/Updated_plots/Tcell_Tfh_props-LM.pdf",
height=3.95, width=4.95, useDingbats=FALSE)
```
```{r, warning=FALSE, message=FALSE}
cellfreq.df <- cell.freq.merge %>% group_by(Status_on_day_collection_summary, CellType) %>% summarise(Mean = mean(Freq))
cellfreq.df$Status_on_day_collection_summary <- factor(cellfreq.df$Status_on_day_collection_summary,
levels=c("Healthy", "Asymptomatic", "Mild",
"Moderate", "Severe", "Critical", "LPS"))
```
Create a new colour palette.
```{r}
cell.order <- c("CD4.Naive", "CD4.CM", "CD4.EM", "CD4.IL22", "CD4.Prolif", "CD4.Th1", "CD4.Th2", "CD4.Th17", "CD4.Tfh",
"Treg", "CD8.Naive", "CD8.Activated", "CD8.Prolif", "CD8.CM", "CD8.TE", "CD8.EM", "gdT", "MAIT", "NKT")
all.qual.cols <- brewer.pal.info[brewer.pal.info$category %in% c("qual") & brewer.pal.info$colorblind, ]
col_vector = unlist(mapply(brewer.pal, all.qual.cols$maxcolors, rownames(all.qual.cols)))
# cell.cols <- col_vector[c(1:7, 9:19)]
cell.cols <- col_vector[c(1:17, 19, 20)]
names(cell.cols) <- cell.order
```
```{r}
pie(rep(1, 18), col=cell.cols)
```
```{r, warning=FALSE, message=FALSE, fig.height=4.25, fig.width=6.95}
# n.cells <- length(cell.order)
# cell.cols <- colorRampPalette(pal_futurama()(12))(n.cells)
# names(cell.cols) <- cell.order
cellfreq.df$CellType <- factor(cellfreq.df$CellType,
levels=cell.order)
ggplot(cellfreq.df,
aes(x=Status_on_day_collection_summary, y = Mean, fill=CellType))+
geom_bar(position="fill", stat="identity", width = 0.9, colour = "black")+
scale_fill_manual(values = cell.cols) +
theme(aspect.ratio = 1/1.5,
axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1, size = 14, colour = "black"),
panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
axis.text.y = element_text(size = 14, colour = "black"),
legend.key.size=unit(0.4, "cm"),
axis.title.x=element_blank(), axis.title.y=element_blank(),
legend.text = element_text(size = 14), legend.title=element_text(size=14),
panel.background = element_blank(), axis.line = element_line(colour = "black")) +
guides(fill=guide_legend(title="Cell type", ncol=1)) +
ggsave("~/Dropbox/COVID19/Updated_plots/Tcell_barchart.pdf",
height=4.25, width=6.95, useDingbats=FALSE)
```
I have a misgiving about this, as it masks the distribution of the proportions - maybe a series of filled boxplots, e.g. a waterfall?
```{r, warning=FALSE, message=FALSE, fig.height=5.25, fig.width=12.95}
cell.freq.merge$CellType <- factor(cell.freq.merge$CellType,
levels=cell.order)
ggplot(cell.freq.merge,
aes(x=sample_id, y = Freq, fill=CellType))+
geom_bar(position="fill", stat="identity", width = 0.8, colour = "black")+
scale_fill_manual(values = cell.cols) +
facet_grid(. ~ Status_on_day_collection_summary, scales="free_x", space="free") +
theme(axis.text.x = element_blank(),
panel.spacing=unit(0, "lines"),
panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
axis.text.y = element_text(size = 18),
axis.title.x=element_blank(), axis.title.y=element_blank(),
axis.ticks.x=element_blank(),
strip.background=element_rect(fill='white', colour='white'),
strip.text=element_text(size=14, angle=90, colour='black'),
legend.text = element_text(size = 18), legend.title=element_text(size=18),
panel.background = element_blank(), axis.line = element_line(colour = "black")) +
# guides(fill=guide_legend(title="Cell type")) +
ggsave("~/Dropbox/COVID19/Updated_plots/Tcell_barchart-full.pdf",
height=5.25, width=8.95, useDingbats=FALSE) +
NULL
```
Make dot plots.
```{r, warning=FALSE, message=FALSE}
clust.model.matrix <- model.matrix(~ 0 + Sub.Annotation, data=annot.merge)
rownames(clust.model.matrix) <- annot.merge$CellID
gene.goi.by.cluster <- t(clust.model.matrix) %*% as.matrix(annot.merge[, setdiff(colnames(goi.df), colnames(tcell.merge))])
gene.ave.by.label <- as.data.frame(apply(gene.goi.by.cluster, 2, function(X) X/colSums(clust.model.matrix)))
rownames(gene.ave.by.label) <- gsub(rownames(gene.ave.by.label), pattern="Sub\\.Annotation", replacement="")
gene.goi.by.cluster.bin <- t(clust.model.matrix) %*% as.matrix(annot.merge[, setdiff(colnames(goi.df), colnames(tcell.merge))] > 0)
gene.bin.by.label <- as.data.frame(apply(gene.goi.by.cluster.bin, 2, function(X) X/sum(X)))
gene.bin.by.label[is.na(gene.bin.by.label)] <- 0
rownames(gene.bin.by.label) <- gsub(rownames(gene.bin.by.label), pattern="Sub\\.Annotation", replacement="")
```
```{r, message=FALSE, warning=FALSE}
gene.ave.label.df <- as.data.frame(apply(gene.ave.by.label, 2, FUN=function(XP) XP/max(XP)))
gene.ave.label.df$Sub.Annotation <- rownames(gene.ave.label.df)
gene.ave.label.melt <- melt(gene.ave.label.df, id.vars='Sub.Annotation')
colnames(gene.ave.label.melt) <- c("Sub.Annotation", "Gene", "AveExprs")
gene.bin.label.df <- gene.bin.by.label
gene.bin.label.df$Sub.Annotation <- rownames(gene.bin.label.df)
gene.bin.label.melt <- melt(gene.bin.label.df, id.vars='Sub.Annotation')
colnames(gene.bin.label.melt) <- c("Sub.Annotation", "Gene", "PropExprs")
```
```{r, message=FALSE, warning=FALSE, fig.height=4.65, fig.width=8.95}
dot.plot.df <- merge(gene.ave.label.melt, gene.bin.label.melt, by=c('Sub.Annotation', 'Gene'))
dot.plot.df$Gene <- as.character(dot.plot.df$Gene)
dot.plot.genes <- c("CD3G", "CD8A", "CD4", "CCR7", "CD27", "CD28", "TOX", "CTLA4", "IFNG", "TBX21",
"PDCD1", "GZMB", "RORC", "IKZF2", "GATA3", "TNFRSF9", "IL7R", "KLRG1",
"TRAV1.2", "TRDV2", "TRGV9", "CD40LG", "RORC", "NCAM1", "CCR6",
"CXCR5", "HLA.DRB1", "LAG3", "FOXP3", "MKI67", "NCR1", "GARP", "FCGR3A")
gene.cluster <- hclust(dist(t(gene.ave.by.label[, intersect(colnames(gene.ave.by.label), dot.plot.genes)])))$order
# manually set the gene order
gene.order <- c("CD3G", "CCR7", "CD27", "CD4", "CD28", "TNFRSF9", "CCR6", "CD40LG", "MKI67", "TBX21", "GATA3", "RORC", "IKZF2",
"FOXP3", "CTLA4", "PDCD1", "CXCR5", "CD8A", "GZMB", "HLA.DRB1", "IFNG", "LAG3", "TOX", "IL7R", "KLRG1", "TRDV2",
"TRGV9", "TRAV1.2", "NCAM1", "NCR1", "FCGR3A")
dot.plot.df <- dot.plot.df[dot.plot.df$Gene %in% dot.plot.genes, ]
dot.plot.df$Gene <- factor(dot.plot.df$Gene,
levels=gene.order)
dot.plot.df <- dot.plot.df[dot.plot.df$Sub.Annotation %in% cell.order, ]
dot.plot.df$Sub.Annotation <- factor(dot.plot.df$Sub.Annotation,
levels=rev(cell.order))
ggplot(dot.plot.df,
aes(x=Gene, y=Sub.Annotation, colour=AveExprs, size=PropExprs*100)) +
geom_point() +
theme_cowplot() +
scale_colour_gradient(low='grey70', high='blue') +
theme(axis.text.x=element_text(angle=90, vjust=0.5, hjust=1, colour='black', size=14),
axis.text.y=element_text(colour='black', size=14),
legend.text=element_text(size=14),
legend.title=element_text(size=14),
panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
panel.background = element_blank(), axis.line = element_line(colour = "black")) +
guides(colour=guide_colorbar(title='Mean\nExpression'),
size=guide_legend(title='% Express')) +
labs(x="Gene", y="Cell type") +
ggsave("~/Dropbox/COVID19/Updated_plots/Tcell_ext_dotplot-genes.pdf",
height=4.65, width=8.95, useDingbats=FALSE)
```
Plot the cytokine profiles to figure out if Th1/2/17 dominate (or all of them).
```{r, message=FALSE, warning=FALSE, fig.height=4.15, fig.width=7.0}
cytokine.dot.plot.df <- merge(gene.ave.label.melt, gene.bin.label.melt, by=c('Sub.Annotation', 'Gene'))
cytokine.dot.plot.df$Gene <- as.character(cytokine.dot.plot.df$Gene)
cytokine.genes <- intersect(c("IL2", "IL1A", "IL1B", "IL10", "IL22", "TNF", "IFNG", "LTA",
"IL4", "IL5", "IL13", "IL17A", "IL17F", "IL21",
"IL12A", "IL12B", "IL25"), colnames(annot.merge))
cytokine.dot.plot.df <- cytokine.dot.plot.df[cytokine.dot.plot.df$Gene %in% cytokine.genes, ]
cytokine.dot.plot.df$Gene <- factor(cytokine.dot.plot.df$Gene,
levels=cytokine.genes)
cytokine.dot.plot.df <- cytokine.dot.plot.df[cytokine.dot.plot.df$Sub.Annotation %in% cell.order, ]
cytokine.dot.plot.df$Sub.Annotation <- factor(cytokine.dot.plot.df$Sub.Annotation,
levels=rev(cell.order))
ggplot(cytokine.dot.plot.df,
aes(x=Gene, y=Sub.Annotation, colour=AveExprs, size=PropExprs*100)) +
geom_point() +
theme_cowplot() +
scale_colour_gradient(low='grey70', high='blue') +
theme(axis.text.x=element_text(angle=90, vjust=0.5, hjust=1, colour='black', size=14),
axis.text.y=element_text(colour='black', size=14),
legend.text=element_text(size=14),
legend.title=element_text(size=14),
legend.key.size=unit(0.4, "cm"),
panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
panel.background = element_blank(), axis.line = element_line(colour = "black")) +
guides(colour=guide_colorbar(title='Mean\nExpression'),
size=guide_legend(title='% Express')) +
labs(x="Cytokine", y="Cell type") +
ggsave("~/Dropbox/COVID19/Updated_plots/Tcell_ext_dotplot-cytokines.pdf",
height=4.65, width=7.0, useDingbats=FALSE) +
NULL
```
I want to see how these cytokines change over the course of severity.
```{r, warning=FALSE, message=FALSE}
cytokine.melt <- melt(annot.merge[, c("sample_id", "Status_on_day_collection_summary", "Sub.Annotation", "Collection_Day", "Sex", "Age", "CellID",
cytokine.genes)],
variable.name="Cytokine",
id.vars=c("sample_id", "Status_on_day_collection_summary", "Sub.Annotation", "Collection_Day", "Sex", "Age", "CellID"))
cytokine.means <- cytokine.melt %>% group_by(Status_on_day_collection_summary, Sub.Annotation, Cytokine) %>% summarise("Mean.Exprs"=mean(value))
cytokine.var <- cytokine.melt %>% group_by(Status_on_day_collection_summary, Sub.Annotation, Cytokine) %>% summarise("Var.Exprs"=var(value))
cytokine.se <- cytokine.melt %>% group_by(Status_on_day_collection_summary, Sub.Annotation, Cytokine) %>% summarise("SE.Exprs"=sqrt(var(value))/sqrt(sum(value >= 0)))
cytokine.df <- Reduce(x=list("mean"=cytokine.means, "var"=cytokine.var, "se"=cytokine.se),
f=function(x, y) merge(x, y, by=c('Status_on_day_collection_summary', 'Sub.Annotation', 'Cytokine')))
# cytokine.df <- merge(cytokine.df, covid.meta, by='sample_id')
# remove superflous or malignant samples
# cytokine.df <- cytokine.df[!cytokine.df$patient_id %in% c("CV0198"), ]
# cytokine.df <- cytokine.df[!cytokine.df$sample_id %in% c("BGCV01_CV0902"), ]
# cytokine.df <- cytokine.df[!cytokine.df$Status_on_day_collection_summary %in% c("Non_Covid"), ]
cytokine.df$D0_satus_summary <- factor(cytokine.df$Status_on_day_collection_summary,
levels=c("Healthy", "Asymptomatic", "Mild",
"Moderate", "Severe", "Critical", "LPS"))
```
Maybe visualise this as a heatmap?
```{r, fig.height=7.15, fig.width=7.95}
cytokine.df <- cytokine.df[!cytokine.df$Sub.Annotation %in% c("NK", "ILCs", "Doublets:Bcell", "Doublets", "Doublets:Platelet") &
!is.na(cytokine.df$D0_satus_summary), ]
ggplot(cytokine.df[!cytokine.df$Cytokine %in% c("IL13", "IL22", "IL12A", "IL17A", "IL17F", "IL10", "IL6"), ],
aes(x=D0_satus_summary, y=Cytokine, fill=Mean.Exprs)) +
geom_tile() +
facet_wrap(~Sub.Annotation) +
scale_fill_viridis(option="magma") +
theme_cowplot() +
guides(fill=guide_colourbar(title="Mean Expression")) +
labs(y="Cytokine Gene", x="Severity") +
theme(axis.text.x=element_text(angle=90, vjust=0.5, hjust=1),
strip.background=element_rect(fill='white', colour='white')) +
ggsave("~/Dropbox/COVID19/Updated_plots/Tcell_cytokines_heatmap.pdf",
height=7.15, width=7.95, useDingbats=FALSE) +
NULL
```
What about exhaustion markers?
```{r, warning=FALSE, message=FALSE}
exhaustion.genes <- c("TOX", "TIGIT", "LAG3", "HAVCR2", "PDCD1", "FOXP3", "MKI67")
exhaustion.melt <- melt(annot.merge[, c("sample_id", "Status_on_day_collection_summary", "Sub.Annotation", "Collection_Day", "Sex", "Age", "CellID",
exhaustion.genes)],
variable.name="Gene",
id.vars=c("sample_id", "Status_on_day_collection_summary", "Sub.Annotation", "Collection_Day", "Sex", "Age", "CellID"))
exhaustion.means <- exhaustion.melt %>% group_by(Status_on_day_collection_summary, Sub.Annotation, Gene) %>% summarise("Mean.Exprs"=mean(value))
exhaustion.var <- exhaustion.melt %>% group_by(Status_on_day_collection_summary, Sub.Annotation, Gene) %>% summarise("Var.Exprs"=var(value))
exhaustion.se <- exhaustion.melt %>% group_by(Status_on_day_collection_summary, Sub.Annotation, Gene) %>% summarise("SE.Exprs"=sqrt(var(value))/sqrt(sum(value >= 0)))
exhaustion.df <- Reduce(x=list("mean"=exhaustion.means, "var"=exhaustion.var, "se"=exhaustion.se),
f=function(x, y) merge(x, y, by=c('Status_on_day_collection_summary', 'Sub.Annotation', 'Gene')))
# exhaustion.df <- merge(exhaustion.df, covid.meta, by='sample_id')
# remove superflous or malignant samples
# exhaustion.df <- exhaustion.df[!exhaustion.df$patient_id %in% c("CV0198"), ]
# exhaustion.df <- exhaustion.df[!exhaustion.df$sample_id %in% c("BGCV01_CV0902"), ]
# exhaustion.df <- exhaustion.df[!exhaustion.df$Status_on_day_collection_summary %in% c("Non_Covid"), ]
exhaustion.df$D0_satus_summary <- factor(exhaustion.df$Status_on_day_collection_summary,
levels=c("Healthy", "Asymptomatic", "Mild",
"Moderate", "Severe", "Critical", "LPS"))
```
```{r, fig.height=5.45, fig.width=9.95}
exhaustion.df <- exhaustion.df[!exhaustion.df$Sub.Annotation %in% c("NK", "Prolif.NK", "ILCs", "B_naive", "CD14_mono",
"DC_prolif", "DC3", "HSC", "ILC1_3", "ILC2", "Mono_prolif",
"NK_16hi", "NK_56hi", "NK_prolif", "pDC", "Plasmablast",
"Platelets", "RBC",
"Doublets:Bcell", "Doublets:Platelet", "Doublets") &
!is.na(exhaustion.df$D0_satus_summary), ]
ggplot(exhaustion.df[!exhaustion.df$Gene %in% c("IL13", "IL22", "IL12A", "IL17A", "IL17F"), ],
aes(x=D0_satus_summary, y=Gene, fill=Mean.Exprs)) +
geom_tile() +
facet_wrap(~Sub.Annotation, nrow=3) +
scale_fill_viridis(option="magma") +
theme_cowplot() +
guides(fill=guide_colourbar(title="Mean Expression")) +
labs(x="Severity", "Gene") +
theme(axis.text.x=element_text(angle=90, vjust=0.5, hjust=1),
strip.background=element_rect(fill='white', colour='white')) +
ggsave("~/Dropbox/COVID19/Updated_plots/Tcell_exhaustion_heatmap.pdf",
height=5.45, width=9.95, useDingbats=FALSE) +
NULL
```
Do the same for the proteins.
```{r, warning=FALSE, message=FALSE}
clust.model.matrix <- model.matrix(~ 0 + Sub.Annotation, data=annot.merge)
rownames(clust.model.matrix) <- annot.merge$CellID
goi.by.cluster <- t(clust.model.matrix) %*% as.matrix(annot.merge[, gsub(setdiff(colnames(prot.df), colnames(tcell.merge)),
pattern="\\.x", replacement="")])
ave.by.label <- as.data.frame(apply(goi.by.cluster, 2, function(X) X/colSums(clust.model.matrix)))
rownames(ave.by.label) <- gsub(rownames(ave.by.label), pattern="Sub\\.Annotation", replacement="")
colnames(ave.by.label) <- gsub(colnames(ave.by.label), pattern="Prot\\.", replacement="")
goi.by.cluster.bin <- t(clust.model.matrix) %*% as.matrix(annot.merge[, gsub(setdiff(colnames(prot.df), colnames(tcell.merge)),
pattern="\\.x", replacement="")] > 0.01)
bin.by.label <- as.data.frame(apply(goi.by.cluster.bin, 2, function(X) X/sum(X)))
bin.by.label[is.na(bin.by.label)] <- 0
rownames(bin.by.label) <- gsub(rownames(bin.by.label), pattern="Sub\\.Annotation", replacement="")
colnames(bin.by.label) <- gsub(colnames(bin.by.label), pattern="Prot\\.", replacement="")
```
```{r, message=FALSE, warning=FALSE}
ave.label.df <- as.data.frame(apply(ave.by.label, 2, FUN=function(XP) (XP-min(XP))/max(XP-min(XP))))
ave.label.df$Sub.Annotation <- rownames(ave.label.df)
ave.label.melt <- melt(ave.label.df, id.vars='Sub.Annotation')
colnames(ave.label.melt) <- c("Sub.Annotation", "Protein", "AveExprs")
bin.label.df <- bin.by.label
bin.label.df$Sub.Annotation <- rownames(bin.label.df)
bin.label.melt <- melt(bin.label.df, id.vars='Sub.Annotation')
colnames(bin.label.melt) <- c("Sub.Annotation", "Protein", "PropExprs")
```
```{r, message=FALSE, warning=FALSE, fig.height=5.15, fig.width=8.95}
protein.dot.plot.df <- merge(ave.label.melt, bin.label.melt, by=c('Sub.Annotation', 'Protein'))
protein.dot.plot.df$Protein <- as.character(protein.dot.plot.df$Protein)
dot.plot.prots <- c("CD3", "CD8", "CD4", "CCR7", "CD45RA", "CD45RO", "CD28", "CD27", "CTLA4", "CD56", "CD16",
"CD274", "CXCR3", "CD40LG", "ICOS", "PD1", "CD25", "LAMP1", "CXCR3", "CCR6", "IL7R",
"CXCR5", "TCR_Va7.2", "TCR_Vg2", "TCR_Vg9", "TCR_Va24-Ja18", "CD62L", "TCR_Va24.Ja18")
protein.order <- c("CD3", "CD4", "CCR7", "CD45RA", "CD45RO", "CD27", "CD28", "CD62L", "CD25", "CTLA4",
"LAMP1", "CCR6", "IL7R",
"CXCR5", "CD40LG", "ICOS", "CXCR3", "CD8", "PD1", "CD274", "TCR_Vg9", "TCR_Vg2",
"TCR_Va7.2", "TCR_Va24.Ja18", "CD56", "CD16")
protein.dot.plot.df <- protein.dot.plot.df[protein.dot.plot.df$Protein %in% dot.plot.prots, ]
protein.dot.plot.df$Protein <- factor(protein.dot.plot.df$Protein,
levels=protein.order)
protein.dot.plot.df <- protein.dot.plot.df[protein.dot.plot.df$Sub.Annotation %in% cell.order, ]
protein.dot.plot.df$Sub.Annotation <- factor(protein.dot.plot.df$Sub.Annotation,
levels=rev(cell.order))
ggplot(protein.dot.plot.df,
aes(x=Protein, y=Sub.Annotation, colour=AveExprs, size=PropExprs*100)) +
geom_point() +
theme_cowplot() +
scale_colour_gradient(low='grey70', high=darken('red')) +
theme(axis.text.x=element_text(angle=90, vjust=0.5, hjust=1, colour='black', size=14),
axis.text.y=element_text(colour='black', size=14),
legend.text=element_text(size=14),
legend.title=element_text(size=14),
panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
panel.background = element_blank(), axis.line = element_line(colour = "black")) +
guides(colour=guide_colorbar(title='Mean\nExpression'),
size=guide_legend(title='% Express')) +
labs(x="Protein", y="Cell type") +
ggsave("~/Dropbox/COVID19/Updated_plots/Tcell_ext_dotplot-protein.pdf",
height=5.15, width=8.95, useDingbats=FALSE)
```
I also want to see some of the CD8 sub-phenotype mRNA and protein expression, namely _KLRG1_ and _CD217_ (_IL7R_).
```{r, warning=FALSE, message=FALSE, fig.height=5.95, fig.width=11.95}
repeat.samps <- covid.meta$patient_id[duplicated(covid.meta$patient_id)]
cell.freq.merge$Day <- ordered(cell.freq.merge$Collection_Day,
levels=c("D0", "D7", "D9", "D12", "D13", "D28"))
ggplot(cell.freq.merge[cell.freq.merge$patient_id %in% repeat.samps &
!cell.freq.merge$CellType %in% c("NK", "Prolif.NK", "ILCs", "Doublets:Bcell", "Doublets:Platelet", "Doublets") &
!cell.freq.merge$Status_on_day_collection_summary %in% c("Non_covid", "LPS"), ],
aes(x=Day, y=Freq, fill=Status_on_day_collection_summary, group=patient_id)) +
#geom_boxplot() +
geom_line(aes(colour=Status_on_day_collection_summary)) +
geom_point(shape=21, size=3) +
theme_cowplot() +
facet_wrap(~CellType, scales="free_y") +
scale_fill_manual(values=group.cols) +
scale_colour_manual(values=group.cols) +
expand_limits(y=c(0)) +
theme(strip.background=element_rect(fill='white', colour='white'),
strip.text=element_text(size=14)) +
labs(x="Day", y="Proportion") +
guides(fill=guide_legend(title="Severity"),
colour=FALSE) +
# ggsave("~/Dropbox/COVID19/Updated_plots/Tcells_proportions-boxplot.pdf",
# height=5.95, width=11.15, useDingbats=FALSE) +
NULL
```
Datasets
Models
