APC.Rmd

---
title: Integration of multiple experiments for the ccRCC project
author: 
- name: Nick Borcherding
  email: ncborch@gmail.com
  affiliation: Washington University in St. Louis, School of Medicine, St. Louis, MO, USA
date: "August 1, 2020"
output:
  BiocStyle::html_document:
    toc_float: true

---

```{r, echo=FALSE, results="hide", message=FALSE}
knitr::opts_chunk$set(error=FALSE, message=FALSE, warning=FALSE)
library(BiocStyle)
```


# Preperation

## Loading Libaries

In general I like to load libraries here that we will use universally, and then call other libraries when we need them in the code chunks that are relevant. 

```{r}
suppressPackageStartupMessages(library(Seurat))
suppressPackageStartupMessages(library(ggplot2))
suppressPackageStartupMessages(library(dplyr))
```

## Loading Data

```{r}
integrated <- readRDS("./data/Processed/integrated_Cluster.rds")
load("./data/Processed/completeMeta.rda")
```

## Selecting Color Palette

```{r setup, include=FALSE}
colorblind_vector <- colorRampPalette(c("#FF4B20", "#FFB433", "#C6FDEC", "#7AC5FF", "#0348A6"))
```

***

# Initial Subsetting and Analysis

```{r} 
#Removing cells with clonotype information
cell <- integrated[[]]
cell <- subset(cell, is.na(cloneType))
cell <- rownames(cell)
integrated <- subset(integrated, cells = cell)
#Subsetting only APC cells
APC <- subset(integrated, Major == "APC")
APC <- NormalizeData(APC)
rm(integrated)
```

We are also going to create a sub-directory for all T cell analyses:
```{r}
dir.create("./DataAnalysis/CellType/APC")
```


## Cell Cycle 

