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

