Integration.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)
```


# Loading Libraries

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

I also like to set a color palette before I begin - this way all the colors are consistent throughout the publication figures.

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

## Normal PBMC

Loading the 10x normal healthy PBMCs - the following series with put "N1_P_" as the prefix on the barcodes and remove the -1. The first step is important as we combine the data to form unique barcodes, the second is pure aesthetics. 

```{r}
N1_P<-  Read10X("data/10x_PBMC/")
colnames(x = N1_P) <- paste('N1_P', colnames(x = N1_P), sep = '_')
colnames(x = N1_P) <- stringr::str_remove(colnames(x = N1_P), "-1")
N1_P <- CreateSeuratObject(N1_P)
```

## Peripheral and Tumor Immune Cell Populations from ccRCC

We will perform the same process for the renal tumor samples, loading both the patient peripheral (P) and tumor-infiltrating (T), immune populations.

```{r}
P1_P <- Read10X("data/ccRCC/GU0700/Peripheral/")
colnames(x = P1_P) <- paste('P1_P', colnames(x = P1_P), sep = '_')
colnames(x = P1_P) <- stringr::str_remove(colnames(x = P1_P), "-1")
P1_P <- CreateSeuratObject(P1_P)

P1_T <- Read10X("data/ccRCC/GU0700/Tumor/")
colnames(x = P1_T) <- paste('P1_T', colnames(x = P1_T), sep = '_')
colnames(x = P1_T) <- stringr::str_remove(colnames(x = P1_T), "-1")
P1_T <- CreateSeuratObject(P1_T)

P2_P <- Read10X("data/ccRCC/GU0744/Peripheral/")
colnames(x = P2_P) <- paste('P2_P', colnames(x = P2_P), sep = '_')
colnames(x = P2_P) <- stringr::str_remove(colnames(x = P2_P), "-1")
P2_P <- CreateSeuratObject(P2_P)

P2_T <- Read10X("data/ccRCC/GU0744/Tumor/")
colnames(x = P2_T) <- paste('P2_T', colnames(x = P2_T), sep = '_')
colnames(x = P2_T) <- stringr::str_remove(colnames(x = P2_T), "-1")
P2_T <- CreateSeuratObject(P2_T)

P3_P <- Read10X("data/ccRCC/GU0715/Peripheral/")
colnames(x = P3_P) <- paste('P3_P', colnames(x = P3_P), sep = '_')
colnames(x = P3_P) <- stringr::str_remove(colnames(x = P3_P), "-1")
P3_P <- CreateSeuratObject(P3_P)

P3_T <- Read10X("data/ccRCC/GU0715/Tumor/")
colnames(x = P3_T) <- paste('P3_T', colnames(x = P3_T), sep = '_')
colnames(x = P3_T) <- stringr::str_remove(colnames(x = P3_T), "-1")
P3_T <- CreateSeuratObject(P3_T)
```

***

#Isolate Immune Cells from normal kidney.

```{r}
normalKidney <- readRDS("./data/processed/Science2018_RCC_immune.rds")
normalKidney <- SplitObject(normalKidney, split.by = "orig.ident")
```

***

# Integrating all the samples

Much in the same way we performed the integration of the normal kidney samples above, we will also do this now across the the single-cell immune samples, by first creating a list and then passing that list for SCT transformation and integration.

```{r}
options(future.globals.maxSize= 4194304000) #Need this to transfer transformation so increasing from 500 Mb to 4 Gb - math: 4000*1024^2 bytes

list <- list(P1_P, P1_T, P2_P, P2_T, P3_P, P3_T, N1_P, normalKidney[[1]], normalKidney[[2]], normalKidney[[3]])

for (i in 1:length(list)) {
    list[[i]] <-  suppressMessages(SCTransform(list[[i]], verbose = FALSE))
}

select.features <- SelectIntegrationFeatures(object.list = list, nfeatures = 3000)
list <- PrepSCTIntegration(object.list = list, anchor.features = select.features, 
    verbose = FALSE)



anchors <- FindIntegrationAnchors(object.list = list, normalization.method = "SCT", 
    anchor.features = select.features, verbose = FALSE)
integrated <- IntegrateData(anchorset = anchors, normalization.method = "SCT", 
    verbose = FALSE)
rm(list)
rm(anchors)

dir.create("data/Processed")
saveRDS(integrated, file = "data/Processed/integrated_PreClustering.rds")
```

Calculating the UMAP and finding clusters.

```{r}
integrated <- ScaleData(object = integrated, verbose = FALSE)
integrated <- RunPCA(object = integrated, npcs = 40, verbose = FALSE)
integrated <- RunUMAP(object = integrated, reduction = "pca", 
    dims = 1:30)