Cell cycle regression as described in the [Satija lab website](https://satijalab.org/seurat/v3.1/cell_cycle_vignette.html).  

```{r}
cc.genes <- Seurat::cc.genes.updated.2019
s.genes <- cc.genes$s.genes
g2m.genes <- cc.genes$g2m.genes
```

Now we can perform cell cycle scoring with the genes. For now, I am not going to regress using the assignments or save the integrated object with the calculations.

```{r}
APC <- CellCycleScoring(APC, s.features = s.genes, g2m.features = g2m.genes, set.ident = FALSE)
```

Like above with the contribution by condition (**type**), we can now look at the phases assignments by cluster and condition. We do not need to scale if we separate **type** into separate bar graphs.

```{r}
dir.create("./DataAnalysis/CellType/APC/CellCycle")
freq_table <- APC[[]]
freq_table <- freq_table[,c("type", "Final_clusters", "Phase")]
freq_table <- subset(freq_table, Phase != "Undecided") #removing undecided phases
freq_table <- freq_table %>%
    group_by(type, Final_clusters, Phase) %>%
    summarise(n = n())

freq_table$Phase <- factor(freq_table$Phase, levels = c("G1", "S", "G2M")) #ordering phases
freq_table$Final_clusters <- factor(freq_table$Final_clusters, levels = c(0,4,10,11,15,16,19,20,22))

ggplot(freq_table, aes(x=Final_clusters, y=n, fill=Phase)) + 
  geom_bar(stat="identity", position="fill", color="black", lwd=0.25) + 
  theme(axis.title.x = element_blank())+ 
    facet_grid(type ~.) + 
   scale_fill_manual(values=colorblind_vector(3)) + 
  theme_classic()
ggsave(path = "./DataAnalysis/CellType/APC/CellCycle", file = "CellCycle_byCluster_byType.pdf", height=6, width=6)

ggplot(freq_table, aes(x=Final_clusters, y=n, fill=Phase)) + 
  stat_summary(geom="bar", position="fill", color="black", lwd=0.25) + 
  theme(axis.title.x = element_blank())+ 
   scale_fill_manual(values=colorblind_vector(3)) + 
  theme_classic()
ggsave(path = "./DataAnalysis/CellType/APC/CellCycle", file = "CellCycle_byCluster.pdf", height=2, width=6)
```


```{r}
APC$Phase <- factor(APC$Phase, levels = c("G1", "S", "G2M")) #ordering phases
DimPlot(object = APC, reduction = 'umap', group.by = "Phase", split.by = "type") +
    scale_color_manual(values=colorblind_vector(3)) 


ggsave(path = "./DataAnalysis/CellType/APC/CellCycle", filename="CellCycle_UMAP_byType.eps", width=10.5, height=3)
```

## Reclustering 

Isolating just the T cells, stripping the Seurat object of non-"RNA" values using *DietSeurat()* and then splitting the object to create a list that can be fed into **SCTransform()**
```{r eval=FALSE}
APC@meta.data$Run <- paste0(APC@meta.data$sample, "_", APC@meta.data$type) #Define the sequence run
APC[["percent.mt"]] <- PercentageFeatureSet(APC, pattern = "^MT-") #Calculate %mito

APC <- DietSeurat(APC, assays = "RNA")
T.list <- SplitObject(object = APC , split.by = "sample")

for (i in 1:length(T.list)) {
    T.list[[i]] <- SCTransform(T.list[[i]], verbose = FALSE, conserve.memory=TRUE)
}
rm(APC)
```

Integrating the data as above with the whole data set.

```{r eval = FALSE}
options(future.globals.maxSize= 2621440000)
features <- SelectIntegrationFeatures(object.list = T.list, nfeatures = 3000)
T.list <- PrepSCTIntegration(object.list = T.list, anchor.features = features,
verbose = FALSE)
T.anchors <- FindIntegrationAnchors(object.list = T.list, normalization.method = "SCT",
anchor.features = features, verbose = FALSE, k.filter = 50)
APC <- IntegrateData(anchorset = T.anchors, normalization.method = "SCT", verbose = FALSE)
rm(T.list)
rm(T.anchors)
saveRDS(APC, file = "./data/Processed/APC_Precluster.rds")
```

Now we can process the integrated APC cells and identify clusters. Calculating the UMAP and finding clusters.

```{r eval = FALSE}
APC <- ScaleData(object = APC, verbose = FALSE)
APC <- RunPCA(object = APC, npcs = 40, verbose = FALSE)
APC <- RunUMAP(object = APC, reduction = "pca", 
    dims = 1:40)
APC <- FindNeighbors(object = APC, dims = 1:40, force.recalc = T)
APC <- FindClusters(object = APC, resolution = 0.8, force.recalc=T)
APC <- readRDS("./data/Processed/APC_FullCluster.rds")
```

# Subclustering Analysis

After the analysis below - Cluster 5 is principally composed of CD3+ T cells. I am going to use the cellselector function in Seurat to remove this cluster and rename the identities.
```{r, eval = F}
plot <- DimPlot(APC, reduction = "umap")
sc5.cells <- CellSelector(plot=plot)
APC <- subset(APC, cells = sc5.cells, invert = TRUE)
APC <- subset(APC, idents = 5, invert = TRUE)
APC <- subset(APC, idents = 20, invert = T)
```

Manually changing the names of the cluster based on the removal above. 

```{r}
APC <- RenameIdents(object = APC, '6' = '5', 
                                                '7' = '6', 
                                                '8' = '7', 
                                                '9' = '8',
                                                '10' = '9', 
                                                '11' = '10', 
                                                '12' = '11', 
                                                '13' = '12',
                                                '14' = '13', 
                                                '15' = '14', 
                                                '16' = '15', 
                                                '17' = '16',
                                                '18' = '17', 
                                                '19' = '18', 
                                                '20' = '19')
```

## Visualizing New Clusters

```{r}
APC@meta.data$newCluster <- APC@active.ident
APC@active.ident <- factor(APC@active.ident, levels = c(0:20))
DimPlot(object = APC, reduction = 'umap', label = T) + NoLegend()
ggsave(path = "./DataAnalysis/CellType/APC", filename="IntegratedObject_byCluster.eps", width=3.5, height=3)
DimPlot(object = APC, reduction = 'umap', group.by = "type") 
ggsave(path = "./DataAnalysis/CellType/APC", filename="IntegratedObject_byType.eps", width=4, height=3)

a <- DimPlot(object = APC, reduction = 'umap', split.by = "type", group.by = "type")  + NoLegend() + NoAxes() + facet_wrap(~type)
a2 <- a + stat_density_2d(a$data, mapping = aes(x = a$data[,"UMAP_1"], y = a$data[,"UMAP_2"]), color = "black") 
a2
ggsave(path = "./DataAnalysis/CellType/APC", filename="IntegratedObject_byType_faceted.eps", width=10.5, height=3)
```


## Calculating the relative proportion of types in each cluster

Here we are using the total number of cell in the integrated data to scale the clusters. 

```{r}
freq_table <- table(APC@active.ident, APC$type)

tTotal <- 12239
pTotal <- 21160
kTotal <- 1778
totals <- c(kTotal,pTotal,tTotal)

ggplot(as.data.frame(freq_table), aes(x=as.factor(Var1), y=Freq, fill=Var2)) + 
  geom_bar(stat="identity", position="fill", color="black", lwd=0.25) + 
  theme(axis.title.x = element_blank()) + 
   scale_fill_manual(values=rev(colorblind_vector(3))) + 
  theme_classic() + 
    xlab("Clusters") + 
    coord_flip() 
ggsave(path = "DataAnalysis/CellType/APC/", file = "relativeContribution_byClusterType_unscaled.pdf", height=5, width=5)

for (i in 1:3) { 
    freq_table[,i] <- freq_table[,i]/totals[i]
}

freq_table <- reshape2::melt(freq_table)

ggplot(freq_table, aes(x=as.factor(Var1), y=value, fill=Var2)) + 
  geom_bar(stat="identity", position="fill", color="black", lwd=0.25) + 
  theme(axis.title.x = element_blank()) + 
   scale_fill_manual(values=rev(colorblind_vector(3))) + 
  theme_classic() + 
    xlab("Clusters") + 
    coord_flip() 
ggsave(path = "DataAnalysis/CellType/APC/", file = "relativeContribution_byClusterType_scaled.pdf", height=5, width=5)


ggplot(freq_table, aes(x=as.factor(Var1), y=value, fill=Var2)) + 
  geom_bar(stat="identity", color="black", lwd=0.25) + 
  theme(axis.title.x = element_blank()) + 
   scale_fill_manual(values=rev(colorblind_vector(3))) + 
  theme_classic() + 
    xlab("Clusters") + 
    coord_flip()
ggsave(path = "DataAnalysis/CellType/APC/", file = "totalContribution_byClusterType_scaled.pdf", height=5, width=5)

```


## Expression Markers in APCs

Like above we also want to look at the expression markers that define clusters and T cell subtypes. First we will make the schex-based reductions. 

```{r}
suppressPackageStartupMessages(library(schex))
APC <- make_hexbin(APC, 40, dimension_reduction = "UMAP")
```

Here is a list I use in APC identification, we can run these first and store them in the selected folder.

```{r}
genelist <- read.delim("./data/immune_sign.txt")
genelist <- c(genelist$M1.Macro, genelist$M2.Macro, genelist$DC1, genelist$DC2, genelist$DC3, genelist$DC4, genelist$DC5, genelist$DC6, genelist$Langerhan, "ITGAX", "LY75", "CD8A", "CD4", "ITGAM", "FCGR1A", "BST2", "TLR9", "HLA-A", "HLA-B", "HLA-C", "HLA-E", "HLA-F", "HLA-G", "CD209")
genelist <- toupper(unique(genelist))


dir.create("DataAnalysis/CellType/APC/markers")
dir.create("DataAnalysis/CellType/APC/markers/selected")

DefaultAssay(APC) <- "RNA"

for (i in seq_along(genelist)) {
        if (length(which(rownames(APC@assays$RNA@counts) == genelist[i])) == 0){
            next() #Need to loop here because plot_hexbin_feature() does not have a built-in function to deal with absence of selected gene
        } else {
        plot <- plot_hexbin_feature(APC, feature = genelist[i], type = "counts", action = "prop_0")+ 
             guides(fill=F, color = F) + 
                scale_fill_gradientn(colors = rev(colorblind_vector(13))) 
        ggsave(path = "DataAnalysis/CellType/APC/markers/selected", file = paste0( genelist[i], "_prop.pdf"), plot, height=3, width=3.25)
        }
    }

```

We can also look at the top markers for each cluster by using the *FindAllMarkers()* call.

```{r}
APC <- NormalizeData(APC)
All.markers <- FindAllMarkers(APC, assay = "RNA", pseudocount.use = 0.1, only.pos = T) 
write.table(All.markers, file = "DataAnalysis/CellType/APC/markers/FindAllMarkers_output.txt", col.names=NA, sep="\t",append=F)
```

Graphing the top markers for each cluster onto the UMAP using schex from the above *FindAllMarkers()* call.

```{r}
All.markers <- read.delim("./DataAnalysis/CellType/APC/markers/FindAllMarkers_output.txt")
top10 <- All.markers %>% group_by(cluster) %>% top_n(n = 10, wt = avg_logFC)
top10 <- top10$gene #just want the IDs

dir.create("DataAnalysis/CellType/APC/markers/TopClusterMarkers")
for (i in seq_along(top10)) {
    plot <- plot_hexbin_feature(APC, feature = top10[i], type = "counts", action = "prop_0")+ 
             guides(fill=F, color = F) + 
                scale_fill_gradientn(colors = rev(colorblind_vector(13))) 
        ggsave(path = "DataAnalysis/CellType/APC/markers/TopClusterMarkers", file = paste0("Top10markers", "_", top10[i], "_prop.pdf"), plot, height=3, width=3.25)
}
```

Graphing the top markers for each cluster using the *DotPlot()* function from the above *FindAllMarkers()* call.

```{r}
top10 <- All.markers %>% group_by(cluster) %>% top_n(n = 10, wt = avg_logFC)
DotPlot(APC, features = unique(top10$gene)) + coord_flip() + 
    scale_color_gradientn(colors = rev(colorblind_vector(11))) + 
    guides(color = F, size = F) +
    scale_size(range = c(0.5,3.5))
ggsave(path = "DataAnalysis/CellType/APC/markers", file = "Top10_dotplot.eps", height=22, width=4)
```

## Correlating Markers with SingleR Output

Unlike the the T cell subsets, there is a lot more variability within myeloid cells, thus we are redoing the singler analysis.

```{r}
library(SingleR)
Average <- AverageExpression(APC, assay = "RNA", return.seurat = T)
expr_matrix <- as.matrix(Average@assays$RNA@counts[,names(Average@active.ident)])
gene_annotation <- data.frame(row.names=rownames(expr_matrix), gene_short_name=rownames(expr_matrix))

singler = CreateSinglerObject(expr_matrix, project.name = "Myo", annot = NULL, min.genes = 200,
  technology = "10X", species = "Human", ref.list = list(), normalize.gene.length = F, variable.genes = "de",
  fine.tune = T, do.signatures = F, do.main.types = T, 
  reduce.file.size = T, numCores =4)
  
singler$seurat = Average 
```

```{r}
library(pheatmap)
library(Rfast)
SingleR.DrawHeatmap2 = function(SingleR,cells.use = NULL, types.use = NULL,
                               clusters=NULL,top.n=40,normalize=F,
                               order.by.clusters=F,cells_order=NULL,silent=F,
                               fontsize_row=9,...) {
    scores = SingleR$scores
  if (!is.null(cells.use)) {
    scores = scores[cells.use,]
  }
  if (!is.null(types.use)) {
    scores = scores[,types.use]
  }
  
  m = apply(t(scale(t(scores))),2,max)
  
  thres = sort(m,decreasing=TRUE)[min(top.n,length(m))]
  
  data = as.matrix(scores)
  
  if (normalize==T) {
      #for (i in 1:nrow(data)) {
         # max <- max(data[i,])
         # min <- min(data[i,])
         # data[,i] <- (data[,i]-min)/(max-min)
     # }
    mmax = rowMaxs(data, value = T)
    mmin = rowMins(data, value = T)
    data = (data-mmin)/(mmax-mmin)
    data = data^3
     
  }
  data = data[,m>(thres-1e-6)]
  
  
  data = t(data)
  
  if (!is.null(clusters)) {
    clusters = as.data.frame(clusters)
    colnames(clusters) = 'Clusters'
    rownames(clusters) = colnames(data)
    
  }
  additional_params = list(...)
  if (is.null(additional_params$annotation_colors)) {
    annotation_colors = NA
  } else {
    annotation_colors = additional_params$annotation_colors
  }
  clustering_method = 'ward.D2'
  if (order.by.clusters==T) {
    data = data[,order(clusters$Clusters)]
    clusters = clusters[order(clusters$Clusters),,drop=F]
    pheatmap(data,border_color=NA,show_colnames=T,
             clustering_method=clustering_method,fontsize_row=fontsize_row,
             annotation_col = clusters,cluster_cols = F,silent=silent, 
             annotation_colors=annotation_colors, color = rev(colorblind_vector(50)))
  } else if (!is.null(cells_order)) {
    data = data[,cells_order]
    clusters = clusters[cells_order,,drop=F]
    pheatmap(data,border_color=NA,show_colnames=T,
             clustering_method=clustering_method,fontsize_row=fontsize_row,
             annotation_col = clusters,cluster_cols = F,silent=silent, 
             annotation_colors=annotation_colors, color = rev(colorblind_vector(50)))
  } else {
    if (!is.null(clusters)) {
      pheatmap(data,border_color=NA,show_colnames=T,
               clustering_method=clustering_method,fontsize_row=fontsize_row,
               annotation_col = clusters,silent=silent, 
               annotation_colors=annotation_colors, color = rev(colorblind_vector(50)))
    } else {
      pheatmap(data[,sample(ncol(data))],border_color=NA,show_colnames=T,
               clustering_method=clustering_method,fontsize_row=fontsize_row,
               silent=silent, annotation_colors=annotation_colors, color = rev(colorblind_vector(50)))
      
    }
  }
}
```

Now we can graph the results by cluster using the newer **SingleR.DrawHeatmap2()** function. There are two data sets in singleR for mice - the first, refereed to #####. There are also two major outputs by cohort *SingleR.single.main* refers to results reduced across cell types, while *SingleR.single* offers finer granularity for cell subtypes.

```{r}
dir.create("./DataAnalysis/CellType/APC/SingleR/")
pdf("./DataAnalysis/CellType/APC/SingleR/CellTypes_complex2.pdf")
SingleR.DrawHeatmap2(singler$singler[[2]]$SingleR.single, top.n = 50, clusters = singler$singler[[2]]$SingleR.single$cell.names, order.by.clusters = F, 
color = rev(colorblind_vector(50)), normalize = T)
dev.off()

pdf("./DataAnalysis/CellType/APC/SingleR/CellTypes_complex1.pdf")
SingleR.DrawHeatmap2(singler$singler[[1]]$SingleR.single, top.n = 50, clusters = singler$singler[[1]]$SingleR.single$cell.names, order.by.clusters = F, normalize = T)
dev.off()
```



## Making Cell Type Assignments
```{r}
assign <- read.delim("./data/myeloid_assignments.txt")
meta <- APC[[]]
meta <- meta[,1:25]
meta$barcode <- rownames(meta)
merge <- merge(meta, assign, by.x = "newCluster", by.y = "Cluster")
rownames(merge) <- merge$barcode
merge <- merge[,c("Class", "Assignment", "Macrophage.Specific")]
APC <- AddMetaData(APC, merge)
```

```{r}
merge <- merge(meta, assign, by.x = "newCluster", by.y = "Cluster")

table <- merge %>%
    group_by(orig.ident, type, Assignment) %>%
    summarise(n = n()) %>%
    mutate(freq = n / sum(n)) 

table$type <- factor(table$type, levels = c("P", "K", "T"))
ggplot(table, aes(x=Assignment, y=freq, fill = type)) +
    geom_boxplot() + 
    facet_grid(.~Assignment, scales = "free_x") + 
    theme_classic() + 
    scale_fill_manual(values = rev(colorblind_vector(3)))

ggsave("./DataAnalysis/CellType/APC/CellType_proportion_byType.pdf", height=2, width=6)
```
```{r}
unique <- unique(merge$Assignment)

for (i in seq_along(unique)) {
    tmp <- subset(table, Assignment == unique[i])
    aov <- aov(tmp$freq~ tmp$type)
    print(summary(aov))
    
}
```

# ssGSEA for antigen presentaty

Setting up the directory and loading the escape R package. This is a package we built for single-cell GSEA analysis - can find it on [GitHub](https://github.com/ncborcherding/escape) and has been accepted for Bioconductor.

```{r}
GS <- getGeneSets(library = c("C5"))
MHC <- which(grepl("MHC", names(GS)))
GS2 <- GS[MHC]
ES <- enrichIt(obj = APC, gene.sets = GS2, groups = 1000, cores = 2)

meta <- APC[[]]
meta <- merge(meta, ES, by = "row.names")
meta <- meta[,c("newCluster", names(GS2))]

melted <- reshape2::melt(meta, id.vars = c("newCluster"))
meanvalues <- melted %>%
  group_by(newCluster, variable) %>%
  summarise(median(value))

matrix <- reshape2::dcast(meanvalues, newCluster ~ variable)
rownames(matrix) <- matrix[,1]
matrix <- matrix[,-1]

pdf("./DataAnalysis/CellType/APC/CellType_proportion_byTypeMHC_heatmap.pdf", height=5, width=12)
pheatmap::pheatmap(t(matrix[-16,]), color = rev(colorblind_vector(50)), scale = "row", fontsize_row = 3, cluster_rows = T, cluster_cols = T)
dev.off()
```

***

# Differential Gene Expression

## Comparing tumor vs peripheral blood in just the ccRCC samples

### Overall comparison

```{r}
dir.create("DataAnalysis/CellType/APC/DGE")
TumorOnly <- subset(APC, orig.ident == "P1" | orig.ident ==  "P2" | orig.ident ==  "P3")
TumorOnly  <- NormalizeData(TumorOnly , assay = "RNA")
TumorOnly <- SetIdent(TumorOnly, value = TumorOnly@meta.data$type)
markers <- FindMarkers(object = TumorOnly, 
                only.pos = F, ident.1 = "T", ident.2="P", min.diff.pct = -Inf, 
                logfc.threshold = -Inf, min.pct = -Inf)
markers$names <- rownames(markers)
markers  <- markers  %>%
        mutate(Difference = pct.1 - pct.2)
mat <- data.frame(markers)
  write.table(mat, file="./DataAnalysis/CellType/APC/DGE/Overall_differentialMarkers_TvPB.txt",
                sep="\t",append=F, row.names = FALSE)
```

### Overall comparison tumor vs normal kidney parenchyma

```{r}
dir.create("DataAnalysis/CellType/APC/DGE")
subset <- subset(APC, Run == "P1_T" |  Run == "P2_T" |  Run == "P3_T" |  Run == "N2_K" |  Run == "N3_K" |  Run == "N4_K")
subset <- NormalizeData(subset  , assay = "RNA")
subset <- SetIdent(subset, value = subset@meta.data$type)
markers <- FindMarkers(object = subset , 
                only.pos = F, ident.1 = "T", ident.2="K", min.diff.pct = -Inf, 
                logfc.threshold = -Inf, min.pct = -Inf)
markers$names <- rownames(markers)
markers  <- markers  %>%
        mutate(Difference = pct.1 - pct.2)
mat <- data.frame(markers)
  write.table(mat, file="./DataAnalysis/CellType/APC/DGE/Overall_differentialMarkers_TvK.txt",
                sep="\t",append=F, row.names = FALSE)
```

## Comparing peripheral blood of patients to peripheral blood of normal healthy

```{r}
subset <- subset(APC, Run == "P1_P" |  Run == "P2_P" |  Run == "P3_P" |  Run == "N1_P")
subset@meta.data$ind <- ifelse(subset@meta.data$Run == "N1_P", "Con", "Tum")
subset <- NormalizeData(subset  , assay = "RNA")
subset <- SetIdent(subset, value = subset@meta.data$ind)
markers <- FindMarkers(object = subset , 
                only.pos = F, ident.1 = "Tum", ident.2="Con", min.diff.pct = -Inf, 
                logfc.threshold = -Inf, min.pct = -Inf)
markers$names <- rownames(markers)
markers  <- markers  %>%
        mutate(Difference = pct.1 - pct.2)
mat <- data.frame(markers)
  write.table(mat, file="./DataAnalysis/CellType/APC/DGE/Overall_differentialMarkers_TumvCon.txt",
                sep="\t",append=F, row.names = FALSE)
```

## Visualizing the differential gene expression

Here we are going to load the differential results from above into a list - marker_list and loop through the visualizations. 

```{r}

file_list <- list.files("./DataAnalysis/CellType/APC/DGE/", pattern = ".txt")
files <- file.path(paste0("./DataAnalysis/CellType/APC/DGE/", file_list))

marker_list <- list()
for (i in 1:length(files)) {
    marker_list[[i]] <- read.delim(files[i])
}

prefix <- suffix <- stringr::str_split(file_list, "_", simplify = T)[,1]
suffix <- stringr::str_split(file_list, "_", simplify = T)[,3] #Isolate just the cell type
suffix <- stringr::str_remove(suffix, ".txt") #remove .txt
suffix <- stringr::str_remove(suffix, " ") #remoce spaces
names <- paste0(prefix, "_", suffix)
names(marker_list) <- names #assign the cell type to each element of the list
```

This loop will 
+  filter for significant genes that are upregulated or downregulated (not filtered by log-fold change).
+  select the top 20 genes will be selected based on the weight of log-fold change and 4*percent difference.
+  visualize the differential genes as 1) scatter plot with percent difference vs. log-fold change and 2) traditional volcano plot.

