---

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