integrated <- FindNeighbors(object = integrated, dims = 1:40, force.recalc = T)
integrated <- FindClusters(object = integrated, resolution = 0.8, force.recalc=T)


dir.create("DataAnalysis/UMAP")
update_geom_defaults("point", list(stroke=0.1))
DimPlot(object = integrated, reduction = 'umap', label = T) + NoLegend()
ggsave(path = "DataAnalysis/UMAP", filename="IntegratedObject_byCluster.eps", width=3.5, height=3)
DimPlot(object = integrated, reduction = 'umap', group.by = "orig.ident") 
ggsave(path = "DataAnalysis/UMAP", filename="IntegratedObject_byOrig.Ident.eps", width=3.75, height=3)
```


Adding the type of cell (or the origin) where K is normal kidney parenychma, P is peripheral blood and T is Tumor.

```{r}
x <- rownames(integrated[[]])
x <- as.data.frame(stringr::str_split(x, "_", simplify = T))

x <- x[,1:2]
colnames(x) <- c("sample", "type")
rownames(x) <- rownames(integrated[[]])
integrated <- AddMetaData(integrated, x)

DimPlot(object = integrated, reduction = 'umap', group.by = "type") 
ggsave(path = "DataAnalysis/UMAP", filename="IntegratedObject_byType.eps", width=3.75, height=3)