```{r}
library(ggrepel)
dir.create("DataAnalysis/CellType/APC/DGE/visualizations")

for (i in seq_along(marker_list)) {
tmp <- marker_list[[i]]
tmp <- tmp %>%
    mutate(Trend = ifelse(p_val_adj <= 0.05 & avg_logFC > 0, "Up",
                    ifelse(p_val_adj <= 0.05 & avg_logFC < 0, "Down", "None")))
filter <- subset(tmp, p_val_adj <= 0.05 & avg_logFC > 0)
top10 <- filter %>% top_n(n =10, wt = avg_logFC + 4*Difference)

ggplot(tmp, aes(x=Difference, y=avg_logFC)) + 
  geom_point(size=0.5, color="#999999") + 
    geom_point(data=subset(tmp, Trend == "Up" | Trend == "Down"), aes(color = Trend), size=0.75) + 
  theme_classic() + 
    geom_hline(yintercept = 0, lty = 2) + 
    geom_vline(xintercept = 0, lty = 2)  + 
    geom_text_repel(data=subset(tmp, names %in% top10$names), aes(label=names), segment.size = 0.25, size=2.5) + 
    scale_color_manual(values = rev(colorblind_vector(2)))+
    guides(color=F)
ggsave(path = "DataAnalysis/CellType/APC/DGE/visualizations", file = paste0("pctDifference_vs_foldchange_", names(marker_list[i]), ".pdf"), height = 2.5, width=2.5) 

tmp <- tmp %>%
    mutate(p_val_adj = ifelse(p_val_adj == 0, min(p_val_adj), p_val_adj)) #mutate the p-values that are lower than detection

ggplot(tmp, aes(x=avg_logFC, y=-log10(p_val_adj))) + 
    geom_point(size=0.5, color="#999999") + 
    geom_point(data=subset(tmp, Trend == "Up" | Trend == "Down"), aes(color = Trend), size=0.75) + 
    theme_classic() + 
    scale_y_sqrt() +
    geom_vline(xintercept = 0, lty = 2) + 
    geom_hline(yintercept = 1.3, lty = 2) + 
    geom_text_repel(data=subset(tmp, names %in% top10$names), aes(label=names), segment.size = 0.25, size=2.5) + 
    scale_color_manual(values = rev(colorblind_vector(2))) + 
    guides(color = F)
ggsave(path = "DataAnalysis/CellType/APC/DGE/visualizations", file = paste0("VolcanoPlot_", names(marker_list[i]), ".pdf"), height = 2.5, width=2.5)
}
```

