---

title: "CD3+ cells"

author: "Ming Tang"

date: '2023-07-26'

output: html_document

---

```{r}

library(devtools)

library(Seurat)

library(tidyverse)

library(presto)

library(here)

library(tictoc)

library(ComplexHeatmap)

# presto is much faster than Seurat::FindAllMarkers

# install_github('immunogenomics/presto')

```

Read in the Seurat object

```{r}

cd3_plus_annotated<- readRDS(here("data/cd3_plus_annotated_seurat_subcluster.rds"))

```

### whole CD3+ subset

```{r}

annotation_colors<-c("#3366cc", "#7d5378", "#c83e27", "#e7580c", "#f68104", "#cc9803", "#649710", "#1d8723", "#58455b", "#950395", "#5740ab", "#1882bf", "#427fae", "#a4598a", "#c55660", "#92842c", "#6ba103", "#8e6a17", "#855499", "#814256", "#465a84", "#4e5a96", "#7b4c97", "#7b4c97", "#29aa92", "#b4332b", "#a0aa1b", "#907d56")

# Figure 1a

DimPlot(cd3_plus_annotated, group.by = "clusters_anno_sub", cols = annotation_colors)

# supp Figure 2

set.seed(123)

library(ComplexHeatmap)

clusters<- cd3_plus_annotated@meta.data$clusters_anno_sub%>%unique()

markers<- presto::wilcoxauc(cd3_plus_annotated, 'clusters_anno_sub', assay = 'data',

                               groups_use = clusters)

markers %>% 

  arrange(desc(auc)) %>% 

  filter(pct_in> 20, pct_out < 30) %>%

  head()

top_markers_unique<- filter(markers, padj <=0.05, abs(logFC) >=0.8) %>% 

  pull(feature) %>%

  unique()

## for plotting heatmap, subsample the seurat object, 100000 or 5000 cells

cd3_plus_annotated_sub<- cd3_plus_annotated[, sample.int(ncol(cd3_plus_annotated), size = 10000)]

all_cells<- cd3_plus_annotated_sub@meta.data$clusters_anno_sub %in% clusters

all_mat<- cd3_plus_annotated_sub[["RNA"]]@data[top_markers_unique, all_cells ] %>% as.matrix()

## scale the rows

all_mat<- t(scale(t(all_mat)))

cluster_anno<- cd3_plus_annotated_sub@meta.data$clusters_anno_sub[all_cells]

column_ha2<- HeatmapAnnotation(

    cluster_anno = cluster_anno,

    col = list(cluster_anno = setNames(annotation_colors, sort(clusters))),

    na_col = "grey",

    show_annotation_name = TRUE

)

quantile(all_mat, c(0.1, 0.5, 0.95, 0.99))

PurpleAndYellow()

#col_fun = colorRamp2(c(-1, 0, 2), c("blue", "white", "red"))

#col_fun = circlize::colorRamp2(c(-1, 0, 2), c("#FF00FF", "black", "#FFFF00"))

col_fun<-  circlize::colorRamp2(c(-1, 0, 2), c("#440154FF", "#238A8DFF", "#FDE725FF"))

## for big matrix, raster helps. https://jokergoo.github.io/2020/06/30/rasterization-in-complexheatmap/

## cluster the slice (without specifiying the factor level) and the columns.

# can not directly rotate column titles (cluster annotation) for 45, has to use decorate() function, see

# https://github.com/jokergoo/ComplexHeatmap/issues/189

pdf(here("results/figures/2022-04-14_cd3_plus_all_clusters_marker_gene_heatmap_remake_new_color.pdf"), width = 16, height = 18)

Heatmap(all_mat, name = "Expression",  

        column_split = factor(cluster_anno),

        cluster_columns = TRUE,

        show_column_dend = FALSE,

        cluster_column_slices = TRUE,

        row_dend_reorder = TRUE,

        column_title_gp = gpar(fontsize = 6),

        column_title_rot = 90,

        column_gap = unit(0.5, "mm"),

        cluster_rows = TRUE,

        column_dend_reorder = TRUE,

        show_row_dend = FALSE,

        col = col_fun,

        row_km = 20,

        row_gap = unit(0.5, "mm"),

        row_title = NULL,

        row_names_gp = gpar(fontsize = 6),

        top_annotation = column_ha2,

        show_column_names = FALSE,

        use_raster = TRUE,

        raster_by_magick = TRUE,

        raster_quality = 10)

dev.off()

StackedVlnPlot(cd3_plus_annotated, features = top_markers_unique, idents = c(0:25, 27,28))

ggsave(here("results/figures/cd3_plus_clusters_marker_gene_stacked_violin.pdf"), width = 6, height = 10)

```