saveRDS(integrated, file = "data/Processed/integrated_Cluster.rds")
integrated <- readRDS( "data/Processed/integrated_Cluster.rds")
```

Looking at the density of distribution for cell types.

```{r}
a <- DimPlot(object = integrated, 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/UMAP", filename="IntegratedObject_byType_faceted.eps", width=10.5, height=3)
```


## Proportion of Clusters by Sample and Type

```{r}
meta <- integrated[[]]
freq_table <- meta  %>% 
  group_by(sample, type, seurat_clusters)  %>% 
  summarise(n=n())

ggplot(freq_table, aes(x=seurat_clusters, y=n, fill = type)) + 
  stat_summary(geom="bar", position = "fill") + 
  theme_classic()
ggsave(path = "DataAnalysis/UMAP", filename="ClusterBreakdown_byType_unscaled.pdf", width=3, height=2)

ggplot(freq_table, aes(x=seurat_clusters, y=n, fill = sample)) + 
  stat_summary(geom="bar", position = "fill") + 
  theme_classic()
ggsave(path = "DataAnalysis/UMAP", filename="ClusterBreakdown_bySample_unscaled.pdf", width=3, height=2)

freq_table <- table(meta$seurat_clusters, meta$type)
for (i in 1:ncol(freq_table)) {
  freq_table[,i] <- freq_table[,i]/sum(freq_table[,i])
}
freq_table <- reshape2::melt(freq_table)

ggplot(freq_table, aes(x=Var1, y=value, fill = Var2)) + 
  geom_bar(stat="identity", position = "fill") + 
  theme_classic()
ggsave(path = "DataAnalysis/UMAP", filename="ClusterBreakdown_byType_Scale.pdf", width=3, height=2)

freq_table <- table(meta$seurat_clusters, meta$sample)
for (i in 1:ncol(freq_table)) {
  freq_table[,i] <- freq_table[,i]/sum(freq_table[,i])
}
freq_table <- reshape2::melt(freq_table)

ggplot(freq_table, aes(x=Var1, y=value, fill = Var2)) + 
  geom_bar(stat="identity", position = "fill") + 
  theme_classic()
ggsave(path = "DataAnalysis/UMAP", filename="ClusterBreakdown_byType_Scale.pdf", width=3, height=2)
```


# Identifying subtypes

## Method 1: Differential Gene Expression

The first step for me is to always look at major markers across clusters, so I will make a folder for differential gene expression (DGE), normalize the *RNA* data, and then use the **FindAllMarkers()** function. This is a generalized function for finding positively expressed genes by cluster (large amount of default filtering will remain intact for now). *Importantly*, RNA data needs to be used over the integrated or sct data as this the former is a reflection of true expression and the latter are values to help with the 2D representation in the UMAP.

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

Graphing markers from the differential genes.

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

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

DefaultAssay(integrated) <- "RNA"
dir.create("DataAnalysis/UMAP/TopClusterMarkers")
for (i in seq_along(top10)) {
        if (length(which(rownames(integrated@assays$RNA@counts) == top10[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(integrated, feature = top10[i], type = "counts", action = "prop_0")+ 
             guides(fill=F, color = F) + 
                scale_fill_gradientn(colors = rev(colorblind_vector(13))) 
        ggsave(path = "DataAnalysis/UMAP/TopClusterMarkers", file = paste0("Top5markers", "_", top10[i], "_prop.pdf"), plot, height=3, width=3.25)
        }
  }

```

### Lineage Marker

#### Loading and organizing the markers

```{r}
dir.create("DataAnalysis/UMAP/LineageMarkers")

file_list <- list.files("./data/markers.genes")
file_list <- file_list[grepl(".txt", file_list)]
files <- file.path(paste0("./data/markers.genes/", file_list))

marker_list <- list()
for (i in 1:length(files)) {
    marker_list[[i]] <- read.delim(files[i], col.names = FALSE)
    marker_list[[i]] <- toupper(unlist(marker_list[[i]]))
}
names <- stringr::str_remove(file_list, ".txt")
names(marker_list) <- names
```

#### Graphing all the genes
```{r}
DefaultAssay(integrated) <- "RNA"

for (i in seq_along(marker_list)) {
    tmp <- as.character(unlist(marker_list[i]))
    for (j in seq_along(tmp)) {
        if (length(which(rownames(integrated@assays$RNA@counts) == tmp[j])) == 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(integrated, feature = tmp[j], type = "counts", action = "prop_0")+ 
             guides(fill=F, color = F) + 
                scale_fill_gradientn(colors = rev(colorblind_vector(13))) 
        ggsave(path = "DataAnalysis/UMAP/LineageMarkers", file = paste0(names(marker_list)[i], "_", tmp[j], "_prop.pdf"), plot, height=3, width=3.25)
        }
    }
}
```

## Method 2: Singler

Singler is a very cool package: it uses large cohorts of isolated bulk sequencing to correlate with single-cell data and makes ID-ing the cell type really intuitive. It allegedly works with Seurat, but I am not the largest fan of its work, so I have made some customization to help.

As opposed to calculating the signatures across every single cell, the first step is to calculate mean expression by cluster using the **AverageExpression()** function in Seurat. Then we will make an expression matrix and load that into the **CreateSinglerObject()** function.

```{r}
library(SingleR)
Average <- AverageExpression(integrated, 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))
```

*Warning:* This is not the greatest call function, sometimes when I delete the defaults the function stops working, so some of these seem random, but they're there for sanity

```{r include=FALSE}
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/SingleR")

pdf("./DataAnalysis/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/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()

pdf("./DataAnalysis/SingleR/CellTypes_simple1.pdf")
SingleR.DrawHeatmap2(singler$singler[[1]]$SingleR.single.main, top.n = 15, clusters = singler$singler[[1]]$SingleR.single$cell.names, order.by.clusters = F, normalize = T)
dev.off()

pdf("./DataAnalysis/SingleR/CellTypes_simple2.pdf")
SingleR.DrawHeatmap2(singler$singler[[2]]$SingleR.single.main, top.n = 15, clusters = singler$singler[[2]]$SingleR.single$cell.names, order.by.clusters = F, normalize = T) 
dev.off()
```

## Method 3: Attaching TCR data

A major issue in the differentiation of cell types is the difference between NK cells and T cells, with a lot of crossover between Th1/CTL expression in the latter. One clear way to differentiate is to use our VDJ sequencing data to identify clusters with prominent TCR recovery. From there, we can say these are more definitively T cells. 

```{r}
library(scRepertoire)
```

### Loading the VDJ data

```{r}
P1_P_contigs <- read.csv("./data/VDJ/P_700_contigs.csv")
P1_T_contigs <- read.csv("./data/VDJ/T_700_contigs.csv")
P2_P_contigs <- read.csv("./data/VDJ/P_744_contigs.csv")
P2_T_contigs <- read.csv("./data/VDJ/T_744_contigs.csv")
P3_P_contigs <- read.csv("./data/VDJ/P_715_contigs.csv")
P3_T_contigs <- read.csv("./data/VDJ/T_715_contigs.csv")
N1_P_contigs <- read.csv("./data/VDJ/vdj_v1_hs_pbmc_t_filtered_contig_annotations.csv")
```

### Matching the Seurat and Contig Barcodes

New integration steps for Seurat have made this a little more tricky - the integration adds a _Number to the end of each sample in the Seurat object - we will need to remove this first. 

```{r}
list <- list(P1_P_contigs, P1_T_contigs, P2_P_contigs, P2_T_contigs, P3_P_contigs, P3_T_contigs, N1_P_contigs)

#Remove the -1 from the end of the barcodes
for (i in seq_along(list)) {
  list[[i]][,"barcode"] <- stringr::str_remove(list[[i]][,"barcode"], "-1")
}

#Remove Prefixes of the ccRCC samples
for (i in 1:6) {
  list[[i]][,"barcode"] <- stringr::str_split(list[[i]][,"barcode"], "_", simplify = T)[,3]
}
```

### Organizing TCR data and adding to the Seurat meta data

```{r}
combined <- combineTCR(list, samples = c("P1", "P1", "P2", "P2", "P3", "P3", "N1"), ID = c("P", "T", "P", "T", "P", "T", "P"), cells = "T-AB")
integrated <- combineExpression(combined, integrated)

#Organizing the order of the factor cloneType
integrated@meta.data$cloneType <- factor(integrated@meta.data$cloneType, levels = c("Hyperexpanded (100 < X <= 500)", "Large (20 < X <= 100)", "Medium (5 < X <= 20)", "Small (1 < X <= 5)", "Single (0 < X <= 1)", NA))
update_geom_defaults("point", list(stroke=0.5))
DimPlot(integrated, group.by = "cloneType") + scale_color_manual(values = c(colorblind_vector(5)), na.value="grey")
ggsave(path = "DataAnalysis/UMAP", filename="IntegratedObject_byClonotypeFreq.eps", width=6.5, height=3)
```

```{r}
x <- table(integrated[[]]$seurat_clusters, integrated[[]]$cloneType, useNA = "ifany")

for (i in 1:nrow(x)) {
  x[i,] <- x[i,]/sum(x[i,])
}
x <- data.frame(x)

ggplot(x, aes(x=Var1, y=Freq, fill = Var2)) +
  geom_bar(stat="identity", position = "fill") + 
  scale_fill_manual(values = colorblind_vector(5), na.value="grey") +
  theme_classic()
ggsave(path = "DataAnalysis/UMAP", filename="PercentofClonotype.pdf", width=6, height=3)
```

*** 

# Assigning Major Cell Types

Based on the 3 different methods, there it appears that there needs to be manual edits to the cluster assignments. Specifically two clusters (9 and 10) have subpopulations with TCR recovered. We can edit these manually using the CellSelector() from Seurat. This needs to be run in the console and use the plot function of Rstudio.

```{r, eval=F}
plot <- DimPlot(integrated, reduction = "umap")
sc11.cells <- CellSelector(plot=plot)
meta <- integrated[[]]
meta <- meta[rownames(meta) %in% sc11.cells, ]
meta <- subset(meta, !is.na(cloneType))
sc11.cells <- rownames(meta)

Idents(integrated, cells = sc11.cells) <- 21
```

## Rerun some previous visualizations 

Now we can regraph some of the previous visualization and write over the former versions with the newer cluster assignments. For instance, not only the major UMAP, but also the contributions to each cluster and the SingleR estimates.

```{r}
integrated@meta.data$Final_clusters <- Idents(integrated)
integrated@active.ident <- factor(integrated@active.ident, levels = 0:21)
integrated@meta.data$Final_clusters <- factor(integrated@meta.data$Final_clusters, levels = 0:21)

update_geom_defaults("point", list(size=1, alpha =1, stroke = 0.1))
DimPlot(object = integrated, reduction = 'umap', label = T) + NoLegend()
ggsave(path = "DataAnalysis/UMAP", filename="IntegratedObject_byCluster.eps", width=3.5, height=3)

meta <- integrated[[]]
freq_table <- meta  %>% 
  group_by(sample, type, Final_clusters)  %>% 
  summarise(n=n())

ggplot(freq_table, aes(x=Final_clusters, y=n, fill = type)) + 
  stat_summary(geom="bar", position = "fill") + 
  theme_classic()
ggsave(path = "DataAnalysis/UMAP", filename="ClusterBreakdown_byType_unscaled.pdf", width=4, height=2)

ggplot(freq_table, aes(x=Final_clusters, y=n, fill = sample)) + 
  stat_summary(geom="bar", position = "fill") + 
  theme_classic()
ggsave(path = "DataAnalysis/UMAP", filename="ClusterBreakdown_bySample_unscaled.pdf", width=4, height=2)

freq_table <- table(meta$Final_clusters, meta$type)
for (i in 1:ncol(freq_table)) {
  freq_table[,i] <- freq_table[,i]/sum(freq_table[,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") + 
  theme_classic()
ggsave(path = "DataAnalysis/UMAP", filename="ClusterBreakdown_byType_Scale.pdf", width=4, height=2)

freq_table <- table(meta$Final_clusters, meta$sample)
for (i in 1:ncol(freq_table)) {
  freq_table[,i] <- freq_table[,i]/sum(freq_table[,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") + 
  theme_classic()
ggsave(path = "DataAnalysis/UMAP", filename="ClusterBreakdown_bySample_Scale.pdf", width=4, height=2)

library(SingleR)
Average <- AverageExpression(integrated, 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 

pdf("./DataAnalysis/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/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()

pdf("./DataAnalysis/SingleR/CellTypes_simple1.pdf")
SingleR.DrawHeatmap2(singler$singler[[1]]$SingleR.single.main, top.n = 15, clusters = singler$singler[[1]]$SingleR.single$cell.names, order.by.clusters = F, normalize = T)
dev.off()

pdf("./DataAnalysis/SingleR/CellTypes_simple2.pdf")
SingleR.DrawHeatmap2(singler$singler[[2]]$SingleR.single.main, top.n = 15, clusters = singler$singler[[2]]$SingleR.single$cell.names, order.by.clusters = F, normalize = T) 
dev.off()

a <- DimPlot(object = integrated, reduction = 'umap', split.by = "type", group.by = "Final_clusters")  + 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/UMAP", filename="IntegratedObject_byType_faceted.eps", width=10.5, height=3)
```

## Attaching the assignments

In the data folder, there is a tab-delimited document called ccRCC_cellAssign, which has the 21 clusters and the major and minor assignments by clusters. We cane take this document, merge with the meta data, and then attach is to the Seurat object.

```{r}
cellAssign <- read.delim("./data/ccRCC_cellAssign.txt")
meta <- integrated[[]]
meta$barcode <-  rownames(meta)
meta <- merge(meta, cellAssign, by.x = "Final_clusters", by.y="Cluster")
add <- meta[,c("Major", "Minor")]
rownames(add) <- meta$barcode
integrated <- AddMetaData(integrated, add)
```

Visualization of major and minor assignments
```{r}
DimPlot(object = integrated, reduction = 'umap', label = T, group.by = "Major")
ggsave(path = "DataAnalysis/UMAP", filename="IntegratedObject_byMajorCellType.eps", width=4, height=3)

DimPlot(object = integrated, reduction = 'umap', group.by = "Minor")
ggsave(path = "DataAnalysis/UMAP", filename="IntegratedObject_byMinorCellType.eps", width=4.85, height=3)
```

***

# Saving Final Seurat Object and Meta Data
```{r}
saveRDS(integrated, file = "./data/Processed/integrated_Cluster.rds")
fullMeta <- integrated[[]]
save(fullMeta, file = "./data/Processed/completeMeta.rda")
```

***

# Summary Statistics for data

```{r}
integrated@meta.data$Run <- paste0(integrated@meta.data$sample, "_", integrated@meta.data$type)
table(integrated$Run)
```

Getting the mean features by sequencing run
```{r}
meta <- integrated[[]]
meta %>%
  group_by(Run) %>%
  summarise(meanF = mean(nFeature_RNA)) 
```

Visualizing the mean cell types by condition.

```{r}
table <- meta %>%
    group_by(orig.ident, type, Minor) %>%
    summarise(n = n()) %>%
    mutate(freq = n / sum(n)) 

table$type <- factor(table$type, levels = c("P", "K", "T"))
ggplot(table, aes(x=Minor, y=freq, fill = type)) +
    geom_boxplot() + 
    facet_grid(.~Minor, scales = "free_x") + 
    theme_classic() + 
    scale_fill_manual(values = colorblind_vector(3))
ggsave("DataAnalysis/UMAP/CellType_proportion_byType.pdf", height=2, width=6)


table <- meta %>%
    group_by(Minor) %>%
    summarise(n = n())

ggplot(table, aes(x=reorder(Minor, n), y=n)) +
    stat_summary(geom="bar", aes(fill = Minor)) + 
    theme_classic()  + 
    guides(fill = F) + 
    coord_flip()
ggsave("DataAnalysis/UMAP/CellType_sum_byType.pdf", height=3, width=2)
```

Calculating the significance of proportion cells by tissue type.

```{r}
unique <- unique(meta$Minor)

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