***

# Only Macrophasges

With the elevated tumor-associated macrophages, we can do further analysis just focusing on them.

```{r}
Macro <- subset(APC, Assignment == "Macrophage")
plot <- DimPlot(Macro, reduction = "umap")
outliers <- CellSelector(plot=plot)
Macro <- subset(Macro, cells = outliers, invert = TRUE)
saveRDS(Macro, file = "./data/Processed/Macrophage_FullCluster.rds")
```


```{r}
Macro <- readRDS("./data/Processed/Macrophage_FullCluster.rds")
dir.create("DataAnalysis/CellType/APC/Macrophage")
colors <- scales::hue_pal()(20)
colors <- colors[c(1,4,8,14,15)]
DimPlot(object = Macro, reduction = 'umap', label = T) + NoLegend() +
  scale_color_manual(values = colors)
ggsave("DataAnalysis/CellType/APC/Macrophage/UMAP_byCluster.eps", width=3.5, height=3)
```

We can also look at the top markers for each cluster by using the *FindAllMarkers()* call. Similar to the APC code above, we are going to look at top markers for the clusters.

```{r}
dir.create("DataAnalysis/CellType/APC/Macrophage/markers/")
Macro <- NormalizeData(Macro)
All.markers <- FindAllMarkers(Macro, assay = "RNA", pseudocount.use = 0.1, only.pos = T) 
write.table(All.markers, file = "DataAnalysis/CellType/APC/Macrophage/markers/FindAllMarkers_output.txt", col.names=NA, sep="\t",append=F)

All.markers <- read.delim("./DataAnalysis/CellType/APC/Macrophage/markers/FindAllMarkers_output.txt")

suppressPackageStartupMessages(library(schex))
Macro <- make_hexbin(Macro, 40, dimension_reduction = "UMAP")

top10 <- All.markers %>% group_by(cluster) %>% top_n(n = 10, wt = avg_logFC)
top10 <- top10$gene #just want the IDs

dir.create("DataAnalysis/CellType/APC/Macrophage/markers/TopClusterMarkers")
for (i in seq_along(top10)) {
    plot <- plot_hexbin_feature(Macro, feature = top10[i], type = "counts", action = "prop_0")+ 
             guides(fill=F, color = F) + 
                scale_fill_gradientn(colors = rev(colorblind_vector(13))) 
        ggsave(path = "DataAnalysis/CellType/APC/Macrophage/markers/TopClusterMarkers", file = paste0("Top10markers", "_", top10[i], "_prop.pdf"), plot, height=3, width=3.25)
}

top10 <- All.markers %>% group_by(cluster) %>% top_n(n = 10, wt = avg_logFC)
DotPlot(Macro, features = unique(top10$gene)) + coord_flip() + 
    scale_color_gradientn(colors = rev(colorblind_vector(11))) + 
    guides(color = F, size = F) +
    scale_size(range = c(0.5,3.5))
ggsave(path = "DataAnalysis/CellType/APC/Macrophage/markers", file = "Top10_dotplot.eps", height=10, width=3)

select <- c("FCGR1A", "CD68", "CD36", "TFRC", "CD20", "MRC1", "ANPEP", "MSR1", "IL3RA", "ITGAM", "CD40", "SLAMF7", "CD3D", "PDCD1LG2", "NRP1", "CD82", "PDCD1", "IFNGR1", "CD7", "CD4", "FCGR2", "FCGR3A", "CD14", "CD163", "SIGLEC1", "CD86", "CD81", "C5AR1", "ICAM1", "HLA-A", "HLA-B", "HLA-C", "HLA-DRA", "HLA-DRB1")


filter <- All.markers[All.markers$gene %in% select,]
top <- filter %>% group_by(cluster) %>% top_n(n = 10, wt = avg_logFC)
top <- rbind.data.frame(top,top10)
top <- top %>% group_by(cluster) %>% top_n(n = 11, wt =  0.5*avg_logFC + 12*(pct.1-pct.2))



DotPlot(Macro, features = unique(top$gene)) + coord_flip() + 
    scale_color_gradientn(colors = rev(colorblind_vector(11))) + 
    guides(color = F, size = F) +
    scale_size(range = c(0.5,3.5))
ggsave(path = "DataAnalysis/CellType/APC/Macrophage/markers", file = "Selected_dotplot.eps", height=7.5, width=2.5)
```