### clustered dotplot 

```{r}

source(here("src/cluster_dotplot.R"))

cd3_genes<- c("IFNG", "IL12A", "LAG3", "CD274", "IL4", "PDCD1", "CCR5", "IL32", "IL10", "HAVCR2", "IRF4", "GZMB", "TBX21", "PRF1", "TGFB1", "KLRG1", "TRGC1", "TRGC2", "TRDC", "B3GAT1", "KLRB1", "RORC", "BCL6", "TNF", "GATA3", "CD44", "CCR4", "IL2RA", "FOXP3", "CD28", "CTLA4", "CD27", "SELL", "CCR7", "CD4", "CD8A", "CD8B", "TRBV30", "TRAV8-2", "TRAV13-1", "TRBV21-1", "TRBV24-1", "TRDV2", "TRGV9", "TRBV11-3", "TRAV16", "TRAV12-1", "TRAV1-2", "ISG15", "IFIT2", "IFIT3", "IFI44L","PMAIP1", "OASL", "ZC3HAV1", "DUSP2", "MAP3K8", "SRGN", "CCL3", "NCR3", "CD247", "FCGR3A", "KLRF1", "FCER1G", "TYROBP", "HCST", "CST7", "NKG7", "CCL5", "GZMA", "CTSW", "GZMM", "HOPX", "CCL4", "KLRD1", "GNLY", "FGFBP2", "GZMH", "NFKBIA", "JUNB", "TNFAIP3", "FOS", "DUSP1", "DDIT4", "FTH1", "ZFP36L2", "LINC02446", "PELI1", "SOX4", "GZMK", "HLA-DMA", "HLA-DQB1", "HLA-DPA1", "HLA-DPB1", "HLA-DRA", "HLA-DRB1", "HLA-DQA1", "CD74", "TXN", "COTL1", "CRIP1", "CD52", "S100A4", "B2M", "PCNA","PCLAF", "TYMS", "STMN1", "DUT", "DEK", "TUBA1B", "TUBB", "HMGN1", "H2AFY", "ANP32B", "NASP", "VIM", "TRDV1", "TRDV3")

cd3_genes[!cd3_genes %in% rownames(cd3_plus_annotated)]

clusters<- cd3_plus_annotated$clusters_anno%>%unique()

cd3_cells<- cd3_plus_annotated$clusters_anno %in% clusters

cd3_plus_annotation<- 

  select(cd3_plus_annotated@meta.data, seurat_clusters, annotation) %>%

  distinct() %>%

  mutate(annotation = str_replace_all(annotation, "-", "_")) %>%

  mutate(annotation = str_replace_all(annotation, "\\+", "")) 

library(Polychrome)

annotation_colors<- kelly.colors(length(unique(cd3_plus_annotation$annotation))) %>%

  unname()

pdf(here("results/figures/cd3_plus_all_dotplot.pdf"), width = 10, height = 24)

cluster_dot_plot(cd3_plus_annotated[, cd3_cells], features = unique(cd3_genes), annotation = cd3_plus_annotation, annotation_colors = annotation_colors, color_scale_limits = c(-1.5, 2.5))

dev.off()

pdf(here("results/figures/cd3_plus_healthy_only_all_dotplot.pdf"), width = 10, height = 24)

cd3_cells<- cd3_plus_annotated_health$clusters_anno %in% clusters

cluster_dot_plot(cd3_plus_annotated_health[, cd3_cells], features = unique(cd3_genes), annotation = cd3_plus_annotation, annotation_colors = annotation_colors, color_scale_limits = c(-1.5, 2.5))

dev.off()

pdf(here("results/figures/cd3_plus_no_health_all_dotplot.pdf"), width = 10, height = 24)

cd3_cells<- cd3_plus_annotated_no_health$clusters_anno %in% clusters

cluster_dot_plot(cd3_plus_annotated_no_health[, cd3_cells], features = unique(cd3_genes), annotation = cd3_plus_annotation, annotation_colors = annotation_colors, color_scale_limits = c(-1.5, 2.5))

dev.off()

```