# ssGSEA for TAMs

```{r}
dir.create("./DataAnalysis/CellType/APC/Macrophage/GSEA")
#library(escape)
GS <- getGeneSets(library = c("H"))
ES <- enrichIt(obj = Macro, gene.sets = GS, groups = 1000, cores = 4)
save(ES, file = "./DataAnalysis/CellType/APC/Macrophage/GSEA/H_enrichment.rda")
```

```{r}
sign <- read.delim("./data/immune_sign.txt")[,c(1:20, 63:69)]
full<- as.list(sign)


unique <- names(full) 
list <- list()
for (i in seq_along(unique)) {
    tmp <- full[[i]]
    tmp <- tmp[tmp != ""]
    tmp <- unique(toupper(tmp))
    tmp <- GSEABase::GeneSet(tmp, setName=paste(unique[i]))
    list[[i]] <- tmp
    
}

list <- GSEABase::GeneSetCollection(list)

ES2 <- enrichIt(obj = Macro, gene.sets = list, groups = 1000, cores = 4)
save(ES2, file = "./DataAnalysis/CellType/APC/Macrophage/GSEA/Requested_enrichment.rda")

filter <- cbind(ES,ES2)
meta <- Macro[[]]
meta <- merge(meta, filter, by = "row.names")

heatmap <- meta[, c("newCluster", colnames(filter))]
melted <- reshape2::melt(heatmap, id.vars = c("newCluster"))
meanvalues <- melted %>%
  group_by(newCluster, variable) %>%
  summarise(median(value))

matrix <- reshape2::dcast(meanvalues, newCluster ~ variable)
rownames(matrix) <- matrix[,1]
matrix <- matrix[,-1]

pdf("test.pdf", height=12, width=5)
pheatmap::pheatmap(t(matrix[c(4,5,1,2,3),]), color = rev(colorblind_vector(50)), scale = "row", fontsize_row = 3, cluster_rows = T, cluster_cols = F)
dev.off()

GOI <- c("HALLMARK_INFLAMMATORY_RESPONSE", "HALLMARK_ANGIOGENESIS", "HALLMARK_TNFA_SIGNALING_VIA_NFKB", "M2.Macro", "Proinflammatory", "Glycolysis", "TCA_Cycle", "M1.Macro", "HALLMARK_INTERFERON_GAMMA_RESPONSE", "HALLMARK_IL6_JAK_STAT3_SIGNALING", "HALLMARK_INTERFERON_ALPHA_RESPONSE", "Lipid_mediators", "HALLMARK_HYPOXIA", "HALLMARK_IL2_STAT5_SIGNALING")
sub <- matrix[,colnames(matrix) %in% GOI]

pdf("./DataAnalysis/CellType/APC/Macrophage/GSEA/SelectGSEAheatmap.pdf", height=4, width=4)
pheatmap::pheatmap(t(sub[c(4,5,1,2,3),]), color = rev(colorblind_vector(50)), scale = "row", fontsize_row = 3, cluster_rows = T, cluster_cols = F)
dev.off()
```

Examining the proportion of clusters by ccRCC.
```{r}

freq_table <- meta[,c("orig.ident", "newCluster")]
freq_table <- freq_table %>%
    group_by(orig.ident, newCluster) %>%
    summarise(n = n())
freq_table <- subset(freq_table, orig.ident == "P1" | orig.ident == "P2" | orig.ident == "P3")
freq_table <- freq_table %>%
  mutate(Prop = ifelse(orig.ident =="P1", n/10694, 
                       ifelse(orig.ident == "P2", n/5174, n/9805)))
freq_table$newCluster <- factor(freq_table$newCluster, levels = c(0,3,7,13,14))
ggplot(freq_table, aes(x=newCluster, y=Prop, fill = orig.ident)) + 
         geom_bar(stat="identity", position = "fill") + 
  theme_classic()
ggsave("./DataAnalysis/CellType/APC/Macrophage/patientBreakdown.pdf")

freq_table <- table(Macro@active.ident, Macro$type)

ggplot(as.data.frame(freq_table), aes(x=as.factor(Var1), y=Freq, fill=Var2)) + 
  geom_bar(stat="identity", position="fill", color="black", lwd=0.25) + 
  theme(axis.title.x = element_blank()) + 
   scale_fill_manual(values=rev(colorblind_vector(3))) + 
  theme_classic() + 
    xlab("Clusters") + 
    coord_flip() 
ggsave(path = "DataAnalysis/CellType/APC/", file = "Macro_relativeContribution_byClusterType_unscaled.pdf", height=5, width=5)
```


***

# Trajectory

```{r}
Myeloid <- subset(APC, Class == "Macrophage" | Assignment == "Mo")
library(slingshot)
dir.create("./DataAnalysis/CellType/APC/Macrophage/slingshot")

sce <- as.SingleCellExperiment(Macro) 

sds <- slingshot(Embeddings(Macro, "umap"), clusterLabels = Macro@active.ident, allow.breaks = TRUE, stretch = 0, reducedDim = "UMAP") #Calcualting the trajectory
```

```{r}
library(scales)
#Making plots more siminmar to ggplot outputs of Seurat
cell_pal <- function(cell_vars, pal_fun,...) {
  if (is.numeric(cell_vars)) {
    pal <- pal_fun(100, ...)
    return(pal[cut(cell_vars, breaks = 100)])
  } else {
    categories <- sort(unique(cell_vars))
    pal <- setNames(pal_fun(length(categories), ...), categories)
    return(pal[cell_vars])
  }
}
#We need color palettes Leiden clusters. These would be the same colors seen in the Seurat plots.

cell_colors_clust <- cell_pal(Macro@active.ident, hue_pal())

pdf("./DataAnalysis/CellType/APC/Macrophage/slingshot/Trajectory.pdf", height=4, width=4)
plot(reducedDim(sds), col = cell_colors_clust, pch = 16, cex = 0.25)
lines(sds, lwd = 2, type = 'lineages', col = 'black')
dev.off()
pdf("./DataAnalysis/CellType/APC/Macrophage/slingshot/Trajectory2.pdf", height=4, width=4)
plot(reducedDim(sds), col = cell_colors_clust, pch = 16, cex = 0.5)
lines(sds, lwd = 2, col = 'black')
dev.off()
```

```{r}
nc <- 2
pt <- slingPseudotime(sds)
nms <- colnames(pt)
nr <- ceiling(length(nms)/nc)
pal <- colorblind_vector(100)
par(mfrow = c(nr, nc))
pdf("./DataAnalysis/CellType/APC/Macrophage/slingshot/Trajectory2_pseudotime.pdf", height=4, width=4)
for (i in nms) {
  colors <- pal[cut(pt[,i], breaks = 100)]
  plot(reducedDim(sds), col = colors, pch = 16, cex = 0.5, main = i)
  lines(sds, lwd = 2, col = 'black')
}
dev.off